resulting from the simulations. Each point is obtained considering only the energy packets transmitted within the corresponding time bin. Dashed lines represent the nominal values for the two distributions. (b) Time evolution of the SLD for *n* = 1 for energy packets transmitted at *t* = 10, 20, 30, 40, 50, 60, 70, 80 and 90 ps. Gray and purple curves show respectively the histogram of the step lengths drawn through the PRNG and of the steps taken inside the sample. The two only differ for the last step of

angle is drawn. Actually, since a typical step will not be in general perfectly parallel to the interfaces, a scattering angle of θ ≈ 180° should provide higher chances of staying inside the sample, hence its higher probability. This results in a cos θ distribution with a slightly negative average value (Figure 6.3a, left axis), which also plays a role in determining the

axis), along with its nominal value of 1 (dashed line). At long times, each curve seems to saturate to an asymptotic value, suggesting the existence of a well-defined effective diffusion coefficient. The random-walk based picture of diffusion as expressed by equation (6.1) is qualitatively supported by the fact that also this figure of merit is enhanced for *n* = 1.016 (red curve), in accordance with Figure 6.1a. In principle, the overall diffusion process will be influenced by both the modified step-length and angular statistics, which in the investigated configurations appear to have opposite effects, as also shown in Figure 6.3a. While the latter would indeed tend to slightly slow down diffusion, the predominant effect is coming from the step lengths being substantially increased, leading to the observed enhanced in-plane diffusion especially for *n* = 1.016. Notably, different configurations might lead to a different balance between these two effects, which also appear to saturate to their respective asymptotic values on slightly different time scales, further illustrating the need for additional investigations even for the simple model of a homogeneous and

The asymptotic nature of the effective diffusion coefficient in a thin slab is further highlighted in Figure 6.3b, where the time evolution of the SLD is shown for *n* = 1 (the *n* = 1.016 and 1.1 cases are analogous). The time-resolved distributions seem to converge towards a single asymptotic envelope distribution with a well-defined asymptotic decay rate which seems to be uniquely determined by the properties of the sample. It is interesting to

104

2 /2 [mm]

10−8

10−4

*p*(; *t*)

100

0 10 20

(b)

/2 appearing in equation (6.1) and of cos θ

2

/*l*s

/2 in Figure 6.3a (right

drawn taken nominal

1

1.5

−0.05

each trajectory.

isotropic single slab.

cos θ

0

0.05

Figure 6.4. (a) MSW expansion for the transmittance of a slab with (*n* = 1.5, *g* = 0, OT−1 = 0.1). As already discussed, in-plane propagation is slightly enhanced with respect to the DA 4*D*DA*t* prediction. This deviation is apparently not supported by the asymptotic (b) step-length distribution (exhibiting a slightly lighter long-step tail than expected) and (c) scattering angle distribution (favoring backwards angles), that we approximate here as the effective distributions after a delay *t* = 15τ.

compare the histogram of the actual steps performed inside the sample (blue curves) with the histogram of the ones drawn through the PRNG (gray curves). The two differ only for the last step, whose length is respectively considered either partially (up to the interface) or totally. At late times the two sets of curves become indistinguishable since, as expected, the contribution of the last step to the whole trajectory becomes eventually statistically negligible.

As a result of the transport statistics being directly altered by the sample configuration, an optically thin sample generally appears to be less scattering than it actually is. In other words, once the diffusive regime is reached, energy packets propagate as if scatterers were further apart than they actually are, i.e., with an *e*ff*ective* transport mean free path greater than the one *intrinsic* to the material. Indeed, because of the asymptotic nature of these effects, only a small fraction of the incoming light is actually subject to this effective transport mean free path when studying thin samples. Yet, the effect is largely accessible experimentally [219] and, as we have seen in Section 4.2, similar deviations are can in fact be found even in more turbid media. As a matter of fact, the asymptotic nature of this effective transport regime makes it even more relevant from an experimental point of view, since it is often believed that the standard DA becomes progressively safer to apply as later times becomes accessible. Moreover, other applications can be envisioned where multiple scattering in thin layers or confined geometries, even if limited to a very small fraction of incident light, could play a significant role [220, 221].

Albeit smaller, it is interesting to investigate whether the discrepancies that we discussed in Chapter 4 for the homogeneous slab arise from the same mechanism that we have described, which would suggest that modifications to the transport statistics can be still appreciable in more turbid media. At the same time, it is also interesting to determine how these deviations are gradually suppressed going towards optically thicker media. To this purpose, we performed a new simulation of a slab sample with *l*s = *l* s = 0.1 mm and *L*0 = 1 mm (OT = 10), using a refractive index contrast of *n* = 1.5 similar to that of our experimental

6.3.1. Anisotropic transport in isotropic media

microscopic optical properties such as *g* or *l*

generic step-length component will be therefore given by

*P*(*UV* ≤ *x*|*U* = *u*)*pU*(*u*) d*u* =

1 − exp (−µs *x*/*u*) d*u* = 1 −

= *P V* ≤ *x u*

the two independent random variables and cos θ.

can be written as

*P*(*UV* ≤ *x*) =

*P*(*X* ≤ *x*) =

integration by parts eventually we obtain

 1 0

As we have seen in the previous Section, our investigation on the random-walk statistics arising in the slab geometry revealed a subtle interplay occurring between the actual thickness of the slab, the refractive index contrast and the scattering anisotropy, determining a transport regime that is basically diffusive on long time scales but which cannot be described in terms of the simple diffusive approximation. A different and asymptotic regime naturally emerges from the overall optical and geometric boundary conditions of the sample, and is univocally determined by them through yet unknown relations. In this respect, our findings recall a recently published work where it is analogously demonstrated that the link between microscopic (i.e., the scattering coefficient) and macroscopic (i.e., the diffusion coefficient) transport parameters remains unknown for diffusive anisotropic media [152]. Analogously, our results show that this link should be further investigated even in the isotropic case, especially for weakly scattering media. In particular, concerning

(e.g., tensorial instead of scalar) depending on incidental geometric conditions.

between an intrinsic and an effective counterpart, where the former is the one that we are typically interested in retrieving while the latter might have a very different value and nature

which enables us to derive the probability density function *p*(*i*) as that of the product of

A convenient choice is that of using the cumulative distribution of the exponential distribution (2.60), which is non-null only for > 0, and the PDF of the uniform distribution, which has a limited support. By plugging their expressions into equation (6.4), following

107

 1 0

In general, the cumulative distribution function (CDF) for the product variable *X* = *UV*

Driven by this observation, we performed a new analysis of our dataset, taking into account the translation symmetry of the problem along the *xy*-plane of the slab, collecting separate statistics on the step-length components taken along the different axes. In order to do so, it is convenient to derive the nominal probability density function (PDF) for the step components. This can be easily done in the simple isotropic (*g* = 0) case, where angles are sampled uniformly on a sphere. This allows to consider interchangeably the (fixed) reference frame of the slab and the (rotating) reference of the simulated energy packet. A

s, it seems appropriate to introduce a distinction

*i* = cos θ, with *i* = *x*, *y*, *z* (6.3)

*P*(*uV* ≤ *x*)*pU*(*u*) d*u*

 *PV x u*

*pU*(*u*) d*u* . (6.4)

*pU*(*u*) d*u* =

exp (−µs *x*/*u*) d*u*

Figure 6.5. Time evolution of the ratio 2 /2 and of cos θ for a slab sample with (*n* = 1.5, *g* = 0, OT−1 = 0.1). Both figures converge asymptotically to a value that is lower than the nominal value of the simulation (dashed lines).

slab sample (see Section 4.2). As we have seen in Chapter 5, the MSW expansion rate of such sample is slightly higher than that predicted by the DA. However, inspecting the step-length and scattering angle distribution reveals an apparently contradictory picture. As Figure 6.4 shows, while the MSW grows at a rate that is appreciably faster (of the order of %) of that expected, this does not reflect in the asymptotic (*t* = 15τ) step-length nor in the scattering-angle distributions. As a matter of fact, the former is only slightly modified with respect to the nominal exponential, perhaps even exhibiting a *lower* probability of taking long steps, while the latter is slightly biased towards backwards scattering angles, which also goes in the direction of a smaller diffusion rate. This behavior is further confirmed by looking at the full evolution of the 2 /2 and cos θ descriptors (Figure 6.5), showing clearly that the asymptotic average step length and scattering angle are both smaller than expected. This interestingly shows how subtle deviations are still clearly present even at an optical thickness of OT = 10 where the DA is commonly employed, and that they are persistent at extremely long delays. In addition to this, the simulated data reveal a further layer of complexity in the effective transport regime that arises in bounded media, which we will address in the following.

## 6.3. A walk on the wild side of diffusion

In the following subsections, we briefly discuss a few preliminary results regarding unexpected properties of transport in bounded media, which became apparent thanks to the magnitude of the simulations that we performed to calculate the exact solution to the transport problem.

### 6.3.1. Anisotropic transport in isotropic media

0 5 10 15

(a)

the nominal value of the simulation (dashed lines).

Figure 6.5. Time evolution of the ratio

we will address in the following.

transport problem.

6.3. A walk on the wild side of diffusion

*t*/τ

0 5 10 15

(b)

/2 and of cos θ for a slab sample with

*t*/τ

−5

0

cos θ

2

(*n* = 1.5, *g* = 0, OT−1 = 0.1). Both figures converge asymptotically to a value that is lower than

slab sample (see Section 4.2). As we have seen in Chapter 5, the MSW expansion rate of such sample is slightly higher than that predicted by the DA. However, inspecting the step-length and scattering angle distribution reveals an apparently contradictory picture. As Figure 6.4 shows, while the MSW grows at a rate that is appreciably faster (of the order of %) of that expected, this does not reflect in the asymptotic (*t* = 15τ) step-length nor in the scattering-angle distributions. As a matter of fact, the former is only slightly modified with respect to the nominal exponential, perhaps even exhibiting a *lower* probability of taking long steps, while the latter is slightly biased towards backwards scattering angles, which also goes in the direction of a smaller diffusion rate. This behavior is further confirmed by looking at the full evolution of the 2 /2 and cos θ descriptors (Figure 6.5), showing clearly that the asymptotic average step length and scattering angle are both smaller than expected. This interestingly shows how subtle deviations are still clearly present even at an optical thickness of OT = 10 where the DA is commonly employed, and that they are persistent at extremely long delays. In addition to this, the simulated data reveal a further layer of complexity in the effective transport regime that arises in bounded media, which

In the following subsections, we briefly discuss a few preliminary results regarding unexpected properties of transport in bounded media, which became apparent thanks to the magnitude of the simulations that we performed to calculate the exact solution to the

106

5

·10−3

2 /2 [µm] As we have seen in the previous Section, our investigation on the random-walk statistics arising in the slab geometry revealed a subtle interplay occurring between the actual thickness of the slab, the refractive index contrast and the scattering anisotropy, determining a transport regime that is basically diffusive on long time scales but which cannot be described in terms of the simple diffusive approximation. A different and asymptotic regime naturally emerges from the overall optical and geometric boundary conditions of the sample, and is univocally determined by them through yet unknown relations. In this respect, our findings recall a recently published work where it is analogously demonstrated that the link between microscopic (i.e., the scattering coefficient) and macroscopic (i.e., the diffusion coefficient) transport parameters remains unknown for diffusive anisotropic media [152]. Analogously, our results show that this link should be further investigated even in the isotropic case, especially for weakly scattering media. In particular, concerning microscopic optical properties such as *g* or *l* s, it seems appropriate to introduce a distinction between an intrinsic and an effective counterpart, where the former is the one that we are typically interested in retrieving while the latter might have a very different value and nature (e.g., tensorial instead of scalar) depending on incidental geometric conditions.

Driven by this observation, we performed a new analysis of our dataset, taking into account the translation symmetry of the problem along the *xy*-plane of the slab, collecting separate statistics on the step-length components taken along the different axes. In order to do so, it is convenient to derive the nominal probability density function (PDF) for the step components. This can be easily done in the simple isotropic (*g* = 0) case, where angles are sampled uniformly on a sphere. This allows to consider interchangeably the (fixed) reference frame of the slab and the (rotating) reference of the simulated energy packet. A generic step-length component will be therefore given by

$$\ell\_i = \ell \cos \theta, \quad \text{with } i = \ge \text{, y, z} \tag{6.3}$$

which enables us to derive the probability density function *p*(*i*) as that of the product of the two independent random variables and cos θ.

In general, the cumulative distribution function (CDF) for the product variable *X* = *UV* can be written as

$$P(UV \le \mathbf{x}) = \int P(UV \le \mathbf{x}|U = u) p\_U(u) \, \mathrm{d}u = \int P(uV \le \mathbf{x}) p\_U(u) \, \mathrm{d}u$$

$$= \int P\Big(V \le \frac{\mathbf{x}}{u}\Big) p\_U(u) \, \mathrm{d}u = \int P\_V\Big(\frac{\mathbf{x}}{u}\Big) p\_U(u) \, \mathrm{d}u. \quad (6.4)$$

A convenient choice is that of using the cumulative distribution of the exponential distribution (2.60), which is non-null only for > 0, and the PDF of the uniform distribution, which has a limited support. By plugging their expressions into equation (6.4), following integration by parts eventually we obtain

$$P(X \le x) = \int\_0^1 1 - \exp\left(-\mu\_s \ge \mu\right) \,\mathrm{d}u = 1 - \int\_0^1 \exp\left(-\mu\_s \ge \mu\right) \,\mathrm{d}u$$

Figure 6.6. Comparison between derived and simulated probability density functions of (a) step lengths and (b) squared step lengths and their components as obtained from a single trajectory of 1010 steps in an infinitely extended volume.

$$=1-\exp\left(-\mu\_s\chi\right)+\mu\_s\text{x}\Gamma(0,\mu\_s\chi) \quad (6.5)$$

0 5 10

(a)

6.8 where we have plotted the 2

/*l*s

*x z* nominal

0 50 100 150

(b)

2/*l* 2 s

2 *x* 2 *z* nominal

10−12

10−8

*p*(2

Figure 6.7. (a) Decomposition of the SLD shown in Figure 6.4b into its *x* and *z* components, revealing how they are differently affected by the presence of the plane boundary. The expression (6.6) for an unbounded medium is shown for comparison. Due to internal reflections, at long times also the *z*-component SLD continues smoothly beyond *L*0. (b) Decomposition of the squared-step-length

> *X p*(*x*) *x* exp (−/*l*s) /*l*s *l*s

*i* Γ(0, /*l*s) /*l*s *l*s/2

exp (/*l*s) 2*l*

A comparison between these expressions and the step-length distributions of an

By analyzing the late-time step-length distribution of our OT = 10 simulation in terms of their *x* and *z* components (*y* is statistically equivalent to *x*), an interesting picture emerges. Figure 6.7 illustrates the difference between the in-plane and axial components of the steps of those packets transmitted at long times. Despite the fact that the simulated system is isotropic and isotropically scattering, an *anisotropic* transport regime establishes at long times characterized by a constrained step-length distribution along *z*. However, also the step lengths along the plane are modified, as can be more clearly appreciated from Figure

the decreased 'diffusivity' along *z* is what determines the enhanced MSW growth that we have found in Figures 6.4 and 5.3. Nonetheless, the overall effect is still that of a decreased

the onset of this anisotropic transport regime does not seem to affect the shape of the transmitted profiles, whose excess kurtosis converges normally to 0 (Figure 6.8b).

109

As expected, even more pronounced deviations can be found in the previously presented OT = 1 configurations, yet with a qualitatively similar behavior (Figure 6.9), where the *z* component is particularly suppressed due to the low refractive index contrast. In such a confined configuration, in fact, long trajectories will be composed of steps mainly lying in

/*l*s) 2*l*

2 s

2 s /3

*i* /2 *i* ratio. This unexpected over-compensation for

*z* /2 *z*. Notably,

distribution, compared with equation (6.9), showing a qualitatively similar behavior.

2 1 2*l*s √

isotropic trajectory in an unbounded medium is shown in Figure 6.6.

probability of taking long steps, due to the more marked reduction of 2

1 2*l*s √ Γ(0, √

2 *i*

*i* ; *t* = 15τ)

10−4

10−5

*p*(*i*; *t* = 15τ)

where Γ(0, *x*) is the incomplete gamma function (also reported as the exponential integral function *E*1(*x*)). The probability density function is then easily found as the derivative

$$p(\mathbf{x}) = \frac{\mathrm{d}P\_X(\mathbf{x})}{\mathrm{d}\mathbf{x}} = \mu\_\mathrm{s} \Gamma(0, \mu\_\mathrm{s} \mathbf{x}). \tag{6.6}$$

Once we have obtained the PDF, the probability functions for the squared variables is also simply found by exploiting the fact that the step lengths and the length of its components are positive quantities. By doing so, we can generally write

$$p\_{\xi^2}(\mathbf{x}) = \frac{\mathbf{d}}{\mathbf{d}\mathbf{x}} P(\xi^2 \le \mathbf{x}) = \frac{\mathbf{d}}{\mathbf{d}\mathbf{x}} P(\xi \le \sqrt{\mathbf{x}}) = \frac{\mathbf{d}P\_{\xi}}{\mathbf{d}\sqrt{\mathbf{x}}} = p(\sqrt{\mathbf{x}})\frac{1}{2\sqrt{\mathbf{x}}},\tag{6.7}$$

which gives

$$p\_{\ell^2}(\mathbf{x}) = \frac{\mu\_s}{2\sqrt{\mathbf{x}}} \exp\left(-\mu\_s \sqrt{\mathbf{x}}\right) \tag{6.8}$$

$$p\_{\ell\_i^2}(\mathbf{x}) = \frac{\mu\_s}{2\sqrt{\mathbf{x}}} \Gamma(0, -\mu\_s\sqrt{\mathbf{x}}) \tag{6.9}$$

respectively for 2 and 2 *i* . The relevant distributions that we obtained along with their first moments are listed in the following table:

10−10

10−2

10−10

10−6

0 5 10 15 20

(a)

steps in an infinitely extended volume.

*p*ξ2 (*x*) = d

moments are listed in the following table:

which gives

respectively for 2 and 2

d*x*

*i*/*l*s

10−14

10−4

10−14

10−9

0 200 400

= 1 − exp (−µs*x*) + µs*x*Γ(0, µs*x*) (6.5)

d*x* = µsΓ(0, µs*x*). (6.6)

√ *x*) 1 2 √*x*

(6.8)

*x*) (6.9)

, (6.7)

(b)

2 *i* /*l* 2 s

10−9

*p*(2)

*p*(2

Figure 6.6. Comparison between derived and simulated probability density functions of (a) step lengths and (b) squared step lengths and their components as obtained from a single trajectory of 1010

where Γ(0, *x*) is the incomplete gamma function (also reported as the exponential integral function *E*1(*x*)). The probability density function is then easily found as the derivative

Once we have obtained the PDF, the probability functions for the squared variables is also simply found by exploiting the fact that the step lengths and the length of its components

*P*(ξ ≤ √

Γ(0, −µs

108

*x*) = d*P*ξ d

√

*i* . The relevant distributions that we obtained along with their first

√*x* = *p*(

*p*(*x*) = d*PX*(*x*)

d*x*

*p*2 (*x*) = µs 2 √*x* exp −µs √ *x*

> (*x*) = µs 2 √*x*

are positive quantities. By doing so, we can generally write

*P*(ξ2 ≤ *x*) = d

*p*2 *i*

*i* )

10−4

10−6

*p*()

*p*(*i*)

10−2

Figure 6.7. (a) Decomposition of the SLD shown in Figure 6.4b into its *x* and *z* components, revealing how they are differently affected by the presence of the plane boundary. The expression (6.6) for an unbounded medium is shown for comparison. Due to internal reflections, at long times also the *z*-component SLD continues smoothly beyond *L*0. (b) Decomposition of the squared-step-length distribution, compared with equation (6.9), showing a qualitatively similar behavior.

A comparison between these expressions and the step-length distributions of an isotropic trajectory in an unbounded medium is shown in Figure 6.6.

By analyzing the late-time step-length distribution of our OT = 10 simulation in terms of their *x* and *z* components (*y* is statistically equivalent to *x*), an interesting picture emerges. Figure 6.7 illustrates the difference between the in-plane and axial components of the steps of those packets transmitted at long times. Despite the fact that the simulated system is isotropic and isotropically scattering, an *anisotropic* transport regime establishes at long times characterized by a constrained step-length distribution along *z*. However, also the step lengths along the plane are modified, as can be more clearly appreciated from Figure 6.8 where we have plotted the 2 *i* /2 *i* ratio. This unexpected over-compensation for the decreased 'diffusivity' along *z* is what determines the enhanced MSW growth that we have found in Figures 6.4 and 5.3. Nonetheless, the overall effect is still that of a decreased probability of taking long steps, due to the more marked reduction of 2 *z* /2 *z*. Notably, the onset of this anisotropic transport regime does not seem to affect the shape of the transmitted profiles, whose excess kurtosis converges normally to 0 (Figure 6.8b).

As expected, even more pronounced deviations can be found in the previously presented OT = 1 configurations, yet with a qualitatively similar behavior (Figure 6.9), where the *z* component is particularly suppressed due to the low refractive index contrast. In such a confined configuration, in fact, long trajectories will be composed of steps mainly lying in

0 10 20

(a)

Figure 6.9. Time evolution of the (a) anisotropic 2

*t*/τ

*i* = *x*, *n* = 1 *i* = *z*, *n* = 1 *i* = *x*, *n* = 1.016 *i* = *z*, *n* = 1.016 *i* = *x*, *n* = 1.1 *i* = *z*, *n* = 1.1

case the excess kurtosis seems to converge towards 0 without any finite offset.

0 10 20

(b)

*i* /2 *i* ratio for the different refractive index

, (6.10)

2

. (6.11)

(*t*). Ideally, the

*t*/τ

*n* = 1 *n* = 1.016 *n* = 1.1

−0.5

0

γ2

contrasts of the OT = 1 configurations. (b) Despite the much larger deviations observed, also in this

[222, 223]. As a matter of fact, many transport processes exhibit moments of displacements

and the study of *q*ν(*q*) as a function of *q* provides a more complete characterization of the process than does the single scalar 2ν(2). In standard diffusive and anomalous diffusive processes, the function ν(*q*) is actually a constant (i.e., ν = 1/2 in the normal diffusive case, ν > 1/2 for superdiffusion and ν = 1 for ballistic propagation). However, characterizing the dichotomy between normal and anomalous diffusion based only on the value of ν(2) fails

complete information is represented by the Green's function or propagator of the process. Often this is not possible to obtain exactly, though asymptotic methods can be exploited to provide useful approximations in the form of similarity solutions. These are expressed, in

− 1/ν P

In the case of normal diffusion ν = 2 and P is a Gaussian, while in the superdiffusive case there are examples in which P is a Lévy density [224]. The approximation (6.11) is valid only in the central scaling region (CSR), i.e., excluding the non-scaling tails of the distribution which are populated by 'particles' that have undergone exceptionally large displacements. The concept of *strong* self-similarity can be introduced to identify those systems where ν(*q*) = 1/ν , ∀*q*, as it is the case for normal diffusion and Lévy flights. On the contrary, *weak* self-similarity refers to propagation dynamics giving rise to non-trivial ν(*q*) functions. Typically, the small-*q* part of ν(*q*) (which passes through the origin of the (ν(*q*), *q*) plane) refers to the central scaling region of the asymptotic profile containing the majority of the 'particles', while the large-*q* range is determined by relatively few tail individuals which have traveled exceptionally far. For this reason, as confirmed by analytical, numerical and experimental evidence, in most known systems ν(*q*)*q* is a piece-

111

 *x t* 1/ν

at providing the full picture [223], and one should look further than |*x*|

*p*(*x*, *t*) ≈ *t*

the limit of *t* → ∞ by the self-similar collapse relation

0.5

0.5

with asymptotic behavior

2

*i* /2 *i* [mm]

1

Figure 6.8. Time evolution of the (a) anisotropic 2 *i* /2 *i* ratios, revealing how an enhanced in-plane MSW expansion is compatible with an overall negative anisotropy factor *g* and smaller probability of taking long steps, as shown in Figure 6.5. In-plane and axial components undergo opposite modifications, with a more marked suppression along *z* accompanied by a net increase along *x* and *y*. (b) These modifications do not seem to affect the shape of the transmitted profile, whose excess kurtosis γ2 still vanishes with time.

the plane of the slab, with the extreme case of matched refractive index contrast where no step can be taken having an axial component larger than the thickness of the slab.

At this stage, to our knowledge, there is no analytical model available to describe how this anisotropic transport regime sets in depending on the optical thickness and the refractive index contrast of the slab. In particular, even though this regime is largely determined by the presence and type of boundaries, its effect is fundamentally different from apparently similar boundary effects described in the literature [75, 76], which can be usually taken into account through some refined extrapolated boundary conditions. This cannot be the case here, since extrapolated boundary conditions correct significantly quantities such as the total transmittance which, conversely, would be negligibly affected by asymptotic modifications of the effective diffusion coefficient, especially in thicker media. This does not mean that the effect that we described cannot be accessed experimentally. On the contrary, we can now interpret the small discrepancy measured for our homogeneous sample in Chapter 4 as a direct experimental evidence of the onset of this anisotropic regime. State-of-the-art detectors used in time-resolved configuration are able to access 8 decades of dynamic range [219], corresponding to a suppression greater than ∼ e18, a range within which even our weakly scattering configurations have almost reached their asymptotic regime.

### 6.3.2. Anomalous diffusion in homogeneous media

In the previous Chapters, we have often used the mean-square width (i.e., the second moment of the time-resolved profiles) as a convenient descriptor of transport. However, as we have seen in Section 4.3 for a few extreme cases, only the profile shape conveys the full information for a proper characterization. An interesting trade-off which has been recently considered in the literature of continuous-time random walks (CTRWs), consists of analyzing the spatio-temporal evolution of the profiles at different (continuous) moments

0 5 10 15

(a)

excess kurtosis γ2 still vanishes with time.

*t*/τ

Figure 6.8. Time evolution of the (a) anisotropic 2

*i* = *x i* = *z* nominal

−0.1

0

γ2

in-plane MSW expansion is compatible with an overall negative anisotropy factor *g* and smaller probability of taking long steps, as shown in Figure 6.5. In-plane and axial components undergo opposite modifications, with a more marked suppression along *z* accompanied by a net increase along *x* and *y*. (b) These modifications do not seem to affect the shape of the transmitted profile, whose

the plane of the slab, with the extreme case of matched refractive index contrast where no

At this stage, to our knowledge, there is no analytical model available to describe how this anisotropic transport regime sets in depending on the optical thickness and the refractive index contrast of the slab. In particular, even though this regime is largely determined by the presence and type of boundaries, its effect is fundamentally different from apparently similar boundary effects described in the literature [75, 76], which can be usually taken into account through some refined extrapolated boundary conditions. This cannot be the case here, since extrapolated boundary conditions correct significantly quantities such as the total transmittance which, conversely, would be negligibly affected by asymptotic modifications of the effective diffusion coefficient, especially in thicker media. This does not mean that the effect that we described cannot be accessed experimentally. On the contrary, we can now interpret the small discrepancy measured for our homogeneous sample in Chapter 4 as a direct experimental evidence of the onset of this anisotropic regime. State-of-the-art detectors used in time-resolved configuration are able to access 8 decades of dynamic range [219], corresponding to a suppression greater than ∼ e18, a range within which even our

step can be taken having an axial component larger than the thickness of the slab.

weakly scattering configurations have almost reached their asymptotic regime.

In the previous Chapters, we have often used the mean-square width (i.e., the second moment of the time-resolved profiles) as a convenient descriptor of transport. However, as we have seen in Section 4.3 for a few extreme cases, only the profile shape conveys the full information for a proper characterization. An interesting trade-off which has been recently considered in the literature of continuous-time random walks (CTRWs), consists of analyzing the spatio-temporal evolution of the profiles at different (continuous) moments

110

6.3.2. Anomalous diffusion in homogeneous media

0.1

0 5 10 15

(b)

*i* /2 *i* ratios, revealing how an enhanced

*t*/τ

MC data DA

2

*i* /2 *i* [µm]

Figure 6.9. Time evolution of the (a) anisotropic 2 *i* /2 *i* ratio for the different refractive index contrasts of the OT = 1 configurations. (b) Despite the much larger deviations observed, also in this case the excess kurtosis seems to converge towards 0 without any finite offset.

[222, 223]. As a matter of fact, many transport processes exhibit moments of displacements with asymptotic behavior

$$
\langle |\mathbf{x}|^{q} \rangle \sim t^{q\gamma(q)},\tag{6.10}
$$

and the study of *q*ν(*q*) as a function of *q* provides a more complete characterization of the process than does the single scalar 2ν(2). In standard diffusive and anomalous diffusive processes, the function ν(*q*) is actually a constant (i.e., ν = 1/2 in the normal diffusive case, ν > 1/2 for superdiffusion and ν = 1 for ballistic propagation). However, characterizing the dichotomy between normal and anomalous diffusion based only on the value of ν(2) fails at providing the full picture [223], and one should look further than |*x*| 2 (*t*). Ideally, the complete information is represented by the Green's function or propagator of the process. Often this is not possible to obtain exactly, though asymptotic methods can be exploited to provide useful approximations in the form of similarity solutions. These are expressed, in the limit of *t* → ∞ by the self-similar collapse relation

$$p(\mathbf{x}, t) \approx t^{-1/\nu} \mathcal{P} \left( \frac{\mathbf{x}}{t^{1/\nu}} \right). \tag{6.11}$$

In the case of normal diffusion ν = 2 and P is a Gaussian, while in the superdiffusive case there are examples in which P is a Lévy density [224]. The approximation (6.11) is valid only in the central scaling region (CSR), i.e., excluding the non-scaling tails of the distribution which are populated by 'particles' that have undergone exceptionally large displacements. The concept of *strong* self-similarity can be introduced to identify those systems where ν(*q*) = 1/ν , ∀*q*, as it is the case for normal diffusion and Lévy flights. On the contrary, *weak* self-similarity refers to propagation dynamics giving rise to non-trivial ν(*q*) functions. Typically, the small-*q* part of ν(*q*) (which passes through the origin of the (ν(*q*), *q*) plane) refers to the central scaling region of the asymptotic profile containing the majority of the 'particles', while the large-*q* range is determined by relatively few tail individuals which have traveled exceptionally far. For this reason, as confirmed by analytical, numerical and experimental evidence, in most known systems ν(*q*)*q* is a piece-

0 1 2 3

0 1 2 3 ln(*t*/τ) (b) *n* = 1.016

Figure 6.11. Time evolution of moments *q* ∈ [0.2, 6.2] calculated for the transmitted profiles of the optically thin samples (OT = 1) with refractive index contrast values of (a) *n* = 1, (b) *n* = 1.016 and

nature of light transport in the turbid slab, with all exponents lying on the analytic diffusive

which simulation noise starts to be appreciable (Figure 6.11). In particular, when plotting the full behavior of *q*ν(*q*) we can actually appreciate that none of the sample is truly diffusive, and that slight deviations in the linearity that are barely appreciable by considering just the second moment (Figure 6.1a) actually fit into a weakly superdiffusive pattern, as can be appreciated in Figure 6.12a. This result is quite surprising given that the velocities of packets in the simulation is constant and the variance of the SLD is finite and is probably due to the time-dependent nature of the SLD. We have fitted the fractional exponents *q* ∈ [0.2, 1.2] of Figure 6.12a to characterize the scaling properties of the CSR, obtaining a good linear agreement. Notably, even in the small-*q* range a small superdiffusive behavior is to be found for all refractive index contrasts, with ν(*q* ≤ 1.2) = νCSR = 0.5605(3), 0.5644(3)

By analyzing the moments in terms of their self-similarity, however, differences between the configurations emerge. In Figure 6.12b, we plot the relative deviation between the line *q*νCSR(*n*) extrapolated to higher values of *q* and the simulated data, highlighting a clear qualitative difference arising in the presence of a refractive-index contrast. In particular, while the index-matched configuration seems to exhibit a strongly self-similar (and superdiffusive) behavior, the other configurations show linearly increasing relative deviations. This is not only a hallmark of weak self-similarity, but is also signaling that the commonly studied piece-wise linear type of weak self-similarity does not apply to this case.

113

exponents have been retrieved in a range of delays between ln 16 ≤ ln *t*

and 0.5461(7) respectively for *n* = 1, 1.016 and 1.1.

The situation is more interesting when we turn to our optically thin configurations, which we have simulated with further increased statistics. In this case, plotting the expansion of the profiles for different widths reveals that the asymptotic regime is reached at even later times, as can also be appreciated by the slow convergence of the step-length distributions to their asymptotic form (see, for instance, Figure 6.2). To this purpose, the power-law

0 1 2 3 ln(*t*/τ)

τ is performed in the range between 8 < *t*/τ < 16 where

(c) *n* = 1.1

τ ≤ ln 20, beyond

*p* = 0.2 *p* = 0.4 *p* = 0.6 *p* = 0.8 *p* = 1.0 *p* = 1.2 *p* = 1.4 *p* = 1.6 *p* = 1.8 *p* = 2.0 *p* = 2.3 *p* = 2.6 *p* = 2.9 *p* = 3.2 *p* = 3.5 *p* = 3.8 *p* = 4.1 *p* = 4.4 *p* = 4.7 *p* = 5.0 *p* = 5.3 *p* = 5.6 *p* = 5.9 *p* = 6.2

ln(*t*/τ)

(c) *n* = 1.1. A linear fit of ln ρ*q* versus ln *t*

the asymptotic power-law growth is established.

(a) *n* = 1

*q*/2.

0

line ρ*q* = *t*

20

ln ρ*q* 40

60

Figure 6.10. (a) Time evolution of moments *q* ∈ [0.2, 6.2] calculated for the transmitted profiles of the *n* = 1.5, OT = 10 sample. A linear fit of ln ρ*q* versus ln *t* τ is performed in the range between 8 < *t*/τ < 16 where the asymptotic power-law growth is established. (b) The obtained exponents reveal a strongly self-similar and diffusive regime, as expected for a turbid, homogeneous medium. Ballistic (*q*ν(*q*) = 1) and diffusive (*q*ν(*q*) = 1/2) lines are drawn as guides to the eye.

wise linear function experiencing a change of slope at a certain value of *q* = *q*c [222, 223, 225, 226]. This is the typical behavior exhibited, for instance, by Lévy walks [227, 228]. For most systems studied in the literature, the large-*q* branch is therefore characterized by a ballistic growth (ν(*q*) = 1) which sets in for those high moments that are eventually dominated by the few particles that propagate ballistically. On the contrary, the presence of ballistically expanding tails does not necessarily imply any form of weak self-similarity. The simple telegrapher model is one such example that is both normally diffusive and strongly self-similar: its tail structure is not self-similar at all times, but its contribution to the profile is mild enough that as *t* → ∞ all moments are still determined by the CSR. As a final note, the value of *q*c at which the slope of weakly self-similar systems can change depending on the parameters of the transport process involved. If this occurs before *q* = 2, measuring the expansion rate of the MSW would actually probe the non-scaling tails rather than the density in the CSR. In other words, a system might still exhibit a perfectly linear 2ν(2) = 1 while being both weakly self-similar and anomalously diffusive. Therefore, a full characterization beyond the simple second moment is always recommended to correctly identify the nature of transport.

In the following, we present a preliminary discussion of the numerical results that we obtained based on the previous simulations on slab geometry systems, considering the asymptotic behavior for different moments of the transmitted profiles using |*x*| = ρ. A first, relevant test case is that of the OT = 10 slab that we used as an analogue of our experimental sample of Section 4.2. As we have seen, this system is characterized by an asymptotic transport regime that is slightly anisotropic, but still diffusive. This is confirmed by the late-*t* analysis of the moments as shown in Figure 6.10. In the analysis, the power-law exponents for each moment can be more accurately obtained by performing a linear fit of ln ρ*q* versus ln *t* τ , considering only late times between ln 8 ≤ ln *t* τ ≤ ln 16 where statistical noise is negligible and transport has already reached its asymptotic regime. As expected, analysis of the fractional moments confirms both the diffusive and the strongly self-similar

Figure 6.11. Time evolution of moments *q* ∈ [0.2, 6.2] calculated for the transmitted profiles of the optically thin samples (OT = 1) with refractive index contrast values of (a) *n* = 1, (b) *n* = 1.016 and (c) *n* = 1.1. A linear fit of ln ρ*q* versus ln *t* τ is performed in the range between 8 < *t*/τ < 16 where the asymptotic power-law growth is established.

0 1 2 3

the *n* = 1.5, OT = 10 sample. A linear fit of ln ρ*q* versus ln *t*

*p* = 0.2 *p* = 0.4 *p* = 0.6 *p* = 0.8 *p* = 1.0 *p* = 1.2 *p* = 1.4 *p* = 1.6 *p* = 1.8 *p* = 2.0 *p* = 2.3 *p* = 2.6 *p* = 2.9 *p* = 3.2 *p* = 3.5 *p* = 3.8 *p* = 4.1 *p* = 4.4 *p* = 4.7 *p* = 5.0 *p* = 5.3 *p* = 5.6 *p* = 5.9 *p* = 6.2

Figure 6.10. (a) Time evolution of moments *q* ∈ [0.2, 6.2] calculated for the transmitted profiles of

8 < *t*/τ < 16 where the asymptotic power-law growth is established. (b) The obtained exponents reveal a strongly self-similar and diffusive regime, as expected for a turbid, homogeneous medium.

wise linear function experiencing a change of slope at a certain value of *q* = *q*c [222, 223, 225, 226]. This is the typical behavior exhibited, for instance, by Lévy walks [227, 228]. For most systems studied in the literature, the large-*q* branch is therefore characterized by a ballistic growth (ν(*q*) = 1) which sets in for those high moments that are eventually dominated by the few particles that propagate ballistically. On the contrary, the presence of ballistically expanding tails does not necessarily imply any form of weak self-similarity. The simple telegrapher model is one such example that is both normally diffusive and strongly self-similar: its tail structure is not self-similar at all times, but its contribution to the profile is mild enough that as *t* → ∞ all moments are still determined by the CSR. As a final note, the value of *q*c at which the slope of weakly self-similar systems can change depending on the parameters of the transport process involved. If this occurs before *q* = 2, measuring the expansion rate of the MSW would actually probe the non-scaling tails rather than the density in the CSR. In other words, a system might still exhibit a perfectly linear 2ν(2) = 1 while being both weakly self-similar and anomalously diffusive. Therefore, a full characterization beyond the simple second moment is always recommended to correctly

In the following, we present a preliminary discussion of the numerical results that we obtained based on the previous simulations on slab geometry systems, considering the asymptotic behavior for different moments of the transmitted profiles using |*x*| = ρ. A first, relevant test case is that of the OT = 10 slab that we used as an analogue of our experimental sample of Section 4.2. As we have seen, this system is characterized by an asymptotic transport regime that is slightly anisotropic, but still diffusive. This is confirmed by the late-*t* analysis of the moments as shown in Figure 6.10. In the analysis, the power-law exponents for each moment can be more accurately obtained by performing a linear fit of

τ , considering only late times between ln 8 ≤ ln *t*

noise is negligible and transport has already reached its asymptotic regime. As expected, analysis of the fractional moments confirms both the diffusive and the strongly self-similar

112

Ballistic (*q*ν(*q*) = 1) and diffusive (*q*ν(*q*) = 1/2) lines are drawn as guides to the eye.

0

1

2

*q*ν(*q*)

ballistic

*n* = 1.5, OT = 10

3

4

0 2 4 6

(b) *n* = 1.5, OT = 10

*q*

τ is performed in the range between

τ ≤ ln 16 where statistical

diffusion

ln(*t*/τ)

identify the nature of transport.

ln ρ*q* versus ln *t*

(a) *n* = 1.5, OT = 10

0

20

ln ρ*q* 40

nature of light transport in the turbid slab, with all exponents lying on the analytic diffusive line ρ*q* = *t q*/2.

The situation is more interesting when we turn to our optically thin configurations, which we have simulated with further increased statistics. In this case, plotting the expansion of the profiles for different widths reveals that the asymptotic regime is reached at even later times, as can also be appreciated by the slow convergence of the step-length distributions to their asymptotic form (see, for instance, Figure 6.2). To this purpose, the power-law exponents have been retrieved in a range of delays between ln 16 ≤ ln *t* τ ≤ ln 20, beyond which simulation noise starts to be appreciable (Figure 6.11). In particular, when plotting the full behavior of *q*ν(*q*) we can actually appreciate that none of the sample is truly diffusive, and that slight deviations in the linearity that are barely appreciable by considering just the second moment (Figure 6.1a) actually fit into a weakly superdiffusive pattern, as can be appreciated in Figure 6.12a. This result is quite surprising given that the velocities of packets in the simulation is constant and the variance of the SLD is finite and is probably due to the time-dependent nature of the SLD. We have fitted the fractional exponents *q* ∈ [0.2, 1.2] of Figure 6.12a to characterize the scaling properties of the CSR, obtaining a good linear agreement. Notably, even in the small-*q* range a small superdiffusive behavior is to be found for all refractive index contrasts, with ν(*q* ≤ 1.2) = νCSR = 0.5605(3), 0.5644(3) and 0.5461(7) respectively for *n* = 1, 1.016 and 1.1.

By analyzing the moments in terms of their self-similarity, however, differences between the configurations emerge. In Figure 6.12b, we plot the relative deviation between the line *q*νCSR(*n*) extrapolated to higher values of *q* and the simulated data, highlighting a clear qualitative difference arising in the presence of a refractive-index contrast. In particular, while the index-matched configuration seems to exhibit a strongly self-similar (and superdiffusive) behavior, the other configurations show linearly increasing relative deviations. This is not only a hallmark of weak self-similarity, but is also signaling that the commonly studied piece-wise linear type of weak self-similarity does not apply to this case.

asymptotic transport regime exhibited in optically thin systems represents a challenging

115

Now that the study of transport phenomena has grown to such a mature research field spanning over many different branches of science, it is perhaps surprising that such a rich transport physics, comprising important and complex concepts as anisotropy, weakly self-similar scaling and anomalous diffusion, can all be found and studied in such a simple and explored model as the single plane-parallel homogeneous slab with isotropic, annealed

issue, even from the point of view of numerical simulations.

disorder.

Figure 6.12. (a) Fractional exponents obtained by the linear fits of Figure 6.11, revealing a superdiffusive characteristic of transport in the optically-thin, low-contrast slab configuration. Error bars (not shown), calculated as the 95 % confidence intervals returned by a linear least square method, are approximately the same size of the symbols. A linear fit in the small-*q* range [0.2, 1.2] reveals a linear superdiffusive behavior with ν(*q* ≤ 1.2) = νCSR = 0.5605(3), 0.5644(3) and 0.5461(7) respectively for *n* = 1, 1.016 and 1.1. Fitted lines are extrapolated to higher moments as guides to the eye. (b) Relative residuals between the values of *q*ν(*q*) obtained from the simulations and the extrapolated linear scaling νCSR(*n*), highlighting the presence of weakly self-similar anomalous diffusion in the presence of a mismatch refractive index.

On the contrary, a second order dependence on *q* seems also to be taken into account.

As we discussed in the introduction, the piece-wise model of weakly self-similar diffusion refers just to the simple case where the ballistic tails of the propagator eventually represent the prevailing contribution in the determination of higher moments. However, this regime is unlikely to be relevant for our samples. The transmitted profiles recorded at 16 ≤ *t*/τ ≤ 20 certainly do not exhibit any ballistic peak on the tails, and basically all energy packets transmitted within that time range have experienced at least a few scattering events, given the large but finite magnitude of our simulations. This might explain why the piecewise linear model fails at describing this configuration, which apparently requires the need of a superlinear term. As a matter of fact, it is believed that weak self-similarity might occur in several forms, of which the piece-wise linear model is only one simple case. Nonetheless, the literature on general, non piece-wise linear examples of weakly self-similar diffusion is extremely limited, and focused on very different models with no losses, power-law SLDs and propagation velocities drawn from a distribution with a time-growing variance [229]. It is possible that the time-varying characteristic (which, in our case, regards the SLD rather than the velocities) is the key element giving rise to more general weakly self-similar dynamics. In this respect, light transport in a simple homogeneous, isotropic scattering slab might represent an interesting physical platform for the experimental study of a broader array of transport regimes than expected. To this end, an optimal trade-off should be sought after between the refractive index contrast, the absolute and the relative optical thickness of the sample, in order to design an experiment where a weakly self-similar transport regime is reached compatibly with the instrumental sensitivity and quantitative accuracy available for spatio-temporal imaging techniques. In fact, as of now, the slow convergence to the

asymptotic transport regime exhibited in optically thin systems represents a challenging issue, even from the point of view of numerical simulations.

Now that the study of transport phenomena has grown to such a mature research field spanning over many different branches of science, it is perhaps surprising that such a rich transport physics, comprising important and complex concepts as anisotropy, weakly self-similar scaling and anomalous diffusion, can all be found and studied in such a simple and explored model as the single plane-parallel homogeneous slab with isotropic, annealed disorder.

0 2 4 6

(a)

presence of a mismatch refractive index.

*q*

diffusion

0 2 4 6

(b)

*q*

1

1.02

ν(*q*)/νCSR

Figure 6.12. (a) Fractional exponents obtained by the linear fits of Figure 6.11, revealing a superdiffusive characteristic of transport in the optically-thin, low-contrast slab configuration. Error bars (not shown), calculated as the 95 % confidence intervals returned by a linear least square method, are approximately the same size of the symbols. A linear fit in the small-*q* range [0.2, 1.2] reveals a linear superdiffusive behavior with ν(*q* ≤ 1.2) = νCSR = 0.5605(3), 0.5644(3) and 0.5461(7) respectively for *n* = 1, 1.016 and 1.1. Fitted lines are extrapolated to higher moments as guides to the eye. (b) Relative residuals between the values of *q*ν(*q*) obtained from the simulations and the extrapolated linear scaling νCSR(*n*), highlighting the presence of weakly self-similar anomalous diffusion in the

On the contrary, a second order dependence on *q* seems also to be taken into account.

As we discussed in the introduction, the piece-wise model of weakly self-similar diffusion refers just to the simple case where the ballistic tails of the propagator eventually represent the prevailing contribution in the determination of higher moments. However, this regime is unlikely to be relevant for our samples. The transmitted profiles recorded at 16 ≤ *t*/τ ≤ 20 certainly do not exhibit any ballistic peak on the tails, and basically all energy packets transmitted within that time range have experienced at least a few scattering events, given the large but finite magnitude of our simulations. This might explain why the piecewise linear model fails at describing this configuration, which apparently requires the need of a superlinear term. As a matter of fact, it is believed that weak self-similarity might occur in several forms, of which the piece-wise linear model is only one simple case. Nonetheless, the literature on general, non piece-wise linear examples of weakly self-similar diffusion is extremely limited, and focused on very different models with no losses, power-law SLDs and propagation velocities drawn from a distribution with a time-growing variance [229]. It is possible that the time-varying characteristic (which, in our case, regards the SLD rather than the velocities) is the key element giving rise to more general weakly self-similar dynamics. In this respect, light transport in a simple homogeneous, isotropic scattering slab might represent an interesting physical platform for the experimental study of a broader array of transport regimes than expected. To this end, an optimal trade-off should be sought after between the refractive index contrast, the absolute and the relative optical thickness of the sample, in order to design an experiment where a weakly self-similar transport regime is reached compatibly with the instrumental sensitivity and quantitative accuracy available for spatio-temporal imaging techniques. In fact, as of now, the slow convergence to the

114

1.04

1.06

*n* = 1 *n* = 1.016 *n* = 1.1

0

1

2

*q*ν(*q*)

ballistic

3

4

*n* = 1 *n* = 1.016 *n* = 1.1

Bibliography

549–560 (1905).

*Tappi* 32, 498–501 (1949).

University Press, 2002).

Business Media, 2005).

Bellingham, 2010).

547–549 (2003).

2172 (2008).

*Photonics* 7, 197–204 (2013).

(2008).

(SPIE press Bellingham, 2007).

*heterogeneous Earth* (Springer, 2012).

*Atmospheric Sciences* 31, 1137–1160 (1974).

mobility patterns. *Nature* 453, 779–782 (2008).

1. Einstein, A. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. *Annalen der Physik* 322,

2. Gonzalez, M. C., Hidalgo, C. A. & Barabasi, A.-L. Understanding individual human

3. Van den Akker, J. Scattering and absorption of light in paper and other diffusing media.

4. Thomas, G. E. & Stamnes, K. *Radiative transfer in the atmosphere and ocean* (Cambridge

5. Marshak, A. & Davis, A. *3D radiative transfer in cloudy atmospheres* (Springer Science &

6. Tuchin, V. *Tissue optics: light scattering methods and instruments for medical diagnosis*

7. Martelli, F., Del Bianco, S., Ismaelli, A. & Zaccanti, G. *Light Propagation through Biological Tissue and Other Di*ff*usive Media: Theory, Solutions and Software* (SPIE Press,

8. Sato, H., Fehler, M. C. & Maeda, T. *Seismic wave propagation and scattering in the*

9. Campillo, M. & Paul, A. Long-range correlations in the diffuse seismic coda. *Science* 299,

10. Hansen, J. E. & Hovenier, J. Interpretation of the polarization of Venus. *Journal of the*

11. Porco, C. C. *et al.* Simulations of the dynamical and light-scattering behavior of Saturn's rings and the derivation of ring particle and disk properties. *The Astronomical Journal* 136,

12. Barthelemy, P., Bertolotti, J. & Wiersma, D. S. A Lévy flight for light. *Nature* 453, 495–498

13. Segev, M., Silberberg, Y. & Christodoulides, D. N. Anderson localization of light. *Nature*

14. Caruso, F., Crespi, A., Ciriolo, A. G., Sciarrino, F. & Osellame, R. Fast escape of a quantum walker from an integrated photonic maze. *Nature Communications* 7 (2016).

117

# Bibliography

117

15. Froufe-Pérez, L. S. *et al.* Role of Short-Range Order and Hyperuniformity in the Formation of Band Gaps in Disordered Photonic Materials. *Physical Review Letters* 117, 053902 (2016).

34. Popoff, S., Goetschy, A., Liew, S., Stone, A. D. & Cao, H. Coherent control of total transmission of light through disordered media. *Physical Review Letters* 112, 133903

35. Redding, B., Liew, S. F., Sarma, R. & Cao, H. Compact spectrometer based on a disordered

36. Goorden, S. A., Horstmann, M., Mosk, A. P., Škoric, B. & Pinkse, P. W. Quantum-secure ´

37. Van Albada, M. P. & Lagendijk, A. Observation of weak localization of light in a random

38. Wiersma, D. S., Bartolini, P., Lagendijk, A. & Righini, R. Localization of light in a

39. Prum, R. O., Dufresne, E. R., Quinn, T. & Waters, K. Development of colour-producing fi-keratin nanostructures in avian feather barbs. *Journal of The Royal Society Interface* 6,

40. Muller, N., Haberko, J., Marichy, C. & Scheffold, F. Silicon hyperuniform disordered photonic materials with a pronounced gap in the shortwave infrared. *Advanced Optical*

41. Leseur, O., Pierrat, R. & Carminati, R. High-density hyperuniform materials can be trans-

42. Badon, A. *et al.* Spatio-temporal imaging of light transport in highly scattering media under

43. Wyman, D. R., Patterson, M. S. & Wilson, B. C. Similarity relations for the interaction

44. Contini, D., Martelli, F. & Zaccanti, G. Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory. *Applied Optics* 36, 4587–4599

45. Vynck, K., Pierrat, R. & Carminati, R. Polarization and spatial coherence of electromagnetic waves in uncorrelated disordered media. *Physical Review A* 89, 013842 (2014). 46. Roux, L., Mareschal, P., Vukadinovic, N., Thibaud, J.-B. & Greffet, J.-J. Scattering by a slab containing randomly located cylinders: comparison between radiative transfer and

47. Fraden, S. & Maret, G. Multiple light scattering from concentrated, interacting suspensions.

48. Bressel, L. & Reich, O. Theoretical and experimental study of the diffuse transmission of light through highly concentrated absorbing and scattering materials: Part I: Monte-Carlo simulations. *Journal of Quantitative Spectroscopy and Radiative Transfer* 146, 190–198

119

parameters in radiation transport. *Applied Optics* 28, 5243–5249 (1989).

authentication of a physical unclonable key. *Optica* 1, 421–424 (2014).

photonic chip. *Nature Photonics* 7, 746–751 (2013).

medium. *Physical Review Letters* 55, 2692 (1985).

disordered medium. *Nature* 390, 671–673 (1997).

white light illumination. *Optica* 3, 1160–1166 (2016).

electromagnetic simulation. *JOSA A* 18, 374–384 (2001).

49. Modest, M. F. *Radiative heat transfer* (Academic Press, 2013).

*Physical Review Letters* 65, 512 (1990).

(2014).

S253–S265 (2009).

(1997).

(2014).

*Materials* 2, 115–119 (2014).

parent. *Optica* 3, 763–767 (2016).

15. Froufe-Pérez, L. S. *et al.* Role of Short-Range Order and Hyperuniformity in the Formation of Band Gaps in Disordered Photonic Materials. *Physical Review Letters* 117, 053902

16. Pauly, M., Kollig, T. & Keller, A. *Metropolis light transport for participating media*

17. Schmidt, T.-W., Pellacini, F., Nowrouzezahrai, D., Jarosz, W. & Dachsbacher, C. *State of the art in artistic editing of appearance, lighting and material* in *Computer Graphics*

18. Planck, M. *The theory of heat radiation* (P. Blakiston's Son & Co., Philadelphia, 1914).

20. Schuster, A. Radiation through a foggy atmosphere. *The Astrophysical Journal* 21, 1 (1905). 21. Mishchenko, M. I. Vector radiative transfer equation for arbitrarily shaped and arbitrarily oriented particles: a microphysical derivation from statistical electromagnetics. *Applied*

22. Ripoll, J. Derivation of the scalar radiative transfer equation from energy conservation of

23. Ishimaru, A. *Wave propagation and scattering in random media* (Academic press New

24. Born, M. & Wolf, E. *Principles of optics: electromagnetic theory of propagation, interfer-*

25. Van de Hulst, H. C. *Multiple light scattering: tables, formulas, and applications* (Elsevier,

27. Henyey, L. G. & Greenstein, J. L. Diffuse radiation in the galaxy. *The Astrophysical Journal*

28. Graaff, R., Aarnoudse, J. G., de Mul, F. F. & Jentink, H. W. Similarity relations for anisotropic scattering in absorbing media. *Optical Engineering* 32, 244–252 (1993). 29. Bertolotti, J. *et al.* Engineering disorder in superdiffusive Levy glasses. *Advanced Func-*

30. Mackowski, D. W. & Mishchenko, M. I. Calculation of the T matrix and the scattering

31. Feinberg, J. Self-pumped, continuous-wave phase conjugator using internal reflection.

32. Yaqoob, Z., Psaltis, D., Feld, M. S. & Yang, C. Optical phase conjugation for turbidity

33. Vellekoop, I., Lagendijk, A. & Mosk, A. Exploiting disorder for perfect focusing. *Nature*

118

suppression in biological samples. *Nature Photonics* 2, 110–115 (2008).

matrix for ensembles of spheres. *JOSA A* 13, 2266–2278 (1996).

19. Chandrasekhar, S. *Radiative transfer* (Oxford, Clarendon Press, 1950).

Maxwell's equations in the far field. *JOSA A* 28, 1765–1775 (2011).

*ence and di*ff*raction of light* (CUP Archive, 2000).

(2016).

(Springer, 2000).

*Forum* (2015).

York, 1978).

2012).

93, 70–83 (1941).

*tional Materials* 20, 965–968 (2010).

*Optics Letters* 7, 486–488 (1982).

*Photonics* 4, 320–322 (2010).

*Optics* 41, 7114–7134 (2002).

50. Machida, M., Panasyuk, G. Y., Schotland, J. C. & Markel, V. A. The Green's function for the radiative transport equation in the slab geometry. *Journal of Physics A: Mathematical and Theoretical* 43, 065402 (2010).

66. Martelli, F., Sassaroli, A., Zaccanti, G. & Yamada, Y. Properties of the light emerging from a diffusive medium: angular dependence and flux at the external boundary. *Physics in*

67. MacKintosh, F. & John, S. Diffusing-wave spectroscopy and multiple scattering of light in

68. Yoo, K., Liu, F. & Alfano, R. When does the diffusion approximation fail to describe photon

69. Yoo, K. & Alfano, R. Time-resolved coherent and incoherent components of forward light

70. Lemieux, P.-A., Vera, M. & Durian, D. Diffusing-light spectroscopies beyond the diffusion limit: The role of ballistic transport and anisotropic scattering. *Physical Review E* 57, 4498

71. Elaloufi, R., Carminati, R. & Greffet, J.-J. Time-dependent transport through scattering media: from radiative transfer to diffusion. *Journal of Optics A: Pure and Applied Optics* 4,

72. Elaloufi, R., Carminati, R. & Greffet, J.-J. Diffusive-to-ballistic transition in dynamic light transmission through thin scattering slabs: a radiative transfer approach. *JOSA A* 21, 1430–

73. Kienle, A. & Patterson, M. S. Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium. *JOSA A* 14, 246–254

74. Venugopalan, V., You, J. & Tromberg, B. Radiative transport in the diffusion approximation: an extension for highly absorbing media and small source-detector separations. *Physical*

75. Chen, N. G. & Bai, J. Monte Carlo approach to modeling of boundary conditions for the

76. Popescu, G., Mujat, C. & Dogariu, A. Evidence of scattering anisotropy effects on boundary

77. Garofalakis, A. *et al.* Characterization of the reduced scattering coefficient for optically thin samples: theory and experiments. *Journal of Optics A: Pure and Applied Optics* 6, 725

78. Farrell, T. J., Patterson, M. S. & Wilson, B. A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties

79. Martelli, F., Bassani, M., Alianelli, L., Zangheri, L. & Zaccanti, G. Accuracy of the diffusion equation to describe photon migration through an infinite medium: numerical and

experimental investigation. *Physics in Medicine and Biology* 45, 1359 (2000). 80. Jacques, S. L. Optical properties of biological tissues: a review. *Physics in Medicine and*

121

conditions of the diffusion equation. *Physical Review E* 61, 4523 (2000).

diffusion equation. *Physical Review Letters* 80, 5321 (1998).

in vivo. *Medical Physics* 19, 879–888 (1992).

correlated random media. *Physical Review B* 40, 2383 (1989).

scattering in random media. *Optics Letters* 15, 320–322 (1990).

transport in random media? *Physical Review Letters* 64, 2647 (1990).

*Medicine and Biology* 44, 1257 (1999).

(1998).

S103 (2002).

1437 (2004).

*Review E* 58, 2395 (1998).

*Biology* 58, R37 (2013).

(1997).

(2004).

50. Machida, M., Panasyuk, G. Y., Schotland, J. C. & Markel, V. A. The Green's function for the radiative transport equation in the slab geometry. *Journal of Physics A: Mathematical*

51. Liemert, A. & Kienle, A. Exact and efficient solution of the radiative transport equation for

52. Wilson, B. & Adam, G. A Monte Carlo model for the absorption and flux distributions of

53. Wang, L., Jacques, S. L. & Zheng, L. MCML—Monte Carlo modeling of light transport in multi-layered tissues. *Computer Methods and Programs in Biomedicine* 47, 131–146

54. Alerstam, E., Svensson, T. & Andersson-Engels, S. Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration. *Journal of*

55. Fang, Q. Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates.

56. Patterson, M. S., Chance, B. & Wilson, B. C. Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties. *Applied Optics* 28, 2331–2336

57. Bassani, M., Martelli, F., Zaccanti, G. & Contini, D. Independence of the diffusion coefficient from absorption: experimental and numerical evidence. *Optics Letters* 22, 853–855

58. Durduran, T., Yodh, A., Chance, B. & Boas, D. Does the photon-diffusion coefficient

59. Pierrat, R., Greffet, J.-J. & Carminati, R. Photon diffusion coefficient in scattering and

60. Alerstam, E. *Optical spectroscopy of turbid media: time-domain measurements and accel-*

61. Haskell, R. C. *et al.* Boundary conditions for the diffusion equation in radiative transfer.

62. Diebold, M. P. *Application of Light Scattering to Coatings: A User's Guide* (Springer,

63. Comsa, D., Farrell, T. & Patterson, M. Quantification of bioluminescence images of point source objects using diffusion theory models. *Physics in Medicine and Biology* 51, 3733

64. Martelli, F. *et al. Hybrid heuristic time dependent solution of the radiative transfer equation for the slab* in *European Conference on Biomedical Optics* (2009), 7369\_1B.

65. Morse, P. M. & Feshbach, H. *Methods of theoretical physics* (McGraw-Hill New York,

120

*erated monte carlo modelling.* PhD thesis (Lund University, 2011).

*and Theoretical* 43, 065402 (2010).

(1995).

(1989).

(1997).

2014).

(2006).

1953).

the semi-infinite medium. *Scientific Reports* 3 (2013).

light in tissue. *Medical Physics* 10, 824–830 (1983).

*Biomedical Optics* 13, 060504–060504 (2008).

depend on absorption? *JOSA A* 14, 3358–3365 (1997).

absorbing media. *JOSA A* 23, 1106–1110 (2006).

*JOSA A* 11, 2727–2741 (1994).

*Biomedical Optics Express* 1, 165 (2010).

81. Durian, D. Influence of boundary reflection and refraction on diffusive photon transport. *Physical Review E* 50, 857 (1994).

96. Drezek, R., Dunn, A. & Richards-Kortum, R. Light scattering from cells: finite-difference time-domain simulations and goniometric measurements. *Applied Optics* 38, 3651–3661

97. Kienle, A. *et al.* Spatially resolved absolute diffuse reflectance measurements for noninvasive determination of the optical scattering and absorption coefficients of biological tissue.

98. Cuccia, D. J., Bevilacqua, F., Durkin, A. J., Ayers, F. R. & Tromberg, B. J. Quantitation and mapping of tissue optical properties using modulated imaging. *Journal of Biomedical*

99. Scheffold, F. & Block, I. D. Rapid high resolution imaging of diffusive properties in turbid

100. Fishkin, J. B. & Gratton, E. Propagation of photon-density waves in strongly scattering media containing an absorbing semi-infinite plane bounded by a straight edge. *JOSA A* 10,

101. Tromberg, B. J., Svaasand, L. O., Tsay, T.-T. & Haskell, R. C. Properties of photon density

102. Svensson, T. *et al.* Exploiting breakdown of the similarity relation for diffuse light transport: simultaneous retrieval of scattering anisotropy and diffusion constant. *Optics Letters* 38,

103. Stott, J. & Boas, D. A. *A practical comparison between time-domain and frequency-domain di*ff*usive optical imaging systems* in *Biomedical Topical Meetings* 71 (2002), 626–628. 104. Dam, J. S., Yavari, N., Sørensen, S. & Andersson-Engels, S. Real-time absorption and scattering characterization of slab-shaped turbid samples obtained by a combination of angular and spatially resolved measurements. *Applied Optics* 44, 4281–4290 (2005). 105. Alerstam, E. & Svensson, T. Observation of anisotropic diffusion of light in compacted

106. Torricelli, A. *et al.* Time domain functional NIRS imaging for human brain mapping.

107. Svensson, T., Andersson-Engels, S., Einarsdóttír, M. & Svanberg, K. In vivo optical characterization of human prostate tissue using near-infrared time-resolved spectroscopy.

108. Svensson, T., Alerstam, E., Johansson, J. & Andersson-Engels, S. Optical porosimetry and investigations of the porosity experienced by light interacting with porous media. *Optics*

109. D'Andrea, C. *et al.* Time-resolved optical imaging through turbid media using a fast data acquisition system based on a gated CCD camera. *Journal of Physics D: Applied Physics*

110. Sperling, T., Buehrer, W., Aegerter, C. M. & Maret, G. Direct determination of the transition to localization of light in three dimensions. *Nature Photonics* 7, 48–52 (2012).

123

waves in multiple-scattering media. *Applied Optics* 32, 607–616 (1993).

granular porous materials. *Physical Review E* 85, 040301 (2012).

*Journal of Biomedical Optics* 12, 014022–014022 (2007).

(1999).

127–140 (1993).

437–439 (2013).

*Neuroimage* 85, 28–50 (2014).

*Letters* 35, 1740–1742 (2010).

36, 1675 (2003).

*Applied Optics* 35, 2304–2314 (1996).

*Optics* 14, 024012–024012 (2009).

media. *Optics Express* 20, 192–200 (2012).

81. Durian, D. Influence of boundary reflection and refraction on diffusive photon transport.

82. Furutsu, K. & Yamada, Y. Diffusion approximation for a dissipative random medium and

83. Ripoll, J., Yessayan, D., Zacharakis, G. & Ntziachristos, V. Experimental determination of photon propagation in highly absorbing and scattering media. *JOSA A* 22, 546–551 (2005).

84. Zhang, Z. *et al.* Wave transport in random media: The ballistic to diffusive transition.

85. Zhang, X. & Zhang, Z.-Q. Wave transport through thin slabs of random media with internal reflection: Ballistic to diffusive transition. *Physical Review E* 66, 016612 (2002). 86. Yaroshevsky, A., Glasser, Z., Granot, E. & Sternklar, S. Transition from the ballistic to the diffusive regime in a turbid medium. *Optics Letters* 36, 1395–1397 (2011).

87. Kop, R. H., de Vries, P., Sprik, R. & Lagendijk, A. Observation of anomalous transport of strongly multiple scattered light in thin disordered slabs. *Physical Review Letters* 79, 4369

88. Swartling, J. *Biomedical and atmospheric applications of optical spectroscopy in scattering*

89. Grabtchak, S. & Whelan, W. M. Separation of absorption and scattering properties of turbid media using relative spectrally resolved cw radiance measurements. *Biomedical Optics*

90. Groenhuis, R., Ferwerda, H. A. & Ten Bosch, J. Scattering and absorption of turbid materials determined from reflection measurements. 1: Theory. *Applied Optics* 22, 2456–

91. Hammer, M., Roggan, A., Schweitzer, D. & Muller, G. Optical properties of ocular fundus tissues-an in vitro study using the double-integrating-sphere technique and inverse Monte

92. Pickering, J. W. *et al.* Double-integrating-sphere system for measuring the optical properties

93. Dam, J. S., Dalgaard, T., Fabricius, P. E. & Andersson-Engels, S. Multiple polynomial regression method for determination of biomedical optical properties from integrating

94. Dimofte, A., Finlay, J. C. & Zhu, T. C. A method for determination of the absorption and scattering properties interstitially in turbid media. *Physics in Medicine and Biology* 50,

95. Hall, G., Jacques, S. L., Eliceiri, K. W. & Campagnola, P. J. Goniometric measurements of thick tissue using Monte Carlo simulations to obtain the single scattering anisotropy

122

*media.* PhD thesis (Department of Physics, Lund University, 2002).

Carlo simulation. *Physics in Medicine and Biology* 40, 963 (1995).

sphere measurements. *Applied Optics* 39, 1202–1209 (2000).

coefficient. *Biomedical Optics Express* 3, 2707–2719 (2012).

of tissue. *Applied Optics* 32, 399–410 (1993).

*Physical Review E* 50, 857 (1994).

*Physical Review E* 60, 4843 (1999).

*Express* 3, 2371–2380 (2012).

(1997).

2462 (1983).

2291 (2005).

the applications. *Physical Review E* 50, 3634 (1994).

111. Delpy, D. T. *et al.* Estimation of optical pathlength through tissue from direct time of flight measurement. *Physics in Medicine and Biology* 33, 1433 (1988).

127. Haus, H. A. Mode-locking of lasers. *IEEE Journal of Selected Topics in Quantum Electron-*

129. Savo, R. *Transport in complex heterogeneous photonic structures.* PhD thesis (Università

130. Trebino, R. *Frequency-Resolved Optical Gating: The Measurement of Ultrashort Laser*

131. Burresi, M. *et al.* Bright-White Beetle Scales Optimise Multiple Scattering of Light.

132. Potenza, M. A. *et al.* Three dimensional imaging of short pulses. *Optics Communications*

133. Firester, A. Image upconversion: part III. *Journal of Applied Physics* 41, 703–709 (1970). 134. Chiou, W. Geometric optics theory of parametric image upconversion. *Journal of Applied*

135. Devaux, F. & Lantz, E. Ultrahigh-speed imaging by parametric image amplification. *Optics*

136. Abraham, E. *et al.* Real-time two-dimensional imaging in scattering media by use of a

137. Lantz, E. & Devaux, F. Parametric amplification of images: from time gating to noiseless amplification. *Selected Topics in Quantum Electronics, IEEE Journal of* 14, 635–647

138. Pedersen, C., Karamehmedovic, E., Dam, J. S. & Tidemand-Lichtenberg, P. Enhanced ´ 2D-image upconversion using solid-state lasers. *Optics Express* 17, 20885–20890 (2009).

139. Bassi, A. *et al.* Time-gated optical projection tomography. *Optics Letters* 35, 2732–2734

140. Dam, J. S., Pedersen, C. & Tidemand-Lichtenberg, P. High-resolution two-dimensional image upconversion of incoherent light. *Optics Letters* 35, 3796–3798 (2010). 142. Vos, W. L., Tukker, T. W., Mosk, A. P., Lagendijk, A. & IJzerman, W. L. Broadband mean free path of diffuse light in polydisperse ensembles of scatterers for white light-emitting

144. Alerstam, E., Svensson, T. & Andersson-Engels, S. *CUDAMCML-User manual and imple-*

146. Milsom, P. A ray-optic, Monte Carlo, description of a Gaussian beam waist–applied to

147. Alerstam, E., Andersson-Engels, S. & Svensson, T. Improved accuracy in time-resolved diffuse reflectance spectroscopy. *Optics Express* 16, 10440–10454 (2008).

148. Hu, H., Strybulevych, A., Page, J. H., Skipetrov, S. E. & van Tiggelen, B. A. Localization of ultrasound in a three-dimensional elastic network. *Nature Physics* 4, 945–948 (2008).

125

reverse saturable absorption. *Applied Physics B* 70, 593–599 (2000).

diode lighting. *Applied Optics* 52, 2602–2609 (2013).

femtosecond Cr4+: forsterite laser. *Optics Letters* 25, 929–931 (2000).

*ics* 6, 1173–1185 (2000).

*Pulses* (Springer, 2000).

*Scientific Reports* 4 (2014).

*Physics* 42, 1985–1993 (1971).

*Communications* 118, 25–27 (1995).

229, 381–390 (2004).

(2008).

(2010).

*mentation notes* (2009).

di Firenze, 2013).

111. Delpy, D. T. *et al.* Estimation of optical pathlength through tissue from direct time of flight

112. Hebden, J. C., Kruger, R. A. & Wong, K. Time resolved imaging through a highly scattering

113. Hebden, J. C., Arridge, S. R. & Delpy, D. T. Optical imaging in medicine: I. Experimental

114. Dunsby, C. & French, P. Techniques for depth-resolved imaging through turbid media including coherence-gated imaging. *Journal of Physics D: Applied Physics* 36, R207

115. Shimizu, K., Ishimaru, A., Reynolds, L. & Bruckner, A. P. Backscattering of a picosecond pulse from densely distributed scatterers. *Applied Optics* 18, 3484–3488 (1979). 116. Wang, L., Ho, P., Liu, C., Zhang, G. & Alfano, R. Ballistic 2-D imaging through scattering

117. Duncan, M. D., Mahon, R., Tankersley, L. L. & Reintjes, J. Time-gated imaging through scattering media using stimulated Raman amplification. *Optics Letters* 16, 1868–1870

118. Shah, J. Ultrafast luminescence spectroscopy using sum frequency generation. *IEEE Jour-*

119. Yodh, A., Kaplan, P. & Pine, D. Pulsed diffusing-wave spectroscopy: High resolution

120. Faris, G. W. & Banks, M. Upconverting time gate for imaging through highly scattering

121. Le Tolguenec, G., Lantz, E. & Devaux, F. Imaging through scattering media by parametric amplification of images: study of the resolution and the signal-to-noise ratio. *Applied Optics*

122. Yoo, K., Xing, Q. & Alfano, R. Imaging objects hidden in highly scattering media using femtosecond second-harmonic-generation cross-correlation time gating. *Optics Letters* 16,

123. Wiersma, D. S., Muzzi, A., Colocci, M. & Righini, R. Time-resolved experiments on light diffusion in anisotropic random media. *Physical Review E* 62, 6681 (2000).

124. Trebino, R. Second-harmonic-generation rings and refractive-index measurement in uniax-

125. Zhang, D., Kong, Y. & Zhang, J.-y. Optical parametric properties of 532-nm-pumped betabarium-borate near the infrared absorption edge. *Optics Communications* 184, 485–491

126. Wasylczyk, P. Ultracompact autocorrelator for femtosecond laser pulses. *Review of Scien-*

124

walls using an ultrafast optical Kerr gate. *Science* 253, 769–771 (1991).

through nonlinear optical gating. *Physical Review B* 42, 4744 (1990).

measurement. *Physics in Medicine and Biology* 33, 1433 (1988).

techniques. *Physics in Medicine and Biology* 42, 825 (1997).

medium. *Applied Optics* 30, 788–794 (1991).

*nal of Quantum Electronics* 24, 276–288 (1988).

media. *Optics Letters* 19, 1813–1815 (1994).

ial crystals. *Applied Optics* 20, 2090–2096 (1981).

*tific Instruments* 72, 2221–2223 (2001).

36, 8292–8297 (1997).

1019–1021 (1991).

(2000).

(2003).

(1991).

149. Cherroret, N., Skipetrov, S. & Van Tiggelen, B. Transverse confinement of waves in three-dimensional random media. *Physical Review E* 82, 056603 (2010).

166. Kienle, A. & Patterson, M. S. Determination of the optical properties of turbid media from a single Monte Carlo simulation. *Physics in Medicine and Biology* 41, 2221 (1996). 167. Pifferi, A., Taroni, P., Valentini, G. & Andersson-Engels, S. Real-time method for fitting time-resolved reflectance and transmittance measurements with a Monte Carlo model.

168. Thueler, P. *et al.* In vivo endoscopic tissue diagnostics based on spectroscopic absorption, scattering, and phase function properties. *Journal of Biomedical Optics* 8, 495–503 (2003).

169. Palmer, G. M. & Ramanujam, N. Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms. *Applied Optics* 45,

170. Xu, H., Farrell, T. J. & Patterson, M. S. Investigation of light propagation models to determine the optical properties of tissue from interstitial frequency domain fluence mea-

171. Alerstam, E., Andersson-Engels, S. & Svensson, T. White Monte Carlo for time-resolved photon migration. *Journal of Biomedical Optics* 13, 041304–041304 (2008). 172. Martinelli, M. *et al.* Analysis of single Monte Carlo methods for prediction of reflectance

173. Leonardi, L. & Burns, D. H. Quantitative measurements in scattering media: Photon timeof-flight analysis with analytical descriptors. *Applied Spectroscopy* 53, 628–636 (1999). 174. Pfefer, T. J. *et al.* Reflectance-based determination of optical properties in highly attenuating

175. Chauchard, F. *et al.* Least-squares support vector machines modelization for time-resolved

176. Alerstam, E. *et al.* Single-fiber diffuse optical time-of-flight spectroscopy. *Optics Letters*

177. Jäger, M., Foschum, F. & Kienle, A. Application of multiple artificial neural networks for the determination of the optical properties of turbid media. *Journal of Biomedical Optics*

178. Rajaram, N., Nguyen, T. H. & Tunnell, J. W. Lookup table–based inverse model for determining optical properties of turbid media. *Journal of Biomedical Optics* 13, 050501–

179. Nilsson, A. M., Berg, R. & Andersson-Engels, S. Measurements of the optical properties of tissue in conjunction with photodynamic therapy. *Applied Optics* 34, 4609–4619 (1995).

180. Roggan, A., Albrecht, H. J., Doerschel, K., Minet, O. & Mueller, G. J. *Experimental set-up and Monte-Carlo model for the determination of optical tissue properties in the wavelength range 330 to 1100 nm* in *International Symposium on Biomedical Optics Europe'94* (1995),

127

surements. *Journal of Biomedical Optics* 11, 041104–041104 (2006).

from turbid media. *Optics Express* 19, 19627–19642 (2011).

tissue. *Journal of Biomedical Optics* 8, 206–215 (2003).

spectroscopy. *Applied Optics* 44, 7091–7097 (2005).

*Applied Optics* 37, 2774–2780 (1998).

1062–1071 (2006).

37, 2877–2879 (2012).

050501 (2008).

21–36.

18, 057005–057005 (2013).

149. Cherroret, N., Skipetrov, S. & Van Tiggelen, B. Transverse confinement of waves in

150. Bouchard, J.-P. *et al.* Reference optical phantoms for diffuse optical spectroscopy. Part 1–Error analysis of a time resolved transmittance characterization method. *Optics Express*

151. Comelli, D. *et al.* In vivo time-resolved reflectance spectroscopy of the human forehead.

152. Alerstam, E. Anisotropic diffusive transport: Connecting microscopic scattering and macro-

153. Marchesini, R., Bertoni, A., Andreola, S., Melloni, E. & Sichirollo, A. Extinction and absorption coefficients and scattering phase functions of human tissues in vitro. *Applied*

154. Premože, S., Ashikhmin, M., Tessendorf, J., Ramamoorthi, R. & Nayar, S. *Practical rendering of multiple scattering e*ff*ects in participating media* in *Proceedings of the Fifteenth*

155. Nilsson, A. M., Sturesson, C., Liu, D. L. & Andersson-Engels, S. Changes in spectral shape of tissue optical properties in conjunction with laser-induced thermotherapy. *Applied*

156. Mourant, J. R. *et al.* Mechanisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics. *Applied Optics* 37, 3586–3593 (1998). 157. Pattelli, L., Savo, R., Burresi, M. & Wiersma, D. Spatio-temporal visualization of light transport in complex photonic structures. *Light: Science* & *Applications* 5, e16090 (2016).

158. Mazzamuto, G., Pattelli, L., Toninelli, C. & Wiersma, D. S. Deducing effective light transport parameters in optically thin systems. *New Journal of Physics* 18, 023036 (2016).

159. Kienle, A. Anisotropic light diffusion: an oxymoron? *Physical Review Letters* 98, 218104

160. Cortese, L. *et al.* Anisotropic Light Transport in White Beetle Scales. *Advanced Optical*

161. Svensson, T. *et al.* Light diffusion in quenched disorder: Role of step correlations. *Physical*

162. Bruck, R. *et al.* Device-level characterization of the flow of light in integrated photonic circuits using ultrafast photomodulation spectroscopy. *Nature Photonics* 9, 54–60 (2015).

163. Liemert, A. & Kienle, A. Analytical solution of the radiative transfer equation for infinite-

164. Liemert, A. & Kienle, A. Infinite space Green's function of the time-dependent radiative

165. Graaff, R. *et al.* Condensed Monte Carlo simulations for the description of light transport.

126

three-dimensional random media. *Physical Review E* 82, 056603 (2010).

scopic transport properties. *Physical Review E* 89, 063202 (2014).

*Eurographics conference on Rendering Techniques* (2004), 363–374.

18, 11495–11507 (2010).

*Optics* 28, 2318–2324 (1989).

*Optics* 37, 1256–1267 (1998).

*Materials* 3, 1337–1341 (2015).

*Review E* 89, 022141 (2014).

*Applied Optics* 32, 426–434 (1993).

space fluence. *Physical Review A* 83, 015804 (2011).

transfer equation. *Biomedical Optics Express* 3, 543–551 (2012).

(2007).

*Applied Optics* 46, 1717–1725 (2007).

181. Karlsson, H., Fredriksson, I., Larsson, M. & Strömberg, T. Inverse Monte Carlo for estimation of scattering and absorption in liquid optical phantoms. *Optics Express* 20, 12233– 12246 (2012).

198. Ugryumova, N., Matcher, S. J. & Attenburrow, D. P. Measurement of bone mineral density

199. Darling, C. L., Huynh, G. D. & Fried, D. Light scattering properties of natural and artificially demineralized dental enamel at 1310nm. *Journal of Biomedical Optics* 11, 034023–034023

200. Pattelli, L., Mazzamuto, G., Wiersma, D. S. & Toninelli, C. Diffusive light transport in

201. Douglass, K. M. & Dogariu, A. Measuring diffusion coefficients independently of boundary

202. Bolin, F. P., Preuss, L. E., Taylor, R. C. & Ference, R. J. Refractive index of some mammalian tissues using a fiber optic cladding method. *Applied Optics* 28, 2297–2303 (1989).

203. Rivas, J. G., Sprik, R., Lagendijk, A., Noordam, L. & Rella, C. Static and dynamic transport of light close to the Anderson localization transition. *Physical Review E* 63, 046613 (2001).

204. Ramakrishna, S. A. & Kumar, N. Diffusion of particles moving with constant speed.

205. Gopal, V., Anantha Ramakrishna, S., Sood, A. & Kumar, N. Photon transport in thin

206. Spinelli, L. *et al.* Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. Time-resolved method. *Optics Express* 15, 6589–6604 (2007).

207. Chance, B. *et al.* Comparison of time-resolved and -unresolved measurements of deoxyhemoglobin in brain. *Proceedings of the National Academy of Sciences* 85, 4971–4975

208. Wilson, B. C. & Jacques, S. L. Optical reflectance and transmittance of tissues: principles and applications. *IEEE Journal of Quantum Electronics* 26, 2186–2199 (1990). 209. Svensson, T. *Pharmaceutical and biomedical applications of spectroscopy in the photon*

210. Fishkin, J. B., Fantini, S., Gratton, E., *et al.* Gigahertz photon density waves in a turbid

211. Ntziachristos, V. & Chance, B. Accuracy limits in the determination of absolute optical properties using time-resolved NIR spectroscopy. *Medical Physics* 28, 1115–1124 (2001).

212. Scheffold, F., Lenke, R., Tweer, R. & Maret, G. Localization or classical diffusion of light?

213. Wiersma, D. S., Rivas, J. G., Bartolini, P., Lagendijk, A. & Righini, R. Reply: Localization

215. Gómez-Medina, R. *et al.* Negative scattering asymmetry parameter for dipolar particles: Unusual reduction of the transport mean free path and radiation pressure. *Physical Review*

129

medium: theory and experiments. *Physical Review E* 53, 2307 (1996).

disordered slabs. *Pramana-Journal of Physics* 56, 767–778 (2001).

*migration regime.* PhD thesis (Lund University, 2008).

or classical diffusion of light? *Nature* 398, 207 (1999).

*Nature* 398, 206–207 (1999).

*A* 85, 035802 (2012).

via light scattering. *Physics in Medicine and Biology* 49, 469 (2004).

semitransparent media. *Physical Review A* 94, 043846 (2016).

conditions. *Optics Letters* 34, 3379–3381 (2009).

*Physical Review E* 60, 1381–1389 (1999).

(2006).

(1988).

181. Karlsson, H., Fredriksson, I., Larsson, M. & Strömberg, T. Inverse Monte Carlo for estimation of scattering and absorption in liquid optical phantoms. *Optics Express* 20, 12233–

182. Hennessy, R., Lim, S. L., Markey, M. K. & Tunnell, J. W. Monte Carlo lookup table-based inverse model for extracting optical properties from tissue-simulating phantoms using diffuse reflectance spectroscopy. *Journal of Biomedical Optics* 18, 037003–037003 (2013).

183. Simpson, C. R., Kohl, M., Essenpreis, M. & Cope, M. Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using the Monte Carlo inversion

185. Alerstam, E. *et al.* Next-generation acceleration and code optimization for light transport in turbid media using GPUs. *Biomedical Optics Express* 1, 658–675 (2010).

186. Ramella-Roman, J., Prahl, S. & Jacques, S. Three Monte Carlo programs of polarized light transport into scattering media: part I. *Optics Express* 13, 4420–4438 (2005). 187. Fang, Q. & Boas, D. A. Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units. *Optics Express* 17, 20178 (2009).

188. Doronin, A. & Meglinski, I. Online object oriented Monte Carlo computational tool for the

189. Doering, C. R., Ray, T. S. & Glasser, M. L. Long transmission times for transport through

190. Boguñá, M., Porra, J. M. & Masoliver, J. Persistent random walk model for transport

191. Duparré, A. & Kassam, S. Relation between light scattering and the microstructure of

192. Perelman, L. *et al.* Observation of periodic fine structure in reflectance from biological tissue: a new technique for measuring nuclear size distribution. *Physical Review Letters* 80,

193. De Bruin, D. M. *et al.* Optical phantoms of varying geometry based on thin building blocks with controlled optical properties. *Journal of Biomedical Optics* 15, 025001–025001

194. Leonetti, M. & López, C. Measurement of transport mean-free path of light in thin systems.

195. Chan, E. K. *et al.* Effects of compression on soft tissue optical properties. *Selected Topics*

196. Backman, V. *et al.* Polarized light scattering spectroscopy for quantitative measurement of epithelial cellular structures in situ. *Selected Topics in Quantum Electronics, IEEE Journal*

197. Du, Y. *et al.* Optical properties of porcine skin dermis between 900 nm and 1500 nm.

128

needs of biomedical optics. *Biomedical Optics Express* 2, 2461 (2011).

a weakly scattering slab. *Physical Review A* 45, 825 (1992).

through thin slabs. *Physical Review E* 59, 6517 (1999).

optical thin films. *Applied Optics* 32, 5475–5480 (1993).

*in Quantum Electronics, IEEE Journal of* 2, 943–950 (1996).

*Physics in Medicine and Biology* 46, 167 (2001).

*Optics Letters* 36, 2824–2826 (2011).

*of* 5, 1019–1026 (1999).

technique. *Physics in Medicine and Biology* 43, 2465 (1998).

12246 (2012).

627 (1998).

(2010).

216. Wald, A. *A Reprint of" A Method of Estimating Plane Vulnerability Based on Damage of Survivors" by Abraham Wald* (Operations Evaluation Group, Center for Naval Analyses, 1980).

@141. Thorlabs Inc. *UV-Curing Optical Adhesives* http://www.thorlabs.de/newgrouppage9.

@143. Mazzamuto, G. & Pattelli, L. *MCPlusPlus Documentation* http://www.lens.unifi.

@145. Wikipedia. *Logistic distribution — Wikipedia, The Free Encyclopedia* https : / / en .

@184. Alerstam, E. *GPU Monte Carlo* http://www.atomic.physics.lu.se/biophotonics/

@214. Mazzamuto, G. & Pattelli, L. *LUT routine for the single slab geometry* http://www.lens.

131

cfm?objectgroup\_id=196.

it/quantum-nanophotonics/mcplusplus/.

wikipedia.org/wiki/Logistic\_distribution.

research/monte-carlo-simulations/gpu-monte-carlo/.

unifi.it/quantum-nanophotonics/mcplusplus/lut.

## Webpages

216. Wald, A. *A Reprint of" A Method of Estimating Plane Vulnerability Based on Damage of Survivors" by Abraham Wald* (Operations Evaluation Group, Center for Naval Analyses,

217. Mangel, M. & Samaniego, F. J. Abraham Wald's work on aircraft survivability. *Journal of*

218. Svensson, T. *et al.* Holey random walks: Optics of heterogeneous turbid composites.

219. Tosi, A. *et al.* Fast-gated single-photon counting technique widens dynamic range and speeds up acquisition time in time-resolved measurements. *Optics Express* 19, 10735–

220. Cao, H., Xu, J., Seelig, E. & Chang, R. P. Microlaser made of disordered media. *Applied*

221. Wu, X. *et al.* Random lasing in weakly scattering systems. *Physical Review A* 74, 053812

222. Castiglione, P., Mazzino, A., Muratore-Ginanneschi, P. & Vulpiani, A. On strong anomalous

223. Ferrari, R., Manfroi, A. & Young, W. Strongly and weakly self-similar diffusion. *Physica*

224. Zumofen, G. & Klafter, J. Scale-invariant motion in intermittent chaotic systems. *Physical*

225. Gal, N. & Weihs, D. Experimental evidence of strong anomalous diffusion in living cells.

226. Metzler, R., Jeon, J.-H., Cherstvy, A. G. & Barkai, E. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. *Physical Chemistry Chemical Physics* 16, 24128–24164 (2014). 227. Allegrini, P. *et al.* Scaling breakdown: a signature of aging. *Physical Review E* 66, 015101

228. Bernabó, P., Burioni, R., Lepri, S. & Vezzani, A. Anomalous transmission and drifts in one-dimensional Lévy structures. *Chaos, Solitons* & *Fractals* 67, 11–19 (2014). 229. Andersen, K., Castiglione, P., Mazzino, A. & Vulpiani, A. Simple stochastic models showing strong anomalous diffusion. *The European Physical Journal B-Condensed Matter*

@26. Prahl, S. *Mie Scattering calculator* http://omlc.org/calc/mie\_calc.html.

http://ocw.mit.edu. License: Creative Commons by-nc-sa.

@128. Kärtner, F. *Ultrafast Optics* Massachusetts Institute of Technology: MIT OpenCourseWare.

130

diffusion. *Physica D: Nonlinear Phenomena* 134, 75–93 (1999).

*the American Statistical Association* 79, 259–267 (1984).

*Physical Review E* 87, 022120 (2013).

*Physics Letters* 76, 2997–2999 (2000).

*D: Nonlinear Phenomena* 154, 111–137 (2001).

*Review E* 47, 851 (1993).

*Physical Review E* 81, 020903 (2010).

*and Complex Systems* 18, 447–452 (2000).

1980).

10746 (2011).

(2006).

(2002).

Webpages

A. Derivations

A.1. Derivation of the radiative transfer equation

+ d*P*abs d*V*

*S* (*r* − *r* ) + d*P*abs d*V*

 δ*V*

where s is the direction of the flow of energy at *r*, as usual.

where *S* (*r*) is the magnitude of the time-averaged Poynting vector at *r*.

*S*inc · s*j* d*V*

such sample volume can be approximated to the volume-averaged incident flow

averaged expression (2.12) that we derived in section 2.1

∂ *S*(*r*) · s*j* ∂*t*

1 *v*

far-field of the particles.

 δ*V* (s · s*j*) 1 *v* ∂ ∂*t*

assume that

In this appendix we review a possible derivation of the radiative transfer equation (RTE) using Poynting's theorem for energy conservation. Let us start considering the time-

which ensures that energy conservation is rotationally invariant and holds in any particular direction s*j*. This equation is still valid also in the presence of a time-varying dependence (i.e., a modulation of the source) on a time-scale that is slower than the electromagnetic oscillation ω. On the other hand, we must remember that it is based on the assumption that the electric and magnetic fields are mutually orthogonal, and therefore holds only in the

In order to obtain the RTE we will apply energy conservation to a small volume δ*V*

Inside the sample volume δ*V* we can decompose the energy flux *S* as the sum of the contribution of a scattered flux *S*sc and an incident flux *S*inc, where the latter is now a general function accounting for the incoming flux from outside δ*V*. If the volume δ*V* is negligible with respect to the volume *V* − δ*V* δ*V* that is responsible for *S*inc, we can

> δ*V*

which is an analogue to the *dipolar approximation* for multiple scattering of small particles, in that it neglects any self-induction term and assumes that the incident field on each particle is simply given by the sum of all scattered fields excluded its own. Analogously, given the small size of δ*V*, we can make the assumption that the incident field at any point inside

Using these approximations, the first and second term of equation (A.2) can be rewritten

133

+ s*j* · ∇*rS* (*r* − *r*

(s · s*j*) + s*j* · ∇(*S*(*r*) · s*j*) = 0, (A.1)

 ) d3 *r*

*S*sc · s*j* d*V* (A.3)

*S*inc(*r*) · s *S*inc(*r*)v (A.4)

= 0 (A.2)

(cfr. Fig. 2.3) containing *N* particles. Integrating equation (A.1) over δ*V* we obtain

## A. Derivations

## A.1. Derivation of the radiative transfer equation

In this appendix we review a possible derivation of the radiative transfer equation (RTE) using Poynting's theorem for energy conservation. Let us start considering the timeaveraged expression (2.12) that we derived in section 2.1

$$\frac{1}{\nu} \frac{\partial \left< \mathbf{S}(\mathbf{r}) \right> \cdot \mathbf{s}\_j}{\partial t} + \left< \frac{\mathrm{d}P\_{\mathrm{abs}}}{\mathrm{d}V} \right> (\mathbf{s} \cdot \mathbf{s}\_j) + \mathbf{s}\_j \cdot \nabla(\left< \mathbf{S}(\mathbf{r}) \right> \cdot \mathbf{s}\_j) = 0,\tag{A.1}$$

which ensures that energy conservation is rotationally invariant and holds in any particular direction s*j*. This equation is still valid also in the presence of a time-varying dependence (i.e., a modulation of the source) on a time-scale that is slower than the electromagnetic oscillation ω. On the other hand, we must remember that it is based on the assumption that the electric and magnetic fields are mutually orthogonal, and therefore holds only in the far-field of the particles.

In order to obtain the RTE we will apply energy conservation to a small volume δ*V* (cfr. Fig. 2.3) containing *N* particles. Integrating equation (A.1) over δ*V* we obtain

$$\int\_{\delta V} (\mathbf{s} \cdot \mathbf{s}\_j) \left[ \frac{1}{\nu} \frac{\partial}{\partial t} S(\mathbf{r} - \mathbf{r}') + \left\langle \frac{\mathrm{d}P\_{\mathrm{abs}}}{\mathrm{d}V} \right\rangle + \mathbf{s}\_j \cdot \nabla\_{\mathbf{r}'} S(\mathbf{r} - \mathbf{r}') \right] \mathrm{d}^3 r' = 0 \tag{A.2}$$

where *S* (*r*) is the magnitude of the time-averaged Poynting vector at *r*.

Inside the sample volume δ*V* we can decompose the energy flux *S* as the sum of the contribution of a scattered flux *S*sc and an incident flux *S*inc, where the latter is now a general function accounting for the incoming flux from outside δ*V*. If the volume δ*V* is negligible with respect to the volume *V* − δ*V* δ*V* that is responsible for *S*inc, we can assume that

$$
\int\_{\delta V} \langle \mathbf{S}\_{\mathrm{inc}} \rangle \cdot \mathbf{s}\_j \, \mathrm{d}V \gg \int\_{\delta V} \langle \mathbf{S}\_{\mathrm{sc}} \rangle \cdot \mathbf{s}\_j \, \mathrm{d}V \tag{A.3}
$$

which is an analogue to the *dipolar approximation* for multiple scattering of small particles, in that it neglects any self-induction term and assumes that the incident field on each particle is simply given by the sum of all scattered fields excluded its own. Analogously, given the small size of δ*V*, we can make the assumption that the incident field at any point inside such sample volume can be approximated to the volume-averaged incident flow

$$
\langle \mathbf{S}\_{\rm inc}(\mathbf{r}) \rangle \cdot \mathbf{s} \simeq \| \mathbf{S}\_{\rm inc}(\mathbf{r}) \|\_{\rm v} \tag{A.4}
$$

where s is the direction of the flow of energy at *r*, as usual.

Using these approximations, the first and second term of equation (A.2) can be rewritten

133

as

$$\frac{1}{\nu} \frac{\partial}{\partial t} \int\_{\delta V} \langle \mathbf{S} \rangle \cdot \mathbf{s}\_j \, \mathrm{d}V \simeq \frac{1}{\nu} \frac{\partial}{\partial t} ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \nu r\_r(\mathbf{s}\_j) \delta V \,\tag{A.5}$$

equivalent. For each point *r* lying on the surface Σ we can then write that

d*S* σtot

flux traversing Σ will not depend on the position of the particles within δ*V*.

*S*incv*wr*(s*j*)δ*V* + *N*σa*S*incv*wr*(s*j*) + s*j* · ∇*S*incv*wr*(s*j*)δ*V*

*S*incv*wr*(s*j*) + µa*S*incv*wr*(s*j*) + s*j* · ∇*S*incv*wr*(s*j*)

+ *N*σs*S*incv*wr*(s*j*) − *N*σtot

+ µs*S*incv*wr*(s*j*) − µtot

*I*(*r*, *t*,s) = 1

+ s · ∇*I*(*r*, *t*,s) − (µs + µa)*I*(*r*, *t*,s) − µtot

4π

It is convenient to summarize the main approximations that we applied in order to

135

and normalized by δ*V* to remove the dependence from the arbitrary averaging volume

Equation (A.15) is the radiative transport equation for the time-averaged incoming energy flux *S*inc flowing in direction s*j*, averaged over a volume δ*V*. It is important to note that the factor µtot from the last term appears through our definition of the phase function (2.27). It is the product µtot*p*(s,s0) that yields the contribution which is solely due to scattering, i.e. the absorption term related to the phase function that appears through µtot = µs + µa has no

(2.36). Making use of the approximation (A.4) we obtain

where in the last step we have identified d*S*

(A.2), which can be rewritten as

d*S*

*N*

*i*=1

*N*σtot*S*inc(*r*)v

 |*r* − *ri*|

At this point we have obtained the expressions for all the terms appearing in equation

 *N*

 Σ

*S*inc(*r*)v

 Σ *wr*(s*i*)

4π

4π

4π

*I*(*r*, *t*,s )*p*(s,s

*S*incv*wr*(s

*S*incv*wr*(s

 )*p*(s

 )*p*(s

*S*incv*wr*(s), (A.16)

 4π *wr*(s )*p*(s ,s*j*) dΩ ,

*S*sc(*r*)*i* · s*i*

 s*j* · s

*p*(s*i*,s*j*) |*r* − *ri*|

2 d*S*

2 as the solid angle dΩ, since the total

d*S* (A.12)

(A.13)

,s*j*) dΩ (A.14)

,s*j*) dΩ . (A.15)

) dΩ = 0. (A.17)

*i*=1

where *S*sc(*r*)*i* is the energy flux scattered by the *i*-th particle inside δ*V*, given by equation

 Σ *S* (in) sc (*r* − *r* )s*j* · s

 Σ *S* (in) sc (*r* − *r* )s*j* · s

1 *v* ∂ ∂*t*

1 *v* ∂ ∂*t*

1 *v* ∂*I*(*r*, *t*,s) ∂*t*

direct physical meaning.

By defining the specific intensity as

equation (A.15) recovers its well-known appearance

$$
\int\_{\delta V} \left< \frac{\mathrm{d}P\_{\mathrm{abs}}}{\mathrm{d}V} \right> (\mathbf{s} \cdot \mathbf{s}\_j) \, \mathrm{d}V \simeq N \sigma\_{\mathrm{a}} \| \mathbf{S}\_{\mathrm{inc}} \| \_{\mathrm{v}} \omega\_r(\mathbf{s}\_j), \tag{A.6}
$$

in terms of the volume-averaged incident flux, where *wr*(s*j*) is the probability distribution of the energy flux to each angle defined in equation (2.39), *N* is the number of particles in δ*V* and σa is their absorption cross-section (2.23).

On the other hand, the third term of equation (A.2) becomes

$$\begin{split} \int\_{\delta V} (\mathbf{s} \cdot \mathbf{s}\_j) \mathbf{s}\_j \cdot \nabla\_{r'} S \,(\mathbf{r} - \mathbf{r'}) \, \mathrm{d}^3 r' &= \left. \mathrm{w}\_r(\mathbf{s}\_j) \right| \left\| \nabla\_{r'} S\_{\mathrm{inc}} \cdot \mathbf{s}\_j \right\|\_{\mathrm{v}} \delta V \\ &+ \left. \mathrm{w}\_r(\mathbf{s}\_j) \int\_{\delta V} \mathbf{s}\_j \cdot \nabla\_{r'} S\_{\mathrm{sc}} (\mathbf{r} - \mathbf{r'}) \, \mathrm{d}^3 r' \right| \end{split} (\mathrm{A.7})$$

where we have introduced the volume-averaged change in *S* as

$$\left\|\left|\nabla\_{\mathbf{r}'}\mathbf{S}\_{\mathrm{inc}}\cdot\mathbf{s}\_{j}\right\|\right\|\_{\mathrm{V}} = \frac{1}{\delta V} \int\_{\delta V} \mathbf{s}\_{j} \cdot \nabla\_{\mathbf{r}'} \mathbf{S}\_{\mathrm{inc}}(\mathbf{r} - \mathbf{r}') \,\mathrm{d}^3 r' \,. \tag{A.8}$$

At this point, assuming that the far-field approximation is valid within δ*V* (which we already used in order for the cross-sections to be additive in equation (A.6)) and yet that δ*V* is small enough to consider *r r* , we can interchange ∇*r* ↔ ∇*r* and write

$$\left\|\left\|\nabla \mathbf{S}\_{\mathrm{inc}} \cdot \mathbf{s}\_{j}\right\|\right\|\_{\mathrm{V}} \simeq \mathbf{s}\_{j} \cdot \nabla \left\|\left\|\mathbf{S}\_{\mathrm{inc}}\right\|\right\|\_{\mathrm{V}}.\tag{A.9}$$

Plugging this substitution into equation (A.7) we obtain

$$\begin{aligned} \int\_{\delta V} (\mathbf{s} \cdot \mathbf{s}\_j) \mathbf{s}\_j \cdot \nabla\_r S(\mathbf{r} - \mathbf{r}') \, \mathrm{d}^3 r' &= \\ \int\_{\delta V} (\mathbf{s}\_j) \mathbf{s}\_j \cdot \nabla \|\mathbf{S}\_{\mathrm{inc}}(\mathbf{r})\|\_{\mathrm{v}} \delta V &+ \mathbf{w}\_r(\mathbf{s}\_j) \int\_{\delta V} \mathbf{s}\_j \cdot \nabla\_{r'} S\_{\mathrm{sc}}(\mathbf{r} - \mathbf{r}') \, \mathrm{d}^3 r' \quad (\text{A.10}) \end{aligned}$$

where the last term is that accounting for the scattered flux, both outgoing and inwards from outside δ*V*. We can apply Gauss' theorem provided that both these terms are properly represented

$$\int\_{\delta V} \mathbf{s}\_j \cdot \nabla\_r S\_{\mathbf{sc}}(\mathbf{r} - \mathbf{r}') \, \mathrm{d}^3 r' = \int\_{\Sigma} S\_{\mathrm{sc}}^{(\mathrm{out})}(\mathbf{r} - \mathbf{r}') \mathbf{s}\_j \cdot \mathbf{s}' \, \mathrm{d} \mathrm{S}' - \int\_{\Sigma} S\_{\mathrm{sc}}^{(\mathrm{in})}(\mathbf{r} - \mathbf{r}') \mathbf{s}\_j \cdot \mathbf{s}' \, \mathrm{d} \mathrm{S}' \tag{A.11}$$

where Σ is the surface enclosing δ*V*. The first term can be written as the sum of the contributions from the *N* scattering particles in δ*V*, while the inward flux represents the contribution from the whole outer volume *V* − δ*V*. To obtain an expression for such inward flux, we can exploit the fact that we have a medium with average optical properties which are homogeneous at every position. We can therefore assume that at each boundary between neighboring volume elements δ*V*, the total inward and outward flux can be considered

equivalent. For each point *r* lying on the surface Σ we can then write that

as

1 *v* ∂ ∂*t* δ*V*

 δ*V* d*P*abs d*V*

δ*V* and σa is their absorption cross-section (2.23).

(s · s*j*)s*j* · ∇*rS* (*r* − *r*

 δ*V*

enough to consider *r r*

(s · s*j*)s*j* · ∇*rS* (*r* − *r*

s*j* · ∇*rS* sc(*r* − *r*

 ) d3 *r* = Σ *S* (out) sc (*r* − *r* )s*j* · s

 δ*V*

represented

 δ*V* *S* · s*j* d*V* 1

On the other hand, the third term of equation (A.2) becomes

where we have introduced the volume-averaged change in *S* as

∇*rS*inc · s*j* v = 1

Plugging this substitution into equation (A.7) we obtain

 ) d3 *r* =  ) d3 *r*

> δ*V* δ*V*

*wr*(s*j*)s*j* · ∇*S*inc(*r*)vδ*V* + *wr*(s*j*)

*v* ∂ ∂*t*

in terms of the volume-averaged incident flux, where *wr*(s*j*) is the probability distribution of the energy flux to each angle defined in equation (2.39), *N* is the number of particles in

+ *wr*(s*j*)

At this point, assuming that the far-field approximation is valid within δ*V* (which we already used in order for the cross-sections to be additive in equation (A.6)) and yet that δ*V* is small

, we can interchange ∇*r* ↔ ∇*r* and write

where the last term is that accounting for the scattered flux, both outgoing and inwards from outside δ*V*. We can apply Gauss' theorem provided that both these terms are properly

where Σ is the surface enclosing δ*V*. The first term can be written as the sum of the contributions from the *N* scattering particles in δ*V*, while the inward flux represents the contribution from the whole outer volume *V* − δ*V*. To obtain an expression for such inward flux, we can exploit the fact that we have a medium with average optical properties which are homogeneous at every position. We can therefore assume that at each boundary between neighboring volume elements δ*V*, the total inward and outward flux can be considered

134

*S*incv*wr*(s*j*)δ*V* (A.5)

(s · s*j*) d*V N*σa*S*incv*wr*(s*j*), (A.6)

s*j* · ∇*rS* sc(*r* − *r*

 ) d3 *r*

∇*S*inc · s*j* v s*j* · ∇*S*incv. (A.9)

s*j* · ∇*rS* sc(*r* − *r*

 ) d3 *r*

(A.10)

d*S* (A.11)

 δ*V*

d*S* −

 Σ *S* (in) sc (*r* − *r* )s*j* · s

 ) d3 *r*

. (A.8)

(A.7)

= *wr*(s*j*)∇*rS*inc · s*j* vδ*V*

 δ*V*

s*j* · ∇*rS* inc(*r* − *r*

$$\int\_{\Sigma} S\_{\rm sc}^{(\rm in)}(\mathbf{r} - \mathbf{r}') \mathbf{s}\_j \cdot \mathbf{s}' \, \mathrm{d}S' \simeq \sum\_{i=1}^{N} \int\_{\Sigma} \left( \langle \mathbf{S}\_{\rm sc}(\mathbf{r}) \rangle\_i \cdot \mathbf{s}\_i \right) \mathbf{s}\_j \cdot \mathbf{s}' \, \mathrm{d}S' \tag{A.12}$$

where *S*sc(*r*)*i* is the energy flux scattered by the *i*-th particle inside δ*V*, given by equation (2.36). Making use of the approximation (A.4) we obtain

$$\begin{split} \int\_{\Sigma} S\_{\text{sc}}^{(\text{in})}(\mathbf{r} - \mathbf{r}') \mathbf{s}\_{j} \cdot \mathbf{s}' \, \mathrm{d}S' &\simeq \sigma\_{\text{tot}} \sum\_{i=1}^{N} \|\mathbf{S}\_{\text{inc}}(\mathbf{r})\|\_{\text{v}} \int\_{\Sigma} w\_{\mathbf{r}}(\mathbf{s}\_{i}) \frac{p(\mathbf{s}\_{i}, \mathbf{s}\_{j})}{\left|\mathbf{r}' - \mathbf{r}\_{i}\right|^{2}} \, \mathrm{d}S' \\ &\simeq N \sigma\_{\text{tot}} \|\mathbf{S}\_{\text{inc}}(\mathbf{r})\|\_{\text{v}} \int\_{4\pi} w\_{\mathbf{r}}(\mathbf{s}') p(\mathbf{s}', \mathbf{s}\_{j}) \, \mathrm{d}\mathcal{Q}', \end{split} \tag{A.13}$$

where in the last step we have identified d*S* |*r* − *ri*| 2 as the solid angle dΩ, since the total flux traversing Σ will not depend on the position of the particles within δ*V*.

At this point we have obtained the expressions for all the terms appearing in equation (A.2), which can be rewritten as

$$\begin{split} \frac{1}{\nu} \frac{\partial}{\partial t} ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}\_{j}) \delta V &+ N \boldsymbol{\sigma}\_{\mathrm{a}} ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}\_{j}) + \mathbf{s}\_{j} \cdot \nabla ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}\_{j}) \delta V \\ &+ N \boldsymbol{\sigma}\_{\mathrm{s}} ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}\_{j}) - N \boldsymbol{\sigma}\_{\mathrm{tot}} \int\_{4\pi} ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}') p(\mathbf{s}', \mathbf{s}\_{j}) \mathrm{d}\mathcal{Q}' \quad (\text{A.14}) \end{split} (4.14)$$

and normalized by δ*V* to remove the dependence from the arbitrary averaging volume

$$\begin{split} \frac{1}{\nu} \frac{\partial}{\partial t} ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}\_{j}) + \mu\_{\mathrm{a}} ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}\_{j}) + \mathbf{s}\_{j} \cdot \nabla ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}\_{j}) \\ &+ \mu\_{\mathrm{s}} ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}\_{j}) - \mu\_{\mathrm{tot}} \int\_{4\pi} ||\mathbf{S}\_{\mathrm{inc}}||\_{\mathrm{v}} \boldsymbol{\nu}\_{\mathrm{r}}(\mathbf{s}') p(\mathbf{s}', \mathbf{s}\_{j}) \, \mathrm{d}\mathcal{Q}' . \end{split} \tag{A.15}$$

Equation (A.15) is the radiative transport equation for the time-averaged incoming energy flux *S*inc flowing in direction s*j*, averaged over a volume δ*V*. It is important to note that the factor µtot from the last term appears through our definition of the phase function (2.27). It is the product µtot*p*(s,s0) that yields the contribution which is solely due to scattering, i.e. the absorption term related to the phase function that appears through µtot = µs + µa has no direct physical meaning.

By defining the specific intensity as

$$I(\mathbf{r}, t, \mathbf{s}) = \frac{1}{4\pi} \|\mathbf{S}\_{\rm inc}\|\_{\rm v} \nu\_{\rm r}(\mathbf{s}),\tag{A.16}$$

equation (A.15) recovers its well-known appearance

$$\frac{1}{\nu} \frac{\partial I(\mathbf{r}, t, \mathbf{s})}{\partial t} + \mathbf{s} \cdot \nabla I(\mathbf{r}, t, \mathbf{s}) - (\mu\_{\mathbf{s}} + \mu\_{\mathbf{a}})I(\mathbf{r}, t, \mathbf{s}) - \mu\_{\text{tot}} \int\_{4\pi} I(\mathbf{r}, t, \mathbf{s}') p(\mathbf{s}, \mathbf{s}') \, \mathrm{d}\mathcal{Q}' = \mathbf{0}. \tag{A.17}$$

It is convenient to summarize the main approximations that we applied in order to

reach equation (A.15). Namely they are

where we have introduced the terms of the *P*1 expansion (2.72). In the second term of this

which are valid for any *A* independent of s. This whole term therefore reduces to ∇*U*(*r*, *t*)/3.

dΩ

3 4π

1 − (s · s)2

*F*(*r*, *t*) Ω *p*(s,s )s*j* · s

dΩ = (s*j* · s) µs

 Ω

*U*(*r*, *t*) + µtot

where in the last term we would need the product s · s rather than s*j* · s to recover the first moment *g* of the phase function. We can reach this expression by recasting s*j* · s in terms

(s · *A*)s dΩ = 4π

3

[s · ∇(*A* · s)]s dΩ = 0 (A.22)

*A* (A.21)

1 − (s*j* · s)2 (A.24)

µtot

dΩ = µs*gF*(*r*, *t*). (A.26)

|*F*(*r*, *t*)|. (A.28)

*q*(*r*, *t*,s)s dΩ (A.27)

dΩ . (A.23)

*g*. (A.25)

 Ω

)*F*(*r*, *t*) · s

)(s*j* · s) +

 Ω *p*(s · s )s · s

3*g*µs 4π

*F*(*r*, *t*) · s

+ (µtot − *g*µs)*F*(*r*, *t*) =

At this stage, we introduce the second main assumption, namely that the temporal

It is worth commenting further this approximation, since by removing the temporal dependence of *F*(*r*, *t*) we are effectively invalidating one of the fundamental similarity relations of the RTE, relating the specific intensity in the presence of absorption with the specific intensity in a non-absorbing medium (2.53). A few authors point out that the diffusion approximation is expected to hold for light that has undergone a multitude of scattering events, and therefore absorption frustrates the diffusive regime in that it selectively extinguishes

137

= µs 4π

 Ω

Regarding the fourth term, it is the integral in dΩ of

expression we have two integrals of the form

µs 4π

of s*j* · s

*U*(*r*, *t*) + µtot

which eventually yields

The fourth term becomes then

1 *v* ∂ ∂*t*

 Ω *p*(s · s )s*j* · s

3 4π Ω *p*(s,s

s*j* · s

 Ω µs 4π

Putting all the pieces together we finally obtain

*F*(*r*, *t*) +

 = (s · s

dΩ = (s*j* · s)

*U*(*r*, *t*) +

∇*U*(*r*, *t*) 3

variation of the flux is negligible with respect to the vector itself

1 µ s*v*  ∂*F*(*r*, *t*) ∂*t*

Together, these assumptions give an idea of the required size of the averaging volume δ*V*. In order to satisfy the conditions that we set throughout this derivation, δ*V* should be a volume much smaller the total system, but still sufficiently large to contain a large number of particles in order to be characterized by statistically meaningful optical properties.

## A.2. Derivation of the diffusion equation

In this appendix we provide more details on the derivation of the diffusion equation following the *P*1 approximation. In order to obtain the expansion (2.72) from the guess (2.71), we have to determine the expressions for the functions *f*0 and *f*1. The isotropic term can be obtained by plugging equation (2.71) into our definition of *U*(*r*, *t*)

$$U(\mathbf{r},t) = f\_0(\mathbf{r},t) \int\_{\mathcal{Q}} \mathbf{d}\mathcal{Q} + f\_1(\mathbf{r},t) \int\_{\mathcal{Q}} \mathbf{s}\_j \cdot \mathbf{s} \,\mathrm{d}\mathcal{Q} \,. \tag{A.18}$$

Using Ω dΩ = 4π and Ω s*j* · s dΩ = 2π 1 −1 cos θ d cos θ = 0 we obtain *f*0(*r*, *t*) = *U*(*r*, *t*) as expected. Analogously, we can use our ansatz in the definition of *F*(*r*, *t*), which we assumed to be pointed in direction s*j*, to obtain

$$F(\mathbf{r},t) \cdot \mathbf{s}\_j = f\_0(\mathbf{r},t) \int\_{\mathcal{Q}} \mathbf{s}\_j \cdot \mathbf{s} \, \mathrm{d}\mathcal{Q} + f\_1(\mathbf{r},t) \int\_{\mathcal{Q}} (\mathbf{s}\_j \cdot \mathbf{s})^2 \, \mathrm{d}\mathcal{Q} = 2\pi \int\_{-1}^1 \cos^2 \theta \, \mathrm{d}\cos\theta = \frac{4\pi}{3} \tag{A.19}$$

and therefore *f*1(*r*, *t*) = *F*(*r*, *t*)4π/3.

The second approximation leading to the diffusive equation concerns the derivation of the time-dependent Fick's law, which is obtained by multiplying the RTE by s and integrating over the whole solid angle. By doing so, we obtain the following terms

$$\begin{split} \frac{1}{\nu} \frac{\partial}{\partial t} F(\mathbf{r}, t) + \int\_{\mathcal{Q}} \Big( \mathbf{s} \cdot \nabla \Big[ \frac{1}{4\pi} U(\mathbf{r}, t) + \frac{3}{4\pi} F(\mathbf{r}, t) \cdot \mathbf{s} \Big] \Big| \mathbf{s} \, \mathrm{d}\mathcal{Q} + \mu\_{\mathrm{tot}} F(\mathbf{r}, t) \\ - \mu\_{\mathrm{tot}} \int\_{\mathcal{Q} \mathcal{Q}'} p(\mathbf{s}, \mathbf{s}') \Big[ \frac{1}{4\pi} U(\mathbf{r}, t) + \frac{3}{4\pi} F(\mathbf{r}, t) \cdot \mathbf{s} \Big] \mathbf{s} \, \mathrm{d}\mathcal{Q}' \, \mathrm{d}\mathcal{Q} = \int\_{\mathcal{Q}} q(\mathbf{r}, t, \mathbf{s}) \mathbf{s} \, \mathrm{d}\mathcal{Q} \quad (\text{A.20}) \end{split}$$

where we have introduced the terms of the *P*1 expansion (2.72). In the second term of this expression we have two integrals of the form

$$\int\_{\mathcal{Q}} (\mathbf{s} \cdot \mathbf{A}) \mathbf{s} \, d\mathcal{Q} = \frac{4\pi}{3} \mathbf{A} \tag{A.21}$$

$$\int\_{\mathcal{Q}} [\mathbf{s} \cdot \nabla (\mathbf{A} \cdot \mathbf{s})] \mathbf{s} \, d\mathcal{Q} = 0 \tag{A.22}$$

which are valid for any *A* independent of s. This whole term therefore reduces to ∇*U*(*r*, *t*)/3.

Regarding the fourth term, it is the integral in dΩ of

$$\begin{split} \frac{\mu\_{s}}{4\pi}U(\mathbf{r},t) + \mu\_{\rm tot}\frac{3}{4\pi} \int\_{\mathcal{Q}^{\prime}} p(\mathbf{s},\mathbf{s}^{\prime})F(\mathbf{r},t) \cdot \mathbf{s}^{\prime} \,\mathrm{d}\mathcal{Q}^{\prime} \\ = \frac{\mu\_{s}}{4\pi}U(\mathbf{r},t) + \mu\_{\rm tot}\frac{3}{4\pi}F(\mathbf{r},t) \int\_{\mathcal{Q}^{\prime}} p(\mathbf{s},\mathbf{s}^{\prime})\mathbf{s}\_{j} \cdot \mathbf{s}^{\prime} \,\mathrm{d}\mathcal{Q}^{\prime} .\end{split} \tag{A.23}$$

where in the last term we would need the product s · s rather than s*j* · s to recover the first moment *g* of the phase function. We can reach this expression by recasting s*j* · s in terms of s*j* · s

$$\mathbf{s}\_{j}\cdot\mathbf{s}'=(\mathbf{s}\cdot\mathbf{s}')(\mathbf{s}\_{j}\cdot\mathbf{s})+\sqrt{1-(\mathbf{s}\cdot\mathbf{s}')^{2}}\sqrt{1-(\mathbf{s}\_{j}\cdot\mathbf{s})^{2}}\tag{A.24}$$

which eventually yields

reach equation (A.15). Namely they are

implies that

Using

1 *v* ∂ ∂*t* Ω dΩ = 4π and

*F*(*r*, *t*)·s*j* = *f*0(*r*, *t*)

*F*(*r*, *t*) +

− µtot

assumed to be pointed in direction s*j*, to obtain

 Ω

and therefore *f*1(*r*, *t*) = *F*(*r*, *t*)4π/3.

 Ω s · ∇ 1 4π

Ω,Ω

*p*(s,s ) 1 4π

• average of the energy flow: *S*inc(*r*) *S*inc(*r*)v = 1

much greater than local average scattered flow: 1

δ*V* ∇*rS* (*r* − *r*

• far-field approximation: *S*sc*i S*inc(*ri*) *p*(s*i*,s)

• incoherent scattering: *S*sc *N*

• neglect depolarization: ∇(∇ · *E*) 0

A.2. Derivation of the diffusion equation

δ*V*

) d3*r*

δ*V*

δ*V S* (*r* − *r*

• average of the incoming energy flow much larger than the local average scattered flow

• average of the gradient of the energy flow: ∇*S*inc · s*j* ∇*S*incv, which directly

Together, these assumptions give an idea of the required size of the averaging volume δ*V*. In order to satisfy the conditions that we set throughout this derivation, δ*V* should be a volume much smaller the total system, but still sufficiently large to contain a large number of particles in order to be characterized by statistically meaningful optical properties.

In this appendix we provide more details on the derivation of the diffusion equation following the *P*1 approximation. In order to obtain the expansion (2.72) from the guess (2.71), we have to determine the expressions for the functions *f*0 and *f*1. The isotropic term

dΩ + *f*1(*r*, *t*)

 Ω

2 dΩ = 2π

 1 −1

s dΩ + µtot*F*(*r*, *t*)

 Ω

s dΩ dΩ =

can be obtained by plugging equation (2.71) into our definition of *U*(*r*, *t*)

 Ω

1

 Ω (s*j* ·s)

3 4π

> 3 4π

as expected. Analogously, we can use our ansatz in the definition of *F*(*r*, *t*), which we

The second approximation leading to the diffusive equation concerns the derivation of the time-dependent Fick's law, which is obtained by multiplying the RTE by s and integrating over the whole solid angle. By doing so, we obtain the following terms

*F*(*r*, *t*) · s

*F*(*r*, *t*) · s

136

*U*(*r*, *t*) = *f*0(*r*, *t*)

Ω s*j* · s dΩ = 2π

s*j* ·s dΩ+ *f*1(*r*, *t*)

*U*(*r*, *t*) +

*U*(*r*, *t*) +

*i*=1 *S*sc*i*

) d3*r* ∇*r*

• statistically equivalent properties throughout the medium: µa,s = *N*

δ*V S* (*r* − *r*

δ*V S*inc·s d*V* 1

2 d*S* (see equation (2.21))

) d3*r*

δ*V*

δ*V* σa,s *N*tot

s*j* · s dΩ . (A.18)

cos2 θ d cos θ = 4π

3 (A.19)

*q*(*r*, *t*,s)s dΩ (A.20)

−1 cos θ d cos θ = 0 we obtain *f*0(*r*, *t*) = *U*(*r*, *t*)

δ*V S*sc·s d*V*

*V* σa,s

$$\int\_{\mathcal{Q}'} p(\mathbf{s} \cdot \mathbf{s}') \mathbf{s}\_j \cdot \mathbf{s}' \, \mathrm{d}\mathcal{Q}' = (\mathbf{s}\_j \cdot \mathbf{s}) \int\_{\mathcal{Q}'} p(\mathbf{s} \cdot \mathbf{s}') \mathbf{s} \cdot \mathbf{s}' \, \mathrm{d}\mathcal{Q}' = (\mathbf{s}\_j \cdot \mathbf{s}) \frac{\mu\_s}{\mu\_{\mathrm{tot}}} \mathrm{g}.\tag{A.25}$$

The fourth term becomes then

$$\int\_{\mathcal{Q}} \left[ \frac{\mu\_{\rm s}}{4\pi} U(\mathbf{r}, t) + \frac{\Im g \mu\_{\rm s}}{4\pi} F(\mathbf{r}, t) \cdot \mathbf{s} \right] d\mathcal{Q} = \mu\_{\rm s} \mathbf{g} F(\mathbf{r}, t). \tag{A.26}$$

Putting all the pieces together we finally obtain

$$\frac{1}{\nu} \frac{\partial}{\partial t} F(\mathbf{r}, t) + \frac{\nabla U(\mathbf{r}, t)}{\Im} + (\mu\_{\text{tot}} - g\mu\_s) F(\mathbf{r}, t) = \int\_{\varOmega} q(\mathbf{r}, t, \mathbf{s}) \mathbf{s} \, d\varOmega \tag{A.27}$$

At this stage, we introduce the second main assumption, namely that the temporal variation of the flux is negligible with respect to the vector itself

$$\frac{1}{\mu\_{\rm s}^{\prime}\nu} \left| \frac{\partial F(\mathbf{r}, t)}{\partial t} \right| \ll |F(\mathbf{r}, t)|. \tag{A.28}$$

It is worth commenting further this approximation, since by removing the temporal dependence of *F*(*r*, *t*) we are effectively invalidating one of the fundamental similarity relations of the RTE, relating the specific intensity in the presence of absorption with the specific intensity in a non-absorbing medium (2.53). A few authors point out that the diffusion approximation is expected to hold for light that has undergone a multitude of scattering events, and therefore absorption frustrates the diffusive regime in that it selectively extinguishes

= *U* 4π

= −*U* 4 + 3 4π *F*u 2π

= −*U* 4 + 3 4π *F*q 2π 3

s·q>0

 2π

coefficient *A* as

= *U* 2

= *U* 2

π/2

0

0

*R*(θi)(u · s)(s · q) dΩ =

diffusive and a transparent media

π/2

and

π

cos θi sin θi dθi +

0

dφ π

*R*(θi)*I*(*r*, *t*,s)(s · q) dΩ =

where we have used the fact that, for symmetry reasons

π/2

3 4π *F*u

2π

s·q>0

3 4π *F*q 2π

3 2 *F*q 2π

Substituting the expressions obtained in (A.31) and (A.32) into (A.30) and denoting the

we finally obtain the partial current boundary condition (PCBC) at the interface between a

0 *R*(θi) cos2 θi sin θi dθi

0 *R*(θi) cos θi sin θi dθi

π/2

π/2

*U*(*r*, *t*) − 2*AF*(*r*, *t*) · q

139

*R*(θi) cos θi sin θi dθi +

*R*(θi) cos θi sin θi dθi +

 2π 0 dφ π/2 0

*A* = 1 + 3

1 − 2

cos φ cos θi sin2 θi dθi +

(u · s)(s · q) dΩ +

3 4π *F*u

0

3 4π *F*q 2π

*R*(θi)[*U* + 3*F* · s](s · q) dΩ

dφ π/2

0

0

0

2π

dφ π

π/2

*R*(θi)(s · q)

*R*(θi) cos φ sin2 θi cos θi dθi = 0. (A.33)

2 dΩ

, (A.34)

*r*∈Σ = 0. (A.35)

*R*(θi) cos2 θi sin θi dθi (A.32)

(s · q)