Theoretical Framework¶

This section provides the mathematical foundation for the homodyne scattering analysis implemented in the package, based on the theoretical framework developed by He et al. (2024).

Mathematical Foundation¶

Core Correlation Equation¶

The homodyne intensity correlation uses per-angle scaling:

\[c_2(\phi, t_1, t_2) = \text{offset} + \text{contrast} \times \left[c_1(\phi, t_1, t_2)\right]^2\]

with a separable field correlation:

\[c_1(\phi, t_1, t_2) = c_1^{\text{diff}}(t_1, t_2) \times c_1^{\text{shear}}(\phi, t_1, t_2)\]

Diffusion Contribution¶

\[c_1^{\text{diff}}(t_1, t_2) = \exp\!\left[-\frac{q^2}{2} \int_{t_{\min}}^{t_{\max}} D(t') \, dt'\right]\]

where \(t_{\min} = \min(t_1, t_2)\) and \(t_{\max} = \max(t_1, t_2)\).

Shear Contribution¶

\[c_1^{\text{shear}}(\phi, t_1, t_2) = \left[\operatorname{sinc}\!\big(\Phi(\phi, t_1, t_2)\big)\right]^2\]

with

\[\Phi(\phi, t_1, t_2) = \frac{1}{2\pi} \, q \, L \, \cos(\phi_0 - \phi) \, \int_{t_{\min}}^{t_{\max}} \dot{\gamma}(t') \, dt'\]

where \(t_{\min} = \min(t_1, t_2)\) and \(t_{\max} = \max(t_1, t_2)\).

Time-Dependent Transport Coefficients¶

\[D(t) = D_0 \, t^{\alpha} + D_{\text{offset}}\]

\[\dot{\gamma}(t) = \dot{\gamma}_0 \, t^{\beta} + \dot{\gamma}_{\text{offset}}\]

Parameter Sets¶

Static Mode (3 parameters)

\(D_0\): Reference diffusion coefficient [Å²/s]
\(\alpha\): Diffusion time-dependence exponent [-]
\(D_{\text{offset}}\): Baseline diffusion [Å²/s]
Shear parameters set to zero, \(\phi_0 = 0\) (irrelevant)

Laminar Flow Mode (7 parameters)

\(D_0\), \(\alpha\), \(D_{\text{offset}}\)
\(\dot{\gamma}_0\): Reference shear rate [s⁻¹]
\(\beta\): Shear rate time-dependence exponent [-]
\(\dot{\gamma}_{\text{offset}}\): Baseline shear rate [s⁻¹]
\(\phi_0\): Angular offset [degrees]

Experimental Parameters¶

\(q\): Scattering wavevector magnitude [Å⁻¹]
\(L\): Characteristic length scale / gap size [Å]
\(\phi\): Scattering angle [degrees]
\(dt\): Time step between frames [s/frame]
\(\alpha\): Diffusion scaling exponent (dimensionless)
- \(\alpha = 0\): Normal diffusion
- \(\alpha < 0\): Subdiffusion (constrained motion)
- \(\alpha > 0\): Superdiffusion (enhanced motion)
\(D_{\text{offset}}\): Additive offset for baseline correction

Shear Rate¶

The time-dependent shear rate is parameterized as:

\[\dot{\gamma}(t) = \dot{\gamma}_0 \cdot t^{\beta} + \dot{\gamma}_{\text{offset}}\]

where:

\(\dot{\gamma}_0\): Baseline shear rate [units depend on \(\beta\)]
\(\beta\): Shear rate scaling exponent (dimensionless)
\(\dot{\gamma}_{\text{offset}}\): Additive offset [s^-1]

Integral Functions¶

Diffusion Integral¶

The diffusion integral over the time window \([t_1, t_2]\):

\[J(t_1, t_2) = \int_{t_1}^{t_2} D(t') \, dt'\]

Shear Integral¶

The shear integral over the time window \([t_1, t_2]\):

\[\Gamma(t_1, t_2) = \int_{t_1}^{t_2} \dot{\gamma}(t') \, dt'\]

Both integrals are evaluated numerically via cumulative trapezoid on the discrete time grid (see Computational Methods for implementation details).

Numerical Evaluation (NLSQ vs CMC)¶

While the analytic expressions above are useful for intuition, the code evaluates the time integrals numerically. The implementation differs between NLSQ and CMC backends, and also within NLSQ depending on array size.

Cumulative Trapezoid Method (Preferred)

The cumulative trapezoid approach builds a running sum on the discrete time grid:

Build trapezoid averages between adjacent samples: trap_avg[k] = 0.5*(f[k] + f[k+1]).
Cumulative sum of those averages (no dt applied here): cumsum[0]=0 and cumsum[i] = Σ_{k=0}^{i-1} trap_avg[k].
The integral between any two frames is the absolute difference of cumulative sums: |cumsum[i2] - cumsum[i1]| (a smooth abs is used for JAX gradients). The actual dt scaling is applied outside via the precomputed physics factors (wavevector_q_squared_half_dt and sinc_prefactor).

Concrete example (uniform grid, no dt scaling shown):

Samples: f = [f0, f1, f2].
Trapezoid averages: trap_avg = [0.5(f0+f1), 0.5(f1+f2)].
Cumulative sums: cumsum = [0, 0.5(f0+f1), 0.5(f0+f1)+0.5(f1+f2)].
Interval [t0, t2] uses all interior trapezoids: |cumsum[2] - cumsum[0]| = 0.5(f0+f1) + 0.5(f1+f2).
Interval [t1, t2] uses only the last trapezoid: |cumsum[2] - cumsum[1]| = 0.5(f1+f2).

Implementation Differences

Integration Method by Backend and Array Size¶
Backend	Array Size	Integration Method
NLSQ (matrix mode)	n ≤ 2000	Cumulative trapezoid via `create_time_integral_matrix`
NLSQ (element-wise)	n > 2000	Cumulative trapezoid via `trapezoid_cumsum` with fixed grid
CMC	All sizes	Cumulative trapezoid via `trapezoid_cumsum` with index lookup

Detailed Implementation:

NLSQ matrix mode (homodyne/core/jax_backend.py, n ≤ 2000): Builds a full cumulative-trapezoid matrix once on the 1D grid, then uses |cumsum[i] - cumsum[j]| for every meshgrid pair via _create_time_integral_matrix_impl_jax.

NLSQ element-wise mode (homodyne/core/jax_backend.py, n > 2000): Builds a cumulative trapezoid on a fixed grid (MAX_GRID_SIZE=10001 points with exact dt spacing), maps each (t1, t2) pair to grid indices with searchsorted, and computes the integral via cumulative sum differences. This matches the CMC physics exactly.

# Fixed grid with exact dt spacing (JAX JIT requires static shapes)
time_grid = jnp.arange(MAX_GRID_SIZE) * dt
D_grid = compute_D(time_grid, D0, alpha, D_offset)
D_cumsum = trapezoid_cumsum(D_grid)

# Map times to indices and compute integral
idx1 = searchsorted(time_grid, t1, side='left')
idx2 = searchsorted(time_grid, t2, side='left')
D_integral = |D_cumsum[idx2] - D_cumsum[idx1]|

CMC (homodyne/core/physics_cmc.py): Builds the cumulative trapezoid on a dense time grid (spacing = dt), maps each pooled (t1, t2) pair to grid indices with searchsorted, and uses the absolute cumulative difference per pair.

Note

NLSQ and CMC Physics Now Match

As of version 2.6.0, both NLSQ element-wise mode and CMC use the same cumulative trapezoid integration method with searchsorted index lookup. This ensures identical physics output for the same parameters, enabling direct comparison between optimization methods.

Smooth Absolute for Gradients

Both backends apply a smooth absolute function (sqrt(diff**2 + eps)) so gradients stay finite for zero-length intervals.

Parameter Space and Priors (CMC)¶

CMC builds priors from the parameter_space section of the YAML using homodyne.config.parameter_space and homodyne.optimization.cmc.priors:

Bounds are required and shared with NLSQ. Missing bounds fall back to package defaults (contrast [0.0, 1.0], offset [0.5, 1.5], others [0.0, 1.0]). Those bounds are enforced twice: NLSQ clamps the optimizer to them, and CMC samples only from distributions truncated to that interval (NumPyro uses the bounds directly in Uniform/TruncatedNormal draws).
If a prior is not specified for a parameter, a TruncatedNormal is built with mu at the interval midpoint and sigma at one-quarter of the width. If no prior spec is found at runtime, the fallback is Uniform(min, max).
Supported prior types: TruncatedNormal, Uniform, LogNormal, HalfNormal, Normal, BetaScaled (scaled Beta on [min, max], with concentrations inferred from mu/sigma).
Per-angle parameters contrast_i / offset_i inherit the base bounds/priors from a single contrast/offset entry unless per-angle overrides are provided.

Scattering Geometry¶

Wavevector Definition¶

The scattering wavevector magnitude is related to the scattering angle:

\[q = \frac{4\pi}{\lambda} \sin\left(\frac{\theta}{2}\right)\]

where \(\lambda\) is the X-ray wavelength and \(\theta\) is the scattering angle.

Angular Dependence¶

The azimuthal angle \(\phi\) represents the orientation between the flow direction and the scattering wavevector in the detector plane:

\[\phi = \arctan\left(\frac{q_y}{q_x}\right)\]

The flow direction offset \(\phi_0\) accounts for any misalignment between the nominal flow axis and the actual shear direction.

Physical Interpretation¶

Static Systems¶

For systems at equilibrium without flow (\(\dot{\gamma} = 0\)):

\[c_2(t_1, t_2) = 1 + \beta \exp\left[-q^2 J(t_1, t_2)\right]\]

This simplified form captures pure diffusive dynamics without advection.

Flowing Systems¶

For systems under laminar shear flow:

\[c_2(\phi, t_1, t_2) = 1 + \beta \exp\left[-q^2 J(t_1, t_2)\right] \times \text{sinc}^2\left[\frac{qh}{2\pi} \Gamma(t_1, t_2) \cos(\phi - \phi_0)\right]\]

The sinc function introduces angular anisotropy aligned with the flow direction, with the decorrelation rate depending on the shear rate amplitude.

Chi-Squared Objective Function¶

Parameter estimation proceeds by minimizing the chi-squared objective:

\[\chi^2(\boldsymbol{\theta}) = \sum_{i,j} \frac{\left[c_2^{\text{exp}}(\phi_i, t_j) - c_2^{\text{model}}(\phi_i, t_j; \boldsymbol{\theta})\right]^2}{\sigma_{ij}^2}\]

where:

\(\boldsymbol{\theta}\): Parameter vector
\(c_2^{\text{exp}}\): Experimental correlation function
\(c_2^{\text{model}}\): Theoretical model prediction
\(\sigma_{ij}\): Measurement uncertainty at each data point

Uncertainty Quantification¶

Parameter Uncertainties¶

Confidence intervals are computed from the covariance matrix:

\[\boldsymbol{\theta}_{\text{CI}} = \boldsymbol{\theta}_{\text{opt}} \pm t_{\alpha/2} \sqrt{\text{diag}(\mathbf{C})}\]

where \(\mathbf{C}\) is the parameter covariance matrix estimated from the inverse Hessian at the optimum.

Goodness of Fit¶

The reduced chi-squared statistic assesses fit quality:

\[\chi^2_{\text{red}} = \frac{\chi^2}{N - p}\]

where \(N\) is the number of data points and \(p\) is the number of parameters. Values near 1.0 indicate good fit quality; values >> 1 suggest underestimated uncertainties or model inadequacy.

Numerical Implementation¶

The package implements these equations using JAX for high-performance computation:

JIT Compilation: All core computational kernels are JIT-compiled for optimal performance:

compute_g2_scaled: Computes the full \(c_2\) correlation with per-angle scaling
compute_chi_squared: Fast chi-squared evaluation (includes residual calculation)

Vectorized Operations: Time integrals and correlation matrices are computed using vectorized JAX operations for efficient memory access and parallelization.

Numerical Stability: Special handling for edge cases:

Sinc function: Taylor expansion near zero argument
Exponential overflow: Argument clamping for numerical stability
Division by zero: Epsilon regularization

See Computational Methods for implementation details.

References¶

See References and Citations for complete citation information.