Diffusion of light in structurally anisotropic media with uniaxial symmetry

Anisotropic light transport is extremely common among scattering materials, yet a comprehensive picture of how macroscopic diffusion is determined by microscopic tensor scattering coefficients is not fully established yet. In this work, we present a theoretical and experimental study of diffusion in structurally anisotropic media with uniaxially symmetric scattering coefficients. Exact analytical relations are derived in the case of index-matched turbid media, unveiling the general relation between microscopic scattering coefficients and the resulting macroscopic diffusion tensor along different directions. Excellent agreement is found against anisotropic Monte Carlo simulations up to high degrees of anisotropy, in contrast with previously proposed approaches. The obtained solutions are used to analyze experimental measurements of anisotropic light transport in polystyrene foam samples under different degrees of uniaxial compression, providing a practical example of their applicability.


I. INTRODUCTION
Anisotropic light transport is a ubiquitous feature of many scattering media, ranging from biological materials as wood [1], tendons [2][3][4], teeth [5], skin [6] and bone [7], to several types of plastics [8] and isotropic materials under mechanical deformations [9].Despite their relevance for many applications, however, anisotropic light transport properties are often overlooked or disregarded altogether, due to our incomplete understanding of how they arise from the structural features of a material.
In principle, for large scattering materials in the diffusive regime, a general Anisotropic Diffusion Equation (ADE) can be easily cast by replacing the scalar diffusion coefficient D with a diffusion tensor D containing the diffusive rates along different directions.In contrast with the isotropic case, however, the exact relation between the anisotropic diffusive coefficients and their corresponding microscopic scattering properties remains unknown, as well as the expression for the effective boundary conditions required to describe diffusion in finite media.
Despite the significance of structurally anisotropic materials and the need for accurate physical models to describe their light transport regime, a surprisingly small number of works exist which attempt to address this issue, with inconclusive results.Up to date, the agreement between solutions to the radiative transfer equation and anisotropic diffusion has been mostly qualitative [10], with large discrepancies between Monte Carlo (MC) simulations and ADE which seemed to imply a fundamental irreducible incompatibility between anisotropic disorder and the core assumptions of diffusion theory [11].Attempts have been made to reconcile the two approaches, by suggesting a different relationship between scattering cross-sections and the components of the diffusion tensor [9], which however did not resolve the discrepancies with numerical Monte Carlo simulations [12].In 2014, Alerstam showed that the observed deviations were due to the use of oversimplified assumptions on the interdependence of both the diffusion rates and the extrapolated boundary condition on the microscopic scattering coefficients [13], proposing a numerical recipe to compute all relevant transport quantities starting from a single anisotropic trajectory.When using the so obtained coefficients, excellent agreement was finally recovered between ADE and a general anisotropic Monte Carlo framework, albeit only up to a moderate degree of anisotropy.After this important milestone, only sporadic attempts were made towards the development of an analytical framework to connect microscopic and macroscopic transport parameters, without success [14].
In this work, we lay the foundation for an analytical description of anisotropic diffusion parameters based on microscopic scattering coefficients for materials characterized by uniaxial structural anisotropy.Closed-form relationships are derived in case of index-matched boundaries and no scattering asymmetry, which are found in excellent agreement with anisotropic Monte Carlo simulations up to very high degrees of anisotropy.Results are presented in the context of radiative transport, however, their derivation is general and could be seamlessly applied to other fields such as neutron transport in pebbled bed reactor cores [15] or particles in fluids [16], to name a few.
The paper is organized as follows: a theoretical background to anisotropic transport theory is provided in Section II, introducing the parameters of interest and the problem at hand.Section III presents the derivation of analytical solutions in the case of an index-matched medium with uniaxially symmetric scattering tensor, allowing to evaluate the diffusive tensor elements and the extrapolated length directly from the scattering tensor coefficients.These solutions are validated against an opensource Monte Carlo implementation of anisotropic light transport in Section IV, and compared to time-resolved experimental measurements of diffuse reflectance from insulating foam samples under different degrees of controlled mechanical compression in Section V. Finally, the results are commented in Section VI.

II. ANISOTROPIC DIFFUSION EQUATION
The three-dimensional time-dependent Anisotropic Diffusion Equation (ADE) for a scattering and absorbing medium can be written as: where Φ(r, t) is the fluence rate at a position r = (x, y, z) and time t, µ a is the absorption coefficient, S(r, t) is the source term, v = c/n is the speed of light in the medium and D is the anisotropic diffusion tensor having units of length.At the microscopic level, we express the scattering anisotropy via a direction-dependent scattering coefficient µ s (ŝ) , which may exhibit a complex dependence on the unit vector ŝ along the incoming direction of the energy packet before a scattering event.On the other hand, in the following we will assume that the direction dependence of the scattering coefficient can be expressed using three scalar components of a scattering coefficient tensor µ s , through the relation µ s (ŝ) = ŝµ s ŝT .Similar tensors could be defined for the refractive index n or the singlescattering asymmetry factor g. In fact, birefringence and direction-dependent scattering asymmetry may often be associated to the presence of structural anisotropy in tissues and fibrous materials.All these tensors can be considered diagonal (allowing us to use a single subscript index for each component, for brevity), provided that a suitable set of absolute anisotropy axes can be defined in the laboratory reference frame, which in our case we will consider to be aligned with the Cartesian axes.Even in this simplified scenario where each parameter can be defined by just three scalar components, however, the number of free parameters may still grow too large compared to the available experimental observables, making their full retrieval impractical.
In this work, for the sake of simplicity, we will assume that the transport anisotropy information can be condensed in the elements of the scattering coefficient tensor, assuming both n and g as direction-independent scalars, and further taking g = 0 to focus on the role of the scattering coefficient alone in shaping the anisotropic transport behavior.
Similarly, in analogy to the isotropic case, the scattering mean free paths ℓ i and transport mean free paths ℓ * i along i = {x, y, z} will be indicated as: where it should be stressed that, despite setting g = 0, in general ℓ i ̸ = ℓ * i due to the presence of scattering anisotropy.As a matter of fact, in contrast to the isotropic case, no separable "similarity relation" is expected to hold between the scattering and transport mean free paths, which will be mutually connected via a new adimensional anisotropy tensor with elements defined via the relation: In other words, the tensor K represents the transformation connecting the microscopic scattering properties to the macroscopic diffusion rates: In the isotropic case K becomes the 3 × 3 identity matrix.
The diagonal elements of the diffusion tensor are the coefficients appearing in the anisotropic diffusion equation.In this context, a commonly studied configuration is that of an infinitely extended slab with thickness L illuminated by a pencil beam pulse along the perpendicular z axis.In the diffusive regime, a pencil beam source can be more conveniently modeled as an isotropic point source placed at the depth z 0 = ℓ z from the surface S(r, t) = δ(x)δ(y)δ(z − z 0 )δ(t) (note that in the anisotropic case z 0 = ℓ z ̸ = ℓ * z , even when assuming g = 0).The solution of the time-dependent ADE can be written as: with and z e as the extrapolated length.In order to solve the ADE as a function of the microscopic scattering coefficient elements µ i , closed-form expressions for the diffusive rates D i and extrapolated boundary length z e must be therefore derived.

III. ANALYTICAL SOLUTIONS FOR UNIAXIALLY SYMMETRIC ANISOTROPY
A theoretical framework for the definition of a generalized relationship between scattering coefficients and diffusion rates has been presented recently, covering both non-classical propagation regimes and angle-dependent step length distributions [15].This approach provides integral expressions which can be used to describe anisotropic transport [17], even though explicit solutions to these expressions have not been reported nor validated against Monte Carlo simulations.In the following, we present analytical solutions which are derived in the simplified case of media with uniaxially symmetric scattering coefficient tensor, with isotropic phase function (g = 0) and index-matched boundary conditions with the environment.Under these conditions, simple closed-form expressions can be derived, which shed light on the analytical interdependence between microscopic and macroscopic (i.e., experimentally measurable) quantities.
Without loss of generality, let us define z and ϕ as the anisotropy symmetry axis and the angle of rotation around it, respectively.Consequently, we make a distinction between longitudinal (parallel to the symmetry axis) and azimuthal (perpendicular to the symmetry axis) parameters, renaming µ s,z = µ s, ∥ , µ s,x = µ s,y = µ s,⊥ , and D z = D∥, D x = D y = D ⊥ .Given a generic unit vector associated to the walker direction of propagation ŝ = (sin θ cos ϕ, sin θ sin ϕ, cos θ), with θ and ϕ as the polar and azimuthal angles of the absolute (laboratory) reference frame, we find as expected that the directiondependent scattering coefficient µ s (ŝ) is independent of the azimuthal angle ϕ: (9) where χ = cos θ.
A first relevant quantity that can be derived is the expectation value of the mean free path ⟨ℓ⟩ where ξ(ŝ) represents the probability that a walker travels within a solid angle dŝ about ŝ after a scattering event, and ℓ(ŝ) = 1/µ s (χ) is the direction-dependent scattering mean free path along the same direction [17].For an isotropic phase function, we get that ξ(ŝ) = 1/4π, so that the only source of anisotropy is given by ℓ(ŝ) and the integral can be simplified to: The diffusion coefficients, on the other hand, are defined as [17]: which also admit closed-form solutions: From these analytical equations, it is evident that the macroscopic diffusion rate along a certain direction D i does not depend only on the scattering rates relative to the same direction µ s,i .Based on equations ( 7) and ( 11), the coefficients can be rewritten in terms of the scattering mean free paths and the corresponding anisotropy tensor elements: with: Notably, these equation reveal a conservation property for the trace of the anisotropy tensor: which holds independently of the degree of anisotropy between the longitudinal and perpendicular directions, and is verified to hold numerically also for full 3D anisotropy (i.e., when µ s,x ̸ = µ s,y ̸ = µ s,z ).The analytical solution for the extrapolated length z e is derived based on a general approach which accommodates for an anisotropic radiance [13].Taking advantage of the absence of reflections at boundaries for an index-matched scattering medium, the only contributions determining the extrapolated length are given by the integrals z e = C/B [13]: where P (ŝ) describes the angular distribution of the radiance, i.e., the probability density function that a random walker is traveling along a direction ŝ after many scattering events.In the general case, P (ŝ) depends both on the scattering rate tensor and the asymmetry factor g. In our case, we consider g = 0 so that to get P (ŝ) = ℓ(ŝ)/ ⟨ℓ⟩, which further simplifies for uniaxial symmetry to P (χ) = 1/µ s (χ) ⟨ℓ⟩.We note that our equation (22) deviates slightly from that reported by Alerstam due to its dependence on ℓ(ŝ) rather than ℓ * (ŝ) [13].This change is motivated by the fact that we do not expect a simple similarity relation to hold in the case of anisotropic light transport.Moreover, irrespective of the value of the scattering asymmetry (which we have assumed here to be g = 0), a general expression for the extrapolated length should depend on the microscopic scattering parameters rather than on the transport mean free path which, in the context of anisotropic propagation, is instead a direct expression of the macroscopic diffusion rates.A validation of this correction against Monte Carlo simulations is discussed in Section IV B. The previous integrals simplify to: to give an analytical expression for z e ) Equations ( 11), ( 18), ( 19) and ( 25) represent the main theoretical results introduced in this work.As an important consistency check, it can be easily verified that all derived analytical equations converge to the familiar isotropic limit when µ s,⊥ and µ s, ∥ tend to the same value µ s = 1/ℓ.

IV. MONTE CARLO VALIDATION
We validate our analytic results using the open source Monte Carlo package PyXOpto [18], which has been modified to handle anisotropic scattering and asymmetry factors.The analytical ADE solution in the time and space domain is also provided as a set of MATLAB functions in the supplemental material.

A. Transverse intensity spread
Time-resolved reflected and transmitted intensity profiles for a slab illuminated along the z axis can be obtained from the fluence rate in eq. ( 8) using Fick's law [17,19] by setting R(x, y, t) and T (x, y, t) are bi-variate Gaussian distributions: with mean square displacements (MSD) along x and y and integrated time-domain amplitudes R 0 (t) and T 0 (t).
Accessing either R(x, y, t) or T (x, y, t) at different times allows to retrieve directly the diffusion rates (or, equivalently, the transport mean free paths ℓ * x , ℓ * y ) along x and y.Moreover, the evolution of w 2 x (t) and w 2 y (t) does not depend on the amplitude of the profile, so that any effect that may modify its overall intensity (such as a homogeneous absorption coefficient) factors out exactly from the MSD.
We performed time-and space-resolved Monte Carlo simulations for an illustrative slab geometry of thickness L = 1 mm, setting a uniaxial scattering anisotropy along the y axis (µ s,x = µ s,z = µ s,⊥ , µ s,y = µ s, ∥ ).We set ℓ x = ℓ z = 10 µm, and vary ℓ y between 10 µm to 100 µm with 10 µm steps.The macroscopic diffusion coefficients are retrieved by performing a linear regression on the growth rate of the MSD after an initial transient.The resulting values of ℓ * x and ℓ * y are shown in Figure 1 together with the analytical predictions obtained from equations ( 14) and (15).The deviation between theory and MC across this parameter range is about 0.8 % along the longitudinal axis and 0.2 % in the azimuthal plane, which is representative of the approximations of diffusion theory.This discrepancy grows up to a maximum of 1.4 % for the largest value of ℓ y = 100 µm, corresponding to an extreme anisotropy ratio of ℓ∥/ℓ ⊥ = 10.This is also expected in the context of diffusion theory, since a longer scattering mean free path along one axis will also increase the average scattering mean free path (11), hence reducing the "effective" optical thickness of the slab.Numerical results for pencil beam illumination are compared to the analytical solutions using different definitions for the extrapolated length z e and for the equivalent depth z 0 of the isotropic point source.

B. Time-resolved and steady-state profiles
Steady-state and total time-domain transmittance can be derived by direct integration of equation ( 27) as and with analogue expressions for the reflectance.
A first validation test regards the definition of the extrapolated length z e and the depth z 0 of the isotropic point source required to solve the ADE, which is such to match the case of a pencil beam illumination.In the current literature [13], it is suggested to compute z e as a function of ℓ * (ŝ) and to use z 0 = ℓ * z .However, when testing these assumptions against Monte Carlo simulations of anisotropic media with significant transport anisotropy, a much better agreement is obtained when defining z e as a function of ℓ(ŝ) and setting z 0 = ℓ z (Figure 2).
A final set of MC simulations for the time-and spaceresolved transmittance T (x, y, t) was performed for slabs with g = 0, n = 1, µ a = 0 and different scattering coefficients.The results for one representative case are shown in Figure 3, compared to the analytical predictions.Excellent agreement is found at all times and positions, with average absolute discrepancies below 1 % in the considered spatio-temporal range.The prediction from anisotropic theory (dashed) is compared against the results of a Monte Carlo simulation comprising 10 13 trajectories (solid).

V. EXPERIMENTAL MEASUREMENTS
We finally test anisotropic light transport experimentally in samples of extruded polystyrene (XPS) foam under controlled degrees of uniaxial compression, and use the newly developed ADE to retrieve its microscopic scattering coefficients.Due to its very high porosity, the effective refractive index of this material is close to that of air.Moreover, the absence of a well-defined boundary interface further suppresses the occurrence of internal reflections at boundaries, which conforms well with the assumptions under which our analytical results have been derived.
In pristine conditions, XPS foam is characterized by an isotropic closed-cell micro-structure and isotropic scattering properties.Under uniaxial compression, however, the foam cells become increasingly oblate and eventually collapse, introducing different degrees of uniaxial scattering anisotropy [20][21][22].
We realized different samples with increasing levels of compressive strain (defined as the ratio between the height contraction ∆h and the original height h 0 ), spanning across the elastic-plastic deformation transition of XPS (Figure 4) Anisotropic light transport was studied by means of an optical gating imaging apparatus allowing to capture the evolution of intensity profiles with sub-ps resolution [23,24].A collection of frames recorded for each sample at a probe wavelength of 820 nm is shown in Figure 5 for t ≤ 15 ps.Within this time limited range, all samples can be effectively considered as semi-infinite.
Due to the uniaxial compression condition, each sample becomes characterized by a pair of azimuthal (ℓ x = ℓ z ) and longitudinal (ℓ y ) scattering mean free paths.We further assume g = 0 and no Fresnel reflections at boundaries, given the large pore sizes of the samples.Nonetheless, an effective refractive index is calculated for each sample based on the increasing solid fraction under compression, ranging from n = 1.014 for the uncompressed case (solid fraction 2.8 %) to n = 1.043 for the most compressed sample (solid fraction 14 %).These values are used for the correct determination of the energy velocity v inside the sample.
The reflected intensity profiles were fitted at different times with bi-variate Gaussian distributions to retrieve their MSD evolution along the y (vertical) and x (horizontal) axes.The resulting diffusive constants D x , D y , expressed through their corresponding transport mean free paths, can be then obtained via a linear regression of the MSD growth rates.The microscopic scattering coefficients ℓ x and ℓ y along the two anisotropic directions can be finally retrieved using the analytical relations ( 15) and ( 14).The results are shown in Figure 6a.
As intuitively expected, we find that ℓ y is progressively reduced with increasing compression.On the other hand, ℓ x = ℓ z remains mostly constant for small degrees of compression, and then starts to decrease with increasing compressive strain (Fig. 6a).The overall variation is well captured by the anisotropy tensor elements introduced in equation ( 6), which exhibit opposite trends as a function of compressive strain.As suggested by the relation Tr(K)/3 = 1, the increase of K y is compensated by a corresponding decrease of K x (Fig. 6b).This analysis allows to quantitatively estimate the magnitude of the systematic error introduced when inferring microscopic scattering parameters directly from macroscopic observables in structurally anisotropic media, which may directly affects several results previously reported in the literature [25][26][27][28].

VI. CONCLUSIONS
In this work, we analyzed anisotropic light transport induced by a uniaxial tensor scattering coefficient, and derived analytic expressions for both the diffusion rates and the extrapolated boundary length.In the simple The availability of an analytic solution to the anisotropic diffusion equation can streamline the analysis of experimental data, which we demonstrate in the case of a scattering medium subject to a controlled degree of anisotropic deformation.Additionally, it provides a direct insight into the functional interdependence between the scattering properties along different axes, which may be extended to more general anisotropic configurations.More importantly, our study revealed the presence of systematic discrepancies between the observed diffusive rates and the corresponding scattering coefficients, calling for a careful re-evaluation of previously reported scattering mean free path values of anisotropic media.At the same time, it facilitates the application of a more correct diffusive model in several cases where structural anisotropy is evident and yet disregarded due to the lack of established numerical or theoretical approaches.
The analytical solution that we have presented is applicable only to the specific case of index-matched media with an isotropic phase function.This represents a configuration of prominent relevance, especially considering the fact that in the diffusive regime we do not expect to be able to uniquely determine the (tensor) scattering asymmetry from the diffusive rates.Indeed, experimental studies on anisotropic light transport typically consider only one tensor for the scattering coefficients, making our approach directly applicable to the analysis of a large array of previously reported results.Additionally, as we have shown, there is potentially a large class of materials characterized by suppressed reflection at boundaries and uniaxial anisotropy, such as foams or aerogels which are permeated by their host medium, but also biological tissues which are immersed in an aqueous environment decreasing their index contrast with the surrounding media.
For more general anisotropic media, with different scattering and asymmetry coefficients along each axis, we anticipate that closed-form expressions may very likely be impractical to use or not exist altogether.However, even in this case, our results confirm the validity of the adopted theoretical framework for the calculation of the relevant transport parameters, which may still be solved numerically for more complex configurations at a fraction of the computational cost required to run a full anisotropic Monte Carlo simulation, driving the further development of new solutions for the accurate description of direction-dependent diffusion processes.

FIG. 1 :
FIG.1: Mean free paths along x and y obtained with MC simulations and the analytical formulas as a function of ℓ y .Dashed lines represent the simplistic assumption ℓ * i = ℓ i .

TFIG. 2 :
FIG. 2: Total time-resolved transmittance through a 1 mm-thick slab with ℓ x = 30 µm and ℓ y = ℓ z = 10 µm.Numerical results for pencil beam illumination are compared to the analytical solutions using different definitions for the extrapolated length z e and for the equivalent depth z 0 of the isotropic point source.

1 ( 1 (FIG. 3 :
FIG.3:(a) Steady-state transmittance T (x, y) from a 1 mm-thick slab with ℓ x = 30 µm and ℓ y = ℓ z = 10 µm, plotted as log-scale contour lines marking consecutive decades.(b, c, d) Time-resolved transmittance at different locations (x 0 , y 0 ) on the slab's output surface, identified by the markers in panel (a).The prediction from anisotropic theory (dashed) is compared against the results of a Monte Carlo simulation comprising 10 13 trajectories (solid).

FIG. 4 :FIG. 5 :
FIG. 4: Experimental XPS samples under different compressive strain.The uncompressed height of each foam sample is h 0 = 29 mm.A constant level of compression is maintained by placing the samples inside rigid frames.Scanning electron microscopy images of a cross-sectional slice of samples are shown in the bottom row for a few examples with increasing degree of compression.Scale bars correspond to 100 µm.

3 FIG. 6 :
FIG. 6: (a) Mean free path along x, y and (b) anisotropy tensor components as a function of the compressive strain for the 5 samples of insulation foam.

FUNDING
This work was partially funded by the European European Union's NextGenerationEU Programme with the I-PHOQS Research Infrastructure [IR0000016, ID D2B8D520, CUP B53C22001750006] "Integrated infrastructure initiative in Photonic and Quantum Sciences".E.P acknowledges financial support from Sony Europe B.V.. L.P. further acknowledges the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support (ISCRA-C "ARTTESC"), and NVIDIA Corporation for the donation of the Titan X Pascal GPU.F.M. acknowledges financial support by the European Union's -NextGenera-tionEU, National Recovery and Resilience Plan, MNESYS, PE0000006 (DN 572 1553 11.10.2022)and PRIN 2022, DIRS, grant 578 number: 2022EB4B7E.