Light through clumpy media

A beam crossing a medium that alternates between materials: how much more light gets through than the averaged medium lets through, to every order in the correlation length.

The model

A beam crosses a slab of thickness $L$ made of two materials, A and B. They absorb at rates $\Sigma_A$ and $\Sigma_B$ per unit length and do not scatter. Along the beam the material alternates, and the chord lengths through A and B are exponential with means $\lambda_A$ and $\lambda_B$.

The material seen at path length $s$ is then a two-state Markov chain $y_s$, with path length playing the role of time. Its generator is

Its stationary law is the pair of volume fractions:

The chain forgets its starting material over the correlation length

since the nonzero eigenvalue of $Q$ is $-1/\lambda_c$.

The same chain appears in three dimensions. Cut space with isotropic random (Poisson) planes and give each cell material A with probability $p_A$, independently. Any straight line then meets the planes at the points of a Poisson process, and each new cell draws a fresh material.

Along every line this is the chain above, with $\lambda_A = \lambda_c/p_B$ and $\lambda_B = \lambda_c/p_A$. It is the standard benchmark geometry (Larmier et al. 2017).

Physical examples are broken cloud fields, cloud and clear air (Titov 1990), clumpy interstellar dust (Natta and Panagia 1984; Hobson and Scheuer 1993), and pebble-bed reactor cores, where the chord lengths are not exactly exponential (Vasques and Larsen 2014).

There are three length scales: the correlation length $\lambda_c$, the mean free path $1/\langle\Sigma\rangle$, and the slab thickness $L$. The expansion is in $\lambda_c$ against the other two.

Parameters for the worked example: dense clumps $\Sigma_A = 5$ filling $p_A = 10\%$ of a thin gas $\Sigma_B = 0.05$; slab $L = 10$; correlation lengths $\lambda_c$ from $0.16$ down to $0.01$, in the same units of length.

The quantity

The mean transmission of the slab, averaged over realizations of the medium:

with the entry material drawn from the volume fractions. Write $\tau = \int_0^L \Sigma(y_s)\,ds$ for the optical depth along one ray.

This is the regime-switching survival problem with path length in place of time, forcing $g_i = -\Sigma_i$, and no slow variable.

Averaging

When the materials alternate infinitely fast, each ray sees the volume-averaged opacity. This is the atomic-mix model:

By Jensen's inequality $\mathbb{E}\,e^{-\tau} \ge e^{-\mathbb{E}\tau}$, so atomic mix never overstates transmission. Rays that happen to pass mostly through the thin material carry most of the light.

Reduction to a linear system

Let $a_i(L)$ be the transmission of a ray that enters in material $i$. Conditioning on the first short stretch of path gives

with $\pi = (p_A, p_B)$. So the answer is a matrix exponential:

For two materials its eigenvalues give the classical two-exponential law (Levermore, Pomraning, Sanzo and Wong 1986):

The slow rate $r_-$ governs thick slabs. The fast rate $r_+$ is about $1/\lambda_c$, and its term dies within a few correlation lengths of the entry face.

The first-order rule

The derivation is the engine's, written for this problem. It gives the answer as integrals of the opacity's covariance along the ray, which is what carries over to any number of materials.

Step 1: separate the mean from the shape

Split the vector of conditional transmissions into its average and its shape:

Multiply the linear system on the left by $\pi$. The switching term drops out because $\pi Q = 0$:

The first term is atomic mix. The second couples the opacity to the shape: it is nonzero when rays that start in dense material are, as a group, dimmer than the rest.

Step 2: the shape relaxes quickly

The shape obeys

The term $Q\,w$ is of size $1/\lambda_c$ and dominates. To leading order $Q\,w = \tilde\Sigma$, and the solution with $\pi\cdot w = 0$ uses the group inverse $Q^{\#}$:

Component $i$ is minus the extra optical depth a ray collects by starting in material $i$, before the chain forgets where it started.

Step 3: the Green–Kubo term

Put the shape back into the equation for $s$:

$K$ is a Green–Kubo integral (Kubo 1957) of the opacity autocovariance along the ray. It raises transmission. For two materials the covariance is a single exponential,

so the first-order rule is

Step 4: second order and the entry layer

One more pass through Step 2 gives the next term of the shape:

Its contribution to the decay rate is

a triple correlation. For two materials it equals $\lambda_c^2\,\mathbb{E}_\pi\tilde\Sigma^3$: zero when the materials fill equal volumes, and otherwise of the sign of the skew of the opacity.

The shape above does not match the entry condition $a(0) = \mathbf 1$, which requires $w(0) = 0$: a ray enters in whichever material it enters, with no history. The initial layer removes the mismatch within a few correlation lengths,

and leaves a permanent constant in $\log T$:

The same terms as cumulants

The three corrections are the first cumulants of the optical depth. Since $\log\mathbb{E}\,e^{-\tau}$ is the cumulant series of $\tau$,

For a slab many correlation lengths thick, the variance and third cumulant of $\tau$ are

So $K$ is half the variance of optical depth per unit length, and $C_1$ is the variance missing because the first correlation length of the ray has no past. Each higher cumulant is one power of $\lambda_c$ smaller.

Transmission written out

For two materials every term evaluates in the chord lengths and cross sections. Through second order:

The neglected terms are of order $\lambda_c^3L$ and $\lambda_c^3$, plus terms of size $e^{-L/\lambda_c}$ from the fast exponential. The decay rate $\langle\Sigma\rangle - K + M_3$ is the start of the Taylor series of the slow rate $r_-$ in the exact law, and $-C_1$ is the start of the series of the log of its weight.

For broken cloud with extinction $\Sigma$, cloud fraction $c$ and clear air that does not extinguish ($\Sigma_A = \Sigma$, $\Sigma_B = 0$, $p_A = c$), the result reads

Where the series converges

The exact slow rate $r_-$ is an analytic function of $\lambda_c$ at fixed volume fractions. It fails to be analytic only where the square root in the exact law vanishes. Multiplying the discriminant by $\lambda_c^2$ gives

Its two complex roots have product $1/(\Sigma_A - \Sigma_B)^2$ and equal modulus, because $(p_B - p_A)^2 + 4p_Ap_B = 1$. So the series in $\lambda_c$, taken to all orders, converges exactly when

The correlation length must be shorter than the mean free path set by the contrast between the materials. The volume fractions do not enter.

What the high orders look like

The same two branch points determine the individual high-order coefficients, including the oscillations seen in the numerical partial sums. For fixed nondegenerate fractions $0<p_A,p_B<1$, put

where the square-root branch is the continuation of $D(0)=1$. Define the dimensionless coefficients by

Singularity analysis at the conjugate points $\zeta$ and $\bar\zeta$ gives

Indeed, $D(x)=C(1-x/\zeta)^{1/2}\{1+O(1-x/\zeta)\}$ near $\zeta$. Since $(r_--\langle\Sigma\rangle)/\Delta=\delta/2+(1-D)/(2x)$, the square-root transfer law gives the first formula.

Meanwhile $w_-$ has a pole of order $1/2$ at each branch point, so $\log w_-$ has the local singularity $-\tfrac12\log(1-x/\zeta)$, giving the second formula. This is the coefficient-transfer method of Flajolet and Odlyzko (1990).

Thus the slow-rate coefficients have a $k^{-3/2}$ envelope, whereas the boundary-weight coefficients have the larger $k^{-1}$ envelope. The fractions control the oscillatory phase $\theta$, but both exponential growth rates remain one in $x$, or $|\Delta|$ in the original correlation-length variable. The error estimates are for fixed positive fractions and are not uniform as one material fraction tends to zero.

What the all-orders series sums to

The power series describes the slow outer mode, not the entire finite-slab answer. This distinction is exact. Put $\Delta=\Sigma_A-\Sigma_B$, $\delta=p_B-p_A$, $x=\lambda_c\Delta$, and

The spectral gap and the two nonnegative weights in the exact law are

Here $w_-+w_+=1$. Therefore the full finite-slab logarithm has the exact decomposition

The Taylor series of $r_-$ and $\log w_-$ both converge for $|x|<1$ and sum to the first two terms. The last term is the entry layer.

For fixed $L>0$ it is smaller than every power of $\lambda_c$, but it is not zero: the expansion is a convergent outer series plus a beyond-all-orders exponential. Thus an infinite algebraic expansion does not equal the exact finite-$L$ transmission unless the fast mode is retained.

This statement is not uniform in a thin-slab limit. If $L/\lambda_c$ stays bounded, the displayed exponential need not be small, and the exact two-exponential law should be used. This is the finite-distance version of the small-correlation-length/entry-layer distinction in stochastic transport (Malvagi, Levermore and Pomraning 1989).

The thin-slab crossover

The nonuniform regime also has a closed asymptotic law. Fix $0<p_A,p_B<1$ and a bounded dimensionless thickness $0\le\ell\le H$, and scale the physical slab as $L=\ell\lambda_c$.

Put $v=p_Ap_B$, $\Delta=\Sigma_A-\Sigma_B$, and $\bar\Sigma=p_A\Sigma_A+p_B\Sigma_B$. Then, uniformly for $\ell\in[0,H]$,

The three pieces of the second-order coefficient have different origins. The slow eigenvalue contributes $v\Delta^2\ell\lambda_c^2$, the slow-mode weight contributes $-v\Delta^2\lambda_c^2$, and the fast mode contributes $v\Delta^2e^{-\ell}\lambda_c^2$. In particular, at $\ell=0$ the entry layer cancels the nonzero outer boundary shift exactly, as it must because $T(0)=1$.

To prove the formula, expand the exact quantities jointly at $x=\lambda_c\Delta=0$:

All factors are analytic jointly in $(x,\ell)$ on a neighborhood of $\{0\}\times[0,H]$, which makes the remainder uniform on the stated compact interval. The result is not uniform as $H\to\infty$, as a material fraction degenerates, or as the opacity contrast varies with $\lambda_c$.

The certificate evaluates 161 values of $\ell$ in $[0,8]$. Under five halvings of $\lambda_c$, the maximum log-transmission error has observed order $2.99$, while the scaled entry layer $\log(1+\text{relative layer})/(\lambda_c\Delta)^2$ approaches $ve^{-\ell}$ with order $0.98$. Thus the crossover formula quantitatively joins the zero-thickness cancellation to the fixed-slab beyond-all-orders regime.

Any number of materials

Nothing in Steps 1 to 4 used two materials. For a chain with any number of materials, generator $Q$ and volume fractions $\pi$, the same result holds with matrix coefficients:

Products of vectors are taken entry by entry. $K$ is the one-entry case of the Green–Kubo matrix on the any-equation page, with a single switched parameter, the opacity, multiplying the operator $-1$.

For a Poisson tessellation whose cells take material $i$ with probability $\pi_i$, the generator is $Q = (\mathbf 1\pi - I)/\lambda_c$ and $Q^{\#}\tilde\Sigma = -\lambda_c\tilde\Sigma$. Every coefficient is then a moment of the opacity over the volume fractions:

Structured media are not of this form. The check below uses cores nested inside envelopes inside a thin gas, so that a ray passes gas, envelope, core, envelope, gas and never goes from gas straight into a core:

Envelopes and cores have mean chord $\ell$, and an envelope is left towards the core or the gas with equal probability.

Every order

With constant opacities the engine returns every term of the outer series as a constant. Its output through order $n$ has the form

where $\Sigma_{\mathrm{eff}}^{(n)}$ and $C^{(n)}$ are the Taylor polynomials, of degree $n$ in $\lambda_c$, of the slow decay rate and of the log of its weight. The first three terms are the result box above.

For any number of materials the engine computes the rest in the same way. This is an outer expansion for fixed positive $L$; the exact binary error after summing every algebraic order is the positive entry-layer term bounded above.

Dense clumps in thin gas

Clumps with $\Sigma_A = 5$ fill 10% of a gas with $\Sigma_B = 0.05$, arranged as a coloured Poisson tessellation, across a slab of thickness $L = 10$. Atomic mix gives optical depth $5.45$ and transmission $0.0043$ at every correlation length. The clump chord is $\lambda_A = \lambda_c/0.9$.

At $\lambda_c = 0.08$ a clump is crossed with optical depth $0.44$, yet the slab passes 3.7 times the light that atomic mix predicts. At $\lambda_c = 0.16$ the factor is 8. The Green–Kubo term alone overshoots there, and higher orders are needed; the series converges slowly because $\lambda_c\lvert\Sigma_A - \Sigma_B\rvert = 0.79$ is close to 1.

The closed form at $\lambda_c = 0.04$, term by term:

Checks

The two-exponential law agrees with the matrix exponential to $1\times10^{-15}$ on four parameter sets. The closed-form coefficients for two materials, for the $N$-colour tessellation and for broken cloud agree with the matrix expressions of the $N$-material box to $5\times10^{-15}$ relative, and the closed-form $\log T$ at orders 0, 1 and 2 agrees with the engine to $1\times10^{-15}$. The first three cumulants of the optical depth, computed from its moments by a block matrix exponential, agree with $\langle\Sigma\rangle L$, $2KL - 2C_1$ and a growth of $6M_3$ per unit length to $6\times10^{-11}$ relative.

Errors in $\log T$ for the worked example, against the matrix exponential, as the correlation length halves:

Each order gains one power of $\lambda_c$, as the last row shows.

The radius of convergence. Taylor coefficients of the slow rate, computed by a Cauchy integral, grow like $4.947^k$, against $\lvert\Sigma_A - \Sigma_B\rvert = 4.95$.

For the five worked correlation lengths, the exact relative entry-layer bounds are respectively $1.34\times10^{-48}$, $1.84\times10^{-75}$, $9.69\times10^{-130}$, $1.51\times10^{-238}$ and below floating-point range. Hence the visible errors in the table are algebraic truncation errors, not an unresolved fast mode.

Inside the radius the engine's partial sums close in on the answer, with oscillations because the singularities are complex. Outside it they blow up:

Three nested phases, with the $N$-material formulas and the engine, against the matrix exponential:

Errors in $\log T$:

Monte Carlo. Rays of length 10 in random directions are traced through realizations of a coloured isotropic Poisson plane tessellation in three dimensions, 50 rays per realization, with standard errors taken across realizations. Chord-length sampling draws the material sequence along the ray directly, 400,000 rays per case.

Both agree with the matrix exponential within three standard errors, and neither is near atomic mix.

Prior art

The exact transmission of a purely absorbing binary Markovian mixture is classical (Levermore, Pomraning, Sanzo and Wong 1986). So is the Levermore–Pomraning model of coupled transport equations for the material-conditioned intensities, which is exact for pure absorbers (Pomraning 1991), and atomic mix as its limit. Standard benchmarks are Adams, Larsen and Pomraning (1989), Brantley (2011) and Larmier et al. (2017).

Small-correlation-length analyses exist: Malvagi, Levermore and Pomraning (1989), Levermore, Wong and Pomraning (1988) for general chord distributions, Su and Pomraning (1994) with scattering, and Pomraning (1996). Kiedrowski and Vu (2025) give a Gaussian law for the distance travelled in each material, which is the central-limit form of the first correction. The first-order term is one line from the exact law and is not claimed as new.

For clouds, Titov (1990) and Malvagi, Byrne, Pomraning and Somerville (1993) develop Markovian models of broken cloud fields. For interstellar dust, Hobson and Scheuer (1993) solve an $N$-phase Markov slab with scattering. Kostinski (2001) shows that spatial correlation makes attenuation slower than exponential.

The fast-switching expansion adds the form of the answer: the correction to atomic mix as a Green–Kubo integral of the opacity autocovariance for any number of materials, the second-order rate and the entry-layer constant as a triple correlation and a first moment of the covariance, and every order from one recursion. The radius of convergence for two materials follows from the exact law in two lines.

Scope

The rule needs the material along the ray to be a Markov chain in path length. With scattering it is not: a scattered particle returns to regions it has already crossed, and the material it meets there is correlated with its past.

The Levermore–Pomraning model is then also approximate. Non-exponential chord lengths, as in packed beds, likewise leave the Markov setting.

Certificate: verify_clumpy_media.py.