Light scattering as a Poisson process and first passage probability

The Kubelka-Munk equations describe one-dimensional transport with scattering and absorption. The reflectance for a semi-infinite slab is the Laplace transform of the distribution of the photon path length lambda. It is determined by the first passage probability of an alternating random walk after np peaks. The first-passage probability as a function of the number of peaks is a path-length distribution-free combinatoric expression involving Catalan numbers. The conditional probability as a function of lambda and np, is a Poisson process. We present a novel demonstration that the probability of first-passage of a random walk is step-length-distribution-free. These results are verified with two iterative calculations, one using the properties of Volterras composition products and the other via an exponential distribution. A third verification is based on fluctuation theory of sums of random variables. Particle trajectories with scattering and absorption on the real half-line are mapped into a random walk on the integer number line in a lattice model, therefore connecting to path combinatorics. Including a separate forward scattering Poisson process results in a combinatoric expression related to counting Motzkin paths.


Introduction
The simplest solution of the radiative transfer equation is for a one-dimensional flux traveling perpendicular to a plane-parallel layer of absorbing and scattering medium with isotropic radiation intensity over the forward and backward hemispheres.
The solutions of radiation transport equations, obtained during the previous century, have been applied to a broad variety of practical situations from neutron diffusion, optical tomography, spreading of infra-red and visible light in the atmosphere to prints on paper.
One-dimensional radiative transfer can be solved by the well-known two-flux approximation proposed independently by Schuster [14] and Schwarzschild [15] in astronomy, Darwin [4] and Hamilton [6] in crystallography, and Kubelka and Munk [10]. in graphic-art and print quality. More recently it had been used by Youngquist, Carr and Davies [20] in optical coherence tomography and Haney and van Wijk [7] in geology.
Simon and Trachsler [16], in their paper A random walk approach for light scattering in material gave an explicit expression for the reflectance. They show first that the scattering problem can be treated as a Markov chain involving Narayana polynomials. But this Markov chain does not provide a solution of a first-passage time problem that fits the reflectance calculated with the Kubelka-Munk equations. Then using the compositional optical reflectance and transmittance properties for multilayer specimens, they are able to determine a generating function and find a solution as elementary variants of Chebyshev polynomials. This is another way to interpret the hyperbolic functions of the classical solution of Kubelka-Munk equations.
Wuttke [19], in his paper The zigzag walk with scattering and absorption on the real half line and in a lattice model expands the Darwin-Hamilton equations, (identical to the Kubelka-Munk equations) into a recursion and finds that Catalan numbers gives the recurrence probability as function of scattering order. He assumes that the random character of the zigzag walk comes from the exponentially distributed step length between scattering events. As we demonstrate, it is not necessary to specify the distribution since his result is independent of the distribution of step length.
One of our goals is to describe the Kubelka-Munk solution directly in terms of the statistics of the number of peaks n p and the path length λ of rays. We explore the solution using a mathematically equivalent one-dimensional random walk model with random step lengths between reflection events described by a Poisson process with a rate per unit length S. We demonstrate further that the recurrence probability given by Catalan numbers is independent of the step length distribution. We provide another demonstration of the independence of step length distribution using our formulation of the fluctuation theory introduced by Andersen [2].
Additionally, we include an independent Poisson process for forward scattering with a rate S f . Usually forward scattering is not included in one-dimensional scattering because it does not change the ray propagation. We include it here as a step towards analysis of three-dimensional problems. Therefore we obtain a generalization of the Wuttke zigzag walk.

Traditional solution of the Kubelka-Munk model
Let us consider a homogenous layer with thickness d characterized by its absorption coefficient χ and its scattering coefficient S. In this layer, the incident irradiance I propagates in the positive direction and the reflected irradiance J propagates in the negative direction. Both I and J are functions of the depth x in the layer. Depth 0 corresponds to the layer's boundary receiving the incident irradiance I 0 . Depth d indicates the other boundary. We consider, at an arbitrary depth x, a sub-layer with infinitesimal thickness dx. The effect of the material in a thin element dx on I and J is to: • decrease I by I(S + χ)dx (absorption and scattering) • decrease J by J(S + χ)dx (absorption and scattering) • increase I by JSdx (scattered light from J reinforces I) • increase J by ISdx (scattered light from I reinforces J).
On these assumptions we obtain the system of equations: with solution where β = χ (χ + 2S) −1 and κ = χ (χ + 2S). The coefficients A and B are determined by the boundary conditions at the two surfaces. After some elementary calculations reflectance R and transmittance T of a slab of thickness d are given by: The reflectance of a very thick layer (d → ∞) is: The Kubelka-Munk reflectance, plays a major role in elucidating the connection to combinatorics and random walks.

Distribution of path lengths from the reflectance
The fluxes I and J can be interpreted as ensembles of photons moving in the positive and negative directions. The photons are absorbed at a rate χ and scattered in a Poisson process at a rate S. Each photon path is weighted according to the Beer-Lambert law by the absorption factor e −χλ . The reflectance is therefore the Laplace transform of the path length distribution. Conversely the inverse Laplace transform of R ∞ leads to the path distribution.
We calculate the inverse Laplace transform of the reflectance L −1 χ (R ∞ (S, χ)) by first expanding the reflectance in the scattering order where C(x) is the generating function of the Catalan numbers C n = (2n)! (n!(n + 1)!) −1 . We identify, term-by-term, the path-length distribution of a random walk model with the inverse Laplace transform ofR ∞ (S, χ) The distribution of λ and n p derived from R ∞ (S, χ) is We agree with Wuttke [19] that "while scattering occurs at random, its effect is deterministic when we model the Kubelka-Munk equation with a random walk. Scattering always reverses the direction of motion. Therefore, trajectories form a zigzag walk rather than a drunkard's walk." The random character of the zigzag walk does not come necessarily from an exponentially distributed step length between scattering events. We will demonstrate later in the paper that the effect of the distribution of the step length between scattering events is distribution-free.
With no loss of generality we can include n f forward scattering events at a rate S f thus providing a persistent random walk model compatible with Kubelka-Munk. Therefore the joint probability of the three random variables λ, n p and n f is given by and the joint probability for n p and n f is then given by The combinatorial factor is related to the "Triangular array of Motzkin polynomial coefficients" T (n, k), the number of Motzkin paths of length n with k up steps T (n f + 2n p − 2, n p − 1), where according to the OEIS A055151 [17], 4. Probability of first-passage by convolution

Iterations: building the skeleton
The previous results are derived from the inverse Laplace transform of the reflectance. We now want to do a direct calculation of the trajectory distribution. A skeleton trajectory, a.k.a. zigzag or alternating random walk, starts at 0 and moves in the positive direction until it eventually reflects back towards the origin. If it does not again reverse direction before it reaches the negative half plane, we say it x y x > y x < y Figure 1. Trajectories in the x, y plane left the medium. A trajectory with n p peaks is subject to (2n p − 1) reflections before leaving the medium and is labeled by the index n p .
Consider now the flux of all possible trajectories starting at x = 0 and beginning in the positive direction. The distribution of the first valley is a symmetrical function P 1 (z) with the following properties: A trajectory reflects back towards the origin at a height x. It reflects next time in the positive direction after a distance y. If y > x, then we say the trajectory left the medium. The probability distribution for the height of the second reflection z = x − y is This is a particular case of composition products considered by the Italian mathematician Vito Volterra in 1913. The probability distribution for the location of the following reflections n p = 3, 4, · · · are given by iterative convolutions: Only the trajectories that stay in the medium at step n p , are transferred to step n p +1, giving: Using equations (12) to (15), calculating by induction and using the suitable limits of integration, we obtain the probability of first-passage after peak n p : This explains the relation between the Kubelka-Munk reflection and the generating function of the Catalan numbers.

Dressing the skeleton: connection with combinatorics
The enumeration of lattice paths is a topic in combinatorics, closely related to the study of random walks in probability theory. The ubiquitous presence of Catalan numbers in the joint distribution function P (λ, n p , n f ; S, S f ) suggests a connection with combinatorics. One approach to explain this connection is discretization. Since the statistics of first passage is independent of the distribution of the path length, to discretize the path we just have to integrate over λ. Therefore the discretized skeleton (zigzag) walk is mapped onto a path in a two-dimensional lattice. We expect from this transition to a lattice model to reproduce the analytical result obtained in section (6). The joint probability function is the product of the marginal probability times the conditional probability: The goal is to create new steps by randomly filling the 2n p segments of the skeleton with n f forward random scattering events and to calculate the resulting distribution. Ultimately we want to find the related conditional probability P (n p |n f ). There are 2n p − 1 reflections and m s = 2n p + n f steps.
With our notations we have paths from (2, 0) to 2(n p , n f ) with the constraint that n p ≥ 1. Therefore the number of paths is given by: then the number N C of combinations is given by: We can now write the conditional probability of m s at constant n p with r = S This result has been confirmed by extensive Monte Carlo calculations and is obtained by recursion in section (5). Then the joint probability P JG is given by This is the same as equation ( (10)). This confirms the combinatorial nature (lattice paths enumeration) of the discrete form of the Kubelka-Munk equation.

Analytic calculation of the path distribution
In this section we integrate the scattering and attenuation of an ensemble of light rays. We show directly the connection to Eq.(4), as opposed to the two-step calculation in section (4). We use an exponentially distributed step length consistent with a uniform absorption and scattering rate. Explicit integration will serve as a foundation for exploration of higher dimensions and the effect of inhomogeneities in the media that are important for print quality and medical imaging.
Consider a statistical ensemble of rays of light moving in a one-dimensional diffusive medium on the positive axis. Again, the ray starts at the origin moving upward in the positive direction and after multiple reflections, escapes the media to the negative axis. The distribution of the number of scattering events and the path length for the escaping rays is indicative of the statistical behavior in three dimensions.
A ray traveling through the medium is reflected at a rate of S reflections per unit length and is attenuated at a rate of χ per unit length. The Poisson probability density for reflection of a ray after traveling distance d is where Φ(d) is the Heaviside step function. We will derive information about the distribution of the number of scattering events and the path length when the ray ultimately escapes. We choose to describe the current state of a ray by the direction, the total prior upward movement l, the prior number of peaks n p and the position x.
The upward movement is related to the path length by 2l − x = λ. Later we will add the distribution of forward scattering events as an independent Poisson process based on the path length. The probability density for an upward reflection at a valley at height x with total prior upward movement l at the valley after n p peaks is P V np (x, l). The probability density that the ray escapes with total path length λ in the scattering medium, after the n th p peak is E np (λ). The probability densities for the peaks and valleys are not normalized because the ray may escape before the reflection.
There are constraints on the integrals in iterating from one valley to the next. The peaks are higher than the surrounding valleys. The ray must travel at least the height of the peak to reach it. The path length is an increasing function as the ray progresses. Therefore, as shown in Fig.(2), The sequence of probability densities at the peaks and valleys can be tediously calculated recursively. The probability density at the next peak is found by integrating the probability density at the valley over the height of the previous valley. Substituting for P P n and marginalizing over the peak, we obtain the valley-to-valley transfer.
The functional form of the probability densities at the valleys is evaluated by iterating this expression starting from the initial condition that the ray starts with upward motion from x = 0.
The first peak integral simply involves satisfying the delta functions. Similarly, the first valley integral uses the fact that l is the total upward motion. The probability that the ray escapes is given by dropping the constraint that x is positive and integrating over the negative half space: Iterating from the probability distribution at one valley to that at the next valley, we obtain: Here F np is the total volume of the dimensional space of allowed configurations of the 2n p −2 step path from the origin to a valley at the point x with path length λ = 2l −x.
The same form, an exponential times the volume of path configurations, will apply in higher dimensions and more complicated geometries. Monte Carlo methods can be applied to measure the path configuration volume in these situations.
Integrating over the negative half-space gives the escape probability density as a function of upward path length and the total escape probability after the n th p peak: Integrating over λ gives the escape probability after n p peaks in terms of Catalan numbers. The result is independent of the scattering rate and the formula is identical to Eq. (5) with χ = 0: We can add in absorption because we have the distribution of the path length. The attenuation simply adds χ to S in the exponent: Integrating over the path length again yields Eq. (5) The joint probability for n p and n f is obtained by adding an independent Poisson process for forward scattering This result is identical to equations (10) and (21).

First passage events are distribution-free
The location of an alternating walk alternates between peaks and valleys. First passage for an alternating walk occurs only at a valley where the step number m is even. We start with the finite set of alternating walks A n (c n ) generated by permutations of a set of n step sizes c n . The set c n is an element of the set C n of all sets of lengths.
The first passage events are the subsets F m of all allowed walks that are positive before step m and first become negative at step m. It is possible that a walk never becomes negative, so we add the event F 0 of walks that are positive for the first n steps. The set of n + 1 first passage events is a partition of the finite sample space. We will show that the cardinality of F m is independent of the set of n step sizes c n as long as no sub-walk returns to the origin. We define a boundary subset B n of C n of measure zero, where no sub-walk constructed from c n returns exactly to the origin.
We analyze the changes in event membership as we adjust values in c n continuously by modifying step sizes. For a walk constructed from a set of lengths not in the boundary set, let δ be the closest approach to the origin of any non-empty sub-walk. At step m > 0 the absolute value of the position must be at least δ. If δ > 0 then a step size change with magnitude smaller than δ to any length c k ∈ c n will not cause any change in first passage event membership. First passage event cardinality is therefore locally constant, and so is constant within connected subsets of C n − B n . The set C n − B n is not connected, so we examine changes in event cardinality when crossing the boundary.
Any set of lengths can be reached by changing one length at a time, and so if cardinality of an event does not change when crossing B n , then the event cardinality is constant in C n − B n .
The sequences of lengths in B n are the only elements of C n where walks can leave or enter first passage events with small changes to a length. When changing a length in a walk w m and moving the walk through the boundary point w * m the membership of w m in an event can change. The set of walks generated from permutations or an overall sign change of the steps in the sub-walk w * m also return to the origin. The problem is to show cancellation of the changes in event membership among the walks in this set as lengths are changed so that w m moves through the boundary point w * m . The key step in our approach is to consider uniquely defined pairs of walks w * and w ′ * from the boundary set. Walk w * begins with a positive critical sub-walk w * j as in Fig. (3). The paired walk w ′ * is identical after step j but begins with the time-reversed sub-walk w ′ * j . The sets of positions of the paired walks are the same, but they are in the reverse order for the first j steps. The step sizes in w ′ * m are the same as the step sizes in w * m , but the first j steps occur in the reverse order and with the opposite sign. All the critical walks generated from c * n ∈ B n that could change the cardinality of F m can be uniquely paired in this way. For the first passage problem, the boundary of F m consists of walks that are non-negative for the first m − 1 steps and return exactly to zero at one of the first m steps. Walks that return precisely to the origin on step m and are positive before that step for steps 1, · · · (m − 1) will enter or leave the event with a small change in a length. Similarly, walks that pass zero on step m but return precisely to the origin on step j < m change membership in F m with small length changes.
Suppose a walk from the boundary of F m begins with such a j-step sub-walk w * j which is positive for k < j and returns to 0 at step j. Consider changing a length c k ∈ w * j through the critical value c * k . A small increase in a length that is a downward step in w * j will cause the walk to be in F j while decreasing the same step will cause the walk to be in F m . Similarly, a small increase in a length that is an upward step in w * j will cause the walk to be in F m while decreasing the same step will cause the walk be in F j . Alternatively, if the walk touches zero at step m then increasing a downward length will cause it to be in F m while the paired walk will leave F m .
If a step of length c * k in sub-walk w * j has, say, a negative sign, then it is in w ′ * j with a positive sign. When c * k is changed to a lower value then walk w ′ is in F m and w is not. Similarly, when c k is changed to a higher value, then walk w is in F m and w ′ is not. Thus as c k is changed and the set of lengths crosses B n one of each pair of walks leaves the event F m and the other joins the event. Although the pair of walks switch which one is in event F m , the total contribution of the two walks to the cardinality of the event is the same. The cardinality of F m is therefore unchanged by crossing B n . Examining the complete set of pairs of walks with time-reversed sub-walks thus shows that the cardinality of first passage events, and thus the probability of first passage at step m is independent of the set of n real lengths, except for the boundary set B n of measure 0 in C n . Now averaging over C n of this invariant probability gives the result that the distribution of first passage step for alternating walks is step-size distribution-free.
A critical point is that the sub-walk must have even length for alternating walks so that walks beginning with w and w ′ are both alternating walks in our sample space. First passage always occurs on an even number step for alternating walks, so that is not a problem here. In other types of events on A n the requirement that the subwalks must have even length is important. The argument for first passage statistics for symmetric walks carries through the same as for alternating walks with the exception that the location of the first passage, and the length of the relevant sub-walks, need not be even.
To calculate first passage probabilities for an alternating walk, any set of lengths suffices. Make a convenient choice like any subset of integer powers of 2 where each length is larger than the sum of all smaller lengths. First passage of an alternating walk occurs only in a valley which occurs on even numbered steps m = 2m p .
Given a set of n = 2n p lengths selected from integer powers of 2, the fraction of walks with first passage at step m = 2m p < n is and the fraction that remains nonnegative is The calculation will proceed by induction on n p . The theorem is true for n p = 1 because the probability of first passage in the first valley is 1 2 , as is the probability that the walk stays positive. Either the first step upward is larger or smaller than second step.
Suppose the theorem is true for all m p ≤ n p . Consider alternating walks generated from a set c 2(np+1) of 2(n p + 1) lengths selected from integer powers of 2. Divide the alternating walks generated into a complete set of disjoint subsets where all walks in a subset have the same last two steps. For each of these subsets, the first 2n p steps are all the permutation of the same set of lengths. By induction, the fraction of these in F 2mp for m p ≤ n p is given by the theorem. Similarly, the fraction that stay positive until the last step is given by the theorem.
First passage can occur after the peak n p + 1 only if the walk stayed positive for the first n p valleys and step 2n p is the largest element of the set. The probability of this is The walk will stay positive only if it is positive for the first n p valleys and it does not have a first passage at valley n p + 1, so in agreement with (35). The theorem is true for n p + 1, and so is proved by induction.