One of the main determinants of the accuracy of a Monte Carlo simulation at small length scales is the choice of the phase function. Although the Henyey-Greenstein phase function has been described as a sufficiently accurate choice for prediction of backscattering at intermediate length scales (

*r~ls**), accurate modeling of small length scales (

*r*<<

*ls**) requires a more general choice of phase function. For this purpose, we will follow a recently developed model that is based on the Whittle-Matérn correlation function [

12,

13].

The model implements the Born approximation in order to obtain the phase function from the refractive index correlation function. The Born (i.e. weakly scattering) approximation is valid in the regime relevant for soft biological tissue [

14]. In the Born approximation, the differential cross section and, thus, the phase function, are completely defined through the Fourier transform of the refractive index correlation function [

15]. In turn, the refractive index is a linear function of the local density of tissue macromolecules and the refractive index correlation function is proportional to that of the mass density. Thus, one can calculate the phase function if the mass density correlation is known (and vice versa). There have been several hypotheses on the best functional form that can model the refractive index correlation function in tissue, but one convenient expression can encompass nearly all of these possibilities. The refractive index fluctuations of biological tissue can be modeled with the Whittle-Matérn correlation function [

16]:

, where

*Δn*^{2} is the variance of the refractive index fluctuations,

*l*_{c} is the correlation length, and

*m* is a parameter that determines the form of the function. The function

*K*_{m-3/2} denotes the modified Bessel function of the second kind of order

*m-3/2*. When

*m* < 1.5,

*B*_{n}(r) is a power law, thus corresponding to a mass fractal medium with mass fractal dimension

*D*_{mf} = 2m. 1.5<

*m* <2 corresponds to a stretched exponential function,

*m* = 2 corresponds to an exponential function, and as

*m* becomes much larger than 2,

*B*_{n}(r) approaches a Gaussian function. The correlation length

*l*_{c} has different physical meaning depending on the type of the correlation function. For

*m* = 2,

while in a mass fractal case of

*m* < 1.5,

*l*_{c} represents the upper length scale at which the correlation function loses its fractal behavior. The differential scattering cross section can be derived by applying the Born approximation to the Whittle-Matérn correlation function [

13]. For the scalar wave case, the phase function assumes the following form:

where

*k = 2π/λ* and the normalization is such that

. We will refer to this phase function as the Whittle-Matérn phase function. The phase function can also be expressed as a function of

*m* and

*klc* without any change to the normalization:

It is important to note that a special case is observed when the value of

*m* = 1.5. In this special case the correlation function is that of the space filling random field and the phase function simplifies to the commonly used form known as the Henyey-Greenstein phase function. The parameter

*ĝ* then becomes the average cos(θ), also known as the anisotropy factor

*g*. For other values of

*m*,

*g* is given by taking the forward moment:

shows examples of the Whittle-Matérn phase functions with the same value of

*m*, and varying values of

*g*, while shows examples of phase functions with the same value of

*g* and varying values of

*m*. The parameter

*g* influences the width of the phase function while

*m* influences the shape of the phase function independently of the width. There are two cases in the generalized phase function which are removable discontinuities:

*m* = 1 and

*g* = 0. We can evaluate the phase function for these cases by employing L’Hospital’s rule:

In these two cases, *g* becomes:

The reflectance distribution was calculated with the Whittle-Matérn phase function by implementing the Monte Carlo method. An existing Monte Carlo code was modified to implement the generalized phase function [

9,

17]. The code was validated by comparing the results for

*m* = 1.5 (Henyey-Greenstein case) with existing codes that implement the Henyey-Greenstein phase function [

10]. Simulations with values of

*g* varying from 0 to 0.98 and values of

*m* varying from 1.01 to 1.9 were obtained for the backscattering direction (0-10°). We found that the variations of backscattering probability distributions were small within the 10° angular collection range when the backscattering probability distribution was stored as a function of the position of the final scattering event. Therefore, all reflectance distributions were stored as a function of the position of the last scattering event.

*ls** was maintained at 100μm with a scattering slab thickness of 1cm, resulting in a scattering medium that approaches semi-infinite, with less than 2% of the intensity transmitting through the entire thickness of the slab. The infinitely narrow illumination beam was oriented orthogonally to the scattering medium. The boundary at the interface of the scattering medium was assumed to be index-matched and absorption was not present. The scattering angle,

*θ*, in the Monte Carlo simulation was chosen by expressing the probability of a selected angle as a function of the random variable

*ξ*:

where

*ξ* is uniformly distributed between 0 and 1 [

10]. The azimuthal angle was chosen from a uniform random distribution:

*ψ* = 2π

*ξ*. The remaining Monte Carlo simulation elements were identical to previously developed methodology for light propagation in turbid media [

10].