# Application of null-collision backward Monte Carlo algorithm to digital image rendering of sooting flames in the visible range

Mouna El Hafi, Vincent Eymet, Vincent Forest, Richard Fournier

► To cite this version:

Mouna El Hafi, Vincent Eymet, Vincent Forest, Richard Fournier. Application of null-collision backward Monte Carlo algorithm to digital image rendering of sooting flames in the visible range. Journal of Quantitative Spectroscopy and Radiative Transfer, 2025, 347, pp.109661. 10.1016/j.jqsrt.2025.109661. hal-05269590

HAL Id: hal-05269590

https://imt-mines-albi.hal.science/hal-05269590v1

Submitted on 22 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Application of null-collision backward Monte Carlo algorithm to digital image rendering of sooting flames in the visible range

Mouna El Hafi<sup>a, c, e, *,</sup> Vincent Eymet<sup>b</sup>, Vincent Forest<sup>b</sup>, Richard Fournier<sup>c</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Null collision

Hierarchical grids

Digital image

Sooting flame

Monte Carlo

Meshless method

## ABSTRACT

The objective of this work is to present a new tool to produce a digital image computation inside a sooting flame using Monte Carlo methods. The studied case consists of a laser emission corresponding to a spatially localized source of collimated radiation, that partially illuminates the environment. The laser emission consequently causes problems of convergence with a Monte Carlo method due to the difficulty in statistically linking sensors to sources. A new algorithm has been proposed to calculate the radiative intensity at a given probe position. An open source tool has been developed based on a combination of Null Collision Algorithms (NCA) with hierarchical grids technics from the computer graphics community, allowing to perform efficient radiative transfer (and possibly produce images). These tools are generic and are able to accept any input data such as, for instance, temperature and concentration fields produced by a Large Eddy Simulation (LES) approach over a tetrahedral grid. We show in this paper a first monochromatic digital image of a laminar sooting flame in the visible range. Such a new tool can for instance provide a useful support for experimental characterization of the radiative behavior of soot aggregates in flames, such as tracking effects of scattering.

## 1. Introduction

During the last decades, Monte Carlo methods have shown their capability to treat complex geometries and to take into account complex physical phenomena (reflective surfaces, scattering) [1] in radiative heat transfer. New developments conducted these last years for heat transfer, combine Monte Carlo with computer graphics libraries [2–4], which accelerates ray tracing in both surfaces and volumes. These developments are the equivalent of recent advances in the cinema industry that now make possible the use of Monte Carlo Methods to render cloudy scenes [5]. Concerning work at the interface between statistical physics and computer graphics, in addition to the field of combustion presented here, recent research has been developed in atmospheric physics [3] and in the planetary atmospheres of Titan [6]. The tools make it possible to compute relevant radiative quantities and to produce rendered images. Until then, using acceleration structures to efficiently identify the intersections between rays and surfaces was very common in the field of computer graphics. But the use of grids to sample absorption and scattering events within volumetric domains constitutes a recent and significant development in the field of computer graphics, as shown in the work of Kutz [5] and Burley [7]. This is due to the non-linearity of Beer's law: when dealing with spatial

heterogeneities, the optical thickness, which is the main variable in radiative transfer, varies unpredictably along a given line of sight. If statistically sampling the optical thickness is simple, the identification of the corresponding spatial position is a much more difficult task. In 2013, Galtier [8] proposed a trick using the null collision concept [9, 10]. Disney Research then revisited this Null Collision Algorithms (NCA), that was quickly adopted as a usual practice [10,11].

This idea allows to bypass the non-linearity of the Beer's law by artificially introducing a new type of collisioners, that merely scatter radiation (without any absorption), and only in the forward direction. The total extinction field, which is the sum of the real (physical) extinction, and of this null-scattering (or null-collision) field, has to be uniform through the semi-transparent medium. This medium is therefore optically thicker than the original, and in compensation some of the collision events that are named “null” collisions are ignored because they stand for scattering events in the forward direction. All it takes is a known upper-bound value of the real extinction volumetric field. The underlying integral formulation is rigorous: any null-collision field can be added to the real extinction. The resulting algorithm will converge toward the same result. However, the null-collision field will have an effect on the performance of the algorithm; for instance,

<footnote><sup>a</sup> Centre RAPSODEE, Université de Toulouse, Mines Albi, UMR CNRS 5302, Campus Jarlard, Albi, 81013 Cedex 09, France</footnote>

<footnote><sup>b</sup> Méso-Star, 4 rue du général Giraud, 31200 Toulouse, France</footnote>

<footnote><sup>c</sup> Laboratoire LAPLACE, Université de Toulouse, UMR 5213, 118 route de Narbonne, 31062 Toulouse Cedex 9, France</footnote>

<footnote>* Corresponding author at: Centre RAPSODEE, Université de Toulouse, Mines Albi, UMR CNRS 5302, Campus Jarlard, Albi, 81013 Cedex 09, France. E-mail address: mouna.elhafi@mines-albi.fr (M. El Hafi).</footnote>
whenever the upper-bound value of extinction is excessively high, the fraction of null-collisions becomes prohibitively high, which increases the computational time. The principle of the null-collision technique, and its limitations, are presented in Appendix B.

The solution is to partition space in voxels (volume elements) with possibly different values of the upper-bound values. The initial medium is therefore subdivided into a number of homogeneous volume elements. But this is not enough, since a large number of voxels are also a cause for performance degradation : a lot of time is now spent trying to cross the grid. An efficient volume partitioning strategy consists of refining the acceleration grid in regions of strong total extinction gradients. The resulting acceleration procedure quickly identifies a true absorption position, which reduces the computational cost.

Combining NCA and adaptive grids opens the door to making no compromise on the resolution of the underlying data, and to tackle problems involving large scale ratios. N. Villefranque [3] showed that, by using adaptive grids based on an unitary upper-bound optical thickness criterion, the time to render the image becomes insensitive to the resolution of cloud fields. The idea to extend this methodology to the rendering of a sooting flame in combustion applications is the objective of this paper.

The studied case consists of a laser sheet that illuminates the com-
bustion medium. As we use a backward Monte Carlo algorithm, optical
paths are generated from the sensor (cameral pixels) and should be
stopped when reaching the source, in a possible emission direction. This
is a problem when the source has a small spatial extension and/or when
emitted radiation is collimated. In the present case, both difficulties
are present. The associated difficulties are presented in Appendix C.
An extension of NCA recursive integral formulation to handle laser-like
sources that partially illuminates the scene has been proposed in [12].
The methodology, that was initially developed for atmospheric radia-
tion [3] and that was later extended to flames using this new algorithm,
has raised new implementation difficulties that are presented hereafter.
The input data for the radiative transfer code is generated by LES
combustion code AVBP [13] which provides soot volumetric fractions
over a unstructured tetrahedral grid. This required an important work
of algorithmic and computer development that will be explained in
this paper, and that represents a crucial step in the development of
these rendering tools. The following section emphasizes the difficulties
caused by the laser source and recalls the key points of the new
Monte Carlo algorithm to overcome the convergence difficulties and to
produce the acceleration grid. The last section provides results in the
case of a laminar sooting flame in the form of a digital image in the
visible range, generated using Monte Carlo Method.

## 2. Algorithm and implementation in *hrdr-combustion*

The development of the *hrdr* (High-Tune RenDeR<fnref n="1" />) radiative
transfer simulation tool began in 2018. Its initial purpose was to render
visible images of a scene composed of a vertically stratified atmosphere,
a plane gray ground, and clouds whose properties were provided over a
rectangular volumetric grid. It quickly became a versatile tool with the
possibility to produce not only images in the visible and the thermal
infrared range, but also flux maps in any spectral range for radiation
sources both internal and external to the participating medium, and the
possibility to take into account an arbitrary number of solid surfaces of
any shape and of any spectral reflectivity or emissivity. This code is
based on Monte Carlo algorithms that make possible to compute the
intensity at any location of a scene, in any direction, both for radiation
that is emitted by the medium and for radiation that is emitted by a
external source and interacts with the participating medium.

The legacy part of the code that deals with a plane ground and a
vertically stratified atmosphere has since been named *hrdr-atmosphere*,

and a more recent branch named *hrdr-combustion* has been developed
as a diagnostic tool for sooting flames that are partially illuminated by
a laser sheet. In this work, we propose to present a new feature of the
HTRDR code dedicated to combustion applications.

In the visible range, unlike soot, the gaseous species do not radiate
in combustion applications. To account for the complex morphology
of soot aggregates, Rayleigh–Debye–Gans theory for Fractal Aggregates
(RDG-FA) has been implemented in the radiative transfer algorithm
to determine absorption coefficient, scattering coefficient and phase
function [14].

Since the problem consists of simulating radiation that is emitted
by a laser sheet, scattered and absorbed by the combustion medium,
and that finally reaches a detector, the development of *hrdr-combustion*
was focused on the Monte Carlo algorithm that deals with external
radiation. This algorithm is commonly referred as the *shortwave* (SW)
radiation algorithm, even if external radiation has no particular reason
to be limited to the shortwave part of the spectrum. It turns out
that the SW radiation algorithm in *hrdr-atmosphere* [3] would be very
inefficient when applied to a medium that is only very partially directly
illuminated : in order to compute the intensity associated with external
radiation (as presented in *Appendix C*), this algorithm schematically
needs to simulate a number of multi-scattering trajectories within
the medium, starting at the position where the intensity needs to be
computed (the so called probe position). A trajectory is composed of a
number of collision positions. These collisions can be of three different
types: an absorption will end the trajectory; a scattering will result in
a modification in the propagation direction; and so called *null-collisions*
are events that are introduced by the final integral solution to the
radiative transfer equation that is the basis to the Monte Carlo algo-
rithm [8,15]. Using this technique makes possible to take into account
heterogeneous and anisothermal participating media while remaining
perfectly rigorous. Null collision events translate in no modification
to either the propagation direction of the trajectory, or its energetic
weight.

The legacy SW algorithm used in *hrdr-atmosphere*, needs to update
the Monte Carlo weight associated with a trajectory at each scattering
event: in particular, it needs to access the fraction of the source inten-
sity that reaches the current scattering position directly, without any
absorption or scattering. Needless to say, the current scattering position
shall be illuminated by the source; otherwise the local increment to the
Monte Carlo weight is null. This algorithm performs well in the case of
atmospheres, since all positions in the medium are directly illuminated
by the Sun, even at the center of an optically thick cloud. However,
in the case of a combustion medium that is illuminated by a very thin
laser sheet, only a small fraction of scattering positions receive direct
radiation. This SW algorithm would therefore perform very poorly in
such a situation: contributions to the Monte Carlo weight would be null
most of the time, with only very seldom positive contributions. The
convergence qualities (variance associated with the estimated intensity)
of such a algorithm would be catastrophic.

A new SW algorithm had therefore to be designed from the ground
up for media that are only very partially illuminated. It was presented
in [12]. The key point of this algorithm is that most positions in the
medium do not receive any direct radiation: the intensity is generally
due to radiation that was scattered at least once. And starting from
the source, the first scattering position necessarily takes place within
the slab of medium that is illuminated by the laser sheet. This new
SW algorithm is then based on a reformulation of the diffuse part
of intensity: it still needs to build reverse multi-scattering trajectories
(from a probe position), but the Monte Carlo weight, that is updated at
each scattering position, now needs to access the fraction of radiation
that was emitted by the laser sheet and was scattered exactly once by
the medium in the direction of the current scattering position along
the reverse path. This algorithm results in the simulation of two ways
paths between the laser source and the detector (to connect the probe
position to the laser source). Starting from a given probe position for

<footnote n="1">¹ https://www.meso-star.com/projects/hrdr/hrdr.html.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Acceleration grids: in order to compute the intensity at position $\vec{x}$ in direction $\vec{u}$, a line of sight must be initiated from $\vec{x}$ in the $-\vec{u}$ direction. The variation of $k(\vec{x} - \sigma\vec{u})$ along this line of sight is represented (and unpredictable by the algorithm) in the presence of a local maximum of $k$. In situation (a), a single (very high) upper bound value $\hat{k}$ is used over the whole domain. Resulting free paths are therefore very small, and a lot of null collisions occur. In (b), a very fine acceleration grid is surimposed, with relevant values of $\hat{k}$ in each cell, but the cost of crossing this grid is prohibitive. The ideal situation is provided in panel (c): an adaptive acceleration grid provides a trade-off between the number of null collisions and the cost of grid crossing.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** The voxelization step consists of identifying the value of $\hat{k}$ over each cell of a fine regular grid. Radiative properties such as the absorption coefficient $k$, are initially known at the nodes of a input grid (a). In a first step, a bounding-box surrounding each grid element (triangle in 2D, tetrahedron in 3D) is drawn, and is associated with the highest value of $k$ among the nodes of this element. In a second step, a detailed grid (b) is surimposed over the whole domain, and each cell is attributed the highest value of $k$ among all bounding boxes that overlap for the current cell (possibly one).</figcaption></figure>

instance, a pixel on a CCD array, the Monte Carlo weight associated to a given trajectory will successively take into account the contributions of every scattering position along the optical path.

The reason why this algorithm has much better convergence qualities than its legacy counterpart is because the main source of radiation in this configuration is radiation that has been scattered once in the fraction of the medium that is directly illuminated by the laser: this contribution is roughly identical for all positions in the medium, which is why the variance associated with the simulated intensity is small.

As mentioned in [12] (and recalled in appendices), using the null collision technique requires the *a-priori* knowledge of an upper bound of the extinction coefficient field $k(\vec{x})$ within the semi-transparent medium, uniform over known spatial zones called voxels (volume elements). This upper bound field is noted $\hat{k}(\vec{x})$, or more commonly $\hat{k}$. At any given location $\vec{x}$ in the medium, the free path length $\sigma$ between two consecutive collision events is sampled according to $p(\sigma) = \hat{k} \exp(-\hat{k}\sigma)$, which means that the mean free path length is smaller when $\hat{k}$ increases. This is the reason why using a infinite value of $\hat{k}$ is not possible: the values of $\sigma$ would then be null, which would make the Monte Carlo method unpractical. Once a non-null free path length $\sigma$ has been sampled, the collision position $\vec{x} - \sigma\vec{u}$ is known, and the nature of the collision can then be determined: an absorption event according to probability $k_a(\vec{x} - \sigma\vec{u})/\hat{k}$, a scattering event according to probability $k_s(\vec{x} - \sigma\vec{u})/\hat{k}$ and a null collision according to probability $1 - (k_a(\vec{x} - \sigma\vec{u}) + k_s(\vec{x} - \sigma\vec{u}))/\hat{k}$ with $k_a(\vec{x})$ and $k_s(\vec{x})$ respectively the absorption and

scattering coefficient fields in the medium. Arbitrarily increasing the value of $\hat{k}$ in a given voxel has the following consequences: not only the mean free path will be smaller, which means a higher number of collisions will have to be sampled in order to sweep a given line of sight, but also the probability to sample a null collision will be higher (see Fig. 1a).

This induces a higher computation time, because in this case, the value of single $\hat{k}$ over the whole domain will only be adequate in the immediate vicinity of the position where $k(\vec{x})$ reaches its maximum; and in the rest of the volume, $k(\vec{x})$ can possibly be orders of magnitude smaller than $\hat{k}$, resulting in a high number of time consuming null collisions. An optimization process is thus necessary to find a good trade-off between selecting an appropriate upper bound and reducing the number of null collisions. Leveraging null collisions approach, we extended surface-based grid acceleration methods that are commonly used in the computer graphics community [16] to *volumes*. Such techniques allow to improve the efficiency of data access.

In highly inhomogeneous media, it is nonetheless possible to partition the volume into a number of voxels (partitions of homogeneous $\hat{k}$); this will translate into local values of $\hat{k}$ that are locally better suited. The volume partition and associated values of $\hat{k}$ is called an acceleration grid. In addition to the identification of collision positions, the Monte Carlo algorithm also has to evaluate the intersection positions between the current line of sight and this acceleration grid. Finding these positions also has a numerical cost; in the case when
<figure><img src="image_6.png" /><figcaption>(a) Each group of 8 voxels (4 in this 2D representation) is tested for merging against an arbitrary criterion over the optical thickness of the possible resulting voxel of higher order</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) The procedure is repeated until the final adaptive acceleration grid is produced</figcaption></figure>

<caption>**Fig. 3.** The octree is constructed by accepting or rejecting merging between each group of voxels recursively.</caption>

the acceleration grid is very detailed, although the number of null collisions will be minimal, the overall computation time might not be optimal due to the cost of grid crossing (Fig. 1b). A solution was found by producing an acceleration grid refined in zones of maximum gradient of $k(\bar{x})$, with larger voxels everywhere else (see Fig. 1c). This type of adaptative acceleration grid makes possible to reach the better compromise between a low number of null collisions and a fast grid crossing, which results in a smaller computation time (see [3]).

In terms of spectral model, the null-collision algorithm is inherently monochromatic (the reason for this is given in Appendix A). The *htrdr-combustion* code must therefore produce the acceleration grid for each quadrature point the spectral data is provided over. Spectral data is defined over a tetrahedral volumetric mesh (see Fig. 2a), and is used as follows in order to produce the adaptative acceleration grid:

* voxelization: the value of $\hat{k}$ over each cell of a detailed initial grid can easily be found by identifying the tetrahedra the cell intersects with, and retaining the maximum value of $k$ among the nodes of these tetrahedra (see Fig. 2b).

* then a recursive procedure is applied, initially using this detailed grid; at each step of this procedure, voxels of the input grid may be merged by group of 8 when an arbitrary criterion over optical thickness is met (see Fig. 3a). The procedure is repeated until no more merging can take place, and results in a final grid with varying cell dimensions are shown (see Fig. 3b).

The resulting acceleration grid may then be used by the Monte Carlo algorithm that needs to determine, for each statistical realization, the intersection positions between the current line of sight and the acceleration grid, and the value of $\hat{k}$ in each voxel. The choice of the upper bound $\hat{k}$ parameters when constructing an acceleration grid is arbitrary. The acceleration grid does not alter the physical solution of the radiative transfer problem; rather, it serves purely as a computational tool to accelerate simulations. It achieves this by providing an optimal balance between the number of null collisions and the time required to traverse the grid. Moreover, using image synthesis libraries such as Embree [17] often results in highly efficient intersection computations.

### 3. Configuration and results

The *htrdr-combustion* code was employed over a configuration provided by CERFACS; it consists of a laminar diffusion flame of ethylene [18] described by a tetrahedral mesh of 1.1 million nodes that define 6.7 million cells produced by the AVBP combustion solver. The enclosed medium dimensions are 9 × 9 × 35 cm. For each node of

this grid, the following properties are defined: pressure, temperature, concentration of $H_2O$, CO and $CO_2$, volumetric mass, soot primary particles diameter and number density, and soot volumetric fraction (see Fig. 4). In the present case, since *htrdr-combustion* only produces results for the shortwave region (the source of radiation is external to the medium), the gas mixture is considered as non-participating and the only relevant information is about soot. It is used to produce radiative properties of the medium (absorption and scattering coefficients, and discretized phase function) at runtime, for the unique wavelength of the laser source.

The images presented in Fig. 5 have been produced by *htrdr-combustion* in the following configuration: for the vertical rendering (left image), the laser sheet is vertical, and runs through the central axis of the medium. The main direction of observation (the camera) is perpendicular to the plane of the laser sheet, with the Z-axis as the up direction and a field of view of 30 degrees. The camera center is located at an altitude of 17 mm and its distance to the laser sheet is 12 mm. This image (left image) has a definition of 1280 × 500 pixels, and 4096 samples have been used in order to render each pixel, which means that more than 2.6 billion of optical trajectories have been simulated in order to produce this image. Renderings from above (right images): the laser sheet is horizontal, and intersects the medium at various heights. These images (right images) have a 1920 × 1280 resolution, and 4096 samples have been used in order to render each pixel, which makes a total of more than 10 billion optical trajectories per image. Concerning computation costs, for this flame configuration, which contains few heterogeneities, the computational cost ratio between a simulation with and without a grid acceleration, is around 10. This ratio increases and therefore the computational performance improves as the flame becomes more heterogeneous. In [15], for a single level grid we demonstrate the impact of optical thickness on the computation time required for grid traversal. In fact, the criterion that governs the construction of the grids is the optical thickness. When dealing with heterogeneous media, computation time would be prohibitively high if only a single level grid is used. Constructing a hierarchical grid makes it possible to handle such media with significantly reduced computation time.

The color map used to render these images indicates the intensity of the signal that reaches the camera (radiation that was emitted by the laser source, then scattered by soot particles in the medium before reaching a pixel of the camera). The brighter the color, the more intense is this signal. The corresponding spatial zones are therefore characterized by a higher scattering coefficient, which in turn indicates zones containing smaller soot aggregates or a higher density of soot aggregates. The vertical rendering shows the weakening of the laser
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** (a) volumetric fraction field of soot, multiplied by 10⁶. (b) Diameter of primary particles, in nanometers (c) Average number of primary particles per aggregate.</figcaption></figure>

sheet intensity: the laser sheet propagates from the left to the right. And the signal over the right part of the image is clearly reduced compared to the left part: incoming radiation is partly absorbed and scattered as soon as it enters the combustion medium (on the left side), which leaves less radiation to be scattered by the medium on the right side. The same effect can be seen on renderings from above: the laser sheet also

propagates from the left to the right, and the right part of the image always display a weaker signal than the left part.

4. Conclusions

This article presents the use of a new Monte Carlo algorithm that makes possible to render images of a sooting flame that is illuminated
<figure><img src="image_9.png" /><figcaption>Fig. 5. Left: vertical rendering of the medium; Right: rendering from above, at various heights. These images have been produced with a number of samples = 10⁷ to ensure a low variance over each pixel, to reduce the noise in the rendering.</figcaption></figure>

by a monochromatic laser sheet. The algorithm was initially presented in [12], and has been implemented in HTRDR (which was initially an image renderer). We are using an acceleration grid that provides a value of the upper bound of k at each spatial position. In theory, any acceleration grid can be used: it does not change the result of the algorithm. In practice, the acceleration structure that we used achieves the optimal computation time. In this paper, results are presented in terms of simulated images only. Since the underlying Monte Carlo method consists of estimating the intensity field, the code also computes radiative fluxes over a simulated detector of arbitrary size, position and opening angle: results have been extensively validated to this regard. Images can really be seen as a by-product of the radiative transfer simulation. In the present case, the simulated images could help to analyze measurements. The numerical simulation, however, provides a way to investigate dominant phenomena. For instance, the simulation can be stopped for an any given scattering order. This would provide insightful information when trying to identify the number of scattering events that must be taken into account in approximate models for inversion purposes of experimental data for instance. To sum up, two main points have been demonstrated:

- The operability of the algorithm [12] for a real laminar flame calculated by AVBP based on LES.

- A new feature of HTRDR code for combustion applications based on acceleration grids to render flame images in the visible range.

It should be noted that this approach can be extended to other types of flames by modifying the CAD input files and the concentrations of radiative species. Regarding future work, the question of applying the same techniques to radiation emitted by the system can be studied. The Monte Carlo algorithm for computing the intensity due to radiation emitted, absorbed and possibly scattered by the medium already exists, and could be plugged in place of the algorithm used for external radiation (see Appendix B). However, two issues arise when speaking of internal radiation: first, integrating over the whole thermal infrared spectral range would require dealing with a huge amount of spectral data: a full k-distribution set (over several hundred of narrowbands) would have to be produced and manipulated for every node of the unstructured volumetric grid. The required amount of memory can easily reach a hundred gigabytes, even for quite modest volumetric grids (that use a small number of nodes). This question could find a partial answer

by using out-of-core techniques, where spectral data would be loaded in memory only in the vicinity of the current ray position. The second issue concerns the need to pre-compute an acceleration grid for every quadrature point of every spectral narrowband. This would clearly represent a bottleneck in terms of CPU time, without even speaking about the associated additional memory requirements. A strategy for using a totally different acceleration technique has been discussed, but still needs to prove its practicability.

### CRediT authorship contribution statement

**Mouna El Hafi:** Writing – original draft, Supervision, Project administration, Methodology, Funding acquisition. **Vincent Eymet:** Writing – original draft, Validation, Software, Methodology, Conceptualization. **Vincent Forest:** Validation, Software, Methodology. **Richard Fournier:** Validation, Methodology, Formal analysis, Conceptualization.

### Declaration of competing interest

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.

### Acknowledgments

This work was supported by the Agence Nationale de la Recherche for the ASTORIA project under the grant number ANR-18-CE05-0015, and we would like to thank the CERFACS laboratory for providing the LES data and Jérôme Yon for the fruitful discussions regarding soot morphology models.

## Appendix A. General Monte Carlo algorithm

Eq. (A.1) provides the formulation of the radiative flux density $\phi(\vec{x})$ at position $\vec{x}$ over a surface of local normal $\vec{n}$, integrated over a solid angle $\Omega$ (possibly $2\pi$) and over the $[\nu_1, \nu_2]$ spectral interval (possibly $[0, +\infty]$), where $L_\nu(\vec{x}, \vec{u})$ is the monochromatic intensity at position $\vec{x}$ in direction $\vec{u}$.

$$ \phi(\vec{x}) = \int_{\nu_1}^{\nu_2} d\nu \int_{\Omega} d\vec{u} |\vec{u} \cdot \vec{n}| L_{\nu}(\vec{x}, \vec{u}) \quad (A.1) $$

This formulation serves as the basis for our Monte Carlo algorithm: the $[\nu_1, \nu_2]$ spectral interval is first divided into $N_b$ narrowbands of width $\Delta\nu_b$:

$$ \int_{\nu_1}^{\nu_2} d\nu = \sum_{b=1}^{N_b} \int_{\Delta\nu_b} d\nu \quad (A.2) $$

and since the CK spectral model will be used (if only, because input data is provided under the CK formalism) the integral over each narrowband is transformed into a quadrature of order $N_q$, characterized by a set of quadrature weights $\{\omega_q, q = 1, N_q\}$:

$$ \int_{\Delta\nu_b} d\nu \approx \sum_{q=1}^{N_q} \omega_q \quad (A.3) $$

Finally, discrete probability sets $p_B$ and $p_Q$ are introduced in order to sample the narrowband, then the quadrature index, and a probability density function $p_\Omega$ is also introduced in order to sample a direction; Eq. (A.1) then becomes:

$$ \phi(\vec{x}) = \sum_{b=1}^{N_b} p_B(b) \sum_{q=1}^{N_q} p_Q(q) \int_{\Omega} p_{\Omega}(\vec{u}) d\vec{u} \left\{ \frac{|\vec{u} \cdot \vec{n}| L_{b,q}(\vec{x}, \vec{u}) \omega_q}{p_B(b)p_Q(q)p_{\Omega}(\vec{u})} \right\} \quad (A.4) $$

The Monte Carlo algorithm is then straightforward from this formulation, since it now represents the expectancy of a random variable (called the Monte Carlo weight). This algorithm consists of simulating a large number $N$ of optical trajectories as follows:
* a narrowband index b is randomly sampled according to the discrete {$p_B(b), b = 1, N_b$} probability set

* a quadrature order q is sampled according to the discrete {$p_q, q = 1, N_q$} probability set

* a direction $\vec{u}$ is sampled over the $\Omega$ solid angle (defined by to the local normal $\vec{n}$) according to the $p_\Omega$ probability density function

* then the Monte Carlo weight of realization index i is : $\omega_i = \frac{|\vec{u} \cdot \vec{n}| L_{b,q}(\vec{x}, \vec{u}) \omega_q}{p_B(b) p_Q(q) p_{\Omega}(\vec{u})}$, where $L_{b,q}(\vec{x}, \vec{u})$ is the monochromatic radiance (for band index b and quadrature index q) at position $\vec{x}$, propagating in direction $\vec{u}$.

The sum of weights and of squared weights can then be used in order to evaluate an estimation of the radiative flux $\langle \phi \rangle$, but also of its variance $V$:

$$ \langle \phi \rangle = \frac{1}{N} \sum_{i=1}^{N} \omega_i \quad (A.5) $$

$$ V = \frac{1}{N-1} \left[ \frac{1}{N} \sum_{i=1}^{N} \omega_i^2 - \langle \phi \rangle^2 \right] \quad (A.6) $$

The same algorithm makes possible to evaluate the required quantity (here, a surfacic flux density) and its variance ; the standard deviation $s = \sqrt{V}$ can be interpreted as the numerical accuracy attached to the result : we know there is a 99.7% probability that the result is located in the $\langle \phi \rangle \pm 3$ s interval. The various probabilities that have been introduced in Eq. (A.4) do not affect the value of the result (the algorithm will converge toward the solution, no matter what probabilities are used), but only its variance. In other words, the choice of probabilities used in the samplings will condition the convergence speed of the algorithm.

The main task consists of the evaluation of the monochromatic radiance $L_{b,q}(\vec{x}, \vec{u})$. If its value was known, computing Monte Carlo weight would be straightforward. But in all generality, it is unknown. A specific algorithm is associated with the evaluation of the Monte Carlo weight for the monochromatic radiance. The associated algorithm will be detailed in Appendix B for a internal source of radiation (the thermal infrared range) and in Appendix C for a external source.

## Appendix B. Internal source

Developments of the monochromatic radiance ($p$ and $q$ subscripts will be omitted from now) are carried out in the case of an internal source of radiation. This will give the opportunity to explain the necessity for the null-collision technique.

In all generality, monochromatic radiance $L(\vec{x}_0, \vec{u}_0)$ in a infinite medium is the sum of all emitted and scattered radiation along the line of sight $(\vec{x}_0, -\vec{u}_0)$, and transmitted to $\vec{x}_0$:

$$ L(\vec{x}_0, \vec{u}_0) = \int_0^{+\infty} k_a(\vec{x}_0 - \sigma \vec{u}_0) L^0(\vec{x}_0 - \sigma \vec{u}_0) \exp\left(-\int_0^\sigma k(\vec{x}_0 - \sigma' \vec{u}_0) d\sigma'\right) d\sigma \\ + k_s(\vec{x}_0 - \sigma \vec{u}_0) \int_{4\pi} \Phi(\vec{u}_0, \vec{u}_1) L(\vec{x}_0 - \sigma \vec{u}_0, \vec{u}_1) d\vec{u}_1 \\ \times \exp\left(-\int_0^\sigma k(\vec{x}_0 - \sigma' \vec{u}_0) d\sigma'\right) d\sigma \quad (B.1) $$

with $k_a, k_s$ and $\Phi$ respectively the monochromatic absorption coefficient, scattering coefficient, and phase function. $L^0$ is the equilibrium (Planck) intensity, and $k = k_a + k_s$ is the total extinction coefficient. $\sigma$ is the curvilinear abscissa along the line of sight. The two terms on the right hand side respectively account for emitted and scattered radiation.

Eq. (B.1) can be reformulated in a single integral:

$$ L(\vec{x}_0, \vec{u}_0) = \int_0^{+\infty} k(\vec{x}_0 - \sigma \vec{u}_0) \exp\left(-\int_0^\sigma k(\vec{x}_0 - \sigma' \vec{u}_0) d\sigma'\right) d\sigma \\ \times \left[ \frac{k_a(\vec{x}_1)}{k(\vec{x}_1)} L^0(\vec{x}_1) + \frac{k_s(\vec{x}_1)}{k(\vec{x}_1)} \int_{4\pi} \Phi(\vec{u}_0, \vec{u}_1) L(\vec{x}_1, \vec{u}_1) d\vec{u}_1 \right] \quad (B.2) $$

This formulation could be the basis for the following Monte Carlo algorithm:

* since $p_\sigma(\sigma) = k(\vec{x}_0 - \sigma \vec{u}_0) \exp(-\int_0^\sigma k(\vec{x}_0 - \sigma' \vec{u}_0) d\sigma')$ is normalized over $[0, +\infty)$, it can be chosen as the probability density function used to sample a curvilinear abscissa $\sigma$. Once it has been sampled, position $\vec{x}_1 = \vec{x}_0 - \sigma \vec{x}_0$ can be found.

* then the nature of the event can be sampled according to the $[p_a, 1-p_a]$ probability set with $p_a = \frac{k_a(\vec{x}_1)}{k(\vec{x}_1)}$ the probability of absorption and $1-p_a = \frac{k_s(\vec{x}_1)}{k(\vec{x}_1)}$ the probability of scattering.

* in the case a absorption event is sampled, the weight of the realization is $L^0(\vec{x}_1)$ the blackbody radiance at position $\vec{x}_1$.

* in the case a scattering event is sampled, a new propagation direction after scattering $\vec{u}_1$ is sampled over the unit sphere according to the (normalized) phase function, and the Monte Carlo weight is $L(\vec{x}_1, \vec{u}_1)$. Since it is unknown but can still be expressed by Eq. (B.2), the algorithm loops for a new position and a new direction. As a result, a optical path is now constructed from the various scattering positions in the medium.

This algorithm is recursive, and stops only when an absorption event is sampled, or when a limit condition is reached (if radiance $L(\vec{x}_1, \vec{u}_1)$ is known).

The main issue concerns the first step of the previous algorithm: sampling a curvilinear abscissa. Sampling $\sigma$ according to $p_\sigma$ over $[0, +\infty)$ is equivalent to sampling a optical thickness $\tau$ over $[0, +\infty)$ according to the $\exp(-\tau)$ probability density function. However, a problem arises when a relation has to be established between this optical thickness and the corresponding curvilinear abscissa $\sigma$. Since $\tau = \int_0^\sigma k(\vec{x}_0 - \sigma' \vec{u}_0) d\sigma'$ and the evolution of $k$ as a function of $\sigma'$ is unknown, no bijection can be established between $\tau$ and $\sigma$. In other words, it would be very difficult to identify $\sigma$ once $\tau$ is sampled.

A possible solution to this issue is to introduce a null-collision field $k_n$. Null-collisions are defined as scattering events with a purely forward phase function : these events modify neither the weight of the algorithm (the energy carried along the optical path), neither the propagation direction. But the following term should now be added to relation (B.1) :

$$ k_n(\vec{x}_0 - \sigma \vec{u}_0) \int_{4\pi} \Phi^n(\vec{u}_0, \vec{u}_1) L(\vec{x}_0 - \sigma \vec{u}_0, \vec{u}_1) d\vec{u}_1 \exp\left(-\int_0^\sigma k(\vec{x}_0 - \sigma' \vec{u}_0) d\sigma'\right) d\sigma \quad (B.3) $$

Note that when the (positive) $k_n$ null-collision field is introduced, the total extinction coefficient of the medium is now $k_a + j_s + k_n$; and we chose the $k_n$ field so that the total extinction coefficient is uniform throughout the medium. It will be denoted as $\tilde{k} = k_a + k_s + k_n$. Adding the null-collision term from Eq. (B.3) modifies Eq. (B.2) as follows:

$$ L(\vec{x}_0, \vec{u}_0) = \int_0^{+\infty} \tilde{k} \exp(-\tilde{k}\sigma) d\sigma \left[ \frac{k_a(\vec{x}_1)}{\tilde{k}} L^0(\vec{x}_1) \right. \\ + \frac{k_s(\vec{x}_1)}{\tilde{k}} \int_{4\pi} \Phi(\vec{u}_0, \vec{u}_1) L(\vec{x}_1, \vec{u}_1) d\vec{u}_1 \\ \left. + \frac{k_n(\vec{x}_1)}{\tilde{k}} \int_{4\pi} \Phi^n(\vec{u}_0, \vec{u}_1) L(\vec{x}_1, \vec{u}_1) d\vec{u}_1 \right] \quad (B.4) $$
Now, the phase function associated with null-collision events is : $\Phi^n(\vec{u}_0, \vec{u}_1) = \delta(\vec{u}_0 = \vec{u}_1)$ and therefore relation (B.4) becomes:

$$L(\vec{x}_0, \vec{u}_0) = \int_0^{+\infty} \hat{k} \exp(-\hat{k}\sigma) d\sigma \left[ \frac{k_a(\vec{x}_1)}{\hat{k}} L^0(\vec{x}_1) + \frac{k_s(\vec{x}_1)}{\hat{k}} \int_{4\pi} \Phi(\vec{u}_0, \vec{u}_1) L(\vec{x}_1, \vec{u}_1) d\vec{u}_1 + \left(1 - \frac{k(\vec{x}_1)}{\hat{k}}\right) L(\vec{x}_1, \vec{u}_0) \right] \quad (B.5)$$

This last relation is the basis for a null-collision algorithm in which the first step is to sample a length $\sigma$ over $[0, +\infty)$ according to the $\hat{k} \exp(-\hat{k}\sigma)$ probability density function. Since coefficient $\hat{k}$ is uniform in the medium, this sampling is straightforward. Once $\sigma$ has been sampled, position $\vec{x}_1 = \vec{x}_0 - \sigma\vec{u}_0$ is known. Then, the type of collision taking place at position $\vec{x}_1$ has to be sampled either as a absorption, a scattering or a null-collision according to the $[p_a, p_s, 1 - (p_a + p_s)]$ probability set with $p_a = \frac{k_a(\vec{x}_1)}{\hat{k}}$ the probability of absorption, $p_s = \frac{k_s(\vec{x}_1)}{\hat{k}}$ the probability for scattering, and $1 - (p_a + p_s)$ the probability of null-collision. When a absorption takes place, the weight of the realization is $L^0(\vec{x}_1)$ and the realization stops. When a scattering occurs, a propagation direction after scattering $\vec{u}_1$ has to be sampled over the unit sphere according to the phase function, and the algorithm loops for position $\vec{x}_1$ and direction $\vec{u}_1$; when a null-collision occurs, the algorithm loops for position $\vec{x}_1$, but direction $\vec{u}_0$ has not been modified.

The advantage of such a algorithmic variant is that the path length $\sigma$ can easily be sampled, using the value of $\hat{k}$, the upper-bound value of the extinction coefficient over the volume. But it also introduces a new type of collisions : null-collisions do not modify the optical path or its associated energy, but the more often they take place, the more often the absorption and scattering coefficients have to be evaluated for the current collision position. One could argue that using a infinite value of $\hat{k}$ solves the problem of knowing whether or not $\hat{k}$ is a upper-bound value of $k(\vec{x})$. But in that case, path length $\sigma$ would be null, and an infinite number of null-collision events would take place. More generally, the value of $\hat{k}$ must be kept as low as possible in order to reduce the probability of null-collision.

This can be a problem when the extinction coefficient is very heterogeneous, as in Fig. 1(a) : when the extinction coefficient is locally orders of magnitude higher than its value in the bulk of the volume, using a single value of $\hat{k}$ throughout the whole volume is cumbersome, because many null-collision events will need to be sampled before a absorption event occurs. This issue can be leveraged by using different values of $\hat{k}$ in various volume elements (voxels). The range of $\sigma$ is now limited to $\sigma \in [0, \sigma_{int}]$ with $\sigma_{int}$ the curvilinear abscissa of the intersection position $\vec{x}_{int} = \vec{x}_0 - \sigma_{int}\vec{u}_0$ between the line of sight and the boundary of the current voxel $\mathcal{V}$. Including this constraint in formal developments, Eq. (B.5) now becomes:

$$L(\vec{x}_0, \vec{u}_0) = \int_0^{+\infty} \hat{k} \exp(-\hat{k}\sigma) d\sigma \left\{ \begin{aligned} & \left[ \frac{k_a(\vec{x}_1)}{\hat{k}} L^0(\vec{x}_1) + \frac{k_s(\vec{x}_1)}{\hat{k}} \int_{4\pi} \Phi(\vec{u}_0, \vec{u}_1) L(\vec{x}_1, \vec{u}_1) d\vec{u}_1 \right. \\ & \left. + \left(1 - \frac{k(\vec{x}_1)}{\hat{k}}\right) L(\vec{x}_1, \vec{u}_0) \right] \mathcal{H}(\vec{x}_1 \in \mathcal{V}) \\ & + L(\vec{x}_{int}, \vec{u}_0) \mathcal{H}(\vec{x}_1 \notin \mathcal{V}) \end{aligned} \right\} \quad (B.6)$$

The associated null-collision algorithm for the evaluation of a monochromatic probe intensity at position $\vec{x}_0$ in direction $\vec{u}_0$ consists in :

* sampling a path length $\sigma$ according to the $\hat{k} \exp(-\hat{k}\sigma)$ probability density function, with $\hat{k}$ the upper-bound of the extinction coefficient for the current voxel $\mathcal{V}$ (the voxel that contains position $\vec{x}_0$;

* testing whether collision position $\vec{x}_1 = \vec{x}_0 - \sigma\vec{u}_0$ belongs to current voxel $\mathcal{V}$, or not. If $\vec{x}_1$ belongs to $\mathcal{V}$, the algorithm proceeds as previously (selecting a type of collision, and acting accordingly)

* in the case when $\vec{x}_1$ does not belong to $\mathcal{V}$, the Monte Carlo weight is now $L(\vec{x}_{int}, \vec{u}_0)$; since it is generally unknown, the algorithm loops at intersection position $\vec{x}_{int}$ and in direction $\vec{u}_0$. In other words, the optical path continues, but enters a different voxel.

This makes possible to use any upper-bound values spatial structure. As discussed in the main text, we chose to build a structure in which finding whether or not a position belongs to a given voxel is straightforward, and that minimizes the number of null-collision events.

## Appendix C. External source

In the case when the radiation source is external to the medium, emission by the medium is (by definition) not considered. The total intensity has to be decomposed into its direct part $L_{diff}(\vec{x}, \vec{u}) = L^*\mathcal{T}(\vec{x})$ (with $L^*$ the intensity of the source and $\mathcal{T}(\vec{x})$ the transmissivity between the source and position $\vec{x}$) and its diffuse part $L_{diff}(\vec{x}, \vec{u})$ (radiation that has been scattered at least once before reaching position $\vec{x}$ in direction $\vec{u}$). Also, assuming that the radiation source has a limited angular extension $\Omega^*(\vec{x}_1)$, the integral formulation of the probe intensity is now:

$$ \left\{ \begin{aligned} L(\vec{x}_0, \vec{u}_0) &= L^*\mathcal{T}(\vec{x}_0)\mathcal{H}(\ast) + L_{diff}(\vec{x}_0, \vec{u}_0) \\ L_{diff}(\vec{x}_0, \vec{u}_0) &= \int_0^{+\infty} \hat{k} \exp(-\hat{k}\sigma) d\sigma \left\{ \begin{aligned} & \left[ \frac{k_a(\vec{x}_1)}{\hat{k}} \times 0 \right. \\ & \left. + \frac{k_s(\vec{x}_1)}{\hat{k}} \left( \int_{\Omega^*(\vec{x}_1)} \frac{1}{\Omega^*(\vec{x}_1)} d\vec{u}_1 \{\Phi(\vec{u}_0, \vec{u}_1)\Omega^*(\vec{x}_1)L^*\mathcal{T}(\vec{x}_1)\mathcal{H}(\ast)\} \right) \right. \\ & \left. + \int_{4\pi} \Phi(\vec{u}_0, \vec{u}_1) d\vec{u}_1 \{L_{diff}(\vec{x}_1, \vec{u}_1)\} \right) \right. \\ & \left. + \left(1 - \frac{k(\vec{x}_1)}{\hat{k}}\right) L_{diff}(\vec{x}_1, \vec{u}_0) \right] \mathcal{H}(\vec{x}_1 \in \mathcal{V}) \\ & + L_{diff}(\vec{x}_{int}, \vec{u}_0) \mathcal{H}(\vec{x}_1 \notin \mathcal{V}) \end{aligned} \right\} \quad (C.1) $$

The corresponding Monte Carlo algorithm is the following:

* the realization weight is now initialized with the value of the direct intensity, in the case the source is visible along the $(\vec{x}_0, -\vec{u}_0)$ line of sight (which is accounted for by the $\mathcal{H}(\ast)$ symbol). Transmissivity $\mathcal{T}(\vec{x}_0)$ between the source and position $\vec{x}_0$ has to be evaluated (this will be detailed below).

* then the diffuse part of the intensity has to be evaluated ; for this, a path length is sampled over $[0, +\infty)$ according to the $\hat{k} \exp(-\hat{k}\sigma)$ probability density function, with $\hat{k}$ the upper-bound of the extinction coefficient for the current voxel;

* when position $\vec{x}_1 = \vec{x}_0 - \sigma\vec{u}_0$ belongs to current voxel, the type of collision has to be sampled. When a absorption is selected, the realization just stops, with no update to the realization weight. When a null-collision is sampled, the evaluation of the diffuse intensity loops for position $\vec{x}_1$, in direction $\vec{u}_0$; and when scattering is selected, two things occur: firstly, a direction $\vec{u}_1$ is uniformly sampled in the angular angle that defines the source, as seen from $\vec{x}_1$, and the realization weight is increased by the value of $\Phi(\vec{u}_0, \vec{u}_1)\Omega^*(\vec{x}_1)L^*\mathcal{T}(\vec{x}_1)$ when the line of sight between $\vec{x}_1$ and the source is not obstructed. This also requires the evaluation of the transmissivity $\mathcal{T}(\vec{x}_1)$ between the source and local position $\vec{x}_1$. Secondly, a new direction $\vec{u}_1$ has to be sampled over the unit sphere according to the phase function, and the algorithm for the evaluation of the diffuse intensity loops for position $\vec{x}_1$, in direction $\vec{u}_1$.
• when position $\vec{x}_1$ does not belong to the current voxel, the algorithm for the evaluation of the diffuse intensity loops at position $\vec{x}_{int}$, in direction $\vec{u}_0$ (there is a need to change voxel).

Transmissivity $\mathcal{T}(\vec{x}_0)$ between the source and position $\vec{x}_0$ is itself expressed as:

$$ \mathcal{T}(\vec{x}_0) = \int_0^{+\infty} \hat{k} \exp(-\hat{k}\sigma) d\sigma \left\{ \left( \frac{k(\vec{x}_1)}{\hat{k}} \times 0 + \left(1 - \frac{k(\vec{x}_1)}{\hat{k}}\right) \mathcal{T}(\vec{x}_1) \right) \mathcal{H} \times (\vec{x}_1 \in \mathcal{V}) + \mathcal{T}(\vec{x}_{int}) \mathcal{H}(\vec{x}_{int} \notin \mathcal{V}) \right\} \quad (C.2) $$

Starting from position $\vec{x}_0$, a realization of this transmissivity consists in sampling a path length $\sigma$ over $[0, +\infty)$ according to the $\hat{k} \exp(-\hat{k}\sigma)$ probability density function, with $\hat{k}$ the upper-bound of the extinction coefficient for the current voxel; if position $\vec{x}_1 = \vec{x}_0 - \sigma \vec{u}_0$ belongs to current voxel, the collision is tested for an absorption or null-collision, in which case the algorithm loops from position $\vec{x}_1$. If position $\vec{x}_1$ does not belong to the current voxel, the algorithm loops from the intersection position between the line of sight and the boundary of the voxel. The recursive algorithm only stops in two cases : when an absorption event is sampled, in which case the weight for transmissivity takes a null value ; and when the source is reached, in which case the weight for transmissivity takes a value of 1.

This series of nested algorithms are useful when the participating medium is lit by the external source. But when only a small fraction of the volume is being directly illuminated by the source, a new issue arises : the contribution of each successive scattering position gathers most of the time a null contribution to the Monte Carlo weight. When the medium is illuminated by a laser sheet, as is the case in the present article, the variance of such a algorithm would be catastrophic. This is the reason why a new algorithm, for the specific case of partially illuminated media, was presented in [12].

## Data availability

Data will be made available on request.

## References

[1] Hammersley JM, Handscomb DC. Monte-Carlo methods. John Wiley, New York; 1964.

[2] Delatorre J, Bézian JJ, Blanco S, Caliot C, Cornet JF, Dauchet J, ElHafi M, Eymet V, Fournier R, Gautrais J, Gourmel O, Veynandt F, Meilhac N, Pajot A, Paulin M, Perez P, Piaud B, Roger M, Rolland J, Weitz S. Monte-Carlo advances and concentrated solar applications. Sol Energy 2014;(103):653-81. http://dx.doi.org/10.1016/j.solener.2013.02.035.

[3] Villefranque N, Fournier R, Couvreux F, Blanco S, Cornet C, Eymet V, Forest V, Tregan J-M. A path-tracing Monte Carlo library for 3-D radiative transfer in highly resolved cloudy atmospheres. J Adv Model Earth Syst 2019;11(8):2449-73. http://dx.doi.org/10.1029/2018MS001602.

[4] Villefranque Najda, Hourdin Frédéric, d'Alençon Louis, Blanco Stéphane, Boucher Olivier, Caliot Cyril, Coustet Christophe, Dauchet Jérémi, Hafi Mouna El, Eymet Vincent, Farges Olivier, Forest Vincent, Fournier Richard, Gau-trais Jacques, Masson Valéry, Piaud Benjamin, Schoetter Robert. The “teapot in a city”: A paradigm shift in urban climate modeling. Sci Adv 2022;8(27). http://dx.doi.org/10.1126/sciadv.abp8934, arXiv:https://www.science.org/doi/pdf/10.1126/sciadv.abp8934. URL https://www.science.org/doi/abs/10.1126/sciadv.abp8934.

[5] Kutz Peter, Habel Ralf, Li Yining Karl, Novak Jan. Spectral and decomposition tracking for rendering heterogeneous volumes. ACM Trans Graph 2017;36(4). http://dx.doi.org/10.1145/3072959.3073665.

[6] He Z, Vinatier S, Eymet V, Forest V, Bézard B, Rannou P, Rodriguez S, Marcq E, Fournier R, Blanco S, Mourtaday N, Nyffenegger-Péré Y, Lebonnois S, Määt-tänen A. Simultaneous estimation of radiance and its sensitivities to radiative properties in a spherical-heterogeneous atmospheric radiative transfer model by Monte Carlo: Application to Titan. J Quant Spectrosc Radiat Transfer 2025. http://dx.doi.org/10.2139/ssrn.5105650, Under review.

[7] Burley B, Adler D, Chiang MJ-Y, Driskill H, Habel R, Kelly P, et al. The design and evolution of Disney's hyperion renderer. ACM Trans Graph 2018;37(3).

[8] Galtier M, et al. Integral formulation of null-collision Monte-Carlo algorithms. J Quant Spectrosc Radiat Transfer 2013;125:57-68.

[9] E. Woodcock, Murphy T, Hemmings P, Longworth S. Techniques used in the gem code for Monte-Carlo neutronics calculations in reactors and other systems of complex geometry. In: Proc. conf. applications of computing methods to reactor problems. 1965, p. 557.

[10] Novak I, Hanika J, Jarosz W. Monte Carlo methods for volumetric light transport simulation. Comput Graph Forum 2018;37:551-76.

[11] Novák Jan, Selle Andrew, Jarosz Wojciech. Residual ratio tracking for estimating attenuation in participating media. ACM Trans Graph (Proc SIGGRAPH Asia) 2014;33(6). http://dx.doi.org/10.1145/2661229.2661292.

[12] Sans M, El Hafi M, Eymet V, Forest V, Fournier R, Villefranque N. Null-collision meshless Monte-Carlo - a new reverse Monte-Carlo algorithm designed for laser-source emission in absorbing/scattering inhomogeneous media. J Quant Spectrosc Radiat Transfer 2021;271. http://dx.doi.org/10.1016/j.jqsrt.2021.107725.

[13] AVBP code. https://Www.Cerfacs.Fr/Avbp7x/.

[14] Yon J, Liu F, Bescond A, Caumont-Prim C, Rozé C, Ouf F-X, Copalle A. Effects of multiple scattering on radiative properties of soot fractal aggregates. J Quant Spectrosc Radiat Transfer 2014;133:374-81. http://dx.doi.org/10.1016/j.jqsrt.2013.08.022.

[15] Eymet V, Poitou D, Galtier M, El Hafi M, Terrée G, Fournier R. Null-collision meshless Monte-Carlo. Application to the validation of fast radiative transfer solvers embedded in combustion simulators. J Quant Spectrosc Radiat Transfer 2013;129:145-57. http://dx.doi.org/10.1016/j.jqsrt.2013.06.004.

[16] Rubin Steven M, Whitted Turner. A 3-dimensional representation for fast rendering of complex scenes. SIGGRAPH Comput Graph 1980;14(3):110-6.

[17] Wald Ingo, Woop Sven, Benthin Carsten, Johnson Gregory S, Ernst Manfred. Embree: a kernel framework for efficient CPU ray tracing. ACM Trans Graph 2014;33(4):1-8.

[18] Joo PH, Gigone B, Griffin EA, Christensen M, Gülder ÖL. Soot primary particle size dependence on combustion pressure in laminar ethylene diffusion flames. Fuel 2018;220:464-70.