Mouna El-Hafi
email: mouna.elhafi@mines-albi.fr

Stéphane Blanco

Jérémi Dauchet

Richard A Fournier

Mathieu Galtier

Loris Ibarrart

Jean-Marc Tregan

Najda Villefranque

Mouna El Hafi

Stephane Blanco

Mouna El Hafi

Three viewpoints on null-collision Monte Carlo algorithms

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Since the origin of numerical simulations of radiative transfer, it has been claimed that statistical approaches such as Monte Carlo were the only practical "way or path" towards the simultaneous handling of all the spectral and geometric complexity of radiation in 3D realistic systems. This complexity tends to infinity as soon as either 1. the degree of detail that describes the heterogeneous medium of propagation tends to infinity, 2. the spectral resolution of the optical properties of the medium tends to zero, 3. Fredholm integral equations are involved, e.g., in the context of multiple scattering, which theoretically yields an integration domain of infinite dimension.

The reason why "Monte Carlo is the only numerical tool that passes infinite dimension" is the use of the double randomization technique: the expectation of a linear function of another expectation is still just one expectation, or

E[h(E[X])] = E[h(X)] (1)
where h is a linear function and X a random variable. This is trivial and can be easily demonstrated using the linearity of the integral operator. However, the most immediate implication of this property is the key to the renowned power of Monte Carlo: when estimating the expectation of a random variable that is a linear function h of the expectation of other random variables (representing, for example, spectral, spatial or directional variables), then only one sample of each "secondary" random variable is needed to provide one sample of the global random variable. This technique is described in detail in Sabelfeld [START_REF] Sabelfeld | Monte Carlo methods in boundary value problems[END_REF] and identified as a key point in the Monte Carlo community of applied mathematics (see, for example, [START_REF] Sylvain | Stochastic finite differences for elliptic diffusion equations in stratified domains[END_REF][START_REF] Dimov | A new walk on equations monte carlo method for solving systems of linear algebraic equations[END_REF]).

However powerful, this property is strictly limited to cases where the h function is linear. Radiative transfer is, of course, "linear transport" and the exponential extinction of a beam (after Beer's law) is a signature of this linearity. But in the context of Monte Carlo methods, with the necessary writing of the computed quantity as an expectation, the nonlinearity of the exponential function is at the origin of severe difficulties as soon as the medium of propagation is spatially heterogeneous. Indeed, the integral over the heterogeneous extinction field appears within the exponential function, whereas the integration over multiple scattering optical paths appears "outside" the exponential (these paths define the line-of-sight along which Beer's extinction is applied). The integrals are combined nonlinearly; double randomization can no longer be employed (Eq. ( 1) is no longer true); a crucial feature of the Monte Carlo technique is lost (see [START_REF] Dauchet | Addressing nonlinearities in Monte Carlo[END_REF] for more advanced considerations on the matter).

A solution to bypass this nonlinearity without biasing the results or decreasing convergence rates has been commonly used in various parts of the Monte Carlo literature ([31, 23, 2, 16, 12, 20]). The main idea is to make use of "null colliders", fictitious particles that are added to the true extinction field in order to make it homogeneous, which reduces the optical-depth integral to a simple product.

When a path encounters one of these null colliders, it simply continues forward as if the collision had not occurred. From an intuitive or physical point of view, null colliders are pure scatterers characterized by a strictly-forward phase function. In the scope of probability theory, null collisions can be seen as the rejected samples that come from sampling a density function that overestimates the frequency of collisions instead of the true distribution of free-path lengths, as per a standard rejection method.

The null-collision method has long been considered a numerical trick to avoid a heavier, deterministic integration of the extinction field. The fact that the line-of-sight integration was shifted from inside to outside the exponential in the underlying mathematical formulation was only recently made explicit, in the seminal paper of Galtier et al. in 2013 [START_REF] Galtier | Integral formulation of null-collision monte carlo algorithms[END_REF]. Subsequently, the authors and colleagues have highlighted the important implications of this reinterpretation for a diversity of research and applied fields, including combustion, spectroscopy, solar energy, atmospheric radiative transfer and image rendering [START_REF] Eymet | Null-collision meshless monte-carlo-application to the validation of fast radiative transfer solvers embedded in combustion simulators[END_REF][START_REF] Galtier | Radiative transfer and spectroscopic databases: A line-sampling monte carlo approach[END_REF][START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF][START_REF] Dauchet | Addressing nonlinearities in Monte Carlo[END_REF][START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF].

Some examples of the new ideas that could only be explored and implemented thanks to the revisiting of the null-collision method are listed hereafter:

• Villefranque et al. [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF] investigated the question of the acceleration of path tracing in spatially heterogeneous media using null collisions in combination with tools from the computer graphics community. At the junction of physics and computer graphics, they structured the data describing highly heterogeneous extinction fields into octrees to allow both fast traversal of the arbitrarily complex spatial domain and limited time spent in sampling null-collision events. In the same way, the field of image rendering has largely benefited from the formal framework that resulted from Galtier et al.'s work in 2013, bringing them to develop and implement new null-collision algorithms with renewed efficiency and confidence [START_REF] Novák | Residual ratio tracking for estimating attenuation in participating media[END_REF][START_REF] Kutz | Spectral and decomposition tracking for rendering heterogeneous volumes[END_REF].

• Tregan et al. [START_REF] Tregan | Convergence issues in derivatives of monte carlo null-collision integral formulations: A solution[END_REF] propounded a solution to avoid a convergence issue due to the use of null-collision algorithms to compute sensitivity estimates of a radiative quantity. This was only possible by shifting from the intuitive to the probabilistic point of view on the null-collision method, thereby modifying the integral formulation associated with the first null-collision algorithm into a new formulation (and hence, algorithm) yielding better statistical properties.

• Galtier et al. [START_REF] Galtier | Radiative transfer and spectroscopic databases: A line-sampling monte carlo approach[END_REF] have recently applied the idea of using null collision algorithms to spectrally integrate radiative quantities without pre-computing the absorption spectra. In this proposition, the null collisions are again a way to bypass the nonlinearity of the exponential: by shifting the sum over the energetic transition line contributions to the local absorption coefficient from inside to outside the exponential function, the double randomization technique can again be used where it could not have previously been applied without the null-collision method. In the conception of this algorithm, the formal or mathematical point of view was infinitely more helpful than the intuitive or physical one.

Obviously, since their first apparition in the literature, justifications for the use of null-collision algorithms (NCA) have been based on various arguments, depending on the community and context of the application. As the method enjoys growing interest and is used by many scientists from fields as diverse as those cited above, we find it interesting to suggest a classification of the different interpretations of the concept of null collisions. Our intention is not to compare or rank these viewpoints, but rather:

1. to try and reduce the classification to a limited number of ideas (that is, three points of view);

2. to argue that although the second and third viewpoints are less immediate than the similitude viewpoint (the most standard definition of null collisions), they produce quite different theoretical developments toward the same conclusions, offering new formal perspectives;

3. to illustrate the practical benefits that can be expected from these new perspectives (algorithms that would be difficult to establish from only the similitude viewpoint).

In the following sections, we present three readings of the method that we think allow convenient changes in perspective regarding the concept of null collisions: from the physical viewpoint (Section 2, a similitude), to the statistical (Section 3, a rejection method), to the mathematical (Section 4, a Taylor expansion of the exponential) . To illustrate the practical signification of these three viewpoints, we briefly describe recent research that was made possible or facilitated by a shift in the mental representation of the null-collision concept.

The ideas behind these examples and the supporting illustrations are not always original, in the sense that some have been published already. The originality of this work rather resides in the reformulation of the null-collision concept, as well as in the connections that are made between the reported studies. Furthermore, as this research covers a wide range of scientific fields, publications were sometimes addressed to specialists of a particular domain rather than to the community of radiative transfer scientists. With this paper, we hope to facilitate access to this literature.

Adding forward-scatterers: a similitude

Figure 1 illustrates this first viewpoint. It is the most prevalent one in the literature of radiative transfer, and also the most intuitive one, as it can be enunciated entirely within the framework of the physics of transport: fictitious particles are added to the medium to make it homogeneous, without interacting with radiation.

Indeed, designing a Monte Carlo code raises the question of making the best computer-science choices for fast tracking of multiple-scattering and multiple-reflection paths in complex 3D scenes. Path tracing requires finding the location of the next event along the simulated path: it is at the shortest distance between i/ the distance to the next collision in the volume and ii/ the distance to the next surface intersection. Sampling the next volume collision is not straightforward in a heterogeneous medium because the transmittance law is not invertible analytically.

A common solution is to compute the "line-of-sight" optical depth by numerically integrating the extinction along the path, until a sampled value of "free" optical depth is reached. Often, the field of

A similitude viewpoint

Let us write the stationary Radiative Transfer Equation (RTE) in a purely scattering medium, where f is the distribution function, k s the scattering coefficient, Φ the scattering phase function:

ω.∇f = -k s f + 4π k s f Φ(ω|ω ) dω (2)
We assume that several species in the medium are responsible for the total scattering:

• blue scatterers with scattering coefficient k s1 and phase function Φ 1 ,

• red scatterers with scattering coefficient k s2 and phase function Φ 2

such that k s = k s1 + k s2 , k s Φ = k s1 Φ 1 + k s2 Φ 2 : ω.∇f = -(k s1 + k s2 )f + 4π k s1 f Φ 1 (ω |ω) dω + 4π k s2 f Φ 2 (ω |ω) dω (3)
Once a collision is found in the medium, the type of scatterer is sampled according to the local proportion of species: blue scatterers are found with a probability k s1 /k s while red scatterers are found with a probability k s2 /k s .

(a) Scattering by multiple species: considering that different species with different extinction coefficients and phase functions are responsible for the medium total optical thickness, the RTE can be written as a sum of contributions from each scattering species

In the null-collision method, species that has no effect on the radiation but makes the scattering coefficient field homogeneous is added. The second type of scatterer in equation ( 3) is null : it has no physical meaning, hence k s2 = k n can be chosen such that k s1 +k s2 is constant in the medium, and Φ 2 (ω |ω) set to δ(ω-ω ), i.e. the null phase function is such that the path continues in the same direction as if nothing had happened:

ω.∇f = -(k s1 + k n )f + 4π k s1 f Φ 1 (ω |ω) dω + 4π k n f δ(ω -ω ) dω ( 4
)
The RTE is unchanged because k n f cancels out with 4π k n f δ(ωω ) dω . In the corresponding Monte Carlo algorithm, the next-collision location will be sampled from the total extinction field ks = k s1 + k n . A true collision will occur with probability k s1 / ks while a null-collision will occur with probability k n / ks .

(b) NCA viewpoint: in null-collision algorithms, a specie is added that does not impact the radiation, but makes the field of colliders density homogeneous, hence allowing straightforward sampling of the next-collision location extinction coefficient (the input data) is meshed, with a constant value in each gridcell (see Figure 2 for an illustration). Computing the line-of-sight optical depth in a deterministic way then comes back to iterating over the crossed gridcells, and for each one: computing the length of the ray in the gridcell, accessing the data to retrieve the local extinction coefficient, and adding the current-gridcell lengthcoefficient product to the cumulated optical depth. In this regular tracking method, the original field is the support of the path-tracing algorithm, which is by construction intrinsically dependent upon the data structure and size. If the field is highly resolved and heterogeneous with a large range of extinction variations, then many gridcells might be crossed with negligible contribution to the optical depth. For a given field, doubling the spatial resolution of the data implies that finding the next collision in the volume requires at least twice the computational time.

The first benefit of null collisions is undoubtedly that the resulting path-tracing algorithms are strictly independent of the original representation of the field, as in Figure 2-b). This independence was, for example, pointed out in the Appendix of Marshak et al.'s 1995 paper [START_REF] Marshak | Radiative smoothing in fractal clouds[END_REF] in the context of atmospheric radiative transfer, and was at the heart of Eymet et al.'s work [START_REF] Eymet | Null-collision meshless monte-carlo-application to the validation of fast radiative transfer solvers embedded in combustion simulators[END_REF] for combustion applications. More recently, research groups in the main animation studios of the cinema industry have clearly greeted the proposal to introduce null-collision with considerable enthusiasm [START_REF] Novák | Residual ratio tracking for estimating attenuation in participating media[END_REF][START_REF] Kutz | Spectral and decomposition tracking for rendering heterogeneous volumes[END_REF][START_REF] Szirmay-Kalos | Unbiased estimators to render procedurally generated inhomogeneous participating media[END_REF].

Indeed, data-algorithm orthogonality (independence) is a well-known concept in the image-rendering community who have been capitalizing on its benefits for years, as far as surface rendering is concerned [START_REF] Pharr | Physically based rendering: From theory to implementation[END_REF]. Most of the recent advances in computer graphics for the cinema industry are direct results of this orthogonality: they are mainly related to a hierarchical reorganization of the original data, which allows fast crossing of the resulting structure, yielding fast testing of the ray-surface intersection for a large set of surfaces [START_REF] Wald | Interactive rendering with coherent ray tracing[END_REF][START_REF] Ingo | Realtime Ray Tracing and its use for Interactive Global Illumination[END_REF][START_REF] Wald | Embree: A kernel framework for efficient cpu ray tracing[END_REF]. These advances benefit the artists and modellers who can freely design their numerical scenes, regardless of how the renderer will eventually handle the millions of elementary surfaces that they output.

Since null-collision techniques preserve data-algorithm orthogonality, many of the tools developed to accelerate path-tracing in complex surfaces can now be used for path-tracing in complex volumes. In practice, accelerating structures are built prior to path tracing by merging gridcells from the original highly-resolved field into larger voxels that efficiently capture the heterogeneous features of the true field.

An example of such a structure is presented in Figure 3. The structure is filled in the spirit of the null-collision method: in each merged voxel, the field of extinction coefficient is homogeneous, i.e. it is set to the maximum coefficient value found among the merged gridcells. Using a null-collision algorithm means that the paths will travel through this new medium, where null collisionners have been added.

When a collision is found, the algorithm samples the nature (true or null) of the collision and proceeds accordingly, much as in the examples of Figure 1.

A striking example of a complex volume is the cloudy atmosphere, as shown in Figure 4-a). Clouds occupy a small region of a scene, and are made up of thousands of microscopic water droplets. Here, the extinction field is structured into a regular 3D mesh with a uniform extinction coefficient in each gridcell (one gridcell volume is 25×25×25 m 3 and the domain volume is 6.4×6.4×4 km 3 ). A large number of rays will only intersect non-cloudy, empty gridcells, while others will visit the cloudy, optically-deep parts of the scene, where thousands of collisions can occur.

In Villefranque et al. [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF], it is shown that the efficiency of null-collision algorithms relies on a compromise between two extremes:

1. not merging any gridcells, as in the regular tracking method illustrated in Figure 2-a): null collisions will never happen but a lot of time will be spent crossing the regular mesh and accessing the relevant data in memory;

2. merging the whole field into one homogeneous voxel, as illustrated in Figure 2-b): the null-extinction coefficient might be very large in some regions of the scenes, where a lot of time would be spent testing collision nature and re-sampling free optical depths;

Figures 1 and5 of [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF], reproduced here as Figure 5 and Figure 6, present evidence of the dataalgorithm orthogonality: computing time is insensitive to the complexity of the scene description, for both surface and volume.

Number of voxels in raw data

Figure 6: Insensitivity of computing time to scene complexity is achieved when using hierarchical grids with null collisions (full line with squares, τ = 1). τ is the merging criterion for the construction of the hierarchical grids: the voxels are merged while the resulting optical depth is below the value of τ . For τ = 0, (dashed line with circles) the voxels are never merged, yielding hierarchical grids that are the same resolution as the input data. Reproduced from Figure 5 of [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF].

The computer graphics community organises its developments around open-source libraries that are available to our community and can be used on a regular basis [START_REF] Attila T Áfra | Embree ray tracing kernels: overview and new features[END_REF][START_REF] Star | Complexity engineering[END_REF]. With their growing interest in the null-collision method, new libraries will probably be developed and distributed to physicists. An example of this has been described and illustrated in [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF], in the field of atmospheric science. Their implementation is also illustrated in Figure 4 where a library-based renderer is used to compute a physically-correct synthetic image of a cloudy field that is the output of a high-resolution Large Eddy Simulation; and in Figure 4b in which similar tools are used to compute pixel-averaged ground fluxes, each estimated using an independent reverse Monte Carlo simulation involving 10 5 samples. Also displayed is the corresponding statistical uncertainty. These simulations are then used to analyse the radiative effects of clouds at the surface, to improve their representation in climate models. Using null-collision algorithms allows radiative transfer to be simulated in large datasets, in which computations would otherwise be impracticable.

At this stage, we have retained a viewpoint under which all we see is a rigorous similitude between the original radiative-transfer problem and a new one in which additional scatterers with a strictly-forward phase function have been introduced in the medium. The benefits associated with this similitude are essentially in terms of computer implementation.

Rejection sampling

An alternative way of describing the use of null collisions is to consider that when a null collision is encountered, the sampling of the collision's location is rejected. This viewpoint is illustrated in Figure 7.

The problem of the extinction coefficient heterogeneity appearing within the exponential can indeed be reduced to the difficulty of sampling a collision location x according to the probability density function p(x) = k e (x) exp -x 0 k e (y)dy Thus, the algorithm that rejects null collisions until a true collision is found is nothing more than a rejection-sampling algorithm. Under this viewpoint, no similitudeis used, no "forward scatterers" are added to the field. Null collisions are only a statistical method to sample p(x) when it is not analytically invertible.

In our practice, we observe that there are cases in which switching from one viewpoint to the other leads to the designing of quite distinct Monte Carlo algorithms. This is particularly true when we work on the integral form of the algorithm. In this framework, any Monte Carlo algorithm can be translated into its mathematical counterpart (and reciprocally). This integral formulation can be rigorously transformed using simple operations, and translated back into a new algorithm with better properties (usually with reduced variance [START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF]) or into one that computes entirely new quantities (for instance, the partial derivatives (or sensitivities) of the initial quantity with respect to any model parameter [START_REF] De Lataillade | Monte carlo method and sensitivity estimations[END_REF][START_REF] Roger | Monte carlo estimates of domain-deformation sensitivities[END_REF]).

With the first viewpoint (the similitude), each null collision is a sampling event that appears in the integral formulation. With the second one (rejection sampling), the only retained events are the true A rejection-sampling viewpoint Let X and Z be two independent and identically distributed (i.i.d.) random variables with probability density function p X . Let Y be a random variable with probability density function p Y . If there is α > 1 such that ∀x ∈ R, α p Y (x) p X (x), then:

X = P (Y ) Y + (1 -P (Y )) Z where P (y) = pX (y) α pY (y) . x y p X(π X ) αp Y (π Y ) p Y (π Y ) • • • y 1 y 2 y 3 L x 0
(a) Rejection Sampling: {yn} n∈N is defined as a sequence of realizations of the random variable Y . The sequence {rn} n∈N is defined such that r i is sampled uniformly on [0, αp Y (y i )]. The rejection method consists of considering y i as a realization of X when r i < p X (y i ), in the other case y i is rejected. In the above example, y 1 and y 2 are rejected, while y 3 is retained.

Let X N CA be an exponentially distributed random variable of rate parameter k e . Let Z be a random variable i.i.d. with

X N CA .
Let Y be an exponentially distributed random variable of rate parameter k. If ∀x ∈ [x 0 , +∞[, k k e (x), then:

X N CA = P (Y ) Y + (1 -P (Y )) Z (Y )
where P(y) = k(y) k .

x k(y) Indeed, one can notice that X and XNCA have the same mathematical expression. In both cases, one realization of the random variable X or XNCA is seen as a series of realizations of another random variable Y , each of which are successively rejected until one is retained. The difference lies in the random variable XNCA itself, which has to be exponentially distributed, while X may correspond to any random variable. This major difference allows for the use of the NCA, in which each subsequent realization of Y makes use of the previous realizations. In terms of practical use, while a realization yi is only used to reject or retain the event yi in the case of the rejection sampling method, it is also used to evaluate the next realization yi+1 in the case of NCA. This also means that every sub-realization yi of XNCA improves the chances of the following step yi+1 being the final one.

collisions; null collisions do not appear at all in the integral formulation. The difference is very significant as soon as sensitivities are considered.

The principle of Monte Carlo estimation of parametric sensitivities (or Jacobians) in radiative transfer is quite simple [START_REF] De Lataillade | Monte carlo method and sensitivity estimations[END_REF][START_REF] Roger | Monte carlo estimates of domain-deformation sensitivities[END_REF]. Any Monte Carlo algorithm can be expressed under an integral form. This integral can be derived as a function of any parameter to give another integral formulation that can again be evaluated using Monte Carlo.

The sensitivity evaluation is simultaneous as soon as the very same random sampling algorithm can be used for both the addressed quantity and its sensitivity. Roger et al. [START_REF] Roger | Monte carlo estimates of domain-deformation sensitivities[END_REF] have shown that this is always possible, regardless of the parameter type, and even when the parameter affects the integration domain. However, despite being thoroughly general, this proposition can cause practical difficulties in some contexts, in particular as far as numerical convergence is concerned [START_REF] Tregan | Convergence issues in derivatives of monte carlo null-collision integral formulations: A solution[END_REF].

An example of such difficulty can easily be pictured using a simple example. Let us assume that we want to evaluate the direct transmissivity of a semi-transparent column of length L starting at the position x 0 (Beer-Lambert law for heterogeneous medium), as displayed in Figure 8.

T dir (x 0 ) = exp - x0+L x0 k e (y)

dy

Using a NCA, a collision location x will be sampled as if the extinction coefficient k e were uniform and equal to k = k e + k n , and rejected with probability

P N (x) = k n (x) k = k -k e (x) k
Retaining the similitude point of view for the NCA, one would write the following formulation, which is demonstrated in [START_REF] Longo | Direct derivation of skullerud's monte carlo method for charged particle transport from the linear boltzmann equation[END_REF] (see figure 9):

T dir (x 0 ) = +∞ 0 dx k exp -kx H(L -(x -x 0 )) P N (x)T dir (x)+(1 -P N (x)) {w 0 } + H((x -x 0 ) -L) {w 1 } (5)
Where w 0 = 0 represents the null contribution to the transmissivity of an absorption in the medium, and w 1 = 1 the weight associated with a contribution to the transmissivity. To evaluate the sensitivity of the transmissivity to any parameter ζ on which depends k e , the integrand in eq. ( 5) is derived with respect to ζ. Then, in order to be able to use the same samples for T dir and its derivatives, thus making the computation of these quantities simultaneous, the derivative formulation must be rewritten under the same form as eq. ( 5). This yields the logarithmic derivative of P N with respect to ζ in the sensitivity Monte Carlo weight:

∂ ζ P N P N = - ∂ ζ k e k -k e
where ∂ ζ k e is the derivative of k e with respect to ζ (see Lataillade et al. [START_REF] De Lataillade | Monte carlo method and sensitivity estimations[END_REF]). It can be easily seen that -∂ζ ke k-ke tends towards infinity when k tends towards k e (for more details, see Tregan et al. [START_REF] Tregan | Convergence issues in derivatives of monte carlo null-collision integral formulations: A solution[END_REF]), that is, when the probability of null collisions decreases. For this reason, the variance of the sensitivity estimator becomes infinite and convergence is impossible.

This means that when k is well adjusted (which is required in order to reduce the number of useless null collisions), the sensitivity estimation becomes less accurate. Figure 8 illustrates this convergence difficulty.

It also illustrates a solution that could easily be found by changing the viewpoint on null-collisions, which is described in thorough detail in Tregan et al. 2020 [START_REF] Tregan | Convergence issues in derivatives of monte carlo null-collision integral formulations: A solution[END_REF]. In the first case, the null-collision probability needs to be derived, leading to the problematic -∂ζ ke k-ke mentioned above; in the second case, only p(x) needs to be derived, in which k does not appear, and the convergence difficulty associated with -1 k-ke vanishes (see Figure 8).

Another illustration of sensitivity computation can be found on the previously presented maps of ground-fluxes (Figure 4b). In order to better understand and characterize cloud radiation interactions, the ground-flux density was computed along a line on the ground (see Figure 4c, along with the sensitivity to 1ω, where ω is the single scattering albedo (ω is assumed uniform within the clouds). These two quantities were simultaneously computed, using the exact same technique as described above to avoid convergence difficulties.

Transforming the nonlinearity of the exponential into a linear problem of infinite dimension

The principle underlying this third point of view is exposed in detail in Dauchet et al. [START_REF] Dauchet | Addressing nonlinearities in Monte Carlo[END_REF] and described in Figure 9. The departure point is that double randomization can no longer be used as soon as integrals are combined through nonlinear functions. The first idea is that nonlinear functions can be developed into infinite sums of monomials using a Taylor series. The second idea is that a monomial of order n of a random variable X can be written as the product of n independent variables X i identically distributed as X. To compute one realisation of X n , one only needs to compute n samples distributed according to the probability density function of X and retain their product. The entire algorithm corresponding to the evaluation of a Taylor series consists in sampling n, the order of the monomial, and then compute one realisation of X n , thus retrieving the power of double randomization.

In this third viewpoint, the nonlinearity of the exponential is read and treated with this exact idea in Bottom figures In the particular case where the extinction coefficient ke is uniform, number of samples required for a 1% accurate evaluation of the sensitivity to ke of the column transmissivity, using a null-collision algorithm. In the left-hand figure, we start from an integral formulation in which null-collisions are viewed as forward-scattering events, which leads to convergence difficulties when kn is close to 0 (the number of required samples tends to infinity). In the right-hand figure, we start from an integral formulation in which null-collisions are viewed as a rejection-sampling method: the convergence difficulty vanishes.

From a nonlinear to an infinite-dimensional problem

According to the Beer-Lambert law, the direct transmissivity of a semi-transparent column of length L can be expressed as equation 6 in the case of a heterogeneous medium of absorptivity coefficient k e (x). With the introduction of an overestimate k as explained in Fig. 1, and a Taylor expansion of the exponential, this probability can be written as equation 7.

T dir (x 0 ) = exp -

(a) A Taylor expansion approach

According to Beer-Lambert law, the direct transmissivity of a semi-transparent column of length L can be expressed as equation 8 in the case of a heterogeneous medium of absorptivity coefficient k e (x). This transmissivity can be evaluated by making use of the null-collision algorithm presented in Fig. 1, which is written in its integral form in equation 8. By separation of the integrals with respect to the Heaviside functions (equation 9), the transmissivity can be rewritten as equation 10 (see full details in Appendix A)

T dir (x 0 ) = ∞ 0 dy 1 k exp(-ky 1 ) H(y 1 -L){1} + H(L -y 1 ) k e ( x1 x 0 + y 1 ) k {0} + 1 - k e (x 1 ) k T dir (x 1 ) (8) = ∞ L dy 1 k exp(-ky 1 ) + L 0 dy 1 k exp(-ky 1 ) k e (x 1 ) k {0} + 1 - k e (x 1 ) k T dir (x 1 ) (9)
= exp(-kL) + L-y1 0 dy 2 kk e (x 2 ) ...

L•••-yn-1 0 dy n k -k e (x n ) (10)
With the notation P iCN the probability for the null-collision algorithm to do strictly i null-collisions before fulfilling the condition x i > x 0 + L (which means a positive contribution to the transmissivity), both approaches can be expressed as the product of P 0CN (= exp(-kx)) and a series of lookalike integrals.

In addition, both expressions can be further developed to prove a term by term equality (see [START_REF] Longo | Direct derivation of skullerud's monte carlo method for charged particle transport from the linear boltzmann equation[END_REF] for more details), which can be interpreted as

1 n! P1CN (XNCA>x) P0CN (XNCA>x) n = PnCN (XNCA>x) P0CN (XNCA>x) .
(b) The null-collision approach Figure 9: Key ideas on the transformation of the exponential's nonlinearity of the direct transmissivity into an infinite-dimension linear problem. Textbox (a) shows the regular approach to transforming such a non-linearity into an infinite-dimension linear problem by making use of Taylor series. Textbox (b) shows how the ACN for the same transmissivity can also be transformed into an infinite sum of lookalike terms. However, a term-by-term comparison shows that both approaches do not lead to the same expressions, and therefore interpretations. The mathematical proof that both expressions are equal can be found in [START_REF] Longo | Direct derivation of skullerud's monte carlo method for charged particle transport from the linear boltzmann equation[END_REF]. The approach (in textbox (a), eq. ( 7)) is not implemented here for its lack of performance in terms of CPU time and convergence, compared to ACN (in textbox (b), eq. ( 8)).

thickness of a line segment of length x. f can be expanded in its Taylor series around τ = 0:

p(x) = k e (x) +∞ n=0 (-1) n n! x 0 k e (x )dx n = k e (x) 1 - x 0 dx 1 k e (x 1 ) + 1 2 x 0 dx 1 x 0 dx 2 k e (x 1 )k e (x 2 ) - 1 6 x 0 dx 1 x 0 dx 2 x 0 dx 3 k e (x 1
)k e (x 2 )k e (x 3 ) ...

In Monte Carlo terms, the nonlinearity of the exponential has vanished and been replaced by an infinite sum over the order n and an increase up to n of the dimension of the integration domain at each order.

Consequently, a Monte Carlo algorithm evaluating p(x) would simply sample a value of n and then sample n values of x 1 , x 2 ... x n on the [0, x] interval and retain the product of the n values of k e evaluated at

x 1 , x 2 ... x n .
Of course, the purpose of a radiative-transfer Monte Carlo algorithm is seldom to evaluate p(x)

and several modifications to this initial idea are needed before reconstructing a standard null-collision algorithm. These ideas are all detailed Longo [START_REF] Longo | Direct derivation of skullerud's monte carlo method for charged particle transport from the linear boltzmann equation[END_REF] and here we only state the most important: the removal of the sign alternation. If the above decomposition was used as such, the Monte Carlo weight would be positive when sampling even values of n and negative for odd values, and evaluating the average of these positive and negative weights would require a large number of samples (the variance would be large). To avoid this, the form of p(x) is first modified in the following way:

p(x) = k e (x) exp(-kx)g( kx -τ ) with g(u) = exp(u),
where k is an overestimate of k e at all locations. The function f was a negative exponential applied to τ positive. The function g is a positive exponential applied to kxτ , again positive. When expanding g, the Monte Carlo weight remains positive for all values of n.

The fundamental point here is that the integral structure is linearized and that standard Monte Carlo approaches can be used to evaluate it without further reformulation (see Figure 9). Since k e is now outside of the nonlinear exponential, it can in turn be written as an integral or a discrete sum, while double randomization continues to ensure that the algorithm is practicable.

In Galtier et al. (2016) [START_REF] Galtier | Radiative transfer and spectroscopic databases: A line-sampling monte carlo approach[END_REF], this was used to address the question of handling large spectroscopic databases for molecular gases: it was no longer required that the contribution of all lines be first added to construct the absorption coefficient at a given wavenumber; as this sum was shifted out of the exponential, it could be handled by the Monte Carlo algorithm itself, by sampling line transitions and considering the sampled line's contribution to the local absorption only.

In Galtier et al. (2107) [START_REF] Galtier | A symbolic approach for the identification of radiative properties[END_REF] the same approach allowed the application range of so-called "symbolic"

Monte Carlo algorithms to be significantly extended. Consider a radiative transfer problem in which a quantity A needs to be evaluated knowing a set of field parameters π 1 , π 2 ... that determine the temperature and the optical properties. Symbolic algorithms do not only compute A for a given set of parameter values e.g.,

π 1 = π comp 1 , π 2 = π comp 2 .
.., they also evaluate the coefficients of a functional form that makes it possible to later evaluate A for any other value of π 1 , π 2 .... This was claimed to be very difficult for the parameters defining the absorption and scattering coefficients of a heterogeneous field because they appear within the exponential of Beer's law. Galtier et al. [START_REF] Galtier | A symbolic approach for the identification of radiative properties[END_REF] made a successful use of the above-described

Taylor expansion to solve this difficulty.

Conclusion

In this text, ongoing and published research that makes use of the null-collision method has been classified into three families that reflect different readings of the same concept. The similitude viewpoint is the most intuitive one for it is related to physics, and originally opened the door to computational efficiency in heterogeneous media. The rejection approach compares the null-collision method to a classical rejection sampling method and has proved useful in handling variance issues in the estimation of partial derivatives.Finally, Taylor's expansion, initially used to handle nonlinear physics, fundamentally allows the shifting of the optical depth integral out from the nonlinear exponential. This opens up many perspectives, for instance through the combination of transport with other physics underlying the computation of the optical properties: uncertainty and sensitivity estimations can thenceforth be propagated one step further than transport.

The other examples described in Dauchet et al. [START_REF] Dauchet | Addressing nonlinearities in Monte Carlo[END_REF] go at least partially beyond radiative transfer.

One is related to the solution of Maxwell's equations with the objective of evaluating the absorption and scattering properties of an ensemble of complex-shape particles. Another one deals with the extension to nonlinear transport (the Boltzmann equation) of the null-collision algorithms that are described in the present text. Two further examples illustrate the potential of these approaches, which we believe might be widely used in future research: the coupling of radiative transfer with other processes inside a single Monte Carlo algorithm. This idea has already been presented in Fournier et al. [START_REF] Fournier | Radiative, conductive and convective heat-transfers in a single monte carlo algorithm[END_REF] for coupling radiation with other heat-transfer modes, but the coupling was linear. In Dauchet et al. [START_REF] Dauchet | Addressing nonlinearities in Monte Carlo[END_REF] solar radiation within a concentrated solar plant is coupled with a nonlinear photochemical conversion process, and, very similarly, solar radiation within a photobioreactor is coupled with a nonlinear photosynthesis process. In these examples, there is only a one-way coupling but it is nonlinear. Coupling radiation, conduction and convection heat-transfers was effective, but it was in a linear context. Altogether, we do not see any conceptual hurdle forbidding the design of Monte Carlo algorithms that would couple radiative transfer with other physical, chemical or biological processes, including nonlinear ones.

Figure 1 :

 1 Figure1: Null-collision methods seen under the physical viewpoint where a purely forward-scattering species is added to the medium to make it homogeneous.

Figure 2 :

 2 Figure 2: a) Regular tracking, where the original medium is the support of the path-tracing algorithm. b) Null-collision methods, where the field has been virtually homogenized, yielding data -algorithm independence. Adapted from Figure 2 of [27].

Figure 3 :

 3 Figure 3: Examples of hierarchical grids for fast path-tracing. Above: 3D view of a hierarchical grid for a cumulus clouds scene. Below: 2D cross-section of a hierarchical grid for a congestus clouds scene.

  (a) A cloudy scene simulated using the high-resolution Large Eddy Simulation (LES) tool of Météo-France, rendered with a tool based on Monte Carlo path-tracing. The camera is located in the bottom right corner of the field. (b) Solar fluxes at ground level and associated uncertainties (standard deviations), as simulated by a multiple-scattering Monte Carlo algorithm using null-collisions to handle the geometric complexity of the 3D clouds.(c) Cloud optical depth, simulated solar-fluxes and sensitivities to these fluxes to 1ω, along a line across the projection of a given cloud, for two values of ω (left ω = 0.95, right ω = 0.99). The uncertainties are depicted as an envelope of width 3 times the standard deviation. The variance of the sensitivity evaluation increases when ω (the single scattering albedo) gets close to 1, i.e. when reducing absorption.

Figure 4 :

 4 Figure 4: Using null collision algorithms implemented in Monte Carlo codes to simulate radiative transfer in highly heterogeneous cloudy atmospheres.

  rendering time of scenes of increasing complexity Relative rendering time [%] Number of triangles in scene

Figure 5 :

 5 Figure5: Insensitivity of computing time to scene complexity is achieved when using hierarchical grids to organize the surface data. Reproduced from Figure1of[START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF].

  NCA viewpoint: NCA consists of sampling a realization y of the random variable Y , and to sample uniformly r in [0, k]. If r < k(y)/ k, y is retained as a realization of X; otherwise, a new realization ỹ of the random variable Y(y) (defined on [y, +∞[) is sampled, which is therefore dependent on the previous realization. In the above example, y 1 , y 2 and y 3 are rejected, while y 4 is retained.

Figure 7 :

 7 Figure 7: Textboxes (a) and (b) highlight the close link between the rejection sampling method and the NCA.Indeed, one can notice that X and XNCA have the same mathematical expression. In both cases, one realization of the random variable X or XNCA is seen as a series of realizations of another random variable Y , each of which are successively rejected until one is retained. The difference lies in the random variable XNCA itself, which has to be exponentially distributed, while X may correspond to any random variable. This major difference allows for the use of the NCA, in which each subsequent realization of Y makes use of the previous realizations. In terms of practical use, while a realization yi is only used to reject or retain the event yi in the case of the rejection sampling method, it is also used to evaluate the next realization yi+1 in the case of NCA. This also means that every sub-realization yi of XNCA improves the chances of the following step yi+1 being the final one.

keke

  mind. Let us consider p(x) = k e (x) exp -x 0 k e (x )dx and introduce the function f (u) = exp(-u). Then p(x) = k e (x)f (τ ) with τ = x 0 k e (x )dx the optical Number of samples for 1% accuracy Null-Collision then Sensitivities k e l = 0.01 k e l = 0.1 k e l = 0.5 k e l = 1 k e l = 4 k e l = 9 Number of samples for 1% accuracy Sensitivities then Null-Collision k e l = 0.01 k e l = 0.1 k e l = 0.5 k e l = 1 k e l = 4 k e l = 9 M M C

Figure 8 :

 8 Figure 8: Top figure Profile of the extinction coefficient ke within a heterogeneous column of length L.Bottom figures In the particular case where the extinction coefficient ke is uniform, number of samples required for a 1% accurate evaluation of the sensitivity to ke of the column transmissivity, using a null-collision algorithm. In the left-hand figure, we start from an integral formulation in which null-collisions are viewed as forward-scattering events, which leads to convergence difficulties when kn is close to 0 (the number of required samples tends to infinity). In the right-hand figure, we start from an integral formulation in which null-collisions are viewed as a rejection-sampling method: the convergence difficulty vanishes.

  k e (x )n

L 0 dy 1 1 )dy 1 k -k e (x 1 )dy 2 kdy 1 k

 0111121 exp(-ky 1 ) kk e (x 1 ) T dir (xk e (x 2 ) k e (x 1 )

0 dy 1 dy 1 0 dy 1 0 dy 1 0 dy 1 0 dy 2 0 dy 1 dy 2 0 dy 1 k 0 dy 1 dy 2 0 dy 1 k 0 dy 1 dy 2 0 dy 1 k 0 dy 1 0 dy 1 k -k e (x 1 )dy 2 k -k e (x 2 )L-y1-y2 0 dy 3 kdy 1 k -k e (x 1 )dy 2 k 1 0

 011010101020120101201012010101122031121 This work has been sponsored by the French government "Investissements d'avenir" research program, through the ANR SOLSTICE Laboratories of Excellence program (ANR-10-LABX-22-01), the ANR HIGH-TUNE grant (ANR-16-CE01-0010) and the Occitanie Region.T dir (x 0 ) = ∞ k exp(-ky 1 ) H(y 1 -L){1} + H(Ly 1 ) k exp(-ky 1 ) k e (x 1 ) k {0} + 1 -k e (x 1 ) k T dir (x 1 ) = exp(-kL) + L k exp(-kx 1 ) kk e (x 1 ) k T dir (x 1 ) = exp(-kL) + L exp(-ky 1 ) kk e (x 1 ) T dir (x 1 ) = exp(-kL) + L exp(-ky 1 ) kk e (x 1 ) ∞ k exp(-ky 2 ) H(y 1 + y 2 -L){1} + H(Ly 1y 2 ) k e ( x2 x 1 + y 2 ) k {0} + 1 -k e (x 2 ) k T dir (x 2 ) = exp(-kL) + L exp(-ky 1 ) kk e (x 1 ) k exp(-ky 2 ) k e (x 2 ) k {0} + 1 -k e (x 2 ) k T dir (x 2 ) = exp(-kL) + L k e (x 1 ) ∞ L-y1 dy 2 k exp(-k(y 1 + y 2 )) + L exp(-ky 1 ) kk e (x 1 ) k exp(-ky 2 ) k e (x 2 ) k {0} + 1 -k e (x 2 ) k T dir (x 2 )= exp(-kL) + exp(-kL) L k e (x 1 ) + L exp(-ky 1 ) kk e (x 1 ) exp(-ky 2 ) kk e (x 2 ) T dir (x 2 ) = exp(-kL) + exp(-kL) L k e (x 1 ) + exp(-kL) L exp(-ky 1 ) kk e (x 1 ) L-y1 0 dy 2 exp(-ky 2 ) kk e (x 2 ) + L k e (x 3 ) e -k(y1+y2+y3) T dir (x 3 ) = exp(-kL) k e (x 2 ) ... L•••-yndy n kk e (x n ) (A.2)

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: