1. INTRODUCTION

Monte Carlo (MC) simulations represent the tool of choice when describing “classical” photon migration in propagating media [1]. This is especially true when exact solutions (except for statistical noise) of complex media are required. In fact, when dealing with the simultaneous presence of different media types, complex geometries, and varying optical parameters, it becomes extremely difficult to use alternative approaches such us analytical solutions or finite-difference-based methods. Moreover, the continuous increase in the computational power, even for simple desktops, in general makes the MC methods increasingly attractive and easy to use. These are some of the reasons why data generated by MC simulations have nowadays become a “gold standard,” allowing testing of alternative methods.

Recently, a generalization of the classical theory also allowing one to describe anomalous radiative transfer (ART) has been introduced [2,3]. This theory permits the modeling of photon migration in anomalous media, that is, media where the classical Beer–Lambert–Bouguer law is not valid.

The release of the Beer–Lambert–Bouguer law may strongly complexify the definition of the analytical models describing ART and their related solutions. This increase in complexity may be generated, for example, by the need for the introduction of models based on fractional (integro-)differential equations, where the solutions are often difficult to obtain and manipulate. Moreover, when the investigated optical medium is not infinite, and thus necessitates the introduction of boundary conditions, the definition of the analytical models and the obtainment of their solutions may represent a real challenge.

This is why in communities such as the biomedical optics, where the main aim is not specifically to understand the mathematical theory but to obtain solutions related to complex practical problems in photon migration, more efficient and general purpose tools would be particularly appreciated. It is precisely in the context of ART that the MC methods show their extreme power, i.e., thanks to their relative simplicity in describing and solving the considered problems.

In general, the MC methods describe how photons propagate step by step in the media, through scattering, absorption, reflection, and refraction events. The main difference between the classical and the ART MC approaches is given by the specific choice of the function describing the probability $p(s\!)ds$ for a photon to reach a distance $s \in [s - \frac{{ds}}{2},s + \frac{{ds}}{2}]$ in the medium, without interactions, where $p(s\!)$ is a probability density function (pdf).

As is well known, in the classical case

In ART, it is not mandatory for $p(s\!)$ to be an exponential function, and a panoply of choices, depending on the specific optical properties of the medium, is possible.

At this point, the fundamental remark is that to simulate ART models, by means of the MC method, it is not sufficient to substitute the classical $p(s\!)$ [Eq. (1)] with the relative anomalous one in the MC code. In fact, if we proceed in such a simple way, the fundamental reciprocity law (RL) at the basis of optics will not be correctly reproduced by the simulations. As a short reminder, in the present context the RL tells us that if we exchange light sources and detectors (by inverting the direction of the detected photons; see e.g., red line in Fig. 1), we must always measure (detect) the same photon flux. For a more in-depth approach to the RL, see Ref. [4].

 

Fig. 1. Schematic of the geometrical model utilized for the Monte Carlo simulations, representing two 3D spheres centered at the axis origin.

Download Full Size | PDF

Thus, based on the findings derived from ART theory, the aim of the present tutorial is to describe the suitable way to treat this problem. This will be done by directly giving the necessary practical “recipes” allowing the building of MC codes for ART. We will see that the core of the solution is given by the particular choice of the pdf allowing generation of the length ${s_{\!u}}$ of the first step of photons that are starting from a medium boundary (defined by adjacent regions with different optical parameters) or from a fixed light source. The remaining steps, ${s_{\!c}}$, and the photon interactions with absorbers/scatterers, are treated as in a classical MC code.

In this tutorial, it is assumed that the reader already has a basic knowledge of MC simulations applied to classical problems in photon migration, typical, e.g., for biomedical optics. Very simplified numerical examples will be given, allowing intuitive understanding of the influence of a correct/incorrect approach. Note that even if in the following sections we present only simplified examples with continous-wave light sources, the described method can be applied in general to classical MC codes, e.g., to time-domain simulations implying light pulses.

To simplify the explanations, even if ART theory clearly includes the classical case, in the following sections, we will sometimes adopt the convention that the expression ART does not include the classical case. The context easily clarifies this use.

2. BASIC MC TOOLS IN ART

A. Basic Probability Density Functions for Photon Steps

ART theory is a generalization of the classical radiative transfer theory. We will describe here how photon steps are generated in ART and, in particular, how the interactions with the boundaries and light sources are treated [5,6].

In ART, photon step lengths, ${s_{\!c}}$, occurring inside the medium are generated by giving a pdf,

Up to this point, the ART approach appears to be similar to the classical one, i.e., once the suitable ${p_c}({s_{\!c}};a,b, \ldots)$ is defined, we can generate the photon step lengths, and decide if the photons are scattered, absorbed, reflected, or refracted at the different interaction sites. However, the fundamental difference clearly appears when a photon hits a boundary. In this case, if the boundary is encountered when traveling along the step of length ${s_{\!c}}$, then in ART, the photon simply stops at the boundary. From this reached position a new step is immediately generated, without the need to terminate the previous path by applying the well-known procedure usually utilized in the classical case (see e.g., Ref. [7]).

When we generate a new step of length ${s_{\!u}}$ starting from a fixed boundary or a fixed light source position (i.e., these positions obviously do not result from the probability law ${v_k}({s_{\!c}};a,b, \ldots)$), the pdf can no longer be ${p_c}(.)$. In this case, a different pdf must be considered, i.e. [5],

As we mentioned in Section 1, Eq. (3) derives from the fact that, similarly to the classical case, ART must also satisfy the fundamental RL. In other words, if we do not consider Eq. (3), e.g., by simply putting ${p_u}(.) = {p_s}(.)$, the RL is no longer satisfied (except for the classical case; see below) and photon propagation in the medium will be described incorrectly. Note that Eq. (3) permits the RL to be satisfied for any choice of optical parameters (scattering coefficients, phase functions, refractive indexes, etc.) or geometries.

The random steps of lengths ${s_{\!c}}$ and ${s_{\!u}}$, based on laws ${p_c}({s_{\!c}};a,b, \ldots)$ and ${p_u}({s_{\!u}};a,b, \ldots)$, can be numerically generated as usual by analytically solving

This is all we need to know to generate an ART MC simulation, the remaining part of the MC code being the same as in the classical case.

B. Extracting General Medium Properties from ${p_c}({s_{\!c}};a,b, \ldots)$

Once ${p_c}({s_{\!c}};a,b, \ldots)$, the geometry, the optical parameters, and the phase functions are given, it becomes possible to perform MC simulations in the ART domain. However, before running a specific simulation, it is maybe interesting to know that from ${p_c}({s_{\!c}};a,b, \ldots)$, we can derive some very general properties of the investigated anomalous medium. Strictly speaking, the properties presented here hold for an infinite medium with pdf ${p_c}({s_{\!c}};a,b, \ldots)$, but are of fundamental importance in understanding the physics described by ${p_c}({s_{\!c}};a,b, \ldots)$. Thus, for the sake of completeness, in the following paragraphs we will report four functions, directly derived from ${p_c}({s_{\!c}};a,b, \ldots)$, describing some of these properties. The reported equations are derived from well-known results in renewal theory, and the interested reader can refer to Refs. [8,9]. The equations hold for a photon propagation through any homogeneous part of the medium.

3. ART MC SIMULATIONS: EXPLANATORY EXAMPLES

The aim of the following examples is to show the importance of Eq. (3) (the core of the ART MC simulations). This will be done by demonstrating, with easy-to-understand numerical examples, that if we neglect Eq. (3), a fundamental physical property of the system cannot be reproduced, i.e., the invariance property (IP) will not be satisfied (see below). Note that IP represents a very powerful test allowing us to check the reliability and correctness of any MC code [10].

A. Two General Analytical Models

We will define here two possible models for ${p_c}({s_{\!c}};a,b, \ldots)$, i.e., the power law and the constant step models.

1. Power Law

The first model we will consider is called a “power law” and is expressed as

The specific ${s_{\!c}}$ dependence of Eq. (11) is particularly interesting because it implies that there may be an asymptotic diffusion limit of this anomalous model [11] (related to a “classical” or a “fractional” diffusion equation, depending on the choice of $a$; see below). Thus, this simple model [Eq. (10)] represents a kind of unified summary for the classical and anomalous approach.

Then, the pdf ${p_u}({s_{\!u}},\ell ,a)$ can be derived by applying Eq. (3), to obtain

Finally, the algorithms allowing us to generate random step length ${s_{\!c}}$ and ${s_{\!u}}$ values in the MC codes can be obtained by introducing Eqs. (10) and (12) into Eq. (4), and by solving as a function of ${s_{\!c}}$ and ${s_{\!u}}$, respectively. This gives

Note that the parameter $a$ may be experimentally determined (e.g., by fitting of experimental data). However, in the present context, we will just use two representative values (see below).

2. Constant Step

The second model we will consider in this tutorial is called “constant step” because it describes photons that “jump” always with steps of the same length $\ell$, i.e.,

Equation (3) is then written as

The generating functions for ${s_{\!c}}$ and ${s_{\!u}}$ are then expressed as [Eq. (4)]

B. Numerical MC Examples

Let us now define four different examples of MC simulations based on Eqs. (10) and (15). In the following sub-sections, we will give the geometry of the problem and all the necessary optical quantities, together with the light source and detector positions. The general properties (geometry independent) of the four chosen tutorial examples are treated.

1. MC Simulations

Geometry: The common geometrical model chosen in this tutorial consists of two concentric spheres centered at the axis origin, with radius ${r_{{\rm{\rm in}}}} \lt {r_{{\rm{\rm out}}}}$ and matched refractive indexes at all the boundaries, i.e., the refractive index is the same everywhere (Fig. 1). The mean free paths inside the spheres are ${l_{{\rm{\rm in}}}}$ and ${l_{{\rm{\rm out}}}}$. The absorption coefficients and the anisotropy parameters of the spheres are set to zero (mandatory conditions to test the IP in ART; see Section 3.C).

Forward measurements: An isotropic and uniform radiance illumination impinging on the surface of the external sphere (continuous uniform distribution of Lambertian sources) is utilized as a light source. Note that, in Fig. 1, only one representative photon is reported (black path; forward propagation). Photons are detected on the whole surface of the external sphere. The pathlength of each single photon, $s$, is computed, and the mean pathlength $\langle s\rangle _{{\rm{MC}}}^{{\rm{Fwd}}}$ is assessed. One hundred independent simulations (each of them has ${10^7}$ photon trajectories) were performed to evaluate the standard error of the mean. Isotropic scattering was used. Random photon steps for the three examples of power law are generated using Eqs. (13) and (14), for three representative $a$ values, i.e., $a \in \{1,1/2, + \infty \}$ with $\ell \in \{{\ell _{{\rm in}}},{\ell _{{\rm out}}}\}$ (the choice of ${\ell _{{\rm in}}}$ or ${\ell _{{\rm out}}}$ depends on where the photon is situated along the path). For the constant step example, photon steps are obtained using Eqs. (17) and (18) with $\ell \in \{{\ell _{{\rm in}}},{\ell _{{\rm out}}}\}$.

Adjoint measurements: Adjoint measurements are assessed by inverting the direction of each photon reaching the sphere surface (detector). Then, these photons, with their new inverted directions (red dashed path in Fig. 1), are utilized as light sources to generate the adjoint measurements, with the same conditions utilized for the forward measurements. As was the case for the forward measurements, photons are detected on the whole surface of the external sphere. At the end of the simulation, the mean pathlength, $\langle s\rangle _{{\rm{MC}}}^{{\rm{Adj}}}$, is then computed.

The importance of ${p_u}({s_{\!u}},\ell ,a)$: Forward and adjoint MC simulations are repeated two times: (1) by using the rules for the MC ART presented in Section 2.A, giving $\langle s\rangle _{{\rm{MC}}}^{{\rm{Fwd}}}$ and $\langle s\rangle _{{\rm{MC}}}^{{\rm{Adj}}}$, and (2) by setting ${p_u}(.) = {p_c}(.)$, as is usually done in classical MC, resulting in this case in two mean pathlengths that we will denote $\langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Fwd}}}$ and $\langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Adj}}}$.

2. Power Law: $a = 1$

In this case, random step lengths ${s_{\!c}}$ and ${s_{\!u}}$ are generated as [Eqs. (13) and (14)]

To better understand the physical properties of this model, we will compare it with the classical one (see below), through the functions presented in Section 2.B. In this case, Eq. (5) is expressed as

To obtain the remaining Eqs. (7)–(9), let us first compute Eq. (6), i.e.,

 

Fig. 2. Probabilistic functions (see Section 2.B) characterizing the models chosen as tutorial cases (see Sections 3.B.2–3.B.5).

Download Full Size | PDF

Note that the asymptotic limit [Eq. (11)] of this power law with $a = 1$ can be described by a classical diffusion equation [11].

3. Power Law: $a = 1/2$

The random steps ${s_{\!c}}$ and ${s_{\!u}}$ for $a = 1/2$ are generated as [Eqs. (13) and (14)]

This is an interesting case, because the asymptotic limit of this model is related to a fractional diffusion equation with a derivative of order $3/2$ [11]. In this case [Eq. (5)],

Then, Eqs. (7)–(9) can be obtained numerically with the same procedure presented in Section 3.B.2 (see Fig. 2).

4. Power Law: $a = + \infty$

The random step lengths ${s_{\!c}}$ and ${s_{\!u}}$ are generated as

In this case, ${s_{\!c}} = {s_{\!u}}$, because the correlated and uncorrelated pdfs are the same, i.e.,

We immediately recognize here the pdf of the classical model, related to the underlying Beer–Lambert–Bouguer law. In fact [Eq. (5)],

It is in this sense, i.e., when $a = + \infty$, that the classical model is naturally included in the tutorial ART power model of Section 3.A.

In this special case, by deriving before Eq. (6), i.e.,

Equations (7)–(9) can be completely obtained analytically, i.e.,

5. Constant Step

By following the same procedure utilized for the power law, we obtain for the constant step model

C. Invariance Property

As was mentioned in Section 2.A, the introduction of Eq. (3) allows us to suitably take into account the RL in the MC ART simulations. Neglecting this rule, the RL may not be satisfied. In this context, the examples chosen for this tutorial represent a special instructive case. In fact, it is easy to see that, due to the absence of absorption and the particular geometry of the light source and detector, the RL will always be satisfied, independently of the remaining optical parameters or other functions utilized for the MC simulation. The reason for this special behavior is because all the emitted photons always reach the detector (i.e., same photon flux through the detector in forward or adjoint direction). This means that testing the RL is not always a good way to check a MC code and that a more universal method is needed. For this reason, the IP is introduced in the following paragraph and is then exploited to show the importance of the right choice of ${p_u}(.)$ [Eq. (3)].

The IP [13,14] is a very powerful property that in the case of anomalous transport states the following: let us apply a uniform and isotropic radiance illumination at the external boundary, $S$, of a finite, scattering, and non-absorbing medium of any shape, volume $V$, matched refractive indexes, and isotropic phase function. Then, let us apply the rules reported in Section 2.A (to suitably describe ART); the mean pathlength, ${\langle s\rangle _{{\rm{theory}}}}$, spent inside V by photons outgoing from $S$, is an invariant quantity independent of the scattering strength. Note that the volume $V$ can also be the union of any number of sub-volumes with different scattering coefficients, but matched refractive indexes (external medium also matched). The quantity ${\langle s\rangle _{{\rm{theory}}}}$ is expressed as [13,14]

This property holds for anomalous and classical photon propagation and can be utilized to test the MC code validity. It is worth noting that in the case of classical photon transport (Beer–Lambert–Bouguer law), the hypothesis of an isotropic phase function is not necessary [13–15]. A generalization of Eq. (39) for unmatched refractive indexes also exists [15,16], but this goes far from the aim of the present tutorial. However, in this context, we need a general formulation valid for ART, implying the use of an isotropic phase function.

4. Results

A. General Comparisons of Classical and ART Models

Before proceeding with MC simulations, let us compare (semi-)analytically the power law and the constant step models with the classical one [Eq. (28)].

In Fig. 2, the power law [Eqs. (21) and (25)], for $a = 1$ and $a = 1/2$, and the constant step [Eq. (34)] models, are compared with the special case $a = + \infty$ [Eq. (29)], related to the classical Beer–Lambert–Bouguer law. For simplicity, we have chosen the case $\ell = 1$, but the general considerations remain valid for other $\ell$ values.

In Fig. 2(a), the classical model for ${P_c}(.)$ appears as a straight line (logarithmic scale). It is known from renewal theory that this model is representative of a Markovian or memoryless process. This means that a future photon step length only depends on the present state and not on the past history of the photon propagation. The remaining models that deviate from this behavior, as shown in Fig. 2(a), are not memoryless. This in part explains why it may be more complex to express these models through an (integro-)differential equation formalism. In general, ${P_c}(.)$ may be assessed experimentally by the measurement of the so-called ballistic photons, while the remaining functions in Fig. 2 may only be indirectly derived by means of a mathematical procedure. Figure 2(a) clearly shows that the power low curves for ART go slower to zero compared to the classical case. This means that there is a higher probability for a photon to suddenly undergo very long steps. This may be, e.g., the case of photons that pass through a non-homogeneous distribution of clouds.

Figure 2(b) represents the average number of scattering events (or steps), $\langle {m_c}(.)\rangle$, along a path of length ${s_{{cp}}}$. The figure highlights the fact that even if in general the average step (mean free path) is always $\ell$ for the power law [Eq. (10)] and the constant step [Eq. (15)] models, for the specific ART condition, it is not possible to simply divide the path length ${s_{\!c}}$ by $\ell$ to find the $\langle {m_c}(.)\rangle$ value as can be done in the classical case [Eq. (31)] (black color; straight line). In fact, Fig. 2(b) shows that ART curves are not linear, and thus a specific calculation must be performed [Eq. (7)]. This particular behavior is generated by the “memory” of the ART models. Function $\langle {m_c}(.)\rangle$ is extremely important because renewal theory tells us that it uniquely determines the process generating the photon steps.

Figure 2(c) reports the probability ${v_k}(.,.)$ that $k$ scattering events occur on a path of length ${s_{{cp}}} = 5$ mm. As expected, in the classical case, ${v_k}(.,.)$ is a Poisson distribution [Eq. (32)], owing to the exponentially distributed photon steps. The function ${v_k}(.,.)$ can be utilized to reproduce the statistical spatial distribution of the scattering sites inside the medium, and it gives us an idea of how the medium is structured.

Figure 2(d) reports the probability ${F_k}(.)$ that the sum of the first $k$ photon steps does not exceed ${s_{{cp}}} = 5\,\,{\rm{mm}}$. The classical case has smaller values than the ART power law (apart from the cases $k = 1,2$). The exceptions for $k = 1,2$ can be intuitively explained by looking at, e.g., Fig. 2(a), from which we can deduce that the probability to have a long step $\gt$5 mm is higher for the classical case.

Thus, Fig. 2 gives us the main features of a given ART model compared to the classical one. The above (semi-)analytical approach may be somewhat useful for the interpretation of the MC simulations.

B. MC Results

The results of the MC simulations for the four cases and the geometrical model in Fig. 1 are summarized in Table 1. The theoretical mean pathlength that we must reproduce with the MC simulations is ${\langle s\rangle _{{\rm{theory}}}} = \frac{4}{3}{r_{{\rm out}}}$ mm [Eq. (39)].

Table 1. Explanatory Tests of the IP for the ART Theory^a

View Table

From Table 1, it clearly appears that $\langle s\rangle _{{\rm{MC}}}^{{\rm{Fwd}}}$ perfectly reproduces the theoretical ${\langle s\rangle _{{\rm{theory}}}}$ IP value (within the statistical error) for all cases, i.e., the power law, the constant step, and the classical model. Moreover, $\langle s\rangle _{{\rm{MC}}}^{{\rm{Fwd}}} = \langle s\rangle _{{\rm{MC}}}^{{\rm{Adj}}}$ (within the statistical error).

However, if the RL condition is not satisfied, i.e., if we neglect the special behavior for uncorrelated scattering events and set ${p_u}(.)$ equal to ${p_s}(.)$, then the MC simulations generate results in disagreement with the laws of optics, i.e., $\langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Fwd}}} \ne \langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Adj}}}$, $\langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Fwd}}} \ne {\langle s\rangle _{{\rm{theory}}}}$ and $\langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Fwd}}} \ne {\langle s\rangle _{{\rm{theory}}}}$, when $a \in \{1,1/2\}$. The only exception is the classical case ($a = + \infty$), where all results remain valid (last line of Table 1). This is obviously because in this specific case, ${p_u}(.) = {p_s}(.)$ [Eq. (28)].

Considering the problem geometry, the inequality $\langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Fwd}}} \ne \langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Adj}}}$, for the ART, immediately tells us that the light reaching the detector in the forward direction does not represent an isotropic and uniformly distributed radiance. In other words, one of the conditions for the IP is not satisfied for the adjoint case. As expected, in the classical case, this problem does not subsist, i.e., $\langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Fwd}}} = \langle s\rangle _{{{\rm{MC}}_{{\rm{no}} - {\rm{reciprocity}}}}}^{{\rm{Adj}}}$.

In Table 1, the data are exact up to the third decimal, but obviously, if we increase the number of photon trajectories per simulation, we can also increase the precision of the results at any desired level.

5. DISCUSSION AND CONCLUSIONS

Sections 3 and 4 highlighted the importance of distinguishing between correlated and uncorrelated scattering events through the pdfs ${p_c}(.)$ and ${p_u}(.)$. In fact, not only does this approach satisfy (by construction of the theory) the RL for photon propagation, but it is also mandatory if one wants to satisfy the universal IP.

In Table 1, we have also seen that there is an exception for the case $a = + \infty$, corresponding to the classical case where photons propagate under the Beer–Lambert–Bouguer law. In fact, in this special case (the only one) ${p_u}(.) = {p_s}(.)$, and therefore, the results remain exact in all cases (Table 1; last line).

In practice, the ART theory is a generalization of the classical theory, the latter appearing as a particular case. Thus, the ART MC method can be applied to solve classical problems in photon transport, e.g., in biomedical optics. This may greatly simplify the writing of MC codes (the photon step just stops at the boundary, before the next step), because it does not necessitate treatment of the photons reaching the boundaries with more complex algorithms as is classically done [7].

But then, do we have two different algorithms to treat the photon interactions with the boundaries in classical MC? Which is the right one? Actually, the classical and the ART-based MC methods give the same results. Intuitively, it is probably hard to see that the two approaches are equivalent. For this reason, another simplified example might help the reader to get a better idea. Let there be a medium with only one internal boundary and mean free paths ${\ell _1}$ and ${\ell _2}$. Let us generate a photon step length, $s = {s_1} + {s_2}$, as in Fig. 3. It is well known that in classical MC, the random step $s$ is obtained as

 

Fig. 3. Schematic of a photon that crosses a boundary separating two regions with a different mean free path.

Download Full Size | PDF

In the same way, as was explained in this tutorial, in ART, the MC random step $s$ is generated as

At this point, without entering into too much technical detail, it is maybe worth making a general consideration concerning the RL and the IP. In fact, historically, the pdf ${p_u}(.)$ has been derived independently for the RL and for the IP theory. Each theory having its own independent (but compatible) physical assumptions allows us to assess ${p_u}(.)$. In this didactical contribution, we have chosen a point of view that is not historical but that should better highlight the relationship existing between RL and IP. In fact, the RL is a more “general” law than the IP. This is because RL holds also for any light source, detector configuration, optical property values, or phase function choice. Thus, the pdf for a step starting from a fixed position must be ${p_u}(.)$ if we want RL to be satisfied.

In contrast, historically, the pdf ${p_u}(.)$ necessary to satisfy the IP has also been derived independently [14] (the obtained ${p_u}(.)$ equals the one satisfying the RL), but by using a different hypothesis (actually, a subset of conditions of the RL). We realize here that this derivation for the IP is not strictly necessary because ${p_u}(.)$ is already given by the RL. In practice, it is possible to demonstrate the IP by considering ${p_u}(.)$ already known. We invite the interested reader to revisit the original papers on IP (e.g., Ref. [14]) with this point of view. This might help to have a more coherent vision on the strong relationship existing between RL and IP.

The few examples presented in this tutorial had the aim to show how to generate MC simulations for any ${p_c}(.)$ describing a desired medium. However, we must be aware of the fact that to describe the exact physics of the investigated medium, the first moment of ${p_c}(.)$ must be finite [Eq. (3)]. From a practical point of view, this constraint on ${p_c}(.)$ is not a real problem, because when simulating media representing actual physical systems, it is unlikely that photons go to infinity without any interaction. This keeps the average step length finite.

As highlighted in the Introduction section, the proposed method also holds in general for the CW and time domain. In fact, CW and time domain simulations only differ by the way detected photons are classified in the relevant computed quantities such as spatial or temporal histograms. In other words, CW domain and time domain results can be obtained in any MC approach from the same photons’ trajectories.

In conclusion, it may be worth noting that an MC code compatible with the ART theory may also be of a more general interest than a simple technicality for tissue optics. In fact, it has been shown that some biological tissues, such as bone and lung, have a particular fractal-like structure [17–19]. As reported in Ref. [20], (multi-)fractal structures may produce anomalous photon propagation. This observation may open a new domain of exploration in biomedical optics, and adapted MC codes may become one of the important tools to address this interesting topic.