Probability density function for random photon steps in a binary (isotropic-Poisson) statistical mixture

Monte Carlo (MC) simulations allowing to describe photons propagation in statistical mixtures represent an interest that goes way beyond the domain of optics, and can cover, e.g., nuclear reactor physics, image analysis or life science just to name a few. MC simulations are considered a “gold standard” because they give exact solutions (in the statistical sense), however, in the case of statistical mixtures their implementation is often extremely complex. For this reason, the aim of the present contribution is to propose a new approach that should allow us in the future to simplify the MC approach. This is done through an explanatory example, i.e.; by deriving the ‘exact’ analytical expression for the probability density function of photons’ random steps (single step function, SSF) propagating in a medium represented as a binary (isotropic-Poisson) statistical mixture. The use of the SSF reduces the problem to an ‘equivalent’ homogeneous medium behaving exactly as the original binary statistical mixture. This will reduce hundreds MC simulations, allowing to obtain one set of wanted parameters, to only one equivalent simple MC simulation. To the best of our knowledge the analytically ‘exact’ SSF for a binary (isotropic-Poisson) statistical mixture has never been derived before.


Introduction
Two probability density functions (pdf) are at the heart of Monte Carlo simulations describing photon propagation in biomedical optics, optics or, in general, particle transport in diffusive media: the "phase function" (PF) and the pdf allowing to describe the probability for a photon to reach a distance s ∈ [s − ds 2 , s + ds 2 ] in the medium, without interactions.For simplicity, we will call the latter pdf "single step function" (SSF).The knowledge of the SSF is fundamental because it is mandatory for the implementation of a MC simulation, and because from the SSF we can extract the main characteristics of photon propagation through the medium 1 .
The specific shape of the SSF is determined, in a very complex way, by the physical characteristics of the investigated media.For this reason, in the majority of the cases, the SSF cannot be derived theoretically, but is estimated from experimental data.Among the different SSFs appearing in the literature, one in particular seems to be more recurrent and applicable in many situations; i.e., the one derived from the so called Beer-Lambert-Bouguer law, with pdf p LB (s; µ t ) = µ t e −µ t s , where µ t is the extinction coefficient of the medium.Historically, the SSF p LB (.) has been determined experimentally and, nowadays, its applicability covers a panoply of different domains [2][3][4] .Due to the large number of physical systems that can be described by p LB (.), some caution may be natural when describing systems with an SSF that appears to deviate from the Eq.(1).In fact, the question may arise whether SSFs different from Eq. (1)  are not simply the result generated by compound media of immiscible materials; where each bunch of material -considering the subject to come in the manuscript, we will call them "tessels" -always satisfies the classical Beer-Lambert-Bouguer law [Eq.( 1)], with their own µ t (for an intuitive example see the schematic in Fig. 1).With this question in mind, many studies have already been proposed in the literature; going from numerical MC simulations to theoretical approaches [5][6][7][8] .Interestingly, this approach also applies to image rendering 9 .
Due to their ability to describe actual physical systems, a particular attention has been given to "random" media obeying mixing statistics.In intuitive words, random media similar to the one appearing in Fig. 1 were generated with the desired statistical law.Then, the exact or approximated physical quantities of interest (e.g., in the present context the SSF) were obtained (theoretically or numerically) on many random realizations of the medium and, finally, the "average value" of each Figure 1.Schematic of a 3D medium (the cubic shape is just for simplicity, and may also be infinite) composed by random tessels of two kinds of materials (red and sky blue) satisfying both Eq (1), with µ t a and µ t b .The random oriented dashed line crosses different areas of random lengths 1 , 1 , . . ., 6 (along the line), satisfying Eq. ( 2), with parameters σ a and σ b .If Eq. ( 2) with parameters σ a and σ b always holds for any dashed random oriented line, the medium is called isotropic.quantity was derived.Unfortunately, the common characteristic of these studies is that they never propose an exact analytical solution for the SSF.This is why, the aim of the present contribution is to derive an exact analytical solution for a well celebrated model: the binary (isotropic-Poisson) statistical mixture.
But, why the knowledge of an exact analytical expression for the SSF may represent any advantage?One of the reasons is that repeated MC simulations of random media, sometimes composed by thousands of tessels of complex shapes, may be extremely computational expensive; in particular if simulations must be repeated a large number of times.The knowledge of an analytical expression for SSF, describing our random medium with mixing statistics, allows one, in principle, to perform only one simulation on a equivalent homogeneous medium, and to obtain the same desired quantities (e.g., transmitted or reflected photon fluxes) as for the original problem.

Theory Isotropic Poisson tesselation with a binary mixture
The three dimensional random medium that we will consider in the present contribution, is a medium where a photon that propagates in a straight line, crosses, alternatively, two kind of homogeneous material (tessels) (see, e.g., Fig. 2).The length of each piece of path that crosses a tessel (e.g., 1 , 2 , 3 in Fig. 2) is given by a probabilistic law.If the probabilistic law remains the same for any direction of propagation of the photon, then the medium is called isotropic.It has been demonstrated that if we want an isotropic medium with the above characteristics, then the probabilistic law must be exponential; i.e. 6,10 , where is the (random) length of the tessel in the direction of the photon propagation, and the constant σ ∈ {σ a , σ b }, depending on the type of tessel we want to generate (red or sky blue in Fig. 2).Observe that the mean tessel length is Moreover, to have an isotropic medium, σ a and σ b must satisfy the following relationship 6, 10 and Figure 2. Simplified 2D schematic representing a typical case of photon propagation in a (binary) medium.A: Entrance point in the medium; B: The photon is scattered and changes its direction; C: The photon is absorbed and stops its propagation.s 1 : first step length; s 2 : second step length.In the present theoretical context and in MC simulations s 1 and s 2 are decomposed in sub-steps; e.g., s 1 = 1 + 2 + 3 + s 4 .Note that s 4 is shorter than the tessel length because the photon is scattered in B. The parameter is the distance from point B to the next (red) tessel.
where (1 − P) is the probability to have a tessel of type "a", and P the probability to have a tessel of type "b".The constant L is a parameter utilized in the construction of the medium.For such a medium in the literature we speak about "isotropic Poisson tesselation with a binary mixture".

Single step function: exact analytical derivation
In this section we will analytically derive the exact SSF, p s mix (s; µ t a , µ t b , σ a , σ b ), for the isotropic Poisson tesselation with a binary mixture model; where µ t a and µ t b are the extinction coefficients of the tessels of type "a" and "b", respectively.The random step s (ballistic propagation) is defined as a straight line of random length, going from the starting point to the point where the photon is absorbed or scattered.Thus, a step may go across many tessels (e.g., steps s 1 or s 2 in Fig. 2).By following the approach presented e.g. in Refs. 1,11  step can be decomposed as a sum of sub-steps lengths -independently generatedcovered by the photon in the "a" and "b" regions.To obtain p s mix (.) we need before to derive some intermediate functions.
This is what is done in the following sub-sections.Complex analytical calculations were performed using MATHEMATICA ® software.
Probability mass function P N (.) Probability for a photon to make a step larger than s is The probability for a photon to make a step larger than a random tessel of length s is (the photon starts at the tessel boundary) Considered the fact that the tessels crossed by a photon (ballistic propagation) have alternate µ t a and µ t b values, the probability to make a series of consecutive steps larger than the correspondent random tessel lengths must be of the form Therefore, the probability P N (µ t a , µ t b , σ a , σ b , N) for a photon to make N consecutive steps larger than the correspondent N consecutive random tessel lengths (the first tessel always has µ t a ) is where ∑ +∞ n=0 P N (µ t a , µ t b , σ a , σ b , N) = 1, as expected.In other words, if N = 0 (first line Eq.( 9)), the photon makes a step shorter than the random length of the first tessel.If N is even, then the last term in Eq. ( 8) is P 1 (µ t b , σ b ).If N is odd, the last term in Eq. ( 8) is P 1 (µ t a , σ a ).This explains the terms in curly brackets in the second and third line of Eq. ( 9).The last photon step always has probability 1 − P 1 (µ t a , σ a ) or 1 − P 1 (µ t b , σ b ), ∀N [Eq.( 9)], because it is where the photon stops (i.e., the last step must be shorter than the last random tessel length).By inserting Eq. ( 7) in Eq. ( 9) we find We derive here the pdf p s 0 (s; µ t , σ ) for a photon step of random lengh s to remain inside a tessel of random length (the photon always starts on the tessel boundary).To do this, we apply an upper cut-off to Eq. (1) at the distance , i.e., where Θ(.) is the Heaviside function.Thus, the pdf p s 0 (s; µ t , α, σ ) is obtained as where the denominator appearing of Eq. ( 12) is the normalization factor.We need also the probability p Tes N (s; µ t a , µ t b , σ a , σ b , N)ds for a photon to make N > 1 steps, of total random length ; where i is the (random) length of the i th tessel, in the photon direction (the function p Tes N (.) is a pdf).To this aim, we start with the pdf p Tes 1 (s; µ t , σ t ) for a single tessel (case N = 1).

Probability density function p Tes
By applying a lower cut-off at the distance − ∆ 2 and an upper cut-off at the distance + ∆ 2 to Eq. (1) (0 < ∆ 1), we express the fact that we want the photon falls at a distance (tessel length in the direction of the photon propagation), i.e., Thus where µ t ∈ {µ t a , µ t b } and where the denominator is the normalization factor allowing to obtain the pdf.Note that Equation ( 14) allows one to treat the particular case of the pdf p Tes N (s; µ t , µ t , σ , σ , N) of N consecutive tessels with same µ t and σ ; i.e.
that can be explicitly written as This is the expected result since the pdf p Tes (.) is the sum of N positive random variables described by the same exponential pdf p Tes If N > 1 is odd, by following the method proposed in Ref. 12 for the solution of the integral, we obtain where and where Thus, where I n (x) is the modified Bessel function of the first kind.Note that the method proposed in Ref. 12 has been developed by imposing some constraints on the parameters µ t a , µ t b , σ a and σ b .σ a and σ b value.However, it is easy to show that the method remains valid for any value of µ t a , µ t b , σ a and σ b .

5/13
Probability density function p Tes N (.): case N > 1 even If N > 1 is even, in the same vein, by following the method proposed in Ref. 12 for the solution of the integral, we get where and and where Thus, Note that, in general, +∞ 0 p Tes N (s; µ t a , µ t b , σ a , σ b , N)ds = 1, ∀N ≥ 1.
Probability density function p s mix (s; .)(single step function) Let's first express the pdf p s N (s; µ t a , µ t b , σ a , σ b , N) for a photon to jump over N tessels and reach a distance s.In practice, the photon stops inside the (N + 1) th tessel due to absorption or scattering.See, e.g., an intuitive drawing in Fig. 2 for the scattering case (segment AB where N = 3).This is obtained as [Eqs.(15), (18) and (24)] where , depending if N is even or odd, respectively.Thus, the pdf p s (s; µ t a , µ t b , σ a , σ b ) for a photon to make a step of length s, independently of the number of tessels N, is obtained as [Eqs.(10) and (30)] Equation (31) has been derived for a starting tessel with parameters µ t a and σ a .However, Eq. ( 31) also works if we permute µ t a with µ t b , and σ a with σ b in the equation (i.e., we chose the parameters of the starting tessel equal to µ t b and σ b ) Finally, Eq. (31) allows us to express in general the pdf p s mix (s; µ t a , µ t b , σ a , σ b ) (single step function) for a photon step of length s, for a medium with a probability 1 − P to have the first tessel with parameters µ t a and σ a , and probability P to have the first tessel with parameters µ t b and σ b ; i.e., p s mix (s; µ t a , µ t b , σ a , σ b , P) = (1 − P)p s (s; µ t a , µ t b , σ a , σ b ) + Pp s (s; µ t b , µ t a , σ b , σ a ).
(32) Equation ( 32) has been implemented in MATLAB ® language.Note that, to obtain Eq. (32), we did not use the isotropy conditions [Eqs.( 4) and ( 5)], thus the pdf p s mix (.) remains valid even if the parameters change with the propagation direction.

Derivation of the albedo for the random medium with mixing statistics
The albedo for a homogeneous medium is usually expressed as µ s /µ t , where µ s is the scattering coefficient, but where, in probabilistic terms, µ s ds is interpreted as the probability to be scattered inside a length ds.Thus, µ s can be seen as the probability of scattering per unit length.In the same way, µ t = µ a + µ s may be interpreted as the probability for a photon to be scattered or absorbed (extinction) per unit length; where µ a is the absorption coefficient.In the present case, the "probabilities" µ s and µ t must be re-calculated, to take into account of the effects of the mixing statistics.To this aim, the probability, P t unit ∆s, to occur in an extinction interaction, inside a very small ∆s, is calculated as: where we have developed in Taylor series the integral for ∆s 1. Equivalently, the probability, P s unit ∆s, for a photon to occur in a scattering event in ∆s is expressed as where we have set µ a a = µ a b = 0. Finally, the albedo Λ mix , to be utilized with the SSF p s mix (.), is obtained as Note that if µ t a = µ t b and µ s a = µ s b , we obtain the classical albedo for homogeneous media.

Derivation of the single step function for photons starting from "fixed" positions
For the sake of completeness, in this section we give the essential equations allowing the reader to derive the SSP of photons starting from fixed positions (this topic may be subject of further studies).It has been shown that when p s mix (.) [Eq.(32)] is not a pure exponential law, photon steps starting from fixed positions (e.g., from a fixed light source or after the reflection on a medium boudary) must satisfy the following law 1,13 p s fix (s; where p s fix (.) is always a pdf.Using Eq. ( 31) and (32), Eq. ( 36) can be written as where Note that the denominator of Eq. ( 37) is the mean photon step length s mix in the binary statistical mixture (if µ t a = µ t b we obtain s mix = 1 µ ta , as expected).The function F(s; a, b, N) can be expressed as (see the appendix for calculation details) where .
. is the binomial coefficient and 2 F 1 (., .; .; .) the hypergeometric function; and In MC simulations, the use of p s fix (.) is fundamental, because it allows one to preserve the invariance property for light propagation and the related reciprocity law. 1

Single step function: direct Monte Carlo simulation
Numerical generation of Eq. ( 32) Equation ( 32) has been validated numerically by MC simulation; i.e., by explicitly taking into account each single tessel crossed by the photons.Considering that by definition p s mix (s; µ t a , µ t b , σ a , σ b ) takes into account only "ballistic" photons (until they are absorbed or scattered), the code has been implemented -in MATLAB ® language -in the following manner: 1) A uniformly distributed random number ξ ∈ {0, 1} is generated; 2) If ξ < (1 − P) then the parameters used are µ t a and σ a , otherwise µ t b and σ b ; 3) A random tessel of length is generated by means of Eq. ( 2) and the chosen σ t i of point 2 (i ∈ {a, b}); 4) A random photon step s is generated by means of Eq. ( 1) and the chosen µ t i of point 2 (i ∈ {a, b}); 5) If s > l then permute the parameters (i.e., permute a and b parameters) and go to 3, otherwise stop the photon (because absorbed or scattered inside this last tessel) and store the total length traveled and the number of tessels crossed by the photon; 6) Go to 1 until the desired number of photons are propagated.

8/13
Note that, usually, in complex MC simulations many photons are launched for a given random tessel configuration.Then, this procedure is repeated a number of times and the "ensemble" average of the results is taken.This approach is imposed by the intensive computation demand of the MC code.However, the same results can be obtained by changing tessel configuration for each launched photon.This is what is done in the present simpler MC context (points 1 to 6 above).

Single step function for photons starting at any point inside a tessel
Equation (32) has been derived for photons starting at the medium boundary, and thus also from any tessel boundary (see, e.g., point A in Fig. 2).However, in real MC simulations photons propagate through the medium and steps may start at any point inside a tessel (see, e.g., point B in Fig. 2).In the present context this is not a problem, because Eq. ( 2) is a memoryless law (exponential pdf), and thus starting from a tessel boundary or inside the tessel does not change the present findings.
To give a more intuitive view on this point, some tutorial MC simulations have been performed.These simulations have the aim show that Eq. (32) remains valid even in the case where the photons start from points situated inside the tessels.(see, e.g., the s 2 step in Fig 2).Considered that the medium is isotropic, the MC code has been implemented in the following manner (see also Fig. 3): 1) A random s cte ∈ [0, +∞], representing the distance from any point on the medium boundary (e.g., point A in Fig. 3) along the considered direction of photon propagation, is chosen; 2) Random tessel lengths 1 , 2 , . . .are generated using Eq. ( 2), by alternatively using σ = σ a and σ = σ b , until for a given n, 1 + 2 + . . .n > s cte ; 3) = s cte − ( 1 + 2 + . . .n ) is computed; 4) Points 2 and 3 are repeated many times and the different are saved; 5) The pdf for is obtained by computing the histogram of the data saved in point 4.
In practice, we want to demonstrate that the pdf for n (see, Fig. 3) is the same as the pdf for ; i.e., we always have Eq.( 2); where n ∈ {1, 2, . . .} represents any n th tessel where the photon is extinct.This means that, it is not important if the photon starts from the tessel boundary or inside the tessel, because Eq. (32) remains the same (i.e., the previous mathematical derivation of Eq. (32) does not change).

Explanatory examples
In the following subsections we will show the validity of Eq. (32) through some intuitive explanatory examples, and comparisons with MC simulations.

Example 1
One of the simplest cases we can study is when one of the tessels type (e.g., "b") has mean length zero [Eq.(3)].This means that we remain with an homogeneous medium of type "a" (only sky blue tessels).In this case, we must retrieve the classical exponential law [Eq.(1)].Indeed, lim as it is expected.In the case of an isotropic medium, also σ a must go to +∞ because this quantity is linked to σ b by Eqs. ( 4) and ( 5), i.e.; L → 0. In this latter particular case p s mix (.)=0, because the mean length of the tessels is zero, and thus the propagating medium does not exist.

9/13
Example 2 Another simple case is when the mean lengths [Eq.( 3)] of both tessel types are infinite.In this case, we expect that the photon stops inside the first tessel at the entrance of the medium.Considering that tessel types "a" and "b" have probability P and (1 − P) to be the first tessel, we also expect to obtain the sum of two classical exponential laws [Eq.( 1)] weighted by their probability.In fact, lim as it must be.If the medium is isotropic [Eqs.( 4) and ( 5); L → +∞] the result remains the same.Due to the fact that Eq. ( 43) is not a pure exponential law, photon steps that start from a fixed position must be described using Eq.(36), i.e.; lim Example 3 In the case µ t a = µ t b and σ a = σ b [see, Eqs. ( 4) and ( 5)], we are in the presence of a homogeneous medium, and it is possible to obtain the explicit solution for p s mix (.), as p s mix (s; µ t a , µ t a , σ a , σ a , P) = µ t a e −µ ta s ; 0 where Eq. ( 45) represents the expected pdf for the homogeneous medium with extinction coefficient µ t a .The parameter σ a has disappeared because it is no more possible to distinguish the tessels, i.e.; they all have the same optical parameters.It obviously follows, that if the medium is isotropic -i.e.; for µ t a = µ t b and σ a = σ b = 1 2L [see, Eqs. ( 4) and ( 5)], when P = 1 2 -Eq.(45) remains valid.

Example 4
To investigate more complex cases, we need to compare our analytical model for the SSF [Eq.(32)] with the relative "gold standard" MC simulations.Figure 4 shows the SSF p s mix (s; .)for different optical and geometrical values, compared to the MC simulations.It can be seen that the analytical method gives the same results as the reference MC data.

Example 5
Figure 5 (left panels) shows that the pdf of (see, Fig. 3 for the symbols) is equal to the pdf of n [Eq.( 2)].Two representative cases for a small and a large s cte are reported.The data follow a straight lines due to the expected exponential behavior.For a given s cte , we also observe two different lines, depending whether we consider tessels of type "a" or type "b".Hence, MC data We were not able to perform the integral directly, instead we use Laplace transform techniques.We denote by L { f (s)} = +∞ 0 dte −st f (s) the Laplace transform of the function f (s) and by L −1 { f (t)} = c+i∞ c−i∞ dt 2iπ e st f (t) its inverse Laplace transform.For sake of simplicity, we only treat the case where N = 2p is even (the odd case can be handled in the same way).

Figure 3 .
Figure 3.A photon coming from B scatters and goes in direction C. The scattering point is at distance s cte from the boundary, measured along the scattered direction.

Figure 4 .
Figure 4.The SSF p s mix (.) as a function of s [Eq.(32)] for a medium with isotropic Poisson tesselation and binary mixture, and for a set different optical and geometrical parameters.Black lines represent the relative MC data.
1 (.); resulting in a gamma function with parameters N and µ t + σ [Eq.(17)].Equation (17) allows one to consider two cases for p Tes N (s; µ t a , µ t b , σ a , σ b , N); i.e., for N odd and N even.