We do this for every $S_{i},i\in\\{1,2,\dots,n\\}$, so that we now have

${\rm
Corr}_{q}=\sum_{\lambda\in\Lambda}\left(\sum_{i=1}^{n}q_{i}\zeta(M_{i},\lambda)\right)\nu(\lambda).$
(53)

Note that in all this, the response functions satisfy the assumption of
measurement noncontextuality.

Let us denote by $\Lambda_{\rm detP}$ the following set of
$\lambda\in\Lambda$:

$\Lambda_{\rm
detP}\equiv\\{\lambda\in\Lambda\big{|}\mu(s_{i}|S_{i},\lambda)=1\textrm{ for
some }[s_{i}|S_{i}]\\},$ (54)

where $s_{i}\in[d],i\in[n]$ for any $[s_{i}|S_{i}]$.

Next, we choose $\nu(\lambda)$ in such a way that the inequality is saturated:
in particular, we take $\nu(\lambda)$ to be supported entirely over the set

$\Lambda_{\max}\equiv\\{\lambda\in\Lambda\big{|}\sum_{i=1}^{n}q_{i}\zeta(M_{i},\lambda)=\beta(\Gamma,q)\\},$

so that $\nu(\lambda)=0$ for all $\lambda\in\Lambda\backslash\Lambda_{\max}$.
Recalling the constraint of Eq. (52), that

$\displaystyle\forall[s_{i}|S_{i}]:\sum_{\lambda\in\Lambda}\mu(s_{i}|S_{i},\lambda)\nu(\lambda)=p(s_{i}|S_{i}),$
(55)

we have that our choice of $\mu(s_{i}|S_{i},\lambda)$ and $\nu(\lambda)$ above
should satisfy the property that

$\Lambda_{\rm detP}\cap\Lambda_{\max}\neq\varnothing,$ (56)

since otherwise we cannot satisfy the constraint of Eq. (52). If there exist
no $\lambda$ that lie in $\Lambda_{\rm detP}\cap\Lambda_{\max}$, then we
cannot proceed with our construction of the noncontextual ontological model.
If Eq. (56) is satisfied, then our choices of $\mu(s_{i}|S_{i},\lambda)$ and
$\nu(\lambda)$ must further satisfy the constraint of Eq. (52). Again, if this
cannot be done, our construction cannot proceed.

Assuming we have satisfied the constraint of Eq. (52) while choosing
$\mu(s_{i}|S_{i},\lambda)$ and $\nu(\lambda)$ as prescribed above, we then
have

${\rm Corr}_{q}=\beta(\Gamma,q)$ (57)

and the noise-robust noncontextuality inequality is thus saturated.232323What
we haven’t shown above is a concrete example of a scenario where a
noncontextual ontological model following our prescription can actually be
constructed and, to that extent, we haven’t contributed much to proving
saturation. We have merely suggested a possible route towards it.

## References

  * [1] S. Kochen and E. P. Specker, “The Problem of Hidden Variables in Quantum Mechanics”, https://doi.org/10.1512/iumj.1968.17.17004J. Math. Mech. 17, 59 (1967). Also available at https://www.jstor.org/stable/24902153JSTOR.
  * [2] D. A. Meyer, “Finite Precision Measurement Nullifies the Kochen-Specker Theorem”, https://doi.org/10.1103/PhysRevLett.83.3751Phys. Rev. Lett. 83, 3751 (1999).
  * [3] A. Kent, “Noncontextual Hidden Variables and Physical Measurements”, https://doi.org/10.1103/PhysRevLett.83.3755Phys. Rev. Lett. 83, 3755 (1999).
  * [4] R. Clifton and A. Kent, “Simulating quantum mechanics by non-contextual hidden variables”, https://doi.org/10.1098/rspa.2000.0604Proc. R. Soc. Lond. A: Vol. 456, 2101-2114 (2000).
  * [5] J. Barrett and A. Kent, “Non-contextuality, finite precision measurement and the Kochen-Specker theorem”, https://doi.org/10.1016/j.shpsb.2003.10.003Stud. Hist. Philos. Mod. Phys. 35, 151 (2004).
  * [6] R. Kunjwal and R. W. Spekkens, “From the Kochen-Specker Theorem to Noncontextuality Inequalities without Assuming Determinism”, https://doi.org/10.1103/PhysRevLett.115.110403 Phys. Rev. Lett. 115, 110403 (2015).
  * [7] M. D. Mazurek, M. F. Pusey, R. Kunjwal, K. J. Resch, R. W. Spekkens, “An experimental test of noncontextuality without unphysical idealizations”, https://doi.org/10.1038/ncomms11780Nat. Commun. 7, 11780 (2016).
  * [8] M. F. Pusey, “Robust preparation noncontextuality inequalities in the simplest scenario”, https://doi.org/10.1103/PhysRevA.98.022112Phys. Rev. A 98, 022112 (2018).
  * [9] A. Krishna, R. W. Spekkens, and E. Wolfe, “Deriving robust noncontextuality inequalities from algebraic proofs of the Kochen-Specker theorem: the Peres-Mermin square”, https://doi.org/10.1088/1367-2630/aa9168New J. Phys 19, 123031 (2017).
  * [10] R. Kunjwal and R. W. Spekkens, “From statistical proofs of the Kochen-Specker theorem to noise-robust noncontextuality inequalities”, https://doi.org/10.1103/PhysRevA.97.052110Phys. Rev. A 97, 052110 (2018).
  * [11] D. Schmid, R. W. Spekkens, and E. Wolfe, “All the noncontextuality inequalities for arbitrary prepare-and-measure experiments with respect to any fixed set of operational equivalences”, https://doi.org/10.1103/PhysRevA.97.062103Phys. Rev. A 97, 062103 (2018).
  * [12] R. Kunjwal, “Contextuality beyond the Kochen-Specker theorem”, https://arxiv.org/abs/1612.07250arXiv:1612.07250 [quant-ph] (2016).
  * [13] R. W. Spekkens, “Contextuality for preparations, transformations, and unsharp measurements”, https://doi.org/10.1103/PhysRevA.71.052108Phys. Rev. A 71, 052108 (2005).
  * [14] R. Kunjwal, “Beyond the Cabello-Severini-Winter framework: Making sense of contextuality without sharpness of measurements”, https://doi.org/10.22331/q-2019-09-09-184Quantum 3, 184 (2019).
  * [15] R. W. Spekkens, “The Status of Determinism in Proofs of the Impossibility of a Noncontextual Model of Quantum Theory”, https://doi.org/10.1007/s10701-014-9833-xFound. Phys. 44, 1125-1155 (2014).
  * [16] R. Kunjwal, “Fine’s theorem, noncontextuality, and correlations in Specker’s scenario”, https://doi.org/10.1103/PhysRevA.91.022108Phys. Rev. A 91, 022108 (2015).
  * [17] A. Fine, “Hidden Variables, Joint Probability, and the Bell Inequalities”, https://doi.org/10.1103/PhysRevLett.48.291Phys. Rev. Lett. 48, 291 (1982).
  * [18] R. W. Spekkens, “Leibniz’s principle of the identity of indiscernibles as a foundational principle for quantum theory”, Information-Theoretic Interpretations of Quantum Mechanics (2016), University of Western Ontario, available at https://www.youtube.com/watch?v=HWOkjisIxc4.
  * [19] R. W. Spekkens, “The ontological identity of empirical indiscernibles: Leibniz’s methodological principle and its significance in the work of Einstein”, https://arxiv.org/abs/1909.04628arXiv:1909.04628 [physics.hist-ph] (2019).
  * [20] A. Cabello, Adan, J. Estebaranz, and G. Garcia-Alcaine, “Bell-Kochen-Specker theorem: A proof with 18 vectors,” https://doi.org/10.1016/0375-9601(96)00134-XPhys. Lett. A 212, 183 (1996).
  * [21] A. Peres, “Two simple proofs of the Kochen-Specker theorem”, https://doi.org/10.1088/0305-4470/24/4/003J. Phys. A 24, L175 (1991).
  * [22] N. D. Mermin, “Hidden variables and the two theorems of John Bell”, https://doi.org/10.1103/RevModPhys.65.803Rev. Mod. Phys. 65, 803 (1993).
  * [23] A. A. Klyachko, M. A. Can, S. Binicioğlu, and A. S. Shumovsky, “Simple Test for Hidden Variables in Spin-1 Systems”, https://doi.org/10.1103/PhysRevLett.101.020403Phys. Rev. Lett. 101, 020403 (2008).
  * [24] A. Cabello, S. Severini, and A. Winter, “Graph-Theoretic Approach to Quantum Correlations”, https://doi.org/10.1103/PhysRevLett.112.040401Phys. Rev. Lett. 112, 040401 (2014).
  * [25] A. Acín, T. Fritz, A. Leverrier, and A. B. Sainz, A Combinatorial Approach to Nonlocality and Contextuality, https://doi.org/10.1007/s00220-014-2260-1Comm. Math. Phys. 334(2), 533-628 (2015).
  * [26] S. Abramsky and A. Brandenburger, “The sheaf-theoretic structure of non-locality and contextuality”, https://doi.org/10.1088/1367-2630/13/11/113036New J. Phys. 13, 113036 (2011).
  * [27] S. Abramsky, R. S. Barbosa, K. Kishida, R. Lal, and S. Mansfield, “Possibilities determine the combinatorial structure of probability polytopes”, https://doi.org/10.1016/j.jmp.2016.03.006Journal of Mathematical Psychology, Volume 74, Pages 58-65 (2016).
  * [28] S. Mansfield and R. S. Barbosa, “Extendability in the sheaf-theoretic approach: Construction of Bell models from Kochen-Specker models”, https://arxiv.org/abs/1402.4827arXiv:1402.4827 [quant-ph] (2014).
  * [29] R. Kunjwal, C. Heunen, and T. Fritz, “Quantum realization of arbitrary joint measurability structures”, https://doi.org/10.1103/PhysRevA.89.052126Phys. Rev. A 89, 052126 (2014).
  * [30] T. Gonda, R. Kunjwal, D. Schmid, E. Wolfe, and A. B. Sainz, “Almost Quantum Correlations are Inconsistent with Specker’s Principle”, https://doi.org/10.22331/q-2018-08-27-87Quantum 2, 87 (2018).
  * [31] A. Cabello, “Twin inequality for fully contextual quantum correlations”, https://doi.org/10.1103/PhysRevA.87.010104Phys. Rev. A 87, 010104(R) (2013).
  * [32] P. Lisoněk, P. Badziag, J. R. Portillo, and A. Cabello, “Kochen-Specker set with seven contexts”, https://doi.org/10.1103/PhysRevA.89.042101Phys. Rev. A 89, 042101 (2014).
  * [33] M. Pavicic, J-P. Merlet, B. McKay, and N. D. Megill, “Kochen-Specker vectors”, https://doi.org/10.1088/0305-4470/38/7/013 J. Phys. A: Math. Gen. 38 1577 (2005).
  * [34] S. Mansfield, “The Mathematical Structure of Non-locality and Contextuality”, https://ora.ox.ac.uk/objects/uuid:394bb375-db3f-4a12-bdd8-cd1ab5809573PhD thesis, University of Oxford (2013).
  * [35] G. Carù, “Towards a complete cohomology invariant for non-locality and contextuality”, https://arxiv.org/abs/1807.04203arXiv:1807.04203 [quant-ph] (2018).
  * [36] A. Krishna, “Experimentally Testable Noncontextuality Inequalities Via Fourier-Motzkin Elimination”, http://hdl.handle.net/10012/9752M.Sc. Thesis, University of Waterloo (2015).
  * [37] D. Schmid and R. W. Spekkens, “Contextual Advantage for State Discrimination”, https://doi.org/10.1103/PhysRevX.8.011015Phys. Rev. X 8, 011015 (2018).


# Evidence of Early-stage Tidal Structures of Open Clusters Revealed by
Kinematics with Gaia EDR3

Yezhang Li Department of Physics, Xi’an Jiaotong-Liverpool University, 111
Ren’ai Road, Dushu Lake Science and Education Innovation District, Suzhou
215123, Jiangsu Province, P. R. China<EMAIL_ADDRESS> Xiaoying
Pang Department of Physics, Xi’an Jiaotong-Liverpool University, 111 Ren’ai
Road, Dushu Lake Science and Education Innovation District, Suzhou 215123,
Jiangsu Province, P. R. China<EMAIL_ADDRESS> Shanghai Key
Laboratory for Astrophysics, Shanghai Normal University, 100 Guilin Road,
Shanghai 200234, P. R. China  Shih-Yun Tang Lowell Observatory, 1400 W. Mars
Hill Road, Flagstaff, AZ 86001, USA  Department of Astronomy and Planetary
Sciences, Northern Arizona University, Flagstaff, AZ 86011, USA

###### Abstract

Blanco 1, a 100 Myr open cluster in the solar neighborhood, is well known for
its two 50 pc-long tidal tails. Taking Blanco 1 as a reference, we find
evidence of early-stage tidal disruption in two other open clusters of
$\sim$120 Myr: the Pleiades and NGC 2516, via Gaia EDR3 data. These two
clusters have a total mass of $2-6$ times that of Blanco 1. Despite having a
similar age as Blanco 1, the Pleiades and NGC 2516 have a larger fraction of
their members bound: 86% of their mass is inside the tidal radius, versus 63%
for Blanco 1. However, a correlation between Blanco 1’s 50 pc-long tidal tails
and the “kinematic tails” in velocity space is also found for the Pleiades and
NGC 2516. This evidence supports the idea that the modest elongation seen in
the spatial distribution for the Pleiades and NGC 2516 is a result of early-
stage tidal disruption.

##  1 Motivation

Figure 1:  The spatial (a–b) and the tangential velocity (c) distributions of
the Pleiades, NGC 2516 and Blanco 1. The green ellipses in panels (a)
represent the projection of the fitted ellipsoids onto the $X$-$Z$ plane. The
three ellipsoid axes ($a$, $b$, and $c$) are indicated as red, pink, and
orange. The cluster centers are marked with an ($+$) in panels (b). Dashed
black circles indicate the tidal radii. Grey circles in panels (b) and (c)
represent stars located within tidal radius. $r_{\rm t}$ is the tidal radius,
$b_{\rm f}$ the bound mass fraction, and $M_{\rm cl}$ the total mass. The
($+$) symbol in panels (c) indicates the median value of the tangential
velocity $V_{t_{\alpha}}$ and $V_{t_{\delta}}$, obtained from proper motions.
The red and blue open circles in panels (b) and (c) respectively represent
stars outside the tidal radius with $V_{t_{\delta}}$ values larger, or less
than, the median $V_{t_{\delta}}$. The velocity dispersion is indicated in the
bottom-left corner of panels (c).

Spatially-elongated structures outside open clusters’ tidal radius have been
found in local samples of all ages (Pang et al., 2021). However, extended
tails $\gtrapprox$50 pc are only discovered in star clusters older than
$\sim$600 Myr in the solar neighborhood (e.g., Coma Berenice, Hyades,
Praesepe: Tang et al., 2019; Röser et al., 2019; Röser & Schilbach, 2019),
except for Blanco 1 with 100 Myr (Zhang et al., 2020; Pang et al., 2021). Two
other clusters, Pleiades and NGC 2516 ($\sim$120 Myr), with similar ages as
Blanco 1, have structures outside the tidal radius (Lodieu et al., 2019; Pang
et al., 2021). Whether the small number of stars outside the tidal radius in
the Pleiades and NGC 2516 provides evidence of early-stage tidal tail
formation remains unclear. We therefore investigate the morphology of the
Pleiades, NGC 2516, and Blanco 1, to search for the kinematic evidence
supporting the nature of the tidal structures in clusters without significant
elongation.

##  2 Data analysis

Cluster members used in this study are obtained from Pang et al. (2021, NGC
2516 with 2690 members and Blanco 1 with 703 members), except for the
Pleiades, which is not included in Pang et al. (2021). We follow the same
procedure in Pang et al. (2021) to identify members of the Pleiades via the
machine learning algorithm StarGO (Yuan et al., 2018) with the latest Gaia
Early Data Release 3 (EDR 3; Gaia Collaboration et al., 2021). We find 1407
member stars for the Pleiades, with a mean position of
$(X,Y,Z)=(-121.2,+29.1,-54.2)$ pc (in Heliocentric Cartesian coordinates), and
mean proper motions of
$(\mu_{\alpha}\cos\delta,\mu_{\delta})=(+19.908,-45.341)$ mas yr-1. The member
list for the Pleiades is provided
online111https://doi.org/10.5281/zenodo.5067751 (catalog Zenodo:
10.5281/zenodo.5067751). Our identification of members is purely based on 5D
phase-space information, without pre-assumptions about age, cluster size, and
density threshold.

##  3 Tidal Structures versus Kinematic “Tails”

We calculate the Heliocentric ($X,Y,Z$) coordinates for the Pleiades, NGC
2516, and Blanco 1 using distances that have been corrected from the parallax
errors (see Section 4.1 in Pang et al., 2021). Considering a distance range of
135.9–410.5 pc for these clusters, the uncertainty in the corrected distance
is about 0.8–2.6 pc, according to the simulations of Pang et al. (2021, their
Figure 4).

To quantify the morphology of the clusters, we perform an ellipsoid fitting to
their 3D spatial distribution (see Section 4.3 in Pang et al., 2021). Figure 1
panels (1.a–3.a) present the best-fitted ellipsoids (green ellipses) projected
onto the $X$-$Z$ plane. The semi-major axis “$a$” represents the direction of
morphological elongation, and the axes ratios $b/a$ (semi-intermediate/semi-
major) and $c/a$ (semi-minor/semi-major) indicate the significance of the
elongation. Although both the Pleiades and NGC 2516 are elongated, their
elongation is less significant than that of Blanco 1. All three clusters tend
to stretch along the direction of the Galactic plane. Considering the age of
$\sim$100 Myr, such changes in shape are most likely caused by the external
tidal field-stretching (Dinnbier & Kroupa, 2020).

In Figure 1 panels (1.b–3.b), we display the spatial distributions of the
clusters in the $X$-$Z$ plane. Despite the similarity in age, the morphology
of stellar population outside the tidal radius differs significantly. Earlier
studies have already discovered tidal tails in Blanco 1 (Zhang et al., 2020;
Pang et al., 2021). 37% of the mass of Blanco 1 is unbound and is located in
the tidal tails. For the more massive Pleiades and NGC 2516, however, 14% of
the mass (102.5 $\rm M_{\odot}$ and 274.8 $\rm M_{\odot}$, respectively) is
located outside the tidal radii. At a given Galactic position, the cluster
mass is a dominant factor in shaping the 3D spatial morphology. Therefore, the
more massive Pleiades and NGC 2516 can withstand the external tidal stretching
for a longer time. Their structures outside the tidal radius may indicate an
early tidal disruption, while Blanco 1 is already in an advanced stage of
tidal disruption.

To study the tidal structures in Pleiades and NGC 2516, we convert the proper
motions of members into tangential velocities (panels (1.c–3.c) in Figure 1).
All clusters are elongated along $V_{t_{\delta}}$, forming “kinematic tails”.
The bound stars within the tidal radius (grey circles) are concentrated in
kinematic space, whose kinematics is dominated by the internal two-body
relaxation and follows the cluster potential. On the other hand, the
“kinematic tails” are mainly composed of unbound stars influenced by Galactic
tide. We divide the kinematic tails into two subgroups: one with
$V_{t_{\delta}}$ larger (red circles) than the median value (black crosses),
and the other smaller than the median value (blue circles). These two
subgroups are distinct in the 3D spatial distribution (panels b), located on
opposite sides of the cluster. In Blanco 1, the red and blue kinematic tails
correspond to the trailing and leading tails, which resembles the proper
motion tails found by Zhang et al. (2020). Likewise the structure outside the
tidal radius in the Pleiades and NGC 2516 are well separated into two parts by
the red and blue kinematic tails. This correlation supports an early tidal
disruption in these two clusters, with tidal tails just being formed.

The kinematic tails are most apparent in the Pleiades, with a velocity
difference of 3 km s-1, as compared to a difference of 2 km s-1 in Blanco 1
and NGC 2516. This variation is a consequence of the different viewing angles.
The elongated kinematic tails are induced by the Galactic differential
rotation in a group of mainly unbound stars, spanning several tens to 100 pc
in size. The internal velocity dispersion also plays a role in broadening the
kinematic distribution, e.g., the most massive NGC 2516 has the highest
velocity dispersion ($\sigma_{V_{t_{\alpha}}}$, $\sigma_{V_{t_{\delta}}}$).

##  4 Conclusion

The “kinematic tails” in the tangential velocity space successfully reveal
that stars outside the tidal radius in the Pleiades and NGC 2516 are early-
stage tidal structures, compared to the profoundly disrupting Blanco 1.

This study is supported by the Research Development Fund of Xi’an Jiaotong-
Liverpool University (RDF-18–02–32) and the XJTLU Summer Undergraduate
Research Fund (SURF-2021119).

## References

  * Dinnbier & Kroupa (2020) Dinnbier, F., & Kroupa, P. 2020, A&A, 640, A85, doi: 10.1051/0004-6361/201936572
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Lodieu et al. (2019) Lodieu, N., Pérez-Garrido, A., Smart, R. L., & Silvotti, R. 2019, A&A, 628, A66, doi: 10.1051/0004-6361/201935533
  * Pang et al. (2021) Pang, X., Li, Y., Yu, Z., et al. 2021, ApJ, 912, 162, doi: 10.3847/1538-4357/abeaac
  * Röser & Schilbach (2019) Röser, S., & Schilbach, E. 2019, A&A, 627, A4, doi: 10.1051/0004-6361/201935502
  * Röser et al. (2019) Röser, S., Schilbach, E., & Goldman, B. 2019, A&A, 621, L2, doi: 10.1051/0004-6361/201834608
  * Tang et al. (2019) Tang, S.-Y., Pang, X., Yuan, Z., et al. 2019, ApJ, 877, 12, doi: 10.3847/1538-4357/ab13b0
  * Yuan et al. (2018) Yuan, Z., Chang, J., Banerjee, P., et al. 2018, ApJ, 863, 26, doi: 10.3847/1538-4357/aacd0d
  * Zhang et al. (2020) Zhang, Y., Tang, S.-Y., Chen, W. P., Pang, X., & Liu, J. Z. 2020, ApJ, 889, 99, doi: 10.3847/1538-4357/ab63d4


# Diffusion of light in structurally anisotropic media with uniaxial symmetry

Ernesto Pini<EMAIL_ADDRESS>Department of Physics and Astronomy,
Università di Firenze, Sesto Fiorentino, Italy
European Laboratory for Non-linear Spectroscopy (LENS), Sesto Fiorentino,
Italy
Fabrizio Martelli  Department of Physics and Astronomy, Università di Firenze,
Sesto Fiorentino, Italy
Alexander Gatto  Sony Europe B.V., Stuttgart Technology Center, Stuttgart,
Germany
Henrik Schäfer  Sony Europe B.V., Stuttgart Technology Center, Stuttgart,
Germany
Diederik S. Wiersma  Department of Physics and Astronomy, Università di
Firenze, Sesto Fiorentino, Italy
European Laboratory for Non-linear Spectroscopy (LENS), Sesto Fiorentino,
Italy
Istituto Nazionale di Ricerca Metrologica (INRiM), Turin, Italy
Lorenzo Pattelli<EMAIL_ADDRESS>European Laboratory for Non-linear
Spectroscopy (LENS), Sesto Fiorentino, Italy
Istituto Nazionale di Ricerca Metrologica (INRiM), Turin, Italy
National Research Council – National Institute of Optics (CNR-INO), Sesto
Fiorentino, Italy

###### Abstract

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 ${\bf\it 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 open-source 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:

$\left(\frac{1}{v}\partialderivative{t}-\divergence{{\bf\it
D}}\gradient+\mu_{\text{a}}\right)\Phi({\bf\it r},t)=S({\bf\it r},t),$ (1)

where $\Phi({\bf\it r},t)$ is the fluence rate at a position ${\bf\it
r}=(x,y,z)$ and time $t$, $\mu_{\text{a}}$ is the absorption coefficient,
$S({\bf\it r},t)$ is the source term, $v=c/n$ is the speed of light in the
medium and ${\bf\it D}$ is the anisotropic diffusion tensor

${\bf\it D}=\begin{pmatrix}D_{x}&0&0\\\ 0&D_{y}&0\\\ 0&0&D_{z}\end{pmatrix},$
(2)

having units of length.

At the microscopic level, we express the scattering anisotropy via a
direction-dependent scattering coefficient $\mu_{\text{s}}(\hat{\mathbf{s}})$
, which may exhibit a complex dependence on the unit vector $\hat{\mathbf{s}}$
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 ${\bf\it\mu}_{\text{s}}$,

${\bf\it\mu}_{\text{s}}=\begin{pmatrix}\mu_{\text{s},x}&0&0\\\
0&\mu_{\text{s},y}&0\\\ 0&0&\mu_{\text{s},z}\end{pmatrix},$ (3)

through the relation
$\mu_{\text{s}}(\hat{\mathbf{s}})=\hat{\mathbf{s}}{\bf\it\mu}_{\text{s}}\hat{\mathbf{s}}^{\text{T}}$.
Similar tensors could be defined for the refractive index $n$ or the single-
scattering 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
$\ell_{i}$ and transport mean free paths $\ell_{i}^{*}$ along $i=\\{x,y,z\\}$
will be indicated as:

$\mu_{\text{s},i}=\frac{1}{\ell_{i}},\qquad D_{i}=\frac{1}{3}\ell_{i}^{*},$
(4)

where it should be stressed that, despite setting $g=0$, in general
$\ell_{i}\neq\ell_{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 ${\bf\it K}$

${\bf\it K}=\begin{pmatrix}K_{x}&0&0\\\ 0&K_{y}&0\\\ 0&0&K_{z}\end{pmatrix},$
(5)

with elements defined via the relation:

$\ell^{*}_{i}=K_{i}\ell_{i}.$ (6)

In other words, the tensor ${\bf\it K}$ represents the transformation
connecting the microscopic scattering properties to the macroscopic diffusion
rates:

${\bf\it D}=\frac{1}{3}{\bf\it K}{\bf\it\mu}_{\text{s}}^{-1}.$ (7)

In the isotropic case ${\bf\it K}$ becomes the $3\times 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}=\ell_{z}$ from the surface
$S({\bf\it r},t)=\delta(x)\delta(y)\delta(z-z_{0})\delta(t)$ (note that in the
anisotropic case $z_{0}=\ell_{z}\neq\ell_{z}^{*}$, even when assuming $g=0$).
The solution of the time-dependent ADE can be written as:

$\Phi({\bf\it
r},t)=\frac{v\exp(-\mu_{\text{a}}vt)\exp\left(-\frac{x^{2}}{4D_{x}vt}\right)\exp\left(-\frac{y^{2}}{4D_{y}vt}\right)}{(4\pi
vt)^{3/2}\sqrt{D_{x}D_{y}D_{z}}}\left[\sum_{m=-\infty}^{+\infty}\exp\left(-\frac{(z-z_{+,m})^{2}}{4D_{z}vt}\right)-\exp\left(-\frac{(z-z_{-,m})^{2}}{4D_{z}vt}\right)\right],$
(8)

with

$\displaystyle z_{+,m}$ $\displaystyle=2m\left(L+2z_{\text{e}}\right)+z_{0},$
$\displaystyle z_{-,m}$
$\displaystyle=2m\left(L+2z_{\text{e}}\right)-2z_{\text{e}}-z_{0},$

and $z_{\text{e}}$ as the extrapolated length. In order to solve the ADE as a
function of the microscopic scattering coefficient elements $\mu_{i}$, closed-
form expressions for the diffusive rates $D_{i}$ and extrapolated boundary
length $z_{\text{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 $\phi$ 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
$\mu_{\text{s},z}=\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}$,
$\mu_{\text{s},x}=\mu_{\text{s},y}=\mu_{\text{s},{\perp}}$, and
$D_{z}=D_{\stretchrel*{\parallel}{\perp}}$, $D_{x}=D_{y}=D_{\perp}$. Given a
generic unit vector associated to the walker direction of propagation
$\hat{\mathbf{s}}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$, with
$\theta$ and $\phi$ as the polar and azimuthal angles of the absolute
(laboratory) reference frame, we find as expected that the direction-dependent
scattering coefficient $\mu_{\text{s}}(\hat{\mathbf{s}})$ is independent of
the azimuthal angle $\phi$:

$\mu_{\text{s}}(\theta)=\mu_{\text{s},{\perp}}\sin^{2}\theta+\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}\cos^{2}\theta=\\\
=\mu_{\text{s},{\perp}}(1-\chi^{2})+\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}\chi^{2}=\mu_{\text{s}}(\chi),$
(9)

where $\chi=\cos\theta$.

A first relevant quantity that can be derived is the expectation value of the
mean free path $\expectationvalue{\ell}$

$\expectationvalue{\ell}=\int_{4\pi}\xi(\hat{\mathbf{s}})\ell(\hat{\mathbf{s}})\differential{\hat{\mathbf{s}}},$
(10)

where $\xi(\hat{\mathbf{s}})$ represents the probability that a walker travels
within a solid angle $\differential{\hat{\mathbf{s}}}$ about
$\hat{\mathbf{s}}$ after a scattering event, and
$\ell(\hat{\mathbf{s}})=1/\mu_{\text{s}}(\chi)$ is the direction-dependent
scattering mean free path along the same direction [17]. For an isotropic
phase function, we get that $\xi(\hat{\mathbf{s}})=1/4\pi$, so that the only
source of anisotropy is given by $\ell(\hat{\mathbf{s}})$ and the integral can
be simplified to:

$\expectationvalue{\ell}=\frac{1}{2}\int_{-1}^{1}\frac{1}{\mu_{\text{s}}(\chi)}\differential{\chi}=\frac{\tanh^{-1}\sqrt{\frac{\mu_{\text{s},{\perp}}-\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}{\mu_{\text{s},{\perp}}}}}{\sqrt{\mu_{\text{s},{\perp}}}\sqrt{\mu_{\text{s},{\perp}}-\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}}.$
(11)

The diffusion coefficients, on the other hand, are defined as [17]:

$\displaystyle D_{\perp}$
$\displaystyle=\frac{1}{4\expectationvalue{\ell}}\int_{-1}^{1}\frac{1-\chi^{2}}{\mu_{\text{s}}^{2}(\chi)}\differential{\chi},$
(12) $\displaystyle D_{\stretchrel*{\parallel}{\perp}}$
$\displaystyle=\frac{1}{2\expectationvalue{\ell}}\int_{-1}^{1}\frac{\chi^{2}}{\mu_{\text{s}}^{2}(\chi)}\differential{\chi},$
(13)

which also admit closed-form solutions:

$\displaystyle D_{\perp}$
$\displaystyle=\frac{1}{4}\left(\frac{2\mu_{\text{s},{\perp}}-\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}{\mu_{\text{s},{\perp}}^{2}-\mu_{\text{s},{\perp}}\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}-\frac{1}{\sqrt{\mu_{\text{s},{\perp}}}\sqrt{\mu_{\text{s},{\perp}}-\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}\tanh^{-1}\sqrt{\frac{\mu_{\text{s},{\perp}}-\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}{\mu_{\text{s},{\perp}}}}}\right),$
(14) $\displaystyle D_{\stretchrel*{\parallel}{\perp}}$
$\displaystyle=\frac{1}{2}\left(\frac{1}{\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}-\mu_{\text{s},{\perp}}}+\frac{\sqrt{\mu_{\text{s},{\perp}}}}{\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}\sqrt{\mu_{\text{s},{\perp}}-\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}\tanh^{-1}\sqrt{\frac{\mu_{\text{s},{\perp}}-\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}{\mu_{\text{s},{\perp}}}}}\right).$
(15)

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 $\mu_{\text{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:

$\displaystyle D_{\perp}$
$\displaystyle=\frac{\ell_{\perp}}{4}\left[1-\frac{\ell_{\stretchrel*{\parallel}{\perp}}(\expectationvalue{\ell}-\ell_{\perp})}{\expectationvalue{\ell}(\ell_{\perp}-\ell_{\stretchrel*{\parallel}{\perp}})}\right]=\frac{K_{\perp}\ell_{\perp}}{3},$
(16) $\displaystyle D_{\stretchrel*{\parallel}{\perp}}$
$\displaystyle=\frac{\ell_{\stretchrel*{\parallel}{\perp}}}{2}\left[\frac{\ell_{\perp}(\expectationvalue{\ell}-\ell_{\stretchrel*{\parallel}{\perp}})}{\expectationvalue{\ell}(\ell_{\perp}-\ell_{\stretchrel*{\parallel}{\perp}})}\right]=\frac{K_{\stretchrel*{\parallel}{\perp}}\ell_{\stretchrel*{\parallel}{\perp}}}{3},$
(17)

with:

$\displaystyle K_{\perp}$
$\displaystyle=\frac{3}{4}\left[1-\frac{\ell_{\stretchrel*{\parallel}{\perp}}(\expectationvalue{\ell}-\ell_{\perp})}{\expectationvalue{\ell}(\ell_{\perp}-\ell_{\stretchrel*{\parallel}{\perp}})}\right],$
(18) $\displaystyle K_{\stretchrel*{\parallel}{\perp}}$
$\displaystyle=\frac{3}{2}\left[\frac{\ell_{\perp}(\expectationvalue{\ell}-\ell_{\stretchrel*{\parallel}{\perp}})}{\expectationvalue{\ell}(\ell_{\perp}-\ell_{\stretchrel*{\parallel}{\perp}})}\right].$
(19)

Notably, these equation reveal a conservation property for the trace of the
anisotropy tensor:

$\Tr({\bf\it K})=K_{\perp}+K_{\perp}+K_{\stretchrel*{\parallel}{\perp}}=3,$
(20)

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
$\mu_{\text{s},x}\neq\mu_{\text{s},y}\neq\mu_{\text{s},z}$).

The analytical solution for the extrapolated length $z_{\text{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_{\text{e}}=C/B$ [13]:

$\displaystyle B$
$\displaystyle=\int_{0}^{\pi/2}\differential{\theta}\int_{0}^{2\pi}\differential{\phi}P(\hat{\mathbf{s}})\sin\theta\cos\theta,$
(21) $\displaystyle C$
$\displaystyle=\int_{0}^{\pi/2}\differential{\theta}\int_{0}^{2\pi}\differential{\phi}P(\hat{\mathbf{s}})\sin\theta\cos^{2}\theta\,\ell(\hat{\mathbf{s}}),$
(22)

where $P(\hat{\mathbf{s}})$ describes the angular distribution of the
radiance, i.e., the probability density function that a random walker is
traveling along a direction $\hat{\mathbf{s}}$ after many scattering events.
In the general case, $P(\hat{\mathbf{s}})$ depends both on the scattering rate
tensor and the asymmetry factor $g$. In our case, we consider $g=0$ so that to
get $P(\hat{\mathbf{s}})=\ell(\hat{\mathbf{s}})/\expectationvalue{\ell}$,
which further simplifies for uniaxial symmetry to
$P(\chi)=1/\mu_{\text{s}}(\chi)\expectationvalue{\ell}$. We note that our
equation (22) deviates slightly from that reported by Alerstam due to its
dependence on $\ell(\hat{\mathbf{s}})$ rather than
$\ell^{*}(\hat{\mathbf{s}})$ [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.2.

The previous integrals simplify to:

$\displaystyle B$ $\displaystyle=\frac{2\pi
v}{\expectationvalue{\ell}}\int_{0}^{1}\frac{\chi}{\mu_{\text{s}}(\chi)}\differential{\chi},$
(23) $\displaystyle C$ $\displaystyle=\frac{2\pi
v}{\expectationvalue{\ell}}\int_{0}^{1}\frac{\chi^{2}}{\mu^{2}_{\text{s}}(\chi)}\differential{\chi}.$
(24)

to give an analytical expression for $z_{\text{e}}$

$z_{\text{e}}=\frac{C}{B}=\frac{\frac{1}{\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}-\frac{\tanh^{-1}\sqrt{(\mu_{\text{s},{\perp}}-\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}})/\mu_{\text{s},{\perp}}}}{\sqrt{\mu_{\text{s},{\perp}}}\sqrt{\mu_{\text{s},{\perp}}-\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}}}{\ln\frac{\mu_{\text{s},{\perp}}}{\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}}}=\frac{\ell_{\stretchrel*{\parallel}{\perp}}-\expectationvalue{\ell}}{\ln\frac{\ell_{\stretchrel*{\parallel}{\perp}}}{\ell_{\perp}}}.$
(25)

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 $\mu_{\text{s},{\perp}}$ and
$\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}$ tend to the same value
$\mu_{\text{s}}=1/\ell$.

##  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.

###  IV.1 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

$\displaystyle R(x,y,t)$ $\displaystyle=D_{z}\frac{\partial}{\partial
z}\Phi(x,y,z=0,t),$ (26) $\displaystyle T(x,y,t)$
$\displaystyle=-D_{z}\frac{\partial}{\partial z}\Phi(x,y,z=L,t).$ (27)

$R(x,y,t)$ and $T(x,y,t)$ are bi-variate Gaussian distributions:

$\displaystyle R(x,y,t)$
$\displaystyle=R_{0}(t)\exp\left(-\frac{x^{2}}{4D_{x}vt}\right)\exp\left(-\frac{y^{2}}{4D_{y}vt}\right),$
(28) $\displaystyle T(x,y,t)$
$\displaystyle=T_{0}(t)\exp\left(-\frac{x^{2}}{4D_{x}vt}\right)\exp\left(-\frac{y^{2}}{4D_{y}vt}\right),$
(29)

with mean square displacements (MSD) along $x$ and $y$

$w_{x}^{2}(t)=4D_{x}vt,\qquad w_{y}^{2}(t)=4D_{y}vt$ (30)

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 $\ell_{x}^{*}$, $\ell_{y}^{*}$) along $x$ and $y$. Moreover, the
evolution of $w_{x}^{2}(t)$ and $w_{y}^{2}(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\text{\,}\mathrm{mm}$$, setting a
uniaxial scattering anisotropy along the $y$ axis
($\mu_{\text{s},x}=\mu_{\text{s},z}=\mu_{\text{s},{\perp}}$,
$\mu_{\text{s},y}=\mu_{\text{s},{\stretchrel*{\parallel}{\perp}}}$). We set
$\ell_{x}=\ell_{z}=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, and vary
$\ell_{y}$ between $10\text{\,}\mathrm{\SIUnitSymbolMicro
m}100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ with
$10\text{\,}\mathrm{\SIUnitSymbolMicro 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
$\ell_{x}^{*}$ and $\ell_{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\text{\,}\mathrm{\char 37\relax}$ along the longitudinal axis and
$0.2\text{\,}\mathrm{\char 37\relax}$ in the azimuthal plane, which is
representative of the approximations of diffusion theory. This discrepancy
grows up to a maximum of $1.4\text{\,}\mathrm{\char 37\relax}$ for the largest
value of $\ell_{y}=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, corresponding
to an extreme anisotropy ratio of
$\ell_{\stretchrel*{\parallel}{\perp}}/\ell_{\perp}=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.

Figure 1: Mean free paths along $x$ and $y$ obtained with MC simulations and
the analytical formulas as a function of $\ell_{y}$. Dashed lines represent
the simplistic assumption $\ell^{*}_{i}=\ell_{i}$.

###  IV.2 Time-resolved and steady-state profiles

Steady-state and total time-domain transmittance can be derived by direct
integration of equation (27) as

$T(x,y)=\int_{0}^{\infty}T(x,y,t)\differential{t},$ (31)

and

$T(t)=\iint_{-\infty}^{\infty}T(x,y,t)\differential{x}\differential{y},$ (32)

with analogue expressions for the reflectance.

A first validation test regards the definition of the extrapolated length
$z_{\text{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 $\ell^{*}(\hat{\mathbf{s}})$ and to use $z_{0}=\ell_{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
$\ell(\hat{\mathbf{s}})$ and setting $z_{0}=\ell_{z}$ (Figure 2).

Figure 2: Total time-resolved transmittance through a
$1\text{\,}\mathrm{mm}$-thick slab with
$\ell_{x}=$30\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$\ell_{y}=\ell_{z}=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. Numerical
results for pencil beam illumination are compared to the analytical solutions
using different definitions for the extrapolated length $z_{\text{e}}$ and for
the equivalent depth $z_{0}$ of the isotropic point source.

A final set of MC simulations for the time- and space-resolved transmittance
$T(x,y,t)$ was performed for slabs with $g=0$, $n=1$, $\mu_{\text{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\text{\,}\mathrm{\char 37\relax}$ in the considered spatio-temporal range.

Figure 3: (a) Steady-state transmittance $T(x,y)$ from a
$1\text{\,}\mathrm{mm}$-thick slab with
$\ell_{x}=$30\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$\ell_{y}=\ell_{z}=$10\text{\,}\mathrm{\SIUnitSymbolMicro 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).

##  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 $\mathrm{\Delta}h$ and
the original height $h_{0}$), spanning across the elastic-plastic deformation
transition of XPS (Figure 4)

Figure 4: Experimental XPS samples under different compressive strain. The
uncompressed height of each foam sample is $h_{0}=$29\text{\,}\mathrm{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\text{\,}\mathrm{\SIUnitSymbolMicro m}$.

Anisotropic light transport was studied by means of an optical gating imaging
apparatus allowing to capture the evolution of intensity profiles with
sub-$\mathrm{ps}$ resolution [23, 24]. A collection of frames recorded for
each sample at a probe wavelength of $820\text{\,}\mathrm{nm}$ is shown in
Figure 5 for $t\leq$15\text{\,}\mathrm{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 ($\ell_{x}=\ell_{z}$) and longitudinal ($\ell_{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\text{\,}\mathrm{\char 37\relax}$) to
$n=1.043$ for the most compressed sample (solid fraction
$14\text{\,}\mathrm{\char 37\relax}$). These values are used for the correct
determination of the energy velocity $v$ inside the sample.

Figure 5: Time-resolved imaging measurements of the reflected intensity for
the 5 samples of extruded polystyrene. The relative delay between consecutive
frames is $1\text{\,}\mathrm{ps}$. Each frame covers an area of
$5.35\times$5.35\text{\,}{\mathrm{mm}}^{2}$$ and is normalized to its maximum
value.

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 $\ell_{x}$ and $\ell_{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 $\ell_{y}$ is progressively reduced with
increasing compression. On the other hand, $\ell_{x}=\ell_{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({\bf\it 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].

Figure 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.

##  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 case of
isotropic scattering ($g=0$) and index-matched boundary conditions, the
resulting closed-form expression allow to solve directly the anisotropic
diffusion equation both in the spatial and temporal domains. The results have
been validated against a general anisotropic Monte Carlo implementation,
showing excellent agreement up to high degrees of anisotropy between the
longitudinal and azimuthal directions.

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.

## 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 - NextGenerationEU, National Recovery and Resilience Plan,
MNESYS, PE0000006 (DN 572 1553 11.10.2022) and PRIN 2022, DIRS, grant 578
number: 2022EB4B7E.

## Acknowledgments

M. Bürmen and P. Naglič are kindly acknowledged for fruitful discussion.

## References

  * [1] A. Kienle, C. D’Andrea, F. Foschum, _et al._ , “Light propagation in dry and wet softwood,” Optics Express 16, 9895–9906 (2008).
  * [2] A. Kienle, C. Wetzel, A. Bassi, _et al._ , “Determination of the optical properties of anisotropic biological media using an isotropic diffusion model,” Journal of biomedical optics 12, 014026 (2007).
  * [3] E. Simon, P. Krauter, and A. Kienle, “Time-resolved measurements of the optical properties of fibrous media using the anisotropic diffusion equation,” Journal of Biomedical Optics 19, 075006 (2014).
  * [4] A. Kienle, F. Forster, and R. Hibst, “Anisotropy of light propagation in biological tissue,” Optics Letters 29, 2617–2619 (2004).
  * [5] A. Kienle, F. K. Forster, R. Diebolder, and R. Hibst, “Light propagation in dentin: influence of microstructure on anisotropy,” Physics in Medicine & Biology 48, N7 (2002).
  * [6] S. Nickell, M. Hermann, M. Essenpreis, _et al._ , “Anisotropy of light propagation in human skin,” Physics in Medicine & Biology 45, 2873 (2000).
  * [7] A. Sviridov, V. Chernomordik, M. Hassan, _et al._ , “Intensity profiles of linearly polarized light backscattered from skin and tissue-like phantoms,” Journal of Biomedical Optics 10, 014012 (2005).
  * [8] P. Johnson, S. Faez, and A. Lagendijk, “Full characterization of anisotropic diffuse light.” Optics Express 16, 7435–7446 (2008).
  * [9] P. M. Johnson and A. Lagendijk, “Optical anisotropic diffusion: new model systems and theoretical modeling,” Journal of Biomedical Optics 14, 054036–054036 (2009).
  * [10] J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Physical Review E 68, 031908 (2003).
  * [11] A. Kienle, “Anisotropic light diffusion: an oxymoron?” Physical Review Letters 98, 218104 (2007).
  * [12] A. Kienle, F. Foschum, and A. Hohmann, “Light propagation in structural anisotropic media in the steady-state and time domains,” Physics in Medicine & Biology 58, 6205 (2013).
  * [13] E. Alerstam, “Anisotropic diffusive transport: Connecting microscopic scattering and macroscopic transport properties,” Physical Review E 89, 063202 (2014).
  * [14] S. E. Han, “Transport mean free path tensor and anisotropy tensor in anisotropic diffusion equation for optical media,” Journal of Optics 22, 075606 (2020).
  * [15] R. Vasques and E. W. Larsen, “Non-classical particle transport with angular-dependent path-length distributions. II: Application to pebble bed reactor cores,” Annals of Nuclear Energy 70, 301–311 (2014).
  * [16] T. Tran-Duc, E. Bertevas, N. Phan-Thien, and B. C. Khoo, “Simulation of anisotropic diffusion processes in fluids with smoothed particle hydrodynamics,” International Journal for Numerical Methods in Fluids 82, 730–747 (2016).
  * [17] F. Martelli, T. Binzoni, S. Del Bianco, _et al._ , _Light propagation through biological tissue and other diffusive media: theory, solutions, and validations_ (SPIE, 2022).
  * [18] P. Naglič, Y. Zelinskyi, F. Pernuš, _et al._ , “PyXOpto: an open-source python library with utilities for fast light propagation modeling in turbid media,” in _European Conference on Biomedical Optics,_ (Optica Publishing Group, 2021), pp. EM3C–2.
  * [19] R. C. Haskell, L. O. Svaasand, T.-T. Tsay, _et al._ , “Boundary conditions for the diffusion equation in radiative transfer,” JOSA A 11, 2727–2741 (1994).
  * [20] A. Krundaeva, G. De Bruyne, F. Gagliardi, and W. Van Paepegem, “Dynamic compressive strength and crushing properties of expanded polystyrene foam for different strain rates and different temperatures,” Polymer Testing 55, 61–68 (2016).
  * [21] S. Vaitkus, S. Vėjelis, and A. Kairytė, “Analysis of extruded polystyrene short-term compression dependence on exposure time,” Materials Science 19, 471–474 (2013).
  * [22] E. Sadek and N. Fouad, “Finite element modeling of compression behavior of extruded polystyrene foam using X-ray tomography,” Journal of Cellular Plastics 49, 161–191 (2013).
  * [23] L. Pattelli, R. Savo, M. Burresi, and D. S. Wiersma, “Spatio-temporal visualization of light transport in complex photonic structures,” Light: Science & Applications 5, e16090–e16090 (2016).
  * [24] E. Pini, G. Mazzamuto, F. Riboli, _et al._ , “Breakdown of self-similarity in light transport,” arXiv preprint arXiv:2304.02773 (2023).
  * [25] M. Burresi, L. Cortese, L. Pattelli, _et al._ , “Bright-white beetle scales optimise multiple scattering of light,” Scientific Reports 4, 6075 (2014).
  * [26] L. Cortese, L. Pattelli, F. Utel, _et al._ , “Anisotropic light transport in white beetle scales,” Advanced Optical Materials 3, 1337–1341 (2015).
  * [27] F. Utel, L. Cortese, D. S. Wiersma, and L. Pattelli, “Optimized white reflectance in photonic-network structures,” Advanced Optical Materials 7, 1900043 (2019).
  * [28] B. K. Park, S. M. Han, and S. E. Han, “Surpassing Cyphochilus scales in optical scattering strength by well-controlled electrospun nanostructures,” Optical Materials Express 12, 2529–2540 (2022).