# A full Monte-Carlo resolution scheme of a photon-phonon model for conducto-radiative heat transfer in heterogeneous media

Gaël Poëtte, Augustin de la Vauvre, Gérard Vignoles

► To cite this version:

Gaël Poëtte, Augustin de la Vauvre, Gérard Vignoles. A full Monte-Carlo resolution scheme of a photon-phonon model for conducto-radiative heat transfer in heterogeneous media. International Journal of Heat and Mass Transfer, 2025, 239, pp.126603. 10.1016/j.ijheatmasstransfer.2024.126603. hal-04916277

HAL Id: hal-04916277

https://ifp.hal.science/hal-04916277v1

Submitted on 28 Jan 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>

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

# A full Monte-Carlo resolution scheme of a photon-phonon model for conducto-radiative heat transfer in heterogeneous media

Gaël Poëtte<sup>a,b,*</sup>, Augustin De La Vauvre<sup>c,d</sup>, Gérard Vignoles<sup>e</sup>

## ARTICLE INFO

**Keywords:**

Radiative transfer

Diffusion

Monte Carlo

Numerical scheme

Layered media

Stratified media

## ABSTRACT

In this paper, we tackle the problem of solving the system of radiative heat transfer coupled with the diffusion equation for material temperature with a full Monte-Carlo (MC) scheme in heterogeneous media, i.e. possibly containing discontinuities of thermophysical material parameters. MC schemes based on Brownian paths for the diffusion equation in layered media are known to present some issues: corrections at the interface between media are usually mandatory or must be replaced with complex samplings. In this paper, based on the observations that transport models are easily solved by MC methods in layered media, and that state-of-the-art MC schemes capture the diffusion limit, we suggest to build a phonon like model and to solve it with an MC scheme. For capturing the diffusion limit, physical constraints are explicit and enforced within the MC resolution. These constraints can even help revisiting the Brownian paths for efficient corrections. Numerical results are provided in different regimes in order to highlight the relevance of the strategies suggested in this paper.

## 1. Introduction

Heat transfer in heterogeneous and porous media mixing conduction and radiation has attracted much attention, due to many applications in materials development. For instance, solar radiation receivers in concentrated solar power plants are being engineered considering radiative and conductive heat transfer [1,2]. In space technology, the design of thermal protection systems materials for atmospheric re-entry also involves the evaluation of mixed radiative/conductive heat transfer [3–5]. In industry, porous media made of carbon or ceramic foams and/or fibers are used as high-temperature thermal insulation materials [6] or as burners [7,8]. Currently, cellular ceramics attract a large attention as heat exchange enhancers [9]. In all these applications, a crucial point for the design of efficient parts relies on modeling radiative/conductive coupled heat transfer.

Varieties of solvers have been proposed to address this problem [10–14]; they can be decomposed in two main classes, the deterministic ones (based on finite element, finite differences, finite volumes [15–21]) and the stochastic ones, on which we focus in this paper:

- deterministic solvers are generally fast, generate smooth solutions but are sensitive to the dimension of the problem in the sense

that going from spatial dimension 2 to 3 induces a non negligible overcost. Complex geometries may also be quite hard to handle, especially if thin layers of media are intertwined.

- On another hand, stochastic solvers are based on the *central limit theorem* [22]: it typically means that they are less sensitive to the regularity of the solution or to the curse of dimensionality. The latter point makes them particularly interesting when uncertainties have to be taken into account [23–25]. But the convergence rate is slower ($\propto \frac{1}{\sqrt{N_{MC}}}$ where $N_{MC}$ is the number of particles) and the methods produce noisy solutions if the number of particle is not high enough. These solvers are generally costly but they are also less constraining in terms of closure hypothesis. Furthermore, a *posteriori* error estimations are available.

Here, we are mostly interested by Monte-Carlo/Random Walks schemes. This technique has long been one of the favorites for the solution of the radiative part of heat transfer, especially in case of complex geometries [26,27]; here, the random walk algorithm is directly based on the physics of radiative heat transfer, simulating scattering, absorption and emission with Poissonian processes. On the other hand,

<footnote><sup>a</sup> CEA DAM CESTA, F-33114 Le Barp, France</footnote>

<footnote><sup>b</sup> Institut de Mathématiques de Bordeaux, Université de Bordeaux, CNRS, Bordeaux INP, France</footnote>

<footnote><sup>c</sup> IFPEN, France</footnote>

<footnote><sup>d</sup> LTeN, France</footnote>

<footnote><sup>e</sup> Université Bordeaux, LCTS, France</footnote>

<footnote>* Corresponding author at: CEA DAM CESTA, F-33114 Le Barp, France.</footnote>

<footnote>E-mail addresses: gael-poette@cea.fr (G. Poëtte), augustin.de-la-vauvre@univ-nantes.fr (A. De La Vauvre), vinhola@lcts.u-bordeaux.fr (G. Vignoles).</footnote>
in the diffusive limit – that is, for very small free path lengths –, it has been found handy to design algorithms based on the laws of Brownian motion, sometimes called “walk-on-sphere” methods [28–30].

When heat conduction and radiation are mixed in opaque/ semitransparent heterogeneous media, it has been proposed to design “hybrid walks” whereby the walk rules are differing, depending on whether the walker lies in the opaque or the transparent medium, with specific rules for the transitions between both types of media. Successful implementations with Poissonian behavior in the semi- transparent subdomains and Brownian behavior in the opaque ones were reported by Vignoles [31] and Tregan et al. [32]. Another type of coupling has been proposed in the case of semi-transparent media in which both types of transport occur in the same regions of space [33– 35]. Despite these successes, as will be shown in this paper, there exist some uncertainties linked to the solution of Brownian motion problems close to an opaque/semitransparent interface.

Here our goal is to design a fully “kinetic” (Poissonian) scheme in which the essence of the random walk algorithm is the same every- where, but in which discontinuities of material thermophysical prop- erties are handled properly, i.e. without introducing any bias in the results and recovering the correct diffusive limits when appropriate. We will see that analyzing Poissonian paths will also help us have a better understanding of Brownian ones.

We are interested in the Monte-Carlo (MC) resolution of the time- dependent, nonlinear, radiative transfer equations coupled to heat diffusion in layered media. By layered media, we mean configurations in which several layers of media with very different characteristics (conductivity, density, heat capacity etc.) are superimposed. The model has general form (see [36,37]):

$$
\left\{
\begin{array}{l}
\displaystyle \frac{1}{c} \partial_t I + \omega \cdot \nabla I + \sigma_t I = \sigma_a B(T_m) + \sigma_s \int_{4\pi} I \frac{d\omega'}{4\pi}, \\[2ex]
\displaystyle \partial_t E(T_m) - \nabla \cdot (\kappa(T_m) \nabla T_m) = \int_{4\pi} c \sigma_a \left( \frac{I}{4\pi} - B(T_m) \right) d\omega'.
\end{array}
\right.
\qquad (1)
$$

In the above equations, $I = I(t, x, \omega, v)$ and $T_m(t, x)$ are the unknowns of the system and stand respectively for the density of radiation energy (J m⁻³) and the material temperature (K). Variables $t \ge 0$, $x \in \Omega \subset \mathbb{R}^3$, $\omega \in \mathbb{S}^2$ and $v \in \mathbb{R}^+$ are respectively the time (s), space (m), solid angle (sr) and frequency (s⁻¹) variables. The opacities (m⁻¹) $\sigma_t = \sigma_t(x, t, v)$, $\sigma_a = \sigma_a(x, t, v)$ and $\sigma_s = \sigma_s(x, t, v)$ are given functions standing for the total, absorption and scattering opacities. In particular, we have $\sigma_t = \sigma_a + \sigma_s$. The density of internal energy $E$ (J m⁻³) depends on $T_m$ via an equation of state (EoS) $dE = \rho c_v (T_m) dT_m$ with $\rho$ the mass density (kg m⁻³) and $c_v$ the heat capacity (J K⁻¹). It is constant for a perfect gas. Quantity $B(x, t, v)$ is the Planckian distribution (J m⁻³ s⁻¹) and is such that $B(x, t) = B(T_m(x, t)) = \int B(x, t, v) dv = a T_m^4(x, t) / 4\pi$ is the frequency-integrated Planck function (J m⁻³) with $a$ the radiative constant (J m⁻³ K⁻⁴). Here, $c$ denotes the speed of light (m s⁻¹). The conductivity (J m⁻¹ K⁻¹) of the medium is denoted by $\kappa$. In this paper, emphasis is set onto the case of spatially discontinuous $\rho c_v$ and/or $\kappa$ in order to model layered media and the numerical difficulties it triggers. In other words, we focus on some mathematical and numerical problems and for this reason, almost all the test-cases of the paper are dimensionless. Let $(\Omega_i)_{i \in \{1, \dots, M\}}$ denote $M$ different media (or “phases”) involved in our configuration of interest, i.e. such that $\cup_{i=1}^M \Omega_i = \Omega$, then for $\alpha \in \{\rho c_v, \kappa\}$, we have $\alpha(\mathbf{x}) = \sum_{i=1}^M \alpha_i 1_{\Omega_i}(\mathbf{x})$ where $1_{\Omega_i}(\mathbf{x})$ equals 1 if $\mathbf{x} \in \Omega_i$, 0 otherwise. Note that the discontinuities may also come from the temperature dependence of $\kappa$: suppose there is only one material but the discretisation is such that $T_m$ is piecewise constant in each of the $N$ cells, i.e. $T_m(x) = \sum_{i=1}^N T_m^i 1_{\Omega_i}(x)$: this naturally leads to a piecewise constant conductivity given by $\kappa(x) = \sum_{i=1}^N \kappa(T_m^i) 1_{\Omega_i}(x) = \sum_{i=1}^N \kappa^i 1_{\Omega_i}(x)$. The two different cases shall be treated the same way in this paper, without lack of generality. Initial and boundary conditions must be supplemented to system (1):

$$
I(0, x, \omega) = I^0(x, \omega), \quad T_m(0, x) = T_m^0(x), \quad x \in \Omega, \quad \omega \in \mathbb{S}^2, \qquad (2)
$$

$$
I(t, x, \omega) = I_b(t, \omega), \quad t \ge 0, \quad x \in \partial\Omega, \quad \omega \cdot n_s < 0, \qquad (3)
$$

where $n_s$ is the outward normal to $\Omega$ at $x$. System (1) together with initial and boundary conditions (2) + (3) define the well-posed [38] mathematical problem we want to solve.

System (1) is nonlinear. Its resolution consequently needs some linearisation hypothesis. System (1) is often solved by performing some operator splitting on a time step [39] before choosing numerical schemes for the different splitting steps. The choice of the splitting is not unique but a classical and convenient one is the following:

- A radiative transfer phase

$$
\left\{
\begin{array}{l}
\displaystyle \frac{1}{c} \partial_t I + \omega \cdot \nabla I + \sigma_t I = \sigma_a B(T_m) + \sigma_s \int_{4\pi} I \frac{d\omega'}{4\pi}, \\[2ex]
\displaystyle \partial_t E(T_m) = \int_{4\pi} c \sigma_a \left( \frac{I}{4\pi} - B(T_m) \right) d\omega',
\end{array}
\right.
\qquad (4)
$$

- followed by a heat conduction one

$$
\left\{
\begin{array}{l}
\partial_t I = 0, \\
\partial_t E(T_m) - \nabla \cdot (\kappa(T_m) \nabla T_m) = 0,
\end{array}
\right.
\qquad (5)
$$

during which $I$ remains constant.

This splitting is especially convenient when one only needs to focus on radiative transfer for example (i.e. when $\Omega = \cup_{i=1}^M \Omega_i$ has many non conducting phases for which $\kappa_i \ll 1$):

Once a splitting chosen, it remains to choose the resolution scheme for both steps

- due to the high dimensionality of the problem (with $(x, t, \omega, v) \in \mathbb{R}^7$ in (1) and $(x, t, \omega, v, v) \in \mathbb{R}^8$ later on for the photon-phonon model (12)) we aim at tackling, we focus on MC based numerical solvers for the radiative transfer phase. MC schemes are known to be costly but insensitive to the dimension of the problem and to the regularity of the solution [22]. In order to solve (4), we rely on the ISMC (Implicit Semi-analog Monte-Carlo) scheme described in [40] based on both the seminal works of [41,42]: the MC scheme is implicit and teleportation error free. For the spectral dependance, we rely on the material of [43] (even if in this paper, without lack of generality with respect to the numerical difficulties we aim at tackling, the grey approximation is enough). The ISMC scheme discretises both $I$ and $E$ with MC particles which are transported (at velocity $c$ for photons $I$) or not (for matter at rest $E$). The nature of the MC particles eventually changes at a collision at which a photon can become sensible heat or a phonon and sensible heat becomes photons, see [40,43]. At the end of the time step, $I$ and $E$ are approximated in each cell by tallying the contributions of, respectively, the photon MC particles and the sensible heat MC particles. The temperature $T_m$ is deduced from $E$ by inversion of the EoS.

- The second step (5) of the splitting, in lower dimension (with $T_m$ only depending on $(x, t) \in \mathbb{R}^3 \times \mathbb{R}^+$), is usually solved thanks to deterministic solvers. From the temperature computed during the previous phase and the application of a diffusion scheme, the energy and the temperature are updated.

This solution strategy works in practice but one must realise that in order to end the time iteration, the weights of all the sensible heat MC particles must be updated in accordance with the new values of $E$ determined by the diffusion step, before starting the radiative transfer phase of the next time step. This step is sequential (i.e. hard to parallelise [44]) and can lead to unbalanced weights. On the other hand, Eq. (5) can also be solved with an MC scheme: in practice, this is based on the application of Itô's lemma [37,45–53]. In a homogeneous media, the application is quite simple as it only consists in solving

$$
\left\{
\begin{array}{l}
\displaystyle \partial_t T - D \Delta T = 0, \\
\displaystyle T(t=0, x) = T_0(x),
\end{array}
\right.
\qquad (6)
$$
on the time step where we defined $D = \frac{\kappa}{\rho c_v}$ (for a homogeneous isotropic material). Itô's lemma then states that in order to solve (6), the position of the MC particles $X_t$ must be updated by solving the stochastic ordinary equation

$$ \begin{cases} dX_t = \sqrt{2D} dW_t, \\ X_0 = x, \end{cases} \qquad (7) $$

where $W_t$ denotes a Brownian motion (or a Wiener process). Then $T(x, t) = \mathbb{E}[T_0(X_t)]$ is solution of (6). In practice, (7) is often solved thanks to an explicit Euler scheme.

Let us now consider a simple test-case for (6) and illustrate the difficulty we aim at tackling in this paper: it is a one dimensional test-case $x \in \Omega = [0, 1] \subset \mathbb{R}$, with initial condition

$$ T_0(x) = 2 \times \mathbf{1}_{[0.45,0.55]}(x), \qquad (8) $$

where $\mathbf{1}_{[a,b]}$ is the Heaviside function for interval $[a, b]$. The initial condition in terms of temperature is displayed in the left plots of Fig. 1. Quantities $x \to \kappa(x)$ and $x \to \rho(x)c_v(x)$ are also displayed (left vertical axis) together with the reference solution of the test-case. The boundary conditions are Neumann ones but the final time $t = 10^{-2}$ is such that the solution does not have enough time to significantly interact with the edges of the domain. Now, the domain is composed of two layers of material in $\Omega_1 = [0, \frac{1}{2}]$ and $\Omega_2 = [\frac{1}{2}, 1]$ such that $\Omega_1 \cup \Omega_2 = \Omega$ and with diffusivities $D_1$ and $D_2$. If $D_1 = D_2$, then the previously described MC resolution strategy is efficient: on picture top right of Fig. 1, the reference solution obtained with a deterministic code and the MC ones are in very good agreement (the number of cells is $N_x = 100$), especially as $\Delta t \le 10^{-6}$. On the other hand if $D_1 \ne D_2$ ($D_1 = 10$ and $D_2 = 1.0$ in the bottom picture of Fig. 1), the MC results (Brownian paths) do not converge toward the same temperature profile as the one obtained by the deterministic solver (reference). It is well-known that such an MC resolution scheme is not efficient for layered media and that the physical solution for this problem is the one obtained with the deterministic solver [45,46]. In this case, increasing the number of particles, the number of cells or decreasing the time step does not ensure obtaining more physical results [45,46,49]. In this paper, we want to suggest a new way to deal with this phenomenon, i.e. jumps in $x \to \rho(x)c_v(x)$ and/or $x \to \kappa(x)$.

## 2. An overview of the strategies of the literature

As explained in the introduction, the flaw illustrated in Fig. 1 for the Brownian path in layered media is well-known, see [54] for example. Many authors have suggested ways to circumvent it. The suggested solutions can be decomposed in two different families of methods:

- the ones working on adapting the time steps,

- the ones having a fixed time step but with interface corrections.

In order to improve the results for some non homogeneous diffusivity coefficients $x \to D(x)$, some authors relied on linearisations of $x \to D(x) = D(x) + \Delta x \nabla_x D(x) + \mathcal{O}(\Delta x^2)$. In [48], the authors empirically demonstrate the improvements made by correcting the time step of the resolution when $x \to D(x)$ is not constant and $\nabla_x D(x)$ is available. The numerical scheme does improve the results on the benchmarks of the paper but $x \to D(x)$ is quite smooth ($C^2$), which is clearly not the case in layered media.

In [46], the idea is to accept the lack of regularity of $x \to D(x)$ and consider discontinuous diffusivities. The stochastic process that governs the positions of the particles is fully characterised. Two algorithms are suggested. They rely on fixed time steps and the construction of an

interface layer, a numerical vicinity of the interface, within which skew Brownian motions (SBM) are simulated. Several schemes are proposed in the literature in order to treat the samplings within this layer, see [46,55,56]. The SBM one [55–57] seems to outperform the others. But the scheme is quite complex to implement and its efficiency on multi-layered media is not obvious nor, to our knowledge, has been numerically demonstrated [53]. The question of what happens to the numerical interface layer when the mesh is refined remains: remember we aim at solving (1) and the additional equation may require finer meshes for accuracy.

In [51], the samplings needed at the interface between two media are studied: it is explained that introducing reflections at the boundary aims at maintaining some no-flux boundary conditions for the random walk. Some treatments fail to solve the specified problem while it is explained that the methods of [58,59] succeed. One can show that these methods relate to an analytical solution to the problem of one-dimensional diffusion at an interface. But the use of the analytical solution for the samplings at the interface [52,59,60] is quite complex, does not guarantee convergence properties as the number of layers grows and often needs the introduction of some artificial interface layers [58,59]: in other words, it introduces some additional (in plus of $\Delta t$) numerical parameter which must be tuned (and for which, to our knowledge, convergence is not guaranteed).

In [37,47,61], the authors suggest introducing some kind of albedo at the interface between the two domains: the particles have a certain probability of bouncing back into the domain they come from instead of crossing it.

Note that some authors, in combination to the use of rebound probabilities, introduce the need for having negative weights for their particles [61]. This strategy can be efficient in practice [61] or even [62,63] for some applications of negative weights for neutronics applications but they need to be combined with fine algorithm for population control (or the number of negative weight particles and the variance may increase [61,62]).

To sum up, many numerical methods aiming at dealing with discontinuous material properties for the diffusion equation with random walkers are suggested in the literature. But they either

- rely on smoothness assumptions which do not hold for our applications,

- or seem to be efficient but complex to implement,

- or rely on numerical parameters which may need to be tuned without convergence guarantees.

For all these reasons, we suggest a new numerical strategy to tackle the same kind of problems. It is based on what we consider simple modifications of the model (1) we aim at solving.

## 3. The photon-phonon coupled system for conducto-radiative transfer

From the previous simple example (see Fig. 1 and the comments), we can see that the resolution of the diffusion equation with Brownian paths is very efficient in almost homogeneous media but remains complex to treat layered ones. The main idea of this paper comes from two observations:

- the MC resolution of the transport equation (Poissonian paths rather than Brownian ones) does not suffer from the previous problem: layered media can easily and reliably be simulated.

- The diffusion equation is the limit of the transport one for both fast particles and opaque media and the SMC (Semi-analog Monte-Carlo) and ISMC (Implicit Semi-analog Monte-Carlo) schemes capture this limit [40]. By the way, the material of this paper can easily be applied with other (teleportation error free) MC solvers such as the ones in [64–66].

<footnote n="1">¹ $x \to \kappa(x) = 1$ and $x \to \rho(x)c_v(x) = 1$ are displayed on the right vertical axis.</footnote>

<footnote n="2">² $x \to \kappa(x) = 10^{-1}\mathbf{1}_{[0.5,1]}(x) + \mathbf{1}_{[1.5,1]}(x)$ and $x \to \rho(x)c_v(x) = 1$ are displayed on the right vertical axis.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Left: initial conditions. Right: numerical solutions. For the top pictures, the diffusivities are the same ($D_1 = D_2$) in domains $\Omega_1$ and $\Omega_2$ (the background colors are the same on each side of the domain). The plot presents the simulated temperatures at time $t = 10^{-2}$ obtained with a reference deterministic code and an MC one (Brownian path) with $10^5$ particles and time steps $\Delta t_{ref} \in \{10^{-5}, 10^{-6}, 10^{-7}\}$. Bottom: the calculations are carried out in the same conditions except the diffusivities are different ($D_1 = 10.0 \neq D_2 = 1.0$) in domains $\Omega_1$ and $\Omega_2$. This is highlighted by the different background colors. The plot presents the simulated temperatures at time $t = 10^{-2}$ obtained with a reference deterministic code and an MC one (Brownian path) with $10^5$ particles and time steps $\Delta t_{ref} \in \{10^{-5}, 10^{-6}, 10^{-7}\}$.</figcaption></figure>

Hence the idea to solve a transport model, the one for phonons described in [67], in the equilibrium diffusion limit in order to simulate heat conduction. For this, we suggest solving

$$ egin{cases} 
\partial_t I + c\omega \cdot \nabla I + c\sigma_t I = c\sigma_a B(T_m) + c\sigma_s \int_{4\pi} \frac{I}{4\pi} d\omega', \\
\partial_t e + c_m\omega \cdot \nabla e + c_m\sigma_m e = c_m\sigma_m \int_{4\pi} P_m e \frac{d\omega'}{4\pi} d\sigma_m' \\
+ \int_{4\pi} c\sigma_a \left(\frac{I}{4\pi} - B(T_m)\right) d\omega', \\
\text{with } E(T_m(x,t)) = \iint e(x, t, \omega, c_m) d\omega d\sigma_m. 
\end{cases} 
\quad (9) $$

In the above system, $c_m$ may be understood as a kind of speed of sound in sensible heat and $\sigma_m(x, t) = \sigma_m(x, T_m(x, t))$ is a kind of opacity for phonons of density $e(x, t, \omega, c_m)$. But the truth is their relevance will only be numerical (and not physical) in the following sections. We could better call them “Poissonian sensible heat energy carriers” than “phonons”. The quantity $P_m$ describes how the phonons' velocities and angles are scattered at a collision, i.e. it is a scattering kernel. For the moment, the only property we suppose it satisfies is

$$ \int P_m(x \in \Omega_l, c'_m, v) d\sigma'_m = \delta_{c_{m,l}}(v), \forall x \in \Omega, \forall v \in \mathbb{R}^3. \quad (10) $$

In the above expression, $\delta_c(v)$ is the Dirac delta function such that $\int f(v)\delta_c(v)dv = f(c)$ for all continuous compactly supported function $f$. Of course, in order to solve a heat conduction equation for sensible heat, $c_m$, $\sigma_m$ and $P_m$ must be astutely spatially chosen. Their tuning will be the purpose of Section 4. In this section, we focus on model (9), on its properties and on the resolution strategy we suggest in this paper.

First, system (9) is still conservative: indeed, integrating both Eqs. (9) with respect to $\omega$ and summing them up leads to

$$ \partial_t \left( I + E(T_m) \right) + \int c\omega \cdot \nabla I + \int c_m \omega \cdot \nabla e = 0, \quad (11) $$

so that in a closed system (i.e. such that the fluxes cancel), the sum of the photons and phonons' energies are constant $\int I + E(T_m) = \int I^0 + E(T_m^0)$.

Now, system (9) is nonlinear and must be linearised in order to be solved with an MC scheme. Indeed, an MC method looks for $N_{MC}$ particle $p \in \{1, \dots, N_{MC}\}$ solutions of the form $u_p(x, t, v, \omega) = w_p(t)\delta_x(x_p(t))\delta_v(v_p(t))\delta_\omega(\omega_p(t))$ and uses the linearity of the equation so that $\sum_{p=1}^{N_{MC}} u_p(x, t, v, \omega) \approx u(x, t, v, \omega)$ is also an approximation of $u$, the solution of the linear equation we need to solve. For this linearisation choice, we rely either on the SMC discretisation [42] or the ISMC one [40] (in the stiff case or when stationary computations are needed) which we need to modify in order to take the streaming ($c_m$) and the collision ($\sigma_m$) terms of the phonons into account. Of course, other linearisation choices could have been used [41,68]. Let us rewrite (9) as

$$ egin{cases} 
\partial_t I + c\omega \cdot \nabla I = -c\sigma_t I + c\sigma_a \eta(E) \int_{4\pi} \frac{e}{4\pi} d\omega' + c\sigma_s \int_{4\pi} \frac{I}{4\pi} d\omega', \\
\partial_t e + c_m\omega \cdot \nabla e = -c_m\sigma_m e - c\sigma_a \eta(E)e + \int_{4\pi} c\sigma_a \frac{I}{4\pi} d\omega' \\
+ c_m\sigma_m \int_{4\pi} \int_{4\pi} P_m \frac{e}{4\pi} d\omega' d\sigma'_m, 
\end{cases} 
\quad (12) $$

in which we introduced $B(E) = \eta(E)E$. This system is still nonlinear (due to the dependence of $\eta$ with respect to $E$). From now on, let us
assume that $\eta(E(x,t)) = \bar{\eta}(x,t)$ for $t \in [t^n, t^n + \Delta t]$. With the previous hypothesis, (12) becomes

$$ egin{cases} 
\partial_t I + c\omega \cdot \nabla I = -c\sigma_t I + c\sigma_a \bar{\eta} \int_{4\pi} \frac{e}{4\pi} d\omega' + c\sigma_s \int_{4\pi} \frac{I}{4\pi} d\omega', \\ 
\partial_t e + c_m \omega \cdot \nabla e = -c_m \sigma_m e - c\sigma_a \bar{\eta} e + \int_{4\pi} c\sigma_a \frac{I}{4\pi} d\omega', \\ 
\qquad + c_m \sigma_m \int_{4\pi} P_m \frac{e}{4\pi} d\omega' d\sigma_m'. 
\end{cases} 
\quad (13) $$

Independently of the choice of $\bar{\eta}$, system (13) is linear ($\bar{\eta}$ is a function of $x, t$ but not anymore of $E$). An MC scheme can then be applied to solve it on time step $[t^n, t^n + \Delta t]$. Of course, in practice, a choice must be made for $\bar{\eta}$ (explicit as in [42], implicit as in [40]): let us postpone this choice and rather go through a few properties of the above linearisation (which remains independent of the choice of $\bar{\eta}$).

**Remark 1.** At this stage of the discussion it is important noting that if one has already access to an MC solver in order to solve the radiative transfer Eq. (4), a simple modification of the solver can be made in order to solve system (13). This is what we are going to tackle in the next lines.

As explained before, system (13) is now linear and can be solved with an MC scheme. We here insist on the fact that it has the same structure as a multigroup transport equation (for neutronics for example, see [69,70]): in this context, a P-truncated basis $(\phi_0(v), \dots, \phi_P(v))^T$ is introduced and the solution $u$ is decomposed on this basis $\psi(x, t, \omega, v) \approx \sum_{k=0}^P u_k(x, t, \omega) \phi_k(v)$. But in our case, there are $M+1$ groups with $M$ the number of different materials within the layers and the basis functions are analytical, given by $\phi_0(v) = \delta_c(v)$ and $(\phi_l(v) = \delta_{c_m,l}(v))_{l \in \{1, \dots, M\}}$. To clarify this point, let us put forward the expression of the total and scattering opacities of linear system (13). It consists in recovering the equation whose solution is $\psi(x, t, \omega, v) = I(x, t, \omega)\delta_c(v) + \sum_{l=1}^M e(x, t, \omega)\delta_{c_m,l}(v)$.

For the sake of conciseness of the calculations, we assume that $c_m$ can be continuously distributed even if in practice, it will have $M+1$ discrete values specified by the characteristics of the different layers of the media, see Section 4. Hence, $\psi$ simplifies to $\psi(x, t, \omega, v) = I(x, t, \omega)\delta_c(v) + e(x, t, \omega)\delta_{c_m,l}(v)$. In the latter expression, $\delta_{c_m,l}$, $\delta_c$ are such that

$$ \int_{\{V\}} \delta_c(v) dv = \delta_{V,c} \quad \text{and} \quad \int_{\{V\}} \delta_{c_m}(v) dv = \delta_{V,c_m}, \quad (14) $$

where $\delta_{V,c}$ is the Kronecker symbol<fnref n="3" /> and $\{V\}$ denotes the singleton $V$. In fact, $v$ is nothing more than a velocity which can be $c$ for photons or $c_m$ for sensible heat. Let us now build the linear equation satisfied by $\psi$. Expression (13) can be rewritten (we drop the dependences for conciseness):

$$ \begin{aligned} \partial_t (I \delta_c + e \delta_{c_m}) + v \omega \cdot \nabla (I \delta_c + e \delta_{c_m}) = & - (c \sigma_t I \delta_c + c \sigma_a \bar{\eta} e \delta_{c_m} + c_m \sigma_m e \delta_{c_m} \\ & + (c \sigma_a \delta_c \bar{\eta} \int_{4\pi} e + c \sigma_a \delta_c \int_{4\pi} I \\ & + c_m \sigma_m \delta_{c_m} \int_{4\pi} P_m e + c \sigma_s \delta_c \int_{4\pi} I) \end{aligned} $$

It remains to make $\psi$ appear in the collisional part. For the moment the integration is only over the angular distribution. Let us introduce the Kronecker symbols

$$ \delta_{c_m,c}(v', v) = \delta_{c_m}(v') \delta_c(v) \quad \text{and} \quad \delta_{c,c_m}(v', v) = \delta_c(v') \delta_{c_m}(v). $$

We can then rewrite the above equation as:

$$ \begin{aligned} \partial_t \psi + v\omega \cdot \nabla \psi + (c\sigma_t \delta_c(v) + c\sigma_a \bar{\eta} \delta_{c_m}(v) + c_m \sigma_m \delta_{c_m}(v))\psi = & \\ & + \int_V \int_{4\pi} \left[ c\sigma_a \delta_{c,c_m}(v', v) + c\sigma_s \delta_{c,c}(v', v) + c\sigma_a \bar{\eta} \delta_{c_m,c}(v', v) + c_m \sigma_m P_m(v', v) \right] \\ & \times \psi(v', \omega') dv' d\omega', \end{aligned} \quad (15) $$

with $\mathcal{V} = \{c_m, c\}$ a discrete set with only two elements. One can check that performing $\int_{\{c\}}$ (15) $dv$ allows recovering the first line of system (13) and that $\int_{\{c_m\}}$ (15) $dv$ leads to the second equation of system (13). We can identify the scattering $v\Sigma_s(v, v')$ and total $v\Sigma_t(v)$ opacities to rewrite the system above under the general form:

$$ \partial_t \psi + v\omega \nabla \psi + v\Sigma_t(v)\psi = \int_V \int_{4\pi} v\Sigma_s(v', v)\psi(v', \omega')dv'd\omega'. $$

In the above expression, we have:

$$ \begin{aligned} v\Sigma_t(v) &= c\sigma_t \delta_c(v) + c\sigma_a \bar{\eta} \delta_{c_m}(v) + c_m \sigma_m \delta_{c_m}(v), \\ v\Sigma_s(v', v) &= c\sigma_a \delta_{c,c_m}(v', v) + c\sigma_s \delta_{c,c}(v', v) + c\sigma_a \bar{\eta} \delta_{c_m,c}(v', v) \\ &\quad + c_m \sigma_m P_m(v', v). \end{aligned} $$

Let us rewrite the scattering part as $v\Sigma_s(v)P_s(v', v) = v\Sigma_s(v', v)$. By doing so, we introduce $P_s$, the probability of being scattered with velocity $v$, incoming with velocity $v'$: it is mainly used for a backward resolution. This implies

$$ \begin{aligned} v\Sigma_s(v) &= \int v\Sigma_s(v', v)dv', \\ &= \int \left[ c\sigma_a \delta_{c,c_m}(v', v) + c\sigma_s \delta_{c,c}(v', v) + c\sigma_a \bar{\eta} \delta_{c_m,c}(v', v) \right. \\ & \qquad \left. + c_m \sigma_m P_m(v', v) \right] dv', \\ &= (c\sigma_s + c\sigma_a \bar{\eta})\delta_c(v) + (c\sigma_a + c_m \sigma_m)\delta_{c_m}(v). \end{aligned} $$

By definition of $P_s$ we have:

$$ \begin{aligned} P_s(v', v) &= \frac{v\Sigma_s(v', v)}{v\Sigma_s(v)}, \\ &= \frac{[c\sigma_a \delta_{c,c_m}(v', v) + c\sigma_s \delta_{c,c}(v', v) + c\sigma_a \bar{\eta} \delta_{c_m,c}(v', v) + c_m \sigma_m P_m(v', v)]}{(c\sigma_s + c\sigma_a \bar{\eta})\delta_c(v) + (c\sigma_a + c_m \sigma_m)\delta_{c_m}(v)}. \end{aligned} $$

The above expression can be considerably simplified by noticing that

$$ \text{for } v = c_m, \quad P_s(c_m, v') = \frac{c\sigma_a \delta_c(v') + c_m \sigma_m P_m(v', c_m)}{(c\sigma_a + c_m \sigma_m)}, \quad \text{for } v = c, $$

$$ P_s(c, v') = \frac{\sigma_s \delta_c(v') + \sigma_a \bar{\eta} \delta_{c_m}(v')}{\sigma_s + \sigma_a \bar{\eta}}, $$

so that $P_s$ resumes to

$$ P_s(v, v') = \delta_{c_m}(v) \frac{c\sigma_a \delta_c(v') + c_m \sigma_m P_m(v', c_m)}{(c\sigma_a + c_m \sigma_m)} + \delta_c(v) \frac{\sigma_s \delta_c(v') + \sigma_a \bar{\eta} \delta_{c_m}(v')}{\sigma_s + \sigma_a \bar{\eta}}. $$

Now, we are interested in a direct resolution of (13) on time step $[t^n, t^n + \Delta t]$. In other words, we need (cf. [71]) to characterise the total and scattering opacities of the adjoint form of (13). It is given by

$$ - \partial_t \psi - v\omega \cdot \nabla \psi + v\Sigma_t(v)\psi = \int_V \int_{4\pi} v\Sigma_s(v', v)\psi(v', \omega')dv'd\omega', \quad (16) $$

where

$$ v\Sigma_s(v, v') = v\Sigma_s(v)P_s(v, v') = v'\Sigma_s(v')P_s(v, v'). $$

In the previous expression, $P_s$ is the probability of being scattered from velocity $v$ toward velocity $v'$: it is used for a direct resolution. Then we have

$$ \begin{aligned} v\Sigma_s(v) &= \int_V v' \Sigma_s(v') P_s(v, v') dv', \\ &= \int_V \left( (c\sigma_s + c\sigma_a \bar{\eta})\delta_c(v') + (c\sigma_a + c_m \sigma_m)\delta_{c_m}(v') \right) \times \\ & \quad \left( \delta_{c_m}(v') \frac{c\sigma_a \delta_c(v') + c_m \sigma_m P_m(v', c_m)}{(c\sigma_a + c_m \sigma_m)} \right. \\ & \quad \left. + \delta_c(v') \frac{c\sigma_s \delta_c(v) + c\sigma_a \bar{\eta} \delta_{c_m}(v)}{c\sigma_s + c\sigma_a \bar{\eta}} \right) dv' \\ &= (c\sigma_a + c\sigma_s)\delta_c(v) + (c\sigma_a \bar{\eta} + c_m \sigma_m)\delta_{c_m}(v) \\ &= c\sigma_t \delta_c(v) + (c\sigma_a \bar{\eta} + c_m \sigma_m)\delta_{c_m}(v) = \Sigma_t(v). \end{aligned} $$

<footnote n="3">³ i.e. is such that $\delta_{V,k} = 0$ if $V \neq k$ and $\delta_{V,k} = 1$ if $V = k$.</footnote>
Few calculations, similar to the one already performed to identify $P_S$, show that:

$$P_S(v, v') = \delta_{c_m}(v) \frac{c_m \sigma_m P_m(v', c_m) + c \sigma_a \bar{\eta} \delta_c(v')}{c \sigma_a \bar{\eta} + c_m \sigma_m} + \delta_c(v) \frac{\sigma_s \delta_c(v') + \sigma_a \delta_{c_m}(v')}{\sigma_s + \sigma_a}.$$

With the above calculations, we identified the direct and adjoint opacities for (13) to be revisited as a 2-group linear transport Eq. (16). By selecting $\bar{\eta}$ as

- $\bar{\eta} = \eta^n$, we use the SMC for the radiative counterpart as in [42],

- $\bar{\eta} = \eta^{n+1}$, we use the ISMC for the radiative counterpart as in [40].

In practice, in the next calculations, the ISMC solver is used. It natively degenerates toward SMC when the problem is not stiff and allows taking (stably) larger time steps than SMC for stiffer or stationary problems.

### 4. Choosing $c_m, \sigma_m, P_m$ with respect to $\kappa, \rho c_v$

Until now, we mainly explained how the coupled system (9) can be solved consistently, with Poissonian paths, independently of the choice of $c_m, \sigma_m, P_m$. Our aim now is to make some particular choices for the solution of system (9) to coincide with the solution of system (1). For this, we rely on the fact that in the diffusion limit [22] characterised by $c_m \sim \infty$ and $c_m \sigma_m \sim \infty$, system (9) tends towards

$$\left\{ \begin{aligned} & \frac{1}{c} \partial_t I + \omega \cdot \nabla I + \sigma_t I = \sigma_a B(T_m) + \sigma_s \int_{4\pi} I \frac{d\omega'}{4\pi}, \\ & \partial_t E(T_m) - \nabla \cdot \left( \frac{c_m}{3\sigma_m} \nabla E(T_m) \right) = \int_{4\pi} c \sigma_a \left( \frac{I}{4\pi} - B(T_m) \right) d\omega'. \end{aligned} \right. \quad (17)$$

The above system must be compared to

$$\left\{ \begin{aligned} & \frac{1}{c} \partial_t I + \omega \cdot \nabla I + \sigma_t I = \sigma_a B(T_m) + \sigma_s \int_{4\pi} I \frac{d\omega'}{4\pi}, \\ & \partial_t E(T_m) - \nabla \cdot \left( \frac{\kappa}{\rho c_v} \left[ \nabla E(T_m) - E(T_m) \nabla \ln(\rho c_v) \right] \right) = \int_{4\pi} c \sigma_a \left( \frac{I}{4\pi} - B(T_m) \right) d\omega', \end{aligned} \right. \quad (18)$$

which corresponds to system (1) in which the change of variable $T_m = \frac{E}{\rho c_v}$ has been performed in the spatial operator: $\nabla T_m$ in (1) has been replaced by $\nabla \frac{E}{\rho c_v}$.

We now aim at choosing $c_m, \sigma_m, P_m$ with respect to $\kappa, \rho c_v$ for the solution of system (17) to coincide with the solution of system (18). In order to constrain our choices, let us distinguish between two situations, using the fact that we rely on constant-per-cell conductivities and heat capacities:

(*) within a cell, $\frac{c_m}{3\sigma_m}$ must be close to $\frac{\kappa}{\rho c_v} = D$ and a particle “sees” a homogeneous media.

(**) $c_m, \sigma_m, P_m$ must be able to take into account the term $-D \nabla \ln(\rho c_v)$ which is non-zero only when a particle crosses an interface between two media with different heat capacities.

For condition (*) to be satisfied, we need to make sure that, additionally to having

$$\begin{aligned} c_m &\sim \infty, \\ c_m \sigma_m &\sim \infty, \end{aligned} \quad (19)$$

we must satisfy

$$D = \frac{c_m}{3\sigma_m} \quad \text{leading to} \quad \sigma_m = \frac{c_m}{3D}. \quad (20)$$

With the above expression, we can define $c_m$ and $\sigma_m$ up to a scaling factor $f$ such that $c_m = \tilde{c}_m f$ and $\sigma_m = \tilde{\sigma}_m f$ with $f \gg 1$ “large enough”. We will see in Section 6 what “large enough” must be in practice.

**Remark 2.** It is important noting that the Poissonian approach we suggest introduces $f$ as a numerical parameter (which echoes $\Delta t$ for Brownian paths). One important feature of $f$ is that convergence of (12) toward (1) is ensured as $f \to \infty$, see [72]. This means that any user of a simulation code implementing this method needs to perform a convergence study with respect to $f \to \infty$ in order to obtain reliable results. The value of $f$ leading to reliable results is certainly case dependent but as a guideline, with the numerical experiments we performed, as soon as $f > f_0$ with $f_0$ such that the averaged number of collisions per particle reaches $\approx 30-40$ during the time step, the results seems to be satisfactory. A study of the choice of $f$ relative to the choice of $N_{MC}$ is also provided in Section 6.1.2.

Constraints (19) and (20) are enough to satisfy (*). But they are not for the second condition (**) to hold: the difficulty comes from the fact that at the interface between two phases, $D \nabla \ln(\rho c_v)$ is not defined. We suggest relying on the preservation of some physical invariant at the interface, namely the continuity of temperature, in order to impose additional constraints on $c_m, \sigma_m$ with respect to $\rho c_v$. For this, let us consider a configuration at equilibrium (i.e. for $t \to \infty$ for example), without radiative transfer: in other words, in the case we have $\sigma_a = \sigma_s = \sigma_t = 0$. In this case, (9) degenerates towards

$$\left\{ \begin{aligned} & \partial_t I + c\omega \cdot \nabla I = 0, \\ & \partial_t e + c_m\omega \cdot \nabla e + c_m\sigma_m e = c_m\sigma_m \int_{4\pi} P_m e \frac{d\omega'}{4\pi} d\sigma_m', \\ & \text{with } E(T_m(x, t)) = \iint e(x, t, \omega, c_m) d\omega' d\sigma_m. \end{aligned} \right. \quad (21)$$

In this configuration, photons are not anymore coupled with sensible heat. Let us integrate with respect to both $\omega, c_m$ to get

$$\left\{ \begin{aligned} & \partial_t E(x, t) + \iint c_m \omega \cdot \nabla e = 0, \\ & \text{with } E(T_m(x, t)) = \iint e(x, t, \omega, c_m) d\omega' d\sigma_m. \end{aligned} \right. \quad (22)$$

In an infinite medium, we must be able to recover the fact that independently of $\kappa$, as $t \to \infty$, we must have a homogeneous temperature

$$T_m(x, t) = T_m, \forall x \in \Omega, t \to \infty. \quad (23)$$

Let us now consider two<fnref n="4" /> materials next to each other, with different densities and heat capacities. Note that such test-case has already been studied in order to enforce constraints in numerical methods in for similar problems, see [57]. We have

$$E(x) = \rho(x)c_v(x)T_m(x) = \rho(x)c_v(x)T_m, \quad (24)$$

which means that even if $T_m$ is constant as $t \to \infty$, $E$ may suffer a jump just because $x \to \rho(x)c_v(x)$ suffers a jump due to the presence of two layers of different materials. Let us consider we have material 1 and material 2 in presence, then we must have

$$T_m = \frac{E_1}{\rho_1 c_{v,1}} = \frac{E_2}{\rho_2 c_{v,2}}. \quad (25)$$

This means that equilibrium imposes an additional constraint on the flux at the interface between the two materials for equilibrium to remain fulfilled. Indeed, let us come back to (22) with additional hypothesis that $t \to \infty$, we then have

$$\iint c_m \omega \cdot \nabla e(x, \omega, c_m) d\omega' d\sigma_m = 0. \quad (26)$$

<footnote n="4">⁴ The generalisation to M materials is almost straightforward but the intuition is easier to get with only two.</footnote>
Let us integrate the above equation on two half spaces $\Omega = \Omega_1 \cup \Omega_2$

$$ \int_{\Omega} \int_{\Omega} c_m \omega \cdot \nabla e(x, \omega, c_m) d\omega d\omega c_m dx = 0, \quad (27) $$

$$ \int c_m \int \int_{\partial\Omega_1 \cup \partial\Omega_2} \omega \cdot n(x)e(x, \omega, c_m) dxd\omega d\omega c_m = 0, $$

where $n$ is the normal to the boundary $\partial\Omega = \partial\Omega_1 \cup \Omega_2$. We obtain

$$ \int c_m \int \int_{\partial\Omega_1} \omega \cdot n_1(x)e_1(\omega, c_m) dxd\omega d\omega c_m + \int c_m \int \int_{\partial\Omega_2} \omega \cdot n_2(x)e_2(\omega, c_m) dxd\omega d\omega c_m = 0, \quad (28) $$

$$ \int c_m \int \int_{\partial\Omega_1 \cap \partial\Omega_2} \omega \cdot n(x)e_1(\omega, c_m) dxd\omega d\omega c_m = \int c_m \int \int_{\partial\Omega_2 \cap \partial\Omega_1} \omega \cdot n(x)e_2(\omega, c_m) dxd\omega d\omega c_m. $$

Furthermore, in the equilibrium diffusion limit, $e$ is isotropic (i.e. $e(x, t, \omega, c_m) = e(x, t, c_m)$) leading to

$$ \frac{1}{4} \int c_m e_1(c_m) d\omega c_m = \frac{1}{4} \int c_m e_2(c_m) d\omega c_m. \quad (29) $$

The last equation can be reinterpreted as such: for equilibrium to be fulfilled, the mean with respect to the velocity of $c_m e(c_m)$ must be equal in all the materials. There are several ways to make sure this equality holds. Let us introduce the asymptotic expression ($t \to \infty$) into the previous equation to get

$$ \begin{aligned} \frac{1}{4} \int c_m E_1 d\omega c_m &= \frac{1}{4} \int c_m E_2 d\omega c_m, \\ \frac{1}{4} \int c_m \rho_1 c_{v,1} T_m d\omega c_m &= \frac{1}{4} \int c_m T_m \rho_2 c_{v,2} d\omega c_m, \\ \int c_m \rho_1 c_{v,1} d\omega c_m &= \int c_m \rho_2 c_{v,2} d\omega c_m. \end{aligned} \quad (30) $$

In the next lines, we take $dc_m = \delta_{c_{m,1}}(c_m)dc_m$ in material 1 and $dc_m = \delta_{c_{m,2}}(c_m)dc_m$ in material 2, i.e. we enforce the relation quite naively by making sure that the average velocities $c_{m,1}$ and $c_{m,2}$ in each materials verify

$$ c_{m,1} \rho_1 c_{v,1} = c_{m,2} \rho_2 c_{v,2}. \quad (31) $$

More elaborated strategies may present nicer properties but this kind of study is beyond the scope of this paper. In practice, for $M$ materials (or $N$ cells), we will have

$$ \forall i \in \{1, \dots, M\}, c_{m,i} \rho_i c_{v,i} = \bar{c}_m f. \quad (32) $$

It is easy to check that the last equation is compatible with the previous ones given by (19)–(20). Finally, it remains to choose $P_m$. In practice, we choose to take

$$ P_m(x, v, v') = \delta_{c_{m,i}}(v) \mathbf{1}_{\Omega_i}(x). \quad (33) $$

In other words, a sensible heat MC particle incoming in medium $\Omega_i$ has velocity $c_{m,i}$ after its first collision in $\Omega_i$ and if enough collisions are made, its average values within the material will get closer and closer to $c_{m,i}$.

To sum up, given some material physical quantities ($\kappa_i, \rho_i, c_{v,i})_{i \in \{1, \dots, M\}}$, we choose $\forall i \in \{1, \dots, M\}$

$$ \left\{ \begin{array}{l} c_{m,i} = \frac{\bar{c}_m f}{\rho_1 c_{v,i}}, \\ \sigma_i = \frac{c_{m,i}}{3D_i} = \frac{\bar{c}_m f}{3D_i \rho_i c_{v,i}}, \\ P_i(x, v, v') = \delta_{c_{m,i}}(v) \mathbf{1}_{\Omega_i}(x). \end{array} \right. \quad (34) $$

Parameter $f$ is the only numerical parameter and it guarantees the convergence of (12) toward (1) as it tends to infinity. It must be used in order to satisfy condition (19). It must be acknowledged that satisfying condition (19) may lead to computationally intensive calculations as the larger $f$ is, the stiffer $\sigma_m$ is and the more frequent the collisions and the longer the calculations are for an MC resolution. In the numerical

Section 6, we will see that despite the previous con, the strategy remains interesting.

We finally insist on the fact that other choices could be made for $\sigma_m$. For example, one could introduce a dependence of $\sigma_m$ with respect to $v$. For example in order to enforce a collision as soon as possible when a sensible heat MC particle enters medium $\Omega_m$ with velocity $c_{m'}$ with $m' \neq m$ (incoming from another medium $\Omega_{m'})$, we can imagine an opacity of the form

$$ \sigma_m(v) = \sigma_m \delta_{c_m}(v) + \frac{1}{\epsilon}(1 - \delta_{c_m}(v)), \quad \text{with } \epsilon \sim 0. \quad (35) $$

This topic may be the scope of further studies aiming at accelerating our new solver. For the moment, we remain with the simple choices (34) and show that they already allow significant improvements (Section 6).

## 5. Revisiting the Brownian paths resolutions

Having the previous compatibility results in mind, it is possible to revisit the Brownian paths resolution. The Brownian paths ensure condition (*) is satisfied. It remains to work on condition (**) which was ensured via constraints on the velocities of the particles, see (32). The velocity of particles has no sense for Brownian particles. But displacements can be compared and the same analysis at equilibrium holds. With this in mind, let us compare a Brownian displacements and a Poissonian one: a Brownian displacement on time step $\Delta t$ is proportional to $\sqrt{2D\Delta t}$. On this same time step, the equivalent velocity is:

$$ c_m \propto \sqrt{\frac{2D}{\Delta t}}. \quad (36) $$

Let us now apply the same conditions on the velocities as in the previous case, i.e. $c_{m,1} \rho_1 c_{v,1} = c_{m,2} \rho_2 c_{v,2}$, we then get:

$$ \sqrt{\frac{2\kappa_1}{\rho_1 c_{v,1} \Delta t_1}} \rho_1 c_{v,1} = \sqrt{\frac{2\kappa_2}{\rho_2 c_{v,2} \Delta t_2}} \rho_2 c_{v,2} \Rightarrow \frac{\Delta t_1}{\Delta t_2} = \frac{\kappa_1 \rho_1 c_{v,1}}{\kappa_2 \rho_2 c_{v,2}}. \quad (37) $$

In other words, the ratio of the time steps is the ratio of the squared effusivities. Now introduce a reference time step $\Delta t_{ref}$: the velocity condition for continuity at the boundary translates into a condition on the time steps used in each subdomains as

$$ \kappa_i \rho_i c_{v,i} \Delta t_{ref} = \Delta t_i, \quad \forall i \in \{1, \dots, M\}. \quad (38) $$

In the next section, this numerical strategy will be studied as well as the transport (Poissonian) one. We will see that the previous simple analogy leads to efficient corrections of the Brownian path. Convergence as $\Delta t_{ref} \to 0$ will even be observed.

## 6. Numerical results

In this section, we apply the material of the previous sections to several test-cases. In Section 6.1, we begin by revisiting the benchmarks of Section 1, which motivated our study, but with the strategies described in Sections 3–4 and 5. The results are presented in Sections 6.1.1 and 6.1.2. The benchmark of Section 6.1.1 (with constant diffusivities for which we have an analytical solution) studies the convergence speed of the phonon model with respect to $f$. The benchmark of Section 6.1.2 allows focusing on conditions (19)–(20) for the phonon model resolution and on condition (38) for the Brownian one in the case of a discontinuous conductivity $x \to \kappa(x)$. In Section 6.1.3, we consider a problem without conductivity jump but with discontinuous $\rho c_v$. In Section 6.2, we test both of our suggested strategies for the resolution of the full system (1) in the equilibrium diffusion limit for photons and phonons, in layered media (discontinuous conductivities – Section 6.2.1 – volumic heat capacity – Section 6.2.3 – or both at the same time – Section 6.2.3): these test-cases are convenient because obtaining a reference solution in the equilibrium diffusion limit is quite easy with a deterministic solver. Finally, in Section 6.4, we consider
Kinetic scheme of sections 3–4

<figure><img src="image_5.png" /><figcaption>Fig. 2. Results in terms of temperature profiles $x \mapsto T_m(x, t = 10^{-2})$ obtained with the kinetic solver (left column) and the diffusive one (right column) on the constant diffusivity benchmark of Section 1. Several solutions are displayed with different numerical parameters $f \in \{10^1, 10^2, 10^3\}$ for the Kinetic scheme and $\Delta t_{ref} \in \{10^{-5}, 10^{-6}, 10^{-7}\}$ for the Brownian scheme.</figcaption></figure>

Brownian scheme of section 5

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

some fully coupled physics with transparent media in which radiative transfer cannot be approximated with the diffusion model.

Note that all the test-cases are such that the solution does not interact with the boundaries. This is mainly because we suggest ways to overcome a problem which is more related to jumps of physical quantities rather than boundary conditions. As such, the test-cases are simpler and this eases the reproducibility of the results. Note that taking into account the boundary conditions for the kinetic solver is very close to taking into account boundary condition for radiative transfer (4) solved with an MC scheme. For the corrected Brownian path, the treatment of the boundary remains unaffected by the corrections we suggest in Section 5.

In the next sections, the resolution strategy described in Sections 3–4 is referred to as the kinetic solver or the solution with Poissonian paths. The strategy consisting in Brownian paths with the time step corrections described in Section 5 is referred to as the diffusive solver or the solution with Brownian paths.

## 6.1. Uncoupled benchmarks in layered media: with only diffusion for material temperature

In this section, we revisit the benchmarks of Section 1 but with the new resolution strategies described in Sections 3–4 and 5. We focus on the resolution of the diffusion equation for material temperature (not yet coupled to radiative transfer):

$$ \partial_t E(x, t) - \nabla \cdot (\kappa(x) \nabla T_m(x, t)) = 0, \quad (39) $$

with $\partial_t E(x, t) = \partial_t E(T_m(x, t)) = \rho(x)c_v(x)\partial_t T_m(x, t)$. The test-case of Section 6.1.1 considers constant $\kappa$, $\rho c_v$ and studies the sensitivity of $f$ (Poissonian) and to $\Delta t_{ref}$ (Brownian) on this simple benchmark. In Section 6.1.2, we focus on a problem with discontinuous $x \mapsto \kappa(x)$ and in Section 6.1.3, on a problem with discontinuous $x \mapsto \rho(x)c_v(x)$.

### 6.1.1. Constant diffusivities ($D_1 = D_2$) for heat conduction (no coupling with radiative transfer)

In this section, we revisit the simple test-case of Section 1 (simple case $D_1 = D_2$). We quantitatively exploit the fact that an analytical solution is available and we study the convergence rate of the kinetic strategy toward the equilibrium diffusion limit. Fig. 2 presents the results obtained with the kinetic solver (left column) and the diffusive one (right column) on the constant diffusivity benchmark. The results are displayed in terms of temperature profiles $x \mapsto T_m(x, t = 10^{-2})$. For the Brownian paths (right), to reach census, i.e. final time $t = 10^{-2}$, each MC particle encounters

<figure><img src="image_7.png" /><figcaption>Fig. 3. Behavior of the errors with respect to $N_{MC}$ and $f$ for the Poissonian path and to $N_{MC}$ and $\Delta t_{ref}$ for the Brownian one quantified in terms of total scattering events.</figcaption></figure>

- $10^2$ collisions for $\Delta t_{ref} = 10^{-5}$,

- $10^3$ collisions for $\Delta t_{ref} = 10^{-6}$,

- $10^4$ collisions for $\Delta t_{ref} = 10^{-7}$.

As soon as $\Delta t_{ref} \le 10^{-6}$, equivalent to each particle performing more than 1000 collisions, the results show a good agreement. The top left picture of Fig. 2 displays the results obtained with the Poissonian paths for three values of the factor $f \in \{10^1, 10^2, 10^3\}$. For $f = 10^1$, the equilibrium diffusion limit is not reached and the results are similar to the ones we obtain for transport in a semi transparent medium. As soon as $f \ge 100$, the results are in agreement with the reference.

Fig. 3 presents some more quantitative results: it displays the evolutions of the errors obtained with the Poissonian and Brownian methods with respect to the number of scattering events (related to $f$ for the Poissonian strategy and to $\Delta t_{ref}$ for the Brownian one) for two different MC discretisations (respectively $N_{MC} = 10^5$ and $N_{MC} = 10^7$). On Fig. 3, we can see that the error for the Brownian method does not depend on the number of scattering events (related to $\Delta t_{ref}$) but rather more directly on the number of MC particles used to solve the problem. The convergence with respect to $\Delta t_{ref}$ is fast. On another hand, for the Poissonian method, by increasing $f$, the number of scattering events increases, the diffusive regime is better and better captured and the error tends toward the one obtained with the Brownian scheme, dictated by the number of MC particles $N_{MC}$. The convergence seems to be of order 1 with respect to the number of scattering events. As a matter of fact, when the errors match, the two solvers are in exactly the
Kinetic scheme of sections 3–4

<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Results in terms of temperature profiles $x \mapsto T_m(x, t = 10^{-2})$ obtained with the kinetic solver (left column) and the diffusive one (right column) on the discontinuous conductivity $x \mapsto \kappa(x)$ benchmark of Section 1. Several solutions are displayed with different numerical parameters $f \in \{10^1, 10^2, 10^3\}$ for the Kinetic scheme and $\Delta t_{ref} \in \{10^{-4}, 10^{-5}, 10^{-6}\}$ for the Brownian scheme.</figcaption></figure>

same conditions (same $N_{MC}$, same time step, same number of cells and even same accuracies) and their run times are the almost the same. This point will be more developed in the next section but for discontinuous conductivities.

### 6.1.2. Discontinuous conductivities ($\kappa_1 \neq \kappa_2$ and $\rho c_{v,1} = \rho c_{v,2}$) for heat conduction (no coupling with radiative transfer)

In this section, we revisit the test-case of Section 1 for which $x \mapsto \kappa(x)$ is discontinuous. It is given by $x \mapsto \kappa(x) = 10 \times \mathbf{1}_{[0,1]}(x) + 1 \times \mathbf{1}_{[1,1]}(x)$. For this test-case, $x \mapsto \rho(x)c_{v}(x) = 1, \forall x \in [0, 1]$. We tackle it with the material of the previous Sections 3–4 and 5.

Fig. 4 presents the results obtained with the kinetic solver (left column) and the diffusive one (right column) on the discontinuous conductivity $x \mapsto \kappa(x)$ benchmark of Section 1. The results are displayed in terms of temperature profiles $x \mapsto T_m(x, t = 10^{-2})$. Note that the profiles $x \mapsto \kappa(x)$ and $x \mapsto \rho(x)c_{v}(x)$ are also displayed (the scales for them are on the right vertical axis).

The solutions are obtained in the same numerical settings: $N_x = 100$ cells, $N_{MC} = 10^6$. The two resolution strategies give satisfactory results as soon as $f = 10^2$ for the kinetic solver and $\Delta t = 10^{-6}$ for the diffusion one. For both cases of Fig. 4, the new strategies we suggest in this paper are satisfactory in the sense they converge toward the physical solution as their respective numerical parameters ($f$ and $\Delta t_{ref}$) are refined.

The previous results are satisfactory but one may wonder how the fluctuations induced by the MC resolution (i.e. due to the stochasticity of the solvers) impact the results. In the next lines, we aim at proving that the results of Fig. 4 are significantly sound. For this, let us now tackle some performance tests in this discontinuous context: we cannot rely anymore on any analytical reference solution in this case. For this reason, we suggest comparing the two solvers with respect to their mean and variance of the temperature profile obtained over $N_{seed} = 100$ computations (i.e. by changing 100 times the initial seed of the random number generator of the stochastic simulators).

Fig. 5 presents the results obtained with $N_{MC} = 10^5$ for both solvers: with the coarse choices of numerical parameters ($f = 10$) for the kinetic solver (Fig. 5 top left), the reference solution obtained by the deterministic solver is not within the 95% confidence intervals computed thanks to the $N_{seed}$ simulations. For the Brownian scheme with coarse time step $\Delta t_{ref} = 10^{-4}$ (Fig. 5 top right), the reference temperature profile is already within the 95% confidence intervals but the mean of the simulations remains quite far from it.

For the middle line of Fig. 5, the numerical parameters are refined and the reference solution is now within the 95% confidence intervals for both solvers. Note that the averaged simulation obtained with the kinetic solver is in better agreement with the reference solution than the one obtained with the Brownian one.

Brownian scheme of section 5

<figure><img src="image_9.png" /><figcaption>**Table 1**</figcaption></figure>

Execution times averaged over the number of seeds ($N_{seed} = 100$) for both solvers in order to obtain similar accuracies on the mean and the variance.

<table><thead><tr><th>N<sub>MC</sub> = 10<sup>5</sup></th><th>Kinetic</th><th>Brownian</th></tr></thead><tbody><tr><td>f = 10<sup>2</sup> vs. Δt<sub>ref</sub> = 10<sup>-5</sup></td><td>0.296 s.</td><td>0.995 s.</td></tr><tr><td>f = 10<sup>3</sup> vs. Δt<sub>ref</sub> = 10<sup>-6</sup></td><td>2.622 s.</td><td>5.972 s.</td></tr></tbody><thead><tr><th>N<sub>MC</sub> = 10<sup>6</sup></th><th>Kinetic</th><th>Brownian</th></tr></thead><tbody><tr><td>f = 10<sup>2</sup> vs. Δt<sub>ref</sub> = 10<sup>-5</sup></td><td>1.095 s.</td><td>6.712 s.</td></tr><tr><td>f = 10<sup>3</sup> vs. Δt<sub>ref</sub> = 10<sup>-6</sup></td><td>19.90 s.</td><td>63.24 s.</td></tr></tbody></table>

On the last line of Fig. 5, the numerical parameters $f$ and $\Delta t_{ref}$ are once again refined and for these last choices, the means obtained with both solvers are in excellent agreement with the reference solution. Furthermore, the 95% confidence intervals of both solvers are significantly the same. This typically means that we can compare the solvers for similar accuracies in terms of mean and variance for such choices of parameters. Note that such comparisons could not have been obtained without our corrected Brownian path (as the solutions of both solvers must converge toward the same limit for fair comparisons).

The purpose of Table 1 is precisely to exploit the possibility to perform fair comparisons but before commenting on Table 1, let us notice that Fig. 6 displays the same study as Fig. 5 but with $N_{MC} = 10^6$ instead of $N_{MC} = 10^5$. All in all, similar comments as for Fig. 5 can be made. One can notice that the 95% confidence intervals are now tighter so that even for the coarse numerical parameters (top line of Fig. 5), the reference solution is outside these bounds. For the middle plots, the reference solution is within the confidence intervals for both solvers but we can notice that the mean of the kinetic model is in better agreement with the reference solution than the one of the Brownian path (just as in Fig. 5 with $N_{MC} = 10^5$): this testifies that increasing $N_{MC}$ cannot compensate the need for decreasing $\Delta t_{ref}$ in order to eliminate a bias. But Fig. 6 allows another important comment: the optimal choice of $f$ (kinetic) or $\Delta t_{ref}$ (Brownian) ensuring the true solution to be within the confidence bounds strongly depends on the value of $N_{MC}$. Being able to increase the number of walkers (kinetic or Brownian) may impose finer $f$ or $\Delta t_{ref}$ for the bias induced by coarse choices of such parameters to be negligible with respect to the fluctuations induced by $N_{MC}$. For this reason, it is important identifying the faster solver between the kinetic one and the corrected Brownian one.

Table 1 compares the execution times averaged over the number of seeds ($N_{seed} = 100$) for both solvers in order to obtain similar accuracies on the mean and the variance (i.e. for the numerical parameters of the middle and bottom lines of Figs. 5 and 6). The kinetic solver is about 3 times faster than the corrected Brownian one for this benchmark (the
Kinetic scheme of sections 3–4

$$N_{MC} = 10^5$$

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

Brownian scheme of section 5

$$N_{MC} = 10^5$$

<figure><img src="image_11.png" /><figcaption>Fig. 5. Results in terms of temperature profiles $x \mapsto T_m(x, t = 10^{-2})$ (several seeds, mean, confidence intervals) obtained with the kinetic solver (left column) and the diffusive one (right column) on the discontinuous conductivity $x \mapsto \kappa(x)$ benchmark of Section 1 for $N_{seed} = 100$ different seeds and $N_{MC} = 10^5$. Each line corresponds to a numerical parameter choice: $f \in \{10^1, 10^2, 10^3\}$ for the Kinetic scheme and $\Delta t_{ref} \in \{10^{-4}, 10^{-5}, 10^{-6}\}$ for the Brownian scheme.</figcaption></figure>

first line is non informative, the execution times are too small for fair comparisons).

### 6.1.3. Discontinuous density/volumic heat capacity ($\kappa_1 = \kappa_2$ and $\rho c_{v,1} \neq \rho c_{v,2}$) for heat conduction (no coupling with radiative transfer)

In this section, we present some results on a benchmark for which there are discontinuities for $x \mapsto \rho(x)c_v(x)$ but not for $x \mapsto \kappa(x)$. We take $x \mapsto \kappa(x) = 1, \forall x \in [0, 1]$ and $x \mapsto \rho(x)c_v(x) = 10^{-1} \times \mathbf{1}_{[0, \frac{1}{2}]}(x) + 1 \times \mathbf{1}_{[\frac{1}{2}, 1]}(x)$. Once again, as a reference, we rely on the solution provided by a deterministic solver.

The results are displayed in Fig. 7. For the Poissonian numerical strategy (left), as soon as the numerical parameter is well chosen ($f \gg 1$), the exact solution is captured. The factor $f$ must be equal to $10^4$ for accurate restitutions. The new Brownian strategy also allows capturing

efficiently this problem even if a spike at $x = \frac{1}{2}$, i.e. at the interface between the two material is observable for the coarse values of $\Delta t_{ref}$. We insist on the fact that the numerical parameters $N_x = 100$ and $N_{MC} = 10^6$ are the same for both resolution strategies: similar levels of noise are qualitatively observable.

### 6.2. Coupled benchmark (resolution of system (1)) in the equilibrium diffusion limit in layered media

In this section, we solve the full system (1) coupling radiative transfer and the diffusion equation for material temperature in layered media but in the equilibrium diffusion limit, i.e. in opaque materials for photons. In such conditions [73–75], the solution of (1) coincides with
<figure><img src="image_12.png" /><figcaption>Brownian scheme of section 5
$N_{MC} = 10^6$</figcaption></figure>

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

**Fig. 6.** Results in terms of temperature profiles $x \mapsto T_m(x, t = 10^{-2})$ (several seeds, mean, confidence intervals) obtained with the kinetic solver (left column) and the diffusive one (right column) on the discontinuous conductivity $x \mapsto \kappa(x)$ benchmark of Section 1 for $N_{seed} = 100$ different seeds and $N_{MC} = 10^6$. Each line corresponds to a numerical parameter choice: $f \in \{10^1, 10^2, 10^3\}$ for the Kinetic scheme and $\Delta t_{ref} \in \{10^{-4}, 10^{-5}, 10^{-6}\}$ for the Brownian scheme.

<figure><img src="image_14.png" /><figcaption>Brownian scheme of section 5</figcaption></figure>

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

**Fig. 7.** Results in terms of temperature profiles $x \mapsto T_m(x, t = 10^{-2})$ obtained with the kinetic solver (left column) and the diffusive one (right column) on a benchmark with discontinuous $x \mapsto \rho(x)c_v(x)$. Several solutions are displayed with different numerical parameters $f \in \{10^1, 10^2, 10^3\}$ for the Kinetic scheme and $\Delta t_{ref} \in \{10^{-5}, 10^{-6}, 10^{-7}\}$ for the Brownian scheme.

the solution of

$$ \partial_t (E(T_m(x,t)) + aT_m^4(x,t)) - \nabla \cdot \left( \kappa(x) \nabla T_m(x,t) + \frac{c}{3\sigma_a(x)} \nabla (aT_m^4(x,t)) \right) = 0 \quad (40) $$

with Dirichlet boundary conditions $T_m = 0$ and the same Heaviside initial condition as in the example of Section 1. In order to solve the above model, we rely on a deterministic solver.
Kinetic scheme of sections 3–4

<figure><img src="image_16.png" /><figcaption>**Fig. 8.** Spatial profile of the material temperature for two different solvers in the same configuration. The smooth curve comes from a deterministic simulation code solving the equilibrium diffusion limit for both radiative transfer and material conduction at time $t_f = 5 \times 10^{-3}$. The noisy curves are obtained from the phonon solver (left) and the corrected Brownian one (right).</figcaption></figure>

<table><caption>**Table 2**

Characteristics of the benchmark of Section 6.2.1.</caption><thead><tr><th>Ω = [0, 1]</th><th>Ω<sub>1</sub> = [0, 1/2]</th><th>Ω<sub>2</sub> = [1/2, 1]</th></tr></thead><tbody><tr><td>κ(x)</td><td>10<sup>-1</sup></td><td>1</td></tr><tr><td>ρ(x)c<sub>v</sub>(x)</td><td>1</td><td>1</td></tr><tr><td>σ<sub>a</sub>(x)</td><td>5000</td><td>5000</td></tr><tr><td>σ<sub>s</sub>(x)</td><td>0</td><td>0</td></tr></tbody></table>

### 6.2.1. Resolution of system (1) in the equilibrium diffusion limit in layered media (discontinuous κ)

In this section, we suggest considering a test problem with a full coupling between radiative transfer and diffusion in a layered material in the *equilibrium diffusion limit*. The solution of (1) in this regime coincides with the solution of (40). The characteristics of the test-case are summed up in Table 2. We insist on the fact that these test-cases are more numerical than physical: we take $c = 5000$ for the speed of light and $a = 1$ for the radiative constant. In practice, for the numerical parameters, we take $\Delta t = 10^{-4}$, $N_x = 100$ and $N_{MC} = 10^6$ for both computations. The initial condition is given by $T_m(x, t=0) = 2 \times I_{[0,45,0,55]}(x)$. The final time $t_f = 5 \times 10^{-3}$ has been chosen so that the solutions do not interact significantly with the boundary in order to avoid additional difficulties which would be beyond the scope of this paper.

Fig. 8 displays the results in terms of temperature profiles $x \mapsto T_m(x, t_f = 5 \times 10^{-3})$. They are obtained with both strategies suggested in this paper, compared to a reference solution given by a deterministic code solving (40). Once again, even in this coupled context, the two numerical strategies we suggest in this paper allow recovering the reference solution as soon as $f$ and $\Delta t_{ref}$ are fine enough.

### 6.2.2. Resolution of system (1) in the equilibrium diffusion limit in layered media (discontinuous $\rho c_v$)

In this section, we suggest considering a test problem with a full coupling between radiative transfer and diffusion in a layered material in the *equilibrium diffusion limit*. The characteristics of the test-case are summed up in Table 3. Once again, we take $c = 5000$ for the speed of light and $a = 1$ for the radiative constant and $\Delta t = 10^{-4}$, $N_x = 100$ and $N_{MC} = 10^6$ for both computations. The initial condition is given by $T_m(x, t=0) = 2 \times I_{[0,45,0,55]}(x)$. The final time is $t_f = 1 \times 10^{-3}$.

Fig. 9 displays the results in terms of temperature profiles $x \mapsto T_m(x, t_f = 10^{-3})$. They are obtained with both strategies suggested in this paper, compared to a reference solution given by a deterministic code solving (40). Once again, even in this coupled context, the two numerical strategies we suggest in this paper allow recovering the reference solution as soon as $f$ and $\Delta t_{ref}$ are fine enough.

Brownian scheme of section 5

<figure><img src="image_17.png" /><figcaption>**Table 3**

Characteristics of the benchmark of Section 6.2.2.</figcaption></figure>

<table><caption>**Table 4**

Characteristics of the benchmark of Section 6.2.3.</caption><thead><tr><th>Ω = [0, 1]</th><th>Ω<sub>1</sub> = [0, 1/2]</th><th>Ω<sub>2</sub> = [1/2, 1]</th></tr></thead><tbody><tr><td>κ(x)</td><td>1</td><td>1</td></tr><tr><td>ρ(x)c<sub>v</sub>(x)</td><td>10</td><td>1</td></tr><tr><td>σ<sub>a</sub>(x)</td><td>5000</td><td>5000</td></tr><tr><td>σ<sub>s</sub>(x)</td><td>0</td><td>0</td></tr></tbody></table>

<table><thead><tr><th>Ω = [0, 1]</th><th>Ω<sub>1</sub> = [0, 1/2]</th><th>Ω<sub>2</sub> = [1/2, 1]</th></tr></thead><tbody><tr><td>κ(x)</td><td>10<sup>-1</sup></td><td>1</td></tr><tr><td>ρ(x)c<sub>v</sub>(x)</td><td>10</td><td>1</td></tr><tr><td>σ<sub>a</sub>(x)</td><td>5000</td><td>5000</td></tr><tr><td>σ<sub>s</sub>(x)</td><td>0</td><td>0</td></tr></tbody></table>

### 6.2.3. Resolution of system (1) in the equilibrium diffusion limit in layered media (discontinuous $\rho c_v$ and $\kappa$)

In this section, we consider a test-case which is a mix between the two previous ones in the sense that both $x \mapsto \kappa(x)$ and $x \mapsto \rho(x)c_v(x)$ are discontinuous. The characteristics of the test-case are summed up in Table 4.

Fig. 10 displays the results obtained with both strategies suggested in this paper, compared to a reference solution given by a deterministic code the diffusion limit for both radiative transfer and material diffusion. Even in this more complex case, the two strategies are in very good agreement with the reference solution.

### 6.3. Configurations with layers of air, Al₂O₃ and graphite

The previous studies of the paper were mainly dimensionless. In this section, we consider more physical test-cases. The set-up is close to the one of the previous benchmarks, i.e. 1D, with several layers of different materials. For the material, we consider air, graphite and Al₂O₃. The information about the opacities, conductivities, densities and heat capacities are taken from [76,77] for graphite and from [78,79] for Al₂O₃ and are quite classical for air [80]. They are summed-up in Table 5 with their respective units. The speed of light is $c = 3 \times 10^8$ km s⁻¹.

Once again, the spatial configuration is quite simple in the sense that two materials are on the left hand side and the right hand side of the interface at $x = 1/2$. There are three different types of layers: an air-graphite one with final time $t_f = 50$ s, an Al₂O₃-graphite one with final time $t_f = 200$ s and an air-Al₂O₃-graphite one with final time $t_f = 2000$ s. The material temperature profiles $x \mapsto T_m(x, t_f)$ are displayed
Kinetic scheme of sections 3–4

<figure><img src="image_18.png" /><figcaption>**Fig. 9.** Spatial profile of the material temperature for two different solvers in the same configuration. The smooth curve comes from a deterministic simulation code solving the equilibrium diffusion limit for both radiative transfer and material conduction at time $t_f = 10^{-3}$. The noisy curves are obtained from the phonon solver (left) and the corrected Brownian one (right).</figcaption></figure>

Brownian scheme of section 5

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

Kinetic scheme of sections 3–4

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

Brownian scheme of section 5

<figure><img src="image_21.png" /><figcaption>**Fig. 10.** Spatial profile of the material temperature for two different solvers in the same configuration. The smooth curve comes from a deterministic simulation code solving the equilibrium diffusion limit for both radiative transfer and material conduction at time $t_f = 5 \times 10^{-2}$. The noisy curves are obtained from the phonon solver (left) and the corrected Brownian one (right).</figcaption></figure>

<table><caption>**Table 5**

Characteristics of the benchmark of Section 6.3 for air [80], graphite [76,77] and Al₂O₃ [78,79]. The properties are taken at temperature 1000 °C and atmospheric pressure.</caption><thead><tr><th></th><th>Units</th><th>Air [80]</th><th>Graphite [76,77]</th><th>Al<sub>2</sub>O<sub>3</sub> [78,79]</th></tr></thead><tbody><tr><td>κ</td><td>(W m<sup>-1</sup> K<sup>-1</sup>)</td><td>8.0 × 10<sup>-2</sup></td><td>7.5 × 10<sup>+1</sup></td><td>7.0 × 10<sup>+0</sup></td></tr><tr><td>ρ</td><td>(kg m<sup>-3</sup>)</td><td>1.0 × 10<sup>+0</sup></td><td>1.8 × 10<sup>+3</sup></td><td>3.9 × 10<sup>+3</sup></td></tr><tr><td>c<sub>v</sub></td><td>(J kg<sup>-1</sup> K<sup>-1</sup>)</td><td>1.0 × 10<sup>+3</sup></td><td>1.8 × 10<sup>+3</sup></td><td>1.2 × 10<sup>+3</sup></td></tr><tr><td>σ<sub>a</sub></td><td>(m<sup>-1</sup>)</td><td>0.0 × 10<sup>+0</sup></td><td>4.2 × 10<sup>+6</sup></td><td>1.0 × 10<sup>+4</sup></td></tr></tbody></table>

in Fig. 11 for the three configurations. The first line corresponds to the Al₂O₃-graphite layers and the second one to the air-graphite one. The last line corresponds to a set-up with three layers, the first one is air, the second is Al₂O₃ and the third one on the right is graphite. The initial condition is the same for every configurations, the same as in the previous sections.

Each time, two curves are plotted: the reference one is obtained with a deterministic scheme solving (40), the noisy one is obtained with the Poissonian path strategy (left) or the corrected Brownian one (right). Let us comment first on the first line of Fig. 11, for the Al₂O₃-graphite layers: both solvers are in excellent agreement with the reference curve. Both have exactly the same number of particles $N_{MC} = 10^6$, time step $\Delta t = t_f/10$. The corrected Brownian solver uses $\Delta t_{ref} = 10^{-10}$ whereas the Poissonian one uses $f = 10^4$: these values are chosen so that about the same accuracy is obtained for both solvers. To give an idea, in such configuration, the Poissonian strategy is about 100 times faster.

The second line of Fig. 11 presents the curves obtained for the air-graphite configuration. As can be attested, the resolution in air is much noisier than the resolution in graphite. This is mainly due to the important jump in $x \mapsto \rho(x)c_v(x)$: this fact has already been noticed in [64] in a quite different physical context: the ISMC resolution must be incriminated, and not the Poissonian or corrected Brownian strategies we suggest in this paper. The solver DIMC (Discrete IMC) may represent an interesting alternative to ISMC combined to the strategies suggested in this paper. We think the study of this phenomenon is interesting but beyond the scope of this paper. This will certainly be the object of future works. Note that the temperature profile obtained with the corrected Brownian path has 10 times less MC particles and the simulation takes 10 times more time than the Poissonian resolution. This test-case also allows highlighting one important property of the resolutions: even if the solution in air is noisy, the resolution in graphite is very well handled.

Finally, the last line of Fig. 11 presents the results on the set-up with three layers. For this set-up, the interfaces are at $x_{air-Al_2O_3} = \frac{1}{3}$ and $x_{Al_2O_3-graphite} = \frac{2}{3}$. First, the same remark as for the second line can be made: the resolution in air is much noisier than in the other materials. This is true for both the Poissonian and the corrected Brownian solver. Besides, the Poissonian simulation is much less costly than the corrected Brownian one: for their cost to match, we increased the time step $\Delta t_{ref} = 10^{-7}$ for the corrected Brownian path. We can clearly see that $\Delta t_{ref} = 10^{-7}$ is not small enough to damp the spike at the interface $x = \frac{1}{3}$ between air and Al₂O₃. It is much more difficult
Kinetic scheme of sections 3–4
Al₂O₃–graphite

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

Brownian scheme of section 5
Al₂O₃–graphite

<figure><img src="image_23.png" /><figcaption>**Fig. 11.** Spatial profiles of the material temperature for two different solvers (columns) in the same configurations (lines). The smooth curve comes from a deterministic simulation code solving the equilibrium diffusion limit for both radiative transfer and material conduction at time $t_f$. Note that for this benchmark, the values are not anymore dimensionless: data for air, Al₂O₃ and graphite were taken from the literature.</figcaption></figure>

to perform fair comparisons with such costly test-cases, that is why we relied on the choice of time step and factor *f* inspired from the previous simulations. With this test-case, we highlight the fact that the strategy is not limited to two layers and that the generalisation made in the previous sections for *M* materials does work. But it also highlights the fact that material per material tunings of *f* and $\Delta t_{ref}$ may be relevant in order to improve the accuracy or the computation time. This will certainly be the purpose of future studies.

6.4. Resolution of system (1) in semi-transparent layered media

In this last section, we revisit the test-case of Section 6.2.3 but the materials of the two layers are now semitransparent to photons (see Table 6). With this test-case, we want to insist on the importance of resorting to a transport model for photons in semitransparent material, especially in a calibration context. We compare the results obtained by three different models:

- the first model we consider is model (1). It is valid for semi-transparent materials to photons. It needs to be given quantities $\sigma_a, \sigma_s, \kappa, \rho, c_v$ for each materials involved in the computation (c is the speed of light and is fixed). In the next Fig. 12, it is solved with the kinetic solver (left) and with the corrected Brownian path (right). For this model (1), $x \to T_r(x, t_f)$ and $x \to T_m(x, t_f)$

are displayed in Fig. 12 with $t_f = 5 \times 10^{-2}$ the time of interest. The curves of Fig. 12 relative to this model are obtained by using the values of Table 6. Note that for this model, equilibrium may not be fulfilled in some (semitransparent) materials. This means that if we introduce the “radiative temperature” $T_r$, defined by $\int I(x, t, \omega, v) d\omega dv = aT_r^4(x, t)$, we may have $T_r \neq T_m$ for some $x \in D$. This is the case in Fig. 12 as equilibrium is fulfilled only for $x \le 0.4$ and $0.8 \le x$.

- The second model we consider is model (40). To run a computation with model (40), one only needs to feed the code with $\sigma_a, \rho, c_v, \kappa$ (c is the speed of light, a is the radiative constant and they are both fixed) as $\kappa_{rad} = \frac{c}{3\sigma_a}$ depends on the previous data. Note that $\sigma_s$ is not needed because the equilibrium diffusion limit is defined by $\sigma_a \sim \infty$ and in such condition $\sigma_t = \sigma_a + \sigma_s \sim \sigma_a$ and the model is not sensitive to $\sigma_s$. Model (40) is the limit of model (1) in opaque materials. In a sense, it means that in opaque media, the same $\sigma_a$ used to feed model (1) must be used to feed model (40) and the solutions of both models will coincide. In the plots of Fig. 12, the curve $x \to T_m(x, t_f)$ for this model is the one named reference (diffusion limit). In other words, the “reference (diffusion limit)” curve is obtained by solving (40) with the values of Table 6. The temperature profile obtained in the previously described condition is quite far from the temperature
Kinetic scheme of sections 3–4

<figure><img src="image_24.png" /><figcaption>**Fig. 12.** Spatial profile of the material temperature for two different solvers in the same configuration. The smooth curve comes from a deterministic simulation code solving the equilibrium diffusion limit for both radiative transfer and material conduction at time $t = 10^{-2}$. The noisy curves are obtained from the phonon solver (left) and the corrected Brownian one (right) but with radiative transfer in a semitransparent material.</figcaption></figure>

<caption>**Table 6**</caption>

Characteristics of the benchmark of Section 6.4. The speed of the photons is now $c = 5$, $a = 1$. The numerical parameters are $N_x = 100$, $N_{MC} = 10^6$ and $\Delta t = 10^{-4}$.

<table><thead><tr><th>Ω = [0, 1]</th><th>Ω<sub>1</sub> = [0, 1/2]</th><th>Ω<sub>2</sub> = [1/2, 1]</th></tr></thead><tbody><tr><td>κ(x)</td><td>10<sup>-1</sup></td><td>1</td></tr><tr><td>ρ(x)c<sub>v</sub>(x)</td><td>10</td><td>1</td></tr><tr><td>σ<sub>a</sub>(x)</td><td>5</td><td>5</td></tr><tr><td>σ<sub>x</sub>(x)</td><td>0</td><td>0</td></tr></tbody></table>

profiles obtained by solving model (1): the diffusion limit is not valid, the media being semitransparent. For this model, in this configuration, the temperature is over-evaluated in the vicinities of the propagation fronts and under-estimated in the vicinity of the center of the domain.

- Now, model (40) is cheaper to solve than model (1). It is tempting trying to replace $\frac{c}{3\sigma_a}$ by a fitted quantity in order to obtain temperature profiles as close as possible as the ones obtained with model (1) but with model (40). Let us introduce $\kappa_{rad}^{fitted} \neq \kappa_{rad}$. In practice, in Fig. 12, we take $\kappa_{rad}^{fitted} = 3.333 \times 10^{-3}$ instead of $\kappa_{rad} = 3.333 \times 10^{-1}$. In Fig. 12, the result obtained with this fitted model is denoted by diffusion limit with fitted conductivities. This choice of $\kappa_{rad}^{fitted} = 3.333 \times 10^{-3}$ is motivated by the fact that with such value, the material temperature $T_m$ obtained by model (1) is very close to the one obtained with model (40). In fact, $T_m^{(1)}(x, t_f) \neq T_m^{(40)}$, fitted(x, t_f) only in volume $x \in [0.48, 0.58]$, i.e. for about 10% of the simulation domain. In other words, the temperature profiles are in agreement on 90% of the domain.

Surprisingly, model (40) with $\kappa_{rad}^{fitted}$ captures the vicinities of the propagation waves $x \sim 0.70-0.80$, and coincides with $T_r$ for $x \in [0.45, 0.6]$. But it under-estimates the temperature in the vicinity of $x \sim 0.40-0.50$ which is driven by radiative transfer in a semitransparent medium. Of course, other criterion could be considered to fit $\kappa_{rad}^{fitted}$. But the study of this section mainly aims at showing that our model can also be used in semitransparent media. In this semitransparent configuration, it is the reference solution. The study also shows that fitting a diffusion based model in semitransparent medium may be risky or may need some compromise on the accuracy of the predicted temperature and justifies the efforts made by the community to be able to solve (1) with accuracy, independently of the resolution scheme.

## 7. Conclusion

In this paper, we tackle the problem of solving the conducto-radiative equation with a Monte-Carlo scheme (MC) in some layered

Brownian scheme of section 5

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

media. From the observation that the resolution of the transport equation (valid in transparent media) in layered media with MC schemes is easy and the fact that the equilibrium diffusion limit is well captured by state-of-the-art MC schemes, we suggest modeling the sensible heat diffusion with a transport equation in a diffusive regime. The model is close to what is commonly called a phonon model. In order to capture the diffusion limit of material temperature (for the phonons), we suggest some additional model parameters for the velocities and opacities of phonons in the different layers (defined with respect to the conductivities, densities and heat capacities in the different layers). Those parameters are designed to both capture the diffusion limit and to preserve some physical invariants (continuity of temperature at the interface between materials) and ensure convergence. The model is efficient to restitute analytical and reference solutions in the equilibrium diffusion limit and can be safely used coupled with radiative transfer in transparent materials. From the analysis of the phonon model and by performing some analogies between Poissonian paths (which do not suffer jumps in material properties) and Brownian ones, we even suggest some original, quite simple and efficient corrections to the Brownian paths in layered media. In other words, two very different numerical strategies (Poissonian and Brownian) are suggested in this paper and both give very good results. The paper presents many numerical examples testifying of the relevance of the suggested numerical strategies in (transparent and opaque) layered media.

### CRediT authorship contribution statement

**Gaël Poëtte:** Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Methodology, Investigation. **Aujustin De La Vauvre:** Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Investigation, Conceptualization. **Gérard Vignoles:** Writing – review & editing, Writing – original draft, Investigation, 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

The authors would like to thank Bruno Dubroca for valuable discussions. Gael Poette would like to thank Cédric Eaux for old discussions about the diffusion limit which certainly inspired some of the ideas used in this paper. Gael Poette would also like to thank Jérôme Breil and Cédrine Barandon for valuable discussions during Cédrine’s internship.
# Data availability

No data was used for the research described in the article.

# References

[1] T. Fend, B. Hoffschmidt, R. Pitz-Paal, O. Reutter, P. Rietbrock, Porous materials as open volumetric solar receivers: experimental determination of thermophysical and heat transfer properties, Energy 29 (2004) 823-833, http://dx.doi.org/10.1016/S0360-5442(03)00188-9.

[2] A. Kribus, Y. Gray, M. Grijnevich, G. Mittelman, S. Mey-Cloutier, C. Caliot, The promise and challenge of solar volumetric absorbers, Solar Energy 110 (2014) 463-481, http://dx.doi.org/10.1016/j.solener.2014.09.035.

[3] G. Duffa, Ablative Thermal Protection Systems Modeling, in: AIAA Education Series, 2013, http://dx.doi.org/10.2514/4.101717.

[4] L. Ferrari, M. Barbato, B. Esser, I. Petkov, M. Kuhn, S. Gianella, J. Barcena, C. Jimenez, D. Francesconi, V. Liedtke, A. Ortona, Sandwich structured ceramic matrix composites with periodic cellular ceramic cores: an active cooled thermal protection for space vehicles, Compos. Struct. 154 (2016) 61-68, http://dx.doi.org/10.1016/j.comptstruct.2016.07.043.

[5] G. Pulci, J. Tirilló, F. Marra, F. Fossati, C. Bartuli, T. Valente, Carbon- phenolic ablative materials for re-entry space vehicles: manufacturing and properties, Compos. A 41 (10) (2010) 1483-1490, http://dx.doi.org/10.1016/j.compositesita.2010.06.010.

[6] B. Tang, Y. Wang, L. Hu, L. Lin, C. Ma, C. Zhang, Y. Lu, K. Sun, X. Wu, Preparation and properties of lightweight carbon/carbon fiber composite thermal field insulation materials for high-temperature fur- nace, J. Eng. Fibers Fabr. 14 (2019) http://dx.doi.org/10.1177/1558925019884691.

[7] J. Randrianalisoa, Y. Bréchet, D. Baillis, Materials selection for optimal design of a porous radiant burner for environmentally driven requirements, Adv. Eng. Mater. 11 (2009) 1049-1056, http://dx.doi.org/10.1002/adem.200900089.

[8] S. Gianella, D. Gaia, A. Ortona, High temperature applications of Si-SiC cellular ceramics, Adv. Eng. Mater. 14 (12, SI) (2012) 1074-1081, http://dx.doi.org/10.1002/adem.201200012.

[9] M. Pelaconni, M. Barbato, S. Zavattoni, G. Vignoles, A. Ortona, Thermal design, optimization and additive manufacturing of ceramic regular structures to maximize the radiative heat transfer, Mater. Des. 163 (2019) 107539, http://dx.doi.org/10.1016/j.matdes.2018.107539.

[10] F. Asllanaj, Etude et Analyse Numérique des Transferts de Chaleur Couplés par Rayonnement et Conduction Dans Les Milieux Semi-Transparens : Application aux Milieux Fibreux (Ph.D. thesis), 2001, https://theses.fr/2001NAN10208.

[11] V. Leroy, B. Goyeau, J. Taine, Coupled upscaling approaches for conduction, convection, and radiation in porous media: Theoretical developments, Transp. Porous Media 98 (2) (2013) 323-347, http://dx.doi.org/10.1007/s11242-013-0146-x.

[12] Y. Dauvois, Modélisation du transfert thermique couplé conductif et radiatif au sein de milieux fibreux portés à haute température, (Ph.D. thesis), 2016, https://theses.hal.science/tel-01502882.

[13] M. Ghattassi, J. Roche, F. Asllanaj, M. Boutayeb, Galerkin method for solving combined radiative and conductive heat transfer, Int. J. Therm. Sci. 102 (2016) 122-136, http://dx.doi.org/10.1016/j.ijthermalsci.2015.10.011.

[14] M. Badri, Y. Favennec, P. Jolivet, B. Rousseau, Conductive-radiative heat transfer within SiC-based cellular ceramics at high-temperatures: A discrete-scale finite element analysis, Finite Elem. Anal. Des. 178 (2020) 103410, http://dx.doi.org/10.1016/j.finel.2020.103410.

[15] J. Randrianalisoa, D. Baillis, Thermal conductive and radiative properties of solid foams: Traditional and recent advanced modelling approaches, C. R. Phys. 15 (8-9) (2014) 683-695, http://dx.doi.org/10.1016/j.cphy.2014.09.002.

[16] P. Lea, J. Henry, C. Lorreyte, D. Caron, H. Pron, J. Randrianalisoa, Numerical and experimental investigations of heat and mass transfer properties of open-cell ceramic foams, 2018, pp. 8203-8208, http://dx.doi.org/10.1615/IHTC16.pma.024163.

[17] J. Mora-Monteros, C. Suter, S. Haussener, Effective conductivity of porous ceramics in a radiative environment, Ceram. Int. 46 (3) (2020) 2805-2815, http://dx.doi.org/10.1016/j.ceramint.2019.09.272.

[18] M.A. Mendes, P. Talukdar, S. Ray, D. Trimis, Detailed and simplified models for evaluation of effective thermal conductivity of open-cell porous foams at high temperatures in presence of thermal radiation, Int. J. Heat Mass Transfer 68 (2014) 612-624, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.09.071.

[19] P. Jolivet, M. Badri, Y. Favennec, Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning, J. Comput. Phys. 437 (2021) http://dx.doi.org/10.1016/j.jcp.2021.110313.

[20] A. Kumar, J. Vicente, J.-V. Daurelle, Y. Favennec, B. Rousseau, A Numerical Method Based on Domain Decomposition to Solve Coupled Conduction-Radiation Physics using Parallel Computing Within Large Porous Media, Vol. 2116, 2021, http://dx.doi.org/10.1088/1742-6596/2116/1/012057.

[21] A. Kumar, J. Vicente, J.-V. Daurelle, Y. Favennec, B. Rousseau, A modified zonal method to solve coupled conduction-radiation physics within highly porous large scale digitized cellular porous materials, Heat Mass Transfer/Waerme-Stoffuebertrag. (2023) http://dx.doi.org/10.1007/s00231-023-03341-3.

[22] B. Lapeyre, E. Pardoux, R. Sentis, Méthodes de Monte Carlo pour les équations de transport et de diffusion, in: Mathématiques & Applications, no. 29, Springer-Verlag, ISBN: 9783540633938, 1998.

[23] G. Poëtte, Efficient uncertainty propagation for photonics: combining implicit semi-analog Monte Carlo (ismc) and Monte Carlo generalised polynomial chaos (mc-gpc), J. Comput. Phys. (2021) http://dx.doi.org/10.1016/j.jcp.2021.110807, URL https://www.sciencedirect.com/science/article/pii/S0021999121007026.

[24] G. Poëtte, A gPC-intrusive Monte Carlo scheme for the resolution of the uncertain linear Boltzmann equation, J. Comput. Phys. 385 (2019) 135-162, http://dx.doi.org/10.1016/j.jcp.2019.01.052.

[25] G. Poëtte, E. Brun, Efficient uncertain keff computations with the Monte Carlo resolution of generalised polynomial chaos based reduced models, J. Comput. Phys. 456 (2022) 111007, http://dx.doi.org/10.1016/j.jcp.2022.111007.

[26] M. Tancrez, J. Taine, Direct identification of absorption and scattering coefficients and phase function of a porous medium by a Monte Carlo technique, Int. J. Heat Mass Transfer 47 (2) (2004) 373-383, http://dx.doi.org/10.1016/S0017-9310(03)00146-7.

[27] J. Taine, F. Enguehard, Statistical modelling of radiative transfer within non-Beerian effective phases of macroporous media, Int. J. Therm. Sci. 139 (2019) 61-78, http://dx.doi.org/10.1016/j.ijthermalsci.2019.01.012.

[28] M.E. Muller, Some continuous Monte Carlo methods for the Dirichlet problem, Ann. Math. Stat. 27 (3) (1956) 569-589, http://dx.doi.org/10.1214/aoms/1177728169.

[29] I.C. Kim, S. Torquato, Determination of the effective conductivity of heterogeneous media by Brownian motion simulation, J. Appl. Phys. 68 (8) (1990) 3892-3903, http://dx.doi.org/10.1063/1.346276.

[30] I.C. Kim, S. Torquato, Effective conductivity of suspensions of hard spheres by Brownian motion simulation, J. Appl. Phys. 69 (4) (1991) 2280-2289, http://dx.doi.org/10.1063/1.348708.

[31] G. Vignoles, A. Ortona, Numerical study of effective heat conductivities of foams by coupled conduction and radiation, Int. J. Therm. Sci. 109 (2016) 270-278, http://dx.doi.org/10.1016/j.ijthermalsci.2016.06.013.

[32] J.M. Tregan, J.L. Amestoy, M. Bati, J.-J. Bezian, S. Blanco, L. Brunel, C. Caliot, J. Charon, J.-F. Cornet, C. Coustet, L. d'Alençon, J. Dauchet, S. Dutour, S. Eibner, M. El Hafi, V. Eymet, O. Farges, V. Forest, R. Fournier, M. Galtier, V. Gattepaille, J. Gautrais, Z. He, F. Hourdin, L. Ibarrart, J.-L. Joly, P. Lapeyre, P. Lavieille, M.-H. Lecureux, J. Lluc, M. Miscevic, N. Mourtaday, Y. Nyffeneger-Péré, L. Pelissier, L. Penazzi, B. Piaud, C. Rodrigues-Viguier, G. Roques, M. Roger, T. Saez, G. Terrée, N. Villefranque, T. Vourc'h, D. Yaacoub, Coupling radiative, conductive and convective heat-transfers in a single Monte Carlo algorithm: A general theoretical framework for linear situations, PoS ONE 18 (4) (2023) 1-54, http://dx.doi.org/10.1371/journal.pone.0283681.

[33] L. Seyer, V. Gonneau, F. Enguehard, D. Rochais, Improvements of the Brownian walkers method towards the modeling of conduction-radiation coupling, Int. J. Heat Mass Transfer 211 (2023) http://dx.doi.org/10.1016/j.ijheatmasstransfer.2023.124248.

[34] L. Seyer, F. Enguehard, D. Rochais, Improvement of the modeling of the conduction-radiation flash method by a fully stochastic approach. A numerical validation, J. Phys. Conf. Ser. 2766 (1) (2024) http://dx.doi.org/10.1088/1742-6596/2766/1/012035.

[35] L. Seyer, F. Enguehard, D. Rochais, Deterministic and stochastic approaches for the modeling of conduction-radiation coupling within non-Beerian semi-transparent media, J. Quant. Spectrosc. Radiat. Transfer 325 (2024) http://dx.doi.org/10.1016/j.jqsrt.2024.109086.

[36] D. Mihalas, B.W. Mihalas, Foundations of Radiation Hydrodynamics, Dover Publications, 1999.

[37] V. Gonneau, D. Rochais, F. Enguehard, Modelling of heat transfer within heterogeneous media by Brownian walkers, Int. J. Heat Mass Transfer 184 (2022) 122261, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.122261.

[38] B. Mercier, Application of accretive operators theory to the radiative transfer equations, SIAM J. Math. Anal. 18 (2) (1987) 393-408, http://dx.doi.org/10.1137/0518030.

[39] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (3) (1968) 506-517, http://dx.doi.org/10.1137/0705041.

[40] G. Poëtte, X. Valentin, A new implicit Monte-Carlo scheme for photonics (without teleportation error and without tilts), J. Comput. Phys. 412 (2020) 109405, http://dx.doi.org/10.1016/j.jcp.2020.109405.

[41] J.A. Fleck, J.D. Cummings, An implicit Monte-Carlo scheme for calculating time and frequency dependent nonlinear radiation transport, J. Comput. Phys. 8 (1971) 313-342, http://dx.doi.org/10.1016/0021-9991(71)90015-5.

[42] C. Ahrens, E. Larsen, A semi-analog Monte Carlo method for grey radiative transfer problems, in: Proceedings of the ANS Topical Meeting: International Conference on Mathematical Methods to Nuclear Applications, The American Nuclear Society, ISBN: 9780894486616, 2001.

[43] E. Steinberg, S.I. Heizler, Multi-frequency implicit semi-analog Monte-Carlo (ISM C) radiative transfer solver in two-dimensions (without teleportation), 2022, http://dx.doi.org/10.1016/j.jcp.2021.110806.
[44] D. Dureau, G. Poëtte, Hybrid parallel programming models for AMR neutron Monte Carlo transport, in: Joint International Conference on Supercomputing in Nuclear Applications Monte Carlo, 2013, No. 04202 in Parallelism and HPC, EDP Sciences, 2014, http://dx.doi.org/10.1051/snamic/201404202.

[45] A. Lejay, On the constructions of the skew Brownian motion, Probab. Surv. 3 (2006) 413–466, http://dx.doi.org/10.1214/154957807000000013.

[46] A. Lejay, G. Pichot, Simulating diffusion processes in discontinuous media: Benchmark tests, J. Comput. Phys. 314 (2016) 384–413, http://dx.doi.org/10.1016/j.jcp.2016.03.003.

[47] V. Gonneau, Modelling of Heat Transfer within Heterogeneous Media by Brownian Walkers (Ph.D. thesis), Université Paris-Saclay, 2022, https://theses.hal.science/tel-03099750.

[48] L. Farnell, W.G. Gibson, Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: correction to the Gaussian steplength, J. Comput. Phys. 198 (2004) 65–79, http://dx.doi.org/10.1016/j.jcp.2003.12.019.

[49] L. Farnell, W.G. Gibson, Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: A biased random walk on an asymmetrical lattice, J. Comput. Phys. 208 (2005) 253–265, http://dx.doi.org/10.1016/j.jcp.2005.02.013.

[50] P. Étoré, On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients, Electron. J. Probab. 11 (9) (2006) 249–275, URL https://eudml.org/doc/127427.

[51] E.M. LaBolle, J. Quastel, G.E. Fogg, J. Gravner, Diffusion processes in composite porous media and their umerical integration by random walks: Generalized stochastic ipp differential equations with discontinuous coefficients, Water Resour. Res. 36 (3) (2000) 651–662, http://dx.doi.org/10.1029/1999WR900224.

[52] M. Marseguerra, A. Zoia, Normal and anomalous transport across an interface: Monte Carlo and analytical approach, Ann. Nucl. Energy 33 (17) (2006) 1396–1407, http://dx.doi.org/10.1016/j.anucene.2006.09.012.

[53] J.M. Ramirez, Multi-skewed Brownian motion and diffusion in layered media, Proc. Amer. Math. Soc. 139 (10) (2011) 3739–3752, http://dx.doi.org/10.1090/S0002-9939-2011-10766-4.

[54] H.M.V. Ruiz Barlett, M. Hoyuelos, Comparison between fixed and Gaussian steplength in Monte Carlo simulations for diffusion processes, J. Comput. Phys. 230 (2011) 3719–3726, http://dx.doi.org/10.1016/j.jcp.2011.01.041.

[55] A. Lejay, G. Pichot, L. Lenôtre, Diffusion processes in discontinuous media: numerical algorithms and benchmark tests, in: Workshop Validation Approaches for Multiscale Porous Media Models, Nottingham, United Kingdom, 2018, https://inria.hal.science/hal-01900609.

[56] A. Lejay, G. Pichot, Simulating diffusion processes in discontinuous media: Benchmark tests, J. Comput. Phys. 314 (2016) 384–413, http://dx.doi.org/10.1016/j.jcp.2016.03.003.

[57] E. Baioni, A. Lejay, G. Pichot, G.M. Porta, Random walk modeling of conductive heat transport in discontinuous media, 2023, URL https://inria.hal.science/hal-04166562v2.

[58] G.J.M. Uffink, Random walk method for the simulation of macro dispersion in a stratified aquifer, in: Proceedings of the Hamburg Symposium, IUGG 18th General Assembly, August 1983, Vol. 146, 1985, pp. 103–114.

[59] K. Semra, P. Ackerer, R. Mose, Three dimensional groundwater quality modelling in heterogeneous media, Trans. Ecol. Environ. 2 (1993) URL https://www.witpress.com/Secure/elibrary/papers/WP93/WP93001FU.pdf.

[60] E.M. LaBolle, G.E. Fogg, A.F.B. Tompson, Random-walk simulation of transport in heterogeneous porous media: Local mass-conservation problem and implementation methods, Water Resour. Res. 32 (3) (1996) 583–593, http://dx.doi.org/10.1029/95WR03528.

[61] H. Oukili, R. Ababou, G. Debenest, B. Noetinger, Random walks with negative particles for discontinuous diffusion and porosity, J. Comput. Phys. 396 (2019) 687–701, http://dx.doi.org/10.1016/j.jcp.2019.07.006.

[62] H. Belanger, Development of Numerical Methods for Neutronics in Continuous Media (Ph.D. thesis), Université Paris-Saclay, 2022, https://theses.fr/2022UPASP111.

[63] H. Belanger, D. Mancusi, A. Rouchon, A. Zoia, Variance reduction and noise source sampling techniques for Monte Carlo simulations of neutron noise induced by mechanical vibrations, Nucl. Sci. Eng. 197 (4) (2023) 534–557, http://dx.doi.org/10.1080/00295639.2022.2126719.

[64] E. Steinberg, S.I. Heizler, A new discrete implicit Monte Carlo scheme for simulating radiative transfer problems, Astrophys. J. Suppl. Ser. 258 (1) (2022) 14, http://dx.doi.org/10.3847/1538-4365/ac33a3.

[65] E. Steinberg, S.I. Heizler, Frequency-dependent discrete implicit Monte Carlo scheme for the radiative transfer equation, Nucl. Sci. Eng. 197 (9) (2023) 2343–2355, http://dx.doi.org/10.1080/00295639.2023.2190728.

[66] G. Poëtte, X. Valentin, A. Bernede, Canceling teleportation error in legacy IMC code for photonics (without tilts, with simple minimal modifications), J. Comput. Theor. Transp. 49 (4) (2020) 162–194, http://dx.doi.org/10.1080/23324309.2020.1785893.

[67] T. Jeanneau, Etude et validation d'une approche cinétique couplée pour la modélisation du transfert multimodal et multi-échelle de chaleur en milieu hétérogène, (Ph.D. thesis), Université de Bordeaux, 2022, https://theses.hal.science/tel-03991906.

[68] J.A. Fleck, The Calculation of Nonlinear Radiation Transport by a Monte Carlo Method, Doe Tech. Rep. No. w-7405-eng-48, Lawrence Radiation Laboratory, University of California, 1961, http://dx.doi.org/10.2172/4466389, URL https://www.osti.gov/biblio/4466389.

[69] E.E. Lewis, W.F. Miller, Computational methods of neutron transport, 1984, URL https://www.osti.gov/biblio/5538794.

[70] J. Spanier, E.M. Gelbard, Monte Carlo Principles and Neutron Transport Problems, Addison-Wesley, ISBN: 978-0201070897, 1969.

[71] G.C. Papanicolaou, Asymptotic analysis of transport processes, Bull. Amer. Math. Soc. 81 (2) (1975) 330–392, http://dx.doi.org/10.1090/S0002-9904-1975-13744-X.

[72] S. Bianchini, B. Hanouzet, R. Natalini, Asymptotic behavior of smooth solutions for partially dissipative hyperbolic systems with a convex entropy, Comm. Pure Appl. Math. 60 (11) (2007) 1559–1622, http://dx.doi.org/10.1002/cpa.20195.

[73] J.D. Densmore, E.W. Larsen, Asymptotic equilibrium diffusion analysis of time-dependent Monte Carlo methods for grey radiative transfer, J. Comput. Phys. 199 (1) (2004) 175–201, http://dx.doi.org/10.1016/j.jcp.2004.02.004.

[74] A.G. Irvine, I.D. Boyd, N.A. Gentile, Reducing the spatial discretization error of thermal emission in implicit Monte Carlo simulations, J. Comput. Theor. Transp. 45 (1–2) (2016) 99–122, http://dx.doi.org/10.1080/23324309.2015.1060245.

[75] A.B. Wollaber, Four decades of implicit Monte Carlo, J. Comput. Theor. Transp. 45 (1–2) (2016) 1–70, http://dx.doi.org/10.1080/23324309.2016.1138132.

[76] R.J. Papoular, R. Papoular, Some optical properties of graphite from ir to millimetric wavelengths, Mon. Not. R. Astron. Soc. 443 (2014) 2974–2982, http://dx.doi.org/10.1093/mnras/stu1348.

[77] A. Butland, R. Maddison, The specific heat of graphite: An evaluation of measurements, J. Nucl. Mater. 49 (1973) 45–56, http://dx.doi.org/10.1016/0022-3115(73)90060-3.

[78] J. Tiwari, T. Feng, Accurate prediction of thermal conductivity of al 2 o3 at ultrahigh temperatures, 109 (2024) 075–201, http://dx.doi.org/10.1103/PhysRevB.109.075201.

[79] W.J. Tropf, Handbook of Optical Constants of Solids Volume III - Aluminum Oxide (Al2O3) Revisited, Elsevier, 1997, pp. 653–682, http://dx.doi.org/10.1016/B978-012544415-6.50124-2.

[80] M.F. Modest, Radiative Heat Transfer in Turbulent Combustion Systems Theory and Applications, first ed., McGraw-Hill, New York, ISBN: 978-0070426757, 1993.