Onur Tasar

Yuliya Tarabalka

Pierre Alliez

CONTINUAL LEARNING FOR DENSE LABELING OF SATELLITE IMAGES

Keywords: Continual learning, dense labeling, semantic segmentation, convolutional neural networks, catastrophic forgetting

In dense labeling problem, the major drawback of the convolutional neural networks is their inability to learn new classes without affecting performance for the old classes on the data, having no annotations for the previous classes. In this work, we address the issue of adding new classes continually to the already trained network from a stream of data. Our approach comprises two main components: adaptation and remembering. For adaptation, we keep a clone of the previously trained network, which serves as a memory for the old classes in absence of their annotations on the new data. The updated network learns new as well as old classes on the current data using output of the memory network and the new groundtruth. For remembering, we store a little portion of the previous data, from which we systematically feed samples to the updated network during training. Our results prove that segmentation capabilities for the new classes can be added to the already trained network without catastrophically forgetting the previously learned information.

INTRODUCTION

With the great improvements of satellite sensors, it has been possible to collect massive amount of data, which has opened the door to a wide range of remote sensing applications such as navigation, disaster prevention, and mapping. Over the years, because of the growing interest in such applications, generating maps from satellite images automatically has been one of the most popular topics in remote sensing community. Although it has been shown that convolutional neural networks (CNNs), particularly U-net and its variants [START_REF] Huang | Large-scale semantic classification: outcome of the first year of inria aerial image labeling benchmark[END_REF], can generate high quality maps, most of the works in the literature are limited to learning from only static training data. Since new images are collected from all around the world everyday, one can not assume that all the training images available in the beginning of the training procedure. In addition, annotations obtained from different sources usually have distinct classes. Moreover, it may not always possible to store huge amount of data. Because of these reasons, it is of paramount importance to devise a classification system that can learn new classes continually while retaining performance for the old classes without accessing the entire previous data. The continual learning problem is illustrated in Fig. 1. To the best of our knowledge, in remote sensing community, this problem has not been explored yet.

The main drawback of the common CNNs for the continual learning problem is that their performance for the old classes significantly drops when new classes are added to the trained model on the data, having no ground-truth for the previous classes. The aforementioned abrupt performance loss is referred as catastrophic forgetting in the literature [START_REF] Goodfellow | An empirical investigation of catastrophic forgetting in gradient-based neural networks[END_REF]. To overcome this problem, there have been some attempts, which change the network structure. An example methodology that is in this category is described in [START_REF] Sarwar | Incremental learning in deep convolutional neural networks using partial network sharing[END_REF], where the network is expanded as new classes are added. The obvious limitation of this approach is that the network grows every time new classes are taught. Since deep neural networks contain millions of parameters, different configurations of the parameters may produce the same results. Inspired from this idea, several works, which try to determine the important neurons and prevent them from changing drastically, have been proposed [START_REF] Kirkpatrick | Overcoming catastrophic forgetting in neural networks[END_REF][START_REF] Aljundi | Memory aware synapses: Learning what (not) to forget[END_REF]. The methodology proposed in [START_REF] Rebuffi | iCaRL: Incremental classifier and representation learning[END_REF] uses a small fraction of the previous training data to prevent forgetting the already learned information. The knowledge distillation from one network to another one [START_REF] Hinton | Distilling the knowledge in a neural network[END_REF] gave inspiration to several works. In [START_REF] Li | Learning without forgetting[END_REF][START_REF] Shmelkov | Incremental learning of object detectors without catastrophic forgetting[END_REF], the previous knowledge is distilled from the formerly trained network to the current one on the new data. Even though several methodologies have been proposed to overcome catastrophic forgetting, none of those works investigates continual learning for dense labeling problem.

In this paper, we propose a method, which enables to learn dense labeling capabilities for the new classes without forgetting the old classes even if the whole previous training data are not available.

METHODOLOGY

Our network structure is U-net, comprising an encoder that is architecturally the same as encoder of VGG-16 [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF] and a symmetric decoder. Each convolution and deconvolution operation is followed by a rectified liner unit (ReLU). Since batch normalization uses the memory inefficiently, we prefer not to use it. The output of our network consists of binary predicted maps for all the classes.

Let us assume that the current training data are denoted by D curr . The sets of previously learned classes and the classes in D curr are indicated by L prev and L curr , where L prev ∩ L curr = ∅. During adaptation, we teach L curr and L prev to the network on D curr . Since the annotations for L prev are not available in D curr , we transfer the knowledge from the memory network to the updated network. We denote by y curr , the binary target vectors of the pixels in a training batch from D curr for all the classes. We denote the binary predicted probabilities for L curr and L prev from the updated network by ŷup curr and ŷup prev , respectively. The predictions for L prev from the memory network are indicated by ŷmem . To learn L curr , we minimize the classification loss L class sigmCE(ŷup curr ,ycurr) , which is the sigmoid cross entropy loss between ŷup curr and y curr . In order to distill the knowledge for L prev from the memory network to the updated network, we minimize the distillation loss L distil sigmCE(ŷmem,ŷup prev ) that is the sigmoid cross entropy loss between ŷmem and ŷup prev . The loss that is minimized during adaptation is defined as:

L adapt = L class sigmCE(ŷup curr ,ycurr) +L distil sigmCE(ŷmem,ŷup prev ) .
(1) For remembering, we store a little fraction of the previous training data. Because in satellite images number of samples for the classes are usually imbalanced, randomly selecting the training patches that would be stored might result in having no samples for the less frequent classes. Hence, we handle the class imbalance problem. Let us assume that D

w c = median(f c |c ∈ L (j) prev ) f c , (2)
where f c is frequency of the pixels, belonging to class c. We then compute an importance value I (l) for the l th training prev as:

(l) = c∈L (j) prev w c f (l) c , (3)
where

f (l)
c denotes the number of pixels that are labeled as class c in the patch. We store a certain number of patches, having the highest I values. In order to prevent storing only the patches that are located in very close geographic locations, we also store certain number of randomly selected patches. The number of patches that are selected randomly and using I values have to be determined. Let us assume that y (j) prev is the target vector, and ŷ(j) up prev denotes the predictions of the pixels for L (j) prev in the training batch from the stored, small portion of D (j) prev . The remembering loss L rem is defined as:

L rem = L sigmCE (y (j) prev , ŷ(j) up prev ), (4)
where L sigmCE (y

(j)
prev , ŷ(j) up prev ) is the sigmoid cross entropy loss between y (j) prev and ŷ(j) up prev . The end user needs to determine how frequently and on which previous training data, the samples would be fed to the network to optimize L rem . An example optimization sequence is depicted in Fig. 2.

METHODS USED FOR COMPARISON

We compare our methodology with the following approaches:

Multiple learning : In this approach, we train a new network from scratch every time new training data are retrieved. Because of the growing number of classifiers, this learning approach is inefficient in terms of storage. Besides, since the test images have to be segmented by each trained model to generate label maps for all the classes, the running time in the test stage might be very long.

Fixed representation : Whenever new training data are obtained, we increase the number of filters in the last convolutional layer that is responsible for generating a binary predicted map for each class. In the training stage, we optimize only the filters for the new classes, and freeze the rest of the network. As we do not update the filters for the previous classes, the exact performance for them is retained. However, the network is not able to learn the new classes well, as the extracted features are not representative.

Fine-tuning : This approach is similar to fixed representation. The difference is that we freeze only the filters for the previous classes in the last convolutional layer, and optimize the rest of the network. This methodology suffers from catastrophic forgetting.

Continual learning w/o rem.: This is our approach without optimizing L rem . Since the whole network is adapted to the last training data, the information learned from the previous data might be forgotten.

EXPERIMENTS

In our experiments, we use the Luxcarta dataset, which contains satellite images collected from 11 cities in France and 11 cities in Austria. The spatial resolution is 1 m. The full annotations for building, tree, road, railway, and water classes are available. We use 9 cities from each country as the training data, and the rest of the cities as test. We split the training cities into three groups, as reported in Table 1. We assume that the cities in each group have annotations for only the classes that are indicated in the table. When dividing the training cities into groups, we make sure that cities in each group have acceptable number of samples for the indicated classes. We suppose that the training data are streamed in this order: Train1, Train2, and Train3. For our approach, we store 30% of the previous training patches, in which 15% are selected randomly and the other 15% are chosen using the I values. For the approaches described in Sec. 3, we do not use the previous training data.

In the pre-processing stage, we split all the training images into 384×384 patches, having 32×32 pixels of overlap between the neighboring patches. We divide the test images into 2240×2240 patches with an overlap of 64×64 pixels. After classifying each test patch, we merge them back to get the final predicted maps. For the normalization, we subtract 127 from the pixels of each spectral band.

We train 3 models for multiple learning (one from each group) for 500 epochs, where each epoch consists of 100 iterations. For the rest of the approaches, whenever new training data are retrieved, we update the previous model by training for the same number of epochs and iterations as for multiple learning. For continual learning, when we add road and railway classes on Train2 to the model trained on Train1, we op- Roanne (France) water timize L adapt for 4 iterations and L rem for 1 iteration. To add water class on Train3, the optimization sequence is: L adapt for 2 iterations, L rem on Train1 for 1 iteration, L adapt for 2 iterations, and L rem on Train2 for 1 iteration. We use Adam optimizer to optimize parameters of the networks. The learning rate is 0.0001 and exponential decay rate for the first and the second moment estimates are 0.9 and 0.999. Each training batch consists of 12 patches. When sampling a patch, we first randomly select a country. We then choose a random city from the selected country. We finally randomly pick a training patch from the chosen city. In the training stage, we augment the training patches with 0/90/180/270 degrees of random rotation, random horizontal and vertical flips, random contrast change, and gamma correction. F1-scores [START_REF] Audebert | Semantic segmentation of earth observation data using multimodal and multi-scale deep networks[END_REF] of each method for all the classes are reported in Table 2. Close-ups from the generated maps for multiple learning, continual learning w/o L rem , and continual learning are shown in Fig. 3. Although the output of our network architecture is a binary classification map for each class, to visualize the outputs more compactly, we provide multi-class maps in the figure by assigning each pixel to the class, having the highest probability. As expected, fixed representation can not learn new classes very well, and fine-tuning catastrophically forgets the previously learned classes. For our approach, if we do not remind information from the previous data, performance for the already learned classes might decrease over time. For instance, for tree class, performance of continual learning w/o L rem is lower than multiple learning and continual learning. For multiple learning, learning for each class is limited to the images from only one group, which may cause showing a lower performance. E.g., building class is learned from only the images in Train1. Therefore, performance for this class is lower than our approach. As the proposed continual learning approach enables to both learn from the new data and remember from the previous data, it exhibits the best performance for most of the classes.

CONCLUDING REMARKS

We proposed a novel continual learning approach, enabling to learn segmenting new classes without deteriorating performance for the old classes, although a small portion of the previous training data are accessible. As the future work, we plan to incorporate domain adaptation techniques into the proposed continual learning approach so that new classes could be learned successfully, even when there is a large domain shift between the previous and the current training images.

Fig. 1 .

 1 Fig. 1. Continual learning problem. Every time new training data are retrieved, only a small potion of the previous ones are stored. The classes are building (red), tree (green), water (blue), road (white), and railway (yellow).

  prev is the j th previous training data, and L (j) prev denotes the classes in D (j) prev . We compute a weight w c for each class c ∈ L (j) prev in D (j) prev as:

Fig. 2 .

 2 Fig. 2. Example optimization sequence. The sequence is: L rem on data 1, L adapt on the current data, L rem on data 2, and L adapt on the current data again.

Fig. 3 .

 3 Fig. 3. Close-ups from the test cities. Classes: background (black), building (red), road (white), railway (yellow), high vegetation (green), and water (blue).

Table 1 .

 1 Training and test cities.

	Data Type	Class	City (Country)
			Bad Ischl (Austria)
		building	Osttirol (Austria) Voitsberg (Austria)
	Train1		Bayonne-Biarritz (France)
			Bourges (France)
		tree	Draguignan (France)
			Nîmes (France)
			Enns (Austria)
		road	Innsbruck (Austria)
	Train2		Klagenfurt (Austria) Sankt Pölten (Austria)
		railway	Béziers (France)
			Lyon (France)
			Albi (France)
			Villach (Austria)
	Train3	water	Salzburg (Austria)
			Angers (France)
			Douai (France)
	Test	building tree road railway	Amstetten (Austria) Leibnitz (Austria) Lille (France)

Table 2 .

 2 F1 scores on the Luxcarta dataset.

		Method		Epoch Building Tree Road Railway	Water	Overall
		multiple learning		500	71.25	68.88 59.28	55.65	79.83	66.98
		fixed representation	1000 1500	71.25 71.25	68.88 68.88	2.71 2.71	0.00 0.00	-0.11	28.59
		fine-tuning	1000 1500	28.91 27.90	0.17 7.71	59.30 0.14	60.06 0.01	-90.20	25.19
		continual learning	1000	74.19	66.32 56.57	50.87	-
		w/o L rem		1500	74.91	66.87 58.14	51.70	82.32	66.79
		continual learning	1000 1500	75.98 76.78	72.38 57.29 72.06 59.58	50.18 53.07	-78.94	68.09
						Training Set 1	Training Set 2	Training Set 3
	Image	GT	mult. lear.	cont. lear. cont. lear. w/o Lrem
	Amstetten
	Leibnitz
	Lille
	Roanne


Gaël Poëtte
email: gael.poette@cea.fr

Xavier Valentin
email: xavier.valentin@cea.fr

A new Implicit Monte-Carlo scheme for photonics (without teleportation error and without tilts)

Keywords: Transport, Monte-Carlo, Numerical scheme, Photonics, IMC, SMC

In this paper, we present a new implicit Monte-Carlo scheme for photonics. The new solver combines the benefits of both the IMC solver of Fleck & Cummings and the SMC solver of Ahrens & Larsen. It is implicit hence allows taking affordable time steps (as IMC) and has no teleportation error (as SMC). The paper also provides some original analysis of existing schemes (IMC, tilted IMC, SMC), especially with respect to the teleportation error in the equilibrium diffusion regime. In particular, we demonstrate that any small spatial inaccuracies during the sampling of source particles for IMC lead to a competing behaviour between the spatial and time discretisation parameters. The new scheme we suggest is implicit, conservative, has no teleportation error (and as a consequence does not need tilting), does not rely on source sampling for the emission of source particles, captures the equilibrium diffusion limit (provided a small enough time step), can be used with arbitrary equations of state and does not suffer the above competing behaviour. All those properties are either demonstrated or numerically highlighted in the paper.

Introduction

In this article, we are interested in the Monte-Carlo (MC) resolution of the time-dependent, nonlinear, radiative transfer equations. The model has general form (see [START_REF] Mihalas | Foundations of Radiation Hydrodynamics[END_REF]):

       1 c ∂ t I + ω • ∇I + σ t I = σ a B(T m ) + σ s 4π I dω 4π , ∂ t E(T m ) = 4π cσ a I 4π -B(T m ) dω . (1)
In the above equations, I = I(t, x, ω) and T m (t, x) are the unknowns of the system and stand respectively for the density of radiation energy and the material temperature. Variables t ≥ 0, x ∈ Ω ⊂ R 3 and ω ∈ S 2 are respectively the time, space and angle variables. The cross-sections σ t = σ t (x, t), σ a = σ a (x, t) and σ s = σ s (x, t) are given functions of (x, t). They stand for the total, absorption and scattering opacities. In particular, we have σ t = σ a + σ s . The density of internal energy E depends on T m via an equation of state dE = ρC v (T m )dT m with ρ the mass density and C v the heat capacity (constant for a perfect gaz). Quantity B(x, t) = aT 4 m (x, t)/4π is the frequency-integrated Planck function with a the radiative constant. Quantity c denotes the speed of light. Initial and boundary conditions must be supplemented to system [START_REF] Mihalas | Foundations of Radiation Hydrodynamics[END_REF]:

u(0, x, ω) = u 0 (x, ω), T m (0, x) = T 0 m (x), x ∈ Ω, ω ∈ S 2 , (2)
u(t, x, ω) = u b (t, ω), t ≥ 0, x ∈ ∂Ω, ω • n s < 0, (3)
where n s is the outward normal to Ω at x. System (1) together with initial and boundary conditions (2) + (3) define the well-posed mathematical problem we want to solve. To be more precise, in this paper, we are even interested in being able to accurately capture a particular regime: in diffusive media, system (1) behaves, at leading order, like the nonlinear diffusion equation on Φ r (T r ) = aT 4 r = 4π

I 4π dω ∂ t (Φ r (T r ) + E(T r )) -∇ • c 3σt ∇(Φ(T r ) = O(δ), Φ r (T r ) = 4π I 4π dω = B(T m ) + O(δ). (4)
With 4π B(Tm) 4π dω = aT 4 m and Φ r (T r ) = aT 4 r , the second equation is equivalent to T m = T r : the radiative and matter temperatures are at equilibrium. In the above equation, δ ∼ 0 is a small parameter characterising what is commonly called the equilibrium 1 diffusion 2 limit [START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF][START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF][START_REF] Wollaber | Four decades of implicit monte carlo[END_REF]. The limit can be defined by introducing a characteristic length X , a characteristic time T and a characteristic collision rate λ and assuming we have

c T D = O( 1 δ ), cσ T λ = O( 1 δ 2 ), (5)
with δ ∼ 0 small. System (1) and its limit ( 4) are relevant to model photons incoming into cold media [START_REF] Mihalas | Foundations of Radiation Hydrodynamics[END_REF][START_REF] Castor | Radiation hydrodynamics[END_REF][START_REF] Pomraning | The equations of radiation hydrodynamics[END_REF]. At this stage of the discussion, the problem of interest in this paper may seem to lack generality in the sense we here focus on what is commonly called the grey approximation 3 with isotropic scattering 4 and the opacities are independent of the temperature 5 . Still, it is enough to focus on the main numerical difficulties we aim at considering in this paper: a high-dimensional (6 dimensions in 3D) and highly nonlinear (via B ∝ T 4 m ) system (1) degenerating toward (4) in diffusive media. Of course, care has been taken to make sure the new numerical method we present in this paper can be generalised.

In this article, due to the high-dimensional problem we aim at tackling, we only focus on MC based numerical solvers. MC schemes to solve (1) are known to be extremely costly and inaccurate in diffusive media, i.e where absorptions/emissions of photons dominate transport [START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF][START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF][START_REF] Wollaber | Four decades of implicit monte carlo[END_REF]. The aim of this paper is to provide a new and original analysis of some existing schemes, identify some flaws and suggest a new way, i.e. a new MC scheme, to circumvent them.

The paper is organised as follows: in section 2, we describe one of the most common way to solve (1) together with its limit (4), the IMC scheme [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF]. The scheme is described here because one important ingredient is used later on in our new approach. In particular in this section, we identify a simple configuration and illustrate what is commonly called the teleportation error, see [START_REF] Wollaber | Four decades of implicit monte carlo[END_REF][START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF]. Section 3 is devoted to an original analysis of the aforementioned teleportation error and of the tilts, commonly used to mitigate it [START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF][START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF][START_REF] Wollaber | Four decades of implicit monte carlo[END_REF]. In particular, we demonstrate that small inaccuracies on the emission term (during the source sampling phase) can lead to the construction of a solver with competing numerical parameters (namely the time step ∆t and the grid size ∆x or any spatial discrepancy δ x during the source sampling phase). We also put forward the fact that even if mitigating the teleportation error, tilts can not change the competing behaviour of ∆t and ∆x. Finally, in section 4, we suggest a new implicit MC scheme which can not, by construction, produce any teleportation error (hence it does not need the use of a tilt). We call it the ISMC scheme because it is based on both the SMC scheme of [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF] and an adaptation of the implicitation at the basis of the IMC scheme of [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF]. Section 5 is finally devoted to numerical benchmarks.

The Implicit Monte Carlo (IMC) method

We first describe one of the most common methodology to solve system [START_REF] Mihalas | Foundations of Radiation Hydrodynamics[END_REF]. It is denoted by IMC for Implicit Monte-Carlo and has been introduced by Fleck and Cummings in [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF]. Even though the method has been described many times by multiple authors [START_REF]The calculation of nonlinear radiation transport by a monte carlo method[END_REF][START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF][START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF][START_REF] Wollaber | On the stability of the ahrens-larsen or smc equations for thermal radiative transfer[END_REF], we recall the main ideas introduced in the original paper. These ideas play an important role in our new MC scheme we propose in section 4. The IMC method takes its roots in the statement that the standard explicit scheme developed in [START_REF]The calculation of nonlinear radiation transport by a monte carlo method[END_REF] by Fleck exhibited objectionable features (small time steps are needed), especially in collisional media. The idea of Fleck and Cummings in [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF] is to introduce some sort of implicitness in the time discretization scheme in order to unstiff the coupling between radiation and material. Implicit time discretization of (1) yields the following linear system:

       1 c ∂ t I + ω • ∇I + σ n t I = σ n a B n+1 + σ n s 4π I dω 4π , ∂ t E = 4π cσ n a I 4π -B n+1 dω . (6)
To solve system (6) using an MC scheme, one has to sample the source term σ n a B n+1 which is unknown at the beginning of the time step. In [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF], the authors propose a particular estimation of B n+1 that leads to the Fleck and Cummings equations that approximate the original system [START_REF] Mihalas | Foundations of Radiation Hydrodynamics[END_REF]. We sum up the construction of B n+1 in the few next lines for the grey approximation. System (1) can be rewritten with respect to variables (I, Φ = aT 4 m ) instead of variables (I, T m ):

       1 c ∂ t I + ω • ∇I + σ t I = σ a Φ + σ s 4π I dω 4π , ∂ t Φ = 4π cσ a β(Φ) I 4π dω -cσ a β(Φ)Φ, with β = dΦ dE .
Note that β = dΦ dE has been introduced to express ∂ t E with respect to ∂ t Φ. In a sense, this term deals with all the physics relative to the equation of state E(T m ) = E(T m (Φ)) together with making sure the system can be closed using variables (I, Φ).

The above system is still nonlinear. A linearisation is mandatory to apply an MC scheme. Integration of the radiation-material energy balance equation on time step [t n , t n+1 ], using the backward Euler scheme for Φ and the forward Euler scheme for σ a and β yields:

Φ n+1 = Φ n + t n+1 t n 4π cσ n a β n I 4π dω -cσ n a β n ∆tΦ n+1 .
It can be easily inverted in term of Φ n+1 to give:

Φ n+1 = f n Φ n + (1 -f n ) t n+1 t n 4π I 4π dω dt ∆t , with f n = 1 1 + cσ n a β n ∆t .
The quantity f n is the so-called Fleck factor. Replacing the time integral by its instantaneous value and using the approximation of Φ n+1 in the transport equation yields the time-discretized Fleck and Cummings transport equation :

       1 c ∂ t I + ω • ∇I + σ n t I = f n σ n a Φ n + [σ n s + (1 -f n )σ n a ] 4π I dω 4π , ∂ t E = cf n σ n a 4π I 4π dω -Φ n . (7)
Notice that by redefining σ n a ← f n σ n a and σ n s ← σ n s + (1f n )σ n a , one obtains a transport equation similar to the first one of (6) except that the emission term f n σ n a Φ n is now known because evaluated at the beginning of the time step. The time-discretized transport equation above has been supplemented by the proper material energy equation. Choosing to work on E instead of Φ here ensures, by construction, the conservation of the total energy for the system matter+photons provided Φ n+1 can be deduced from E n+1 (eos).

Equations ( 7) are the Fleck and Cummings equations and form a closed linear system on time step [t n , t n+1 ]. This system is a consistant time-discretized version of the original system [START_REF] Mihalas | Foundations of Radiation Hydrodynamics[END_REF]. Remark that the semi-implicit time discretization introduces an artificial scattering term that unstretches the stiff coupling between radiation and material. This numerical trick is important and efficient in practice and will be used in the derivation of our new MC scheme (section 4). Fleck and Cummings' solver for (1) then consists in solving successively the two equations of system [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF] within every time step. Let us introduce a spatial grid of N x non-overlapping cells such that Ω = Nx i=1 Ω i . Then, during time step [t n , t n+1 = t n + ∆t],

-transport equation ( 7) is solved with an MC method. We denote by N M C the number of MC particles. The source term is classically treated with a source sampling strategy within each cell (Ω i ) i∈{1,...,Nx} , see [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF][START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF][START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF][START_REF] Wollaber | Four decades of implicit monte carlo[END_REF].

-The material energy E n+1 is updated in each cell (Ω i ) i∈{1,...,Nx} by tallying the radiation energy deposit in the material with a track length estimator [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF] during the MC resolution. Now, we aim at analysing the asymptotic behavior of the time-discretized system [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF] in the diffusive regime. This analysis has been carried out in [START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF]. One can show that, during time step [t n , t n+1 ], at leading order with respect to δ in the diffusive regime characterised by [START_REF] Castor | Radiation hydrodynamics[END_REF], system (7) behaves like:

         ∂ t Φ r -∇ • c 3σ n t ∇Φ r = 1 β n ∆t (Φ r -Φ n ) + O(δ), E n+1 = E n + 1 β n ∆t t n+1 t n Φ r dt -Φ n + O(δ). (8)
Notice that, from the comparison of ( 4) and [START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF], one can expect the solution to exhibit isotropy for I (given enough MC particles) and recover the diffusion coefficient c/3σ n t (provided a fine enough spatial discretization). On another hand, the equilibrium property Φ r = Φ of ( 4) is only Φ r -Φ = O(β n ∆t) once the IMC linearisation applied, see [START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF]. Let us now present some numerical results on a 'fil rouge' problem and illustrate why we felt the need to build a new MC scheme.

Let us here consider a particular configuration: note that every details (initial conditions, numerical parameter choices and test-case justifications) of the test-problem are presented in Appendix A. We consider a diffusive medium for which (4) is valid. The initial condition is a Heaviside of (relaxed T m = T r ) temperatures in the center of the spatial domain, see figure A.10. The solution at the final time exhibits steep propagation fronts. The (reference) solution obtained by solving (4) is recalled in every pictures of figure 1 together with the radiative T r and material temperatures T m obtained with the IMC solver for -always the same time step ∆t = 10 -12 , -about N M C ≈ 4 × 10 5 MC particles (analog MC scheme + source sampling), -and several mesh sizes as N x ∈ {20, 40, 80, 160, 320, 640, 1280, 2560}.

In other words, figure 1 presents a convergence study on the spatial temperature profiles (T r , T m ) with respect to the spatial discretization N x for fixed ∆t and N M C obtained from an IMC solver.

The top-left picture of figure 1 shows an important difference between the reference solution and the IMC one. On a coarse mesh (N x = 20), we experience what is commonly called the teleportation error [START_REF] Wollaber | Four decades of implicit monte carlo[END_REF][START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF][START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF]: the IMC solution does not behave as a diffusion solution (care will be taken to characterise its behaviour more precisely in section 3). Then, as testifies figure 1, the finer the mesh (as N x grows), the better the equilibrium diffusion solution is captured. But keeping N M C fixed as N x grows also leads to noisier and noisier results. There are less and less particles per cell as N x grows. The number of MC particles N M C must grow with N x to ensure recovering a smooth/accurate solution. Note that the decoupling between radiative and material temperatures, driven by the numerical relaxation time 6 1 β n ∆t , is satisfactory: the authors in [START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF] emphasized the fact that in practice, time steps ensuring a good agreement for this equilibrium limit are affordable for the IMC scheme (meaning we can afford β n ∆t small enough in common applications).

With this simple example, we wanted to express the stakes of being able to deal with the teleportation error. Being able to capture the equilibrium diffusion limit with an MC scheme on a coarse mesh will lead to considerable gains especially because smooth (i.e. not noisy i.e. accurate) solutions will be available with even less MC particles. For this reason in the following section, we focus on the spatial discretization. More precisely, we focus on the spatial discretization of the source term which is critical in the diffusive regime, cause of what is commonly called the teleportation error [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF][START_REF] Mckinley | Comparison of implicit and symbolic implicit monte carlo line transport with frequency weight vector extension[END_REF][START_REF] Wollaber | Four decades of implicit monte carlo[END_REF][START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF]. The next section is devoted to a deepened and original numerical analysis of this teleportation error. Care is taken to characterise precisely its effect on the solution of [START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF].

N x = 20 N x = 40 T m IMC nx=20 T r IMC nx=20 T m = T r ,

Analysis of the teleportation error

The analysis of the limit equation on time step [t n , t n+1 = t n + ∆t] for the IMC scheme has been carried out continuously with respect to the spatial variable whereas many authors (see [START_REF] Wollaber | Four decades of implicit monte carlo[END_REF] and the references therein) admit the spatial discretisation of the source term in the IMC linearisation is closely related to the teleportation error. See [START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF][START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF][START_REF] Wollaber | Four decades of implicit monte carlo[END_REF] for very pedagogical numerical examples. It is even confirmed by the fact that tilts, i.e. spatial interpolary reconstructions of Φ, can lead to significant improvements, see [START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF][START_REF] Wollaber | Four decades of implicit monte carlo[END_REF]. In practice, a spatial discretisation is introduced (i.e. D = Nx i=1 D i ), together with assumptions such as constant Φ in each cell. At every beginning of time step we have an O(∆x)

approximation of Φ(x) = Nx i=1 Φ i 1 Di (x) + O(∆x) where ∆x = max i |D i | and Φ i = 1
|Di| Di Φ(x)dx. This choice affects Φ(x) but also β and above all the source sampling of the MC particles at every beginning of time steps. In the following paragraphs, we suggest -taking into account a spatial O(δ x ) discrepancy in the source term within cells (inaccurate source sampling),

-together with the asymptotic development O(δ) in the equilibrium diffusion limit defined by (5) of the linearized system (1).

Let us now revisit the limit equation under such condition. Let us perform a Taylor development of Φ(x) with respect to a small spatial parameter δ x . Note that we are also going to perform a Hilbert [START_REF] Hilbert | Begrundung der kinetischen gastheorie[END_REF][START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF] development afterward with respect to δ. For this, we write:

Φ(x) = Φ 0 (x) + δ x Φ 1 (x) + O(δ 2 x ), with of course Φ 0 (x) = Nx i=1 Φ 0 i 1 Di (x) and Φ 1 (x) = ∇ x Φ.
The upperscripts refer to the terms in the development with respect to δ x . Assume constant opacities, then for small δ x , we have:

β(Φ) = β(Φ 0 ) + δ x Φ 1 dβ dΦ (Φ 0 ) β (Φ 0 ) + O(δ 2 x ).
Plugging the spatial Taylor development of Φ within four terms appearing in (7) (we drop the n upperscripts for conciseness), namely cσ a f , cσ a (1f ), cσ a f Φ, leads to

cσ a f Φ = cσ a Φ 0 1 + cσ a β(Φ 0 )∆t + δ x cσ a Φ 1 1 + cσ a [β(Φ 0 ) -Φ 0 β (Φ 0 )]∆t (1 + cσ a β(Φ 0 )∆t) 2 + O(δ 2 x ), cσ a (1 -f ) = cσ a cσ a β(Φ 0 )∆t 1 + cσ a β(Φ 0 )∆t + δ x cσ a Φ 1 cσ a β (Φ 0 )∆t (1 + cσ a β(Φ 0 )∆t) 2 + O(δ 2 x ), cσ a f = cσ a 1 + cσ a β(Φ 0 )∆t -δ x cσ a Φ 1 cσ a β (Φ 0 )∆t (1 + cσ a β(Φ 0 )∆t) 2 + O(δ 2 x ). Now use the fact that c T D = O( c * δ ), cσ a T λ = O( c * σ * a δ 2
) as in [START_REF] Castor | Radiation hydrodynamics[END_REF]. Let us perform a Hilbert development [START_REF] Hilbert | Begrundung der kinetischen gastheorie[END_REF][START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF] of the above system with respect to δ. It consists in -expanding

I = I 0 + I 1 δ + I 2 δ 2 + O(δ 3 ), Φ = Φ 0 + Φ 1 δ + Φ 2 δ 2 + O(δ 3 ) and E = E 0 + E 1 δ + E 2 δ 2 + O(δ 3
) as formal power series of δ, -plugging their expressions in [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF],

-and identify the coefficients of δ 0 , δ 1 , δ 2 .

With this additional development, the above equations become (note that we drop the upperscript * for conciseness)

δ 2 cσ a f Φ = Φ 0 0 β(Φ 0 0 )∆t + Φ 0 1 δ x β(Φ 0 0 )∆t -Φ 0 1 δ x β (Φ 0 0 )Φ 0 0 β(Φ 0 0 ) 2 ∆t + O(δ 2 x ) δ 2 + O(δ 3 ), δ 2 (cσ s + cσ a (1 -f )) = cσ t + - 1 β(Φ 0 0 )∆t + Φ 0 1 δ x β (Φ 0 0 ) (β 2 (Φ 0 0 )∆t) + O(δ 2 x ) δ 2 + O(δ 3 ), cσ a f = 1 β(Φ 0 0 )∆t -Φ 0 1 δ x β (Φ 0 0 ) ∆tβ 2 (Φ 0 0 ) + O(δ 2 x ) + O(δ).
By anticipation, the two first term have been multiplied by δ 2 (to be used in the first equation of ( 7)) and not the last one (to be used in the second equation of ( 7)).

The above developments are then plugged in the collisional part of (7) solved on time step [t n , t n + ∆t] assuming δ → 0 and keeping only the O(δ 4 ) orders. We obtain:

                 δ 2 ∂ t I + δcω • ∇I + cσ t I = δ 2 β(Φ 0 0 )∆t Φ 0 0 + δ 2 [β(Φ 0 0 ) -Φ 0 0 β (Φ 0 0 )] β(Φ 0 ) 2 ∆t δ x Φ 1 + cσ t - δ 2 β(Φ 0 0 )∆t + δ 2 β (Φ 0 ) β 2 (Φ 0 0 )∆t δ x Φ 1 4π I 4π dω + O(δ 2 x ) + O(δ 4 ), δ 2 ∂ t E = cσ a f (Φ) 4π I 4π dω -Φ + O(δ 4 ).
Note that we did not develop the last equation with respect to δ x : we only need the δ development in the following to conclude. The leading order terms allow identifying the asymptotic regime for the linearized equation:

                                   I 0 = 4π I 0 4π dω, cω • ∇I 0 + cσ t I 1 = cσ t 4π I 1 4π dω, ∂ t I 0 + cω • ∇I 1 + cσ t I 2 = 1 β(Φ 0 0 )∆t Φ 0 0 - 4π I 0 4π dω + δ x [β(Φ 0 0 ) -Φ 0 0 β (Φ 0 0 )]Φ 1 0 β 2 (Φ 0 0 )∆t + cσ t 4π I 2 4π dω + β (Φ 0 0 ) β 2 (Φ 0 0 )∆t Φ 1 0 δ x 4π I 0 4π dω + O(δ 2 x ), Φ 0 = Φ 0 0 + δ x Φ 1 0 + O(δ 2 x ) = 4π I 0 4π dω = Φ r,0 ,
Integrating the previous relations over angles and combining them yields:

∂ t Φ r,0 -∇ c 3σ t ∇Φ r,0 = 1 β(Φ 0 0 )∆t (Φ 0 0 -Φ r,0 ) + δ x Φ 1 0 β 2 (Φ 0 0 )∆t β(Φ 0 0 ) + β (Φ 0 0 )(Φ r,0 -2Φ 0 0 ) + O(δ 2 x ).
The leading order with respect to δ ensures for Φ 0 :

Φ 0 = Φ 0 0 + δ x Φ 1 0 + O(δ 2 x ) = Φ r,0 .
This finally leads to:

∂ t Φ r,0 -∇ c 3σ t ∇Φ r,0 = -C(Φ 0 0 ) δ x ∆t Φ 1 0 + O(δ 2 x ), = -C(Φ 0 0 ) δ x ∆t ∇Φ 0 + O(δ 2 x ), = -C(Φ 0 0 ) δ x ∆t ∇Φ r,0 + O(δ 2 x ).
The limit equation on time step [t n , t n + ∆t] of the IMC scheme obtained with an infinitely accurate MC approximation but taking into account an O(δ x ) discrepancy in the source term (during the source sampling phase) yields:

         ∂ t Φ r -∇ c 3σ n t ∇Φ r + C δ x ∆t ∇Φ r = O(δ 2 x ) + O(δ), E n+1 = E n + 1 β n ∆t t n+1 t n Φ r dt -Φ n + O(δ) + O(δ x ). (9)
This limit equation is an advection-diffusion one. The velocity of the advection operator depends on the discretisation parameters ∆t and δ x . They even compete during time step [t n , t n + ∆t]. On one hand, taking ∆t the smaller possible ensures recovering the equilibrium limit, but imposes a finer and finer spatial discretisation. In fact, the error can be very important for steep gradients of Φ r (front of a Marshak wave for example for which ∇Φ r 1). It is well-known the mechanism accumulates discrepancies proportionally to the number of time steps (due to the cycle-to-cycle differences of magnitude δx ∆t ∇Φ r , see [START_REF] Wollaber | Four decades of implicit monte carlo[END_REF]). The behaviour of the IMC scheme is close to the behaviour of the Dufort-Frankel scheme for parabolic equations. For this scheme, if Kδ x = ∆t → 0, the scheme is inconsistent. It is possible to force ∆t = K the greater possible (limit of stability) and make sure δ x → 0 to obtain a converging O(δ x ) scheme.

Let us verify and illustrate the above analysis on our 'fil rouge' test-case. The physical configuration is exactly the same as in figure 1 and Appendix A except we here fix N x = 20 and N M C ≈ 4× 10 5 and perform a convergence study with respect to ∆t ∈ {10 -8 , 5×10 -9 , 10 -9 , 10 -10 , 10 -11 , 10 -12 , 10 -13 , 10 -14 }. The results are displayed in figure 2. Note that the common ground for all the pictures of figure 2 is that the grid is coarse: N x = 20. The top-left picture testifies that even if the grid is coarse, the teleportation error is much less important with such big time step. On another hand, the time step is too big to ensure the relaxation of the radiative and material temperature and T m = T r on this same picture. Of course, as ∆t tends to zero on the successive other pictures of figure 2, we recover the fact that T r → T m as predicted by (9) (second line). But we also recover the fact that the solution looks more and more like the solution of an advection-diffusion equation rather than a diffusion one, as predicted by (9) (first line). In fact, with this study, we characterised the teleportation error: it stands for the resolution of an advection-diffusion equation with numerical advection speed depending on ∆t and δ x in a competitive way. Note that the latter numerical experiments show the previous analysis is not optimal: indeed, the results of figure 2 less as ∆t is decreasing. Still, we do not intend to be optimal here, we only wanted to characterise the competing behaviour observed in several papers [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF][START_REF] Wollaber | Four decades of implicit monte carlo[END_REF].

As emphasized in many publications, see [START_REF] Wollaber | Four decades of implicit monte carlo[END_REF] and the references therein, if the emission is non uniform within the cell, a reconstruction method can be introduced to estimate cf σΦ in at every beginning of time steps. This is commonly called a tilt (see [START_REF] Wollaber | Four decades of implicit monte carlo[END_REF][START_REF] Cleveland | Mitigating teleportation error in frequency-dependent hybrid implicit monte carlo diffusion methods[END_REF]) and it has been experimentally observed it reduces teleportation errors. In fact, it corresponds to a first order approximation of Φ(x) = Φ tilt (x) + O(δ 2

x ). The same asymptotic developments (Taylor of order 2 with respect to δ x and Hilbert one with respect to δ) applying an accurate tilt can be performed. The second order reads:

cσ a f Φ = cσ a Φ 0 1 + cσ a β(Φ 0 )∆t + δ 2 x cσ a Φ 2 1 + cσ a [β(Φ 0 ) -Φ 0 β (Φ 0 )]∆t (1 + cσ a β(Φ 0 )∆t) 2 + O(δ 3 x ), cσ a (1 -f ) = cσ a cσ a β(Φ 0 )∆t 1 + cσ a β(Φ 0 )∆t + δ 2 x cσ a Φ 2 cσ a β (Φ 0 )∆t (1 + cσ a β(Φ 0 )∆t) 2 + O(δ 3 x ), cσ a f = cσ a 1 + cσ a β(Φ 0 )∆t -δ 2 x cσ a Φ 2 cσ a β (Φ 0 )∆t (1 + cσ a β(Φ 0 )∆t) 2 + O(δ 3 x ), with Φ 2 = ∇ 2 xx Φ.
Plugged in the linearized transport equation on time step [t n , t n + ∆t] and going through the same steps as before leads to:

∂ t Φ r,0 -∇ c 3σ t ∇Φ r,0 = -D(Φ 0 0 ) δ 2 x 2∆t ∇ 2 Φ 0 + O(δ 3 x ), = -D(Φ 0 0 ) δ 2 x 2∆t ∇ 2 Φ r,0 + O(δ 3 x ).
Re-arranging the terms produces expression

         ∂ t Φ r -∇ c 3σ n t + D δ 2 x 2∆t ∇Φ r = O δ 3 x ∆t + O(δ), E n+1 = E n + 1 β n ∆t t n+1 t n Φ r dt -Φ n + O(δ) + O(δ x ). ( 10
)
It is a diffusion equation (up to order

O(δ 3 x )). The asymptotic diffusion coefficient is c 3σt + D δ 2 x 2∆t
and depends on discretisation parameters δ x and ∆t. Once again, both discretisation parameters compete during time step [t n , t n + ∆t]. This is attenuated by the fact that δ x is squared but it still accumulates discrepancies proportionally to the number of time steps due to the cycle-to-cycle differences of magnitude

δ 2 x 2∆t ∇ 2 Φ r .
Once again, we suggest verifying numerically our analysis. Let us revisit the two previous convergence studies (with respect to N x as in figure 1 and with respect to ∆t as in figure 2) in the same numerical conditions but with a first order tilted IMC scheme. Basically, the tilt we use here builds a linear approximation of Φ within each cell and corrects the position of the sampled source MC particle together with the Fleck factor. The results are displayed in figures 3 and 4. Figure 3 (to be compared with figure 1) presents the convergence study with respect to the spatial discretisation N x . For a coarse grid (top-left picture of figure 3), the teleportation error is still observable. But we can already see that it has a different nature than in figure 1. The propagation fronts are closer to the reference ones and, as expected, one order of magnitude is gained: the grid N x = 320 of the tilted solution gives similar results to the N x = 640 one for the non-tilted one, see figure 1. Note that the equilibrium is sufficiently well resolved with ∆t = 10 -12 , just as it was already the case with the non-tilted IMC scheme: we recover numerically the fact that tilts only ensure a spatial improvement.

N x = 20 N x = 40 T m IMC O1 tilt nx=20 T r IMC O1 tilt nx=20 T m = T r ,
Let us now comment on figure 4 displaying the convergence study with respect to ∆t in the same numerical conditions as in figure 2 but with the first order tilted IMC scheme. Once again, we recover the competing behaviour of ∆t and δ x : if ∆t is not small enough, equilibrium is not fulfilled, see for example the top pictures of figure 4. But making sure ∆t goes to zero increases the teleportation error, even if less pronunced than with the non-tilted IMC scheme. Still, a gain is observable as the spreading of the solution is more contained: the propagation fronts are less far away from the reference solution. In a sense, these two studies validate our numerical analysis: the tilted IMC scheme behaves a little bit more like a diffusion equation. But a first order tilt is not enough to capture the good diffusion coefficient except for a fine tuning between ∆t and ∆x, see [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF]. The same remark as above still applies here: if ∆t = Kδ 2

x with δ x going to zero, the tilted IMC scheme is inconsistent (behaviour to be compared with the Dufort-Frankel scheme for parabolic equations). On another hand, if δx ∆t = K is kept constant and δ x goes to zero, a convergence behaviour can still be observed. The above analysis shows that an accurate second order (tilt) reconstruction (i.e. Φ(x) = Φ tilt (x)+O(δ 3

x )) is mandatory to recover the correct diffusion coefficient for the regime δ → 0.

At this stage of the discussion, we could present results with higher-order tilts (or second order ones, allowing to, at least recover the good diffusion coefficient). But we prefer here taking few lines to list the pros and cons of having resort to them:

-first, let us keep in mind that the above analysis is independent of the kind of tilt one use. This means that polynomial, piecewise constant, splines etc. are all going to behave as above. They may give slightly different behaviours but they all are going to be the same with respect to the competing behaviour between ∆t and ∆x, as soon as a small spatial inaccuracy in the source sampling is introduced. Note that this behaviour is also observable with the teleportation corrector suggested in [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF] (which is probably much more accurate than a tilt but can not guaranty having strictly δ x = 0).

-Second, the analysis suggested the coefficients to build the high order tilts were computed accurately. In practice, they must be evaluated from an MC, i.e. potentially noisy, solution. This may induce some complications.

-Those difficulties may be amplified by the dimension of the problem: here, we only considered 1D spatial problems and only 2 coefficients needed to be evalutated (namely the coefficients Φ 0 , Φ 1 = ∇Φ in each cell). For 3D spatial problems and a first order tilt, at least 4 coefficients are needed. For 3D spatial problems and a second order tilt, 16 coefficients are needed. The reconstruction depends on the dimensionality of the problem which is kind of contradictory with having resort to an MC scheme, classically used when we need to be the less sensitive possible to the dimension.

-Besides, 3D gradient reconstructions are far from being obvious, especially on unstructured grids.

-Add to these points the fact that the previous analysis is only valid if the Fleck factor is also corrected: to obtain ( 9)-( 10 account in the expression of f n . Correcting the emission position is not enough to recover the previous results. There are many ways to correct the Fleck factor but every of them are complex to verify, even in 1D.

-Finally, the competing behaviour between ∆t and δ x , despite being tedious to tune (maximise ∆t to avoid having too much teleportation error but not too much to make sure the relaxation is fulfilled), can become impossible to handle in practice for models where (1) has to be coupled to another physics demanding small time steps for an accurate resolution (see for example radiative hydrodynamics [START_REF] Pomraning | The equations of radiation hydrodynamics[END_REF][START_REF] Mihalas | Foundations of Radiation Hydrodynamics[END_REF][START_REF] Castor | Radiation hydrodynamics[END_REF]).

For these reasons, instead of introducing an additional reconstruction, we prefer building MC schemes which do not suffer teleportation errors. In the following section, we present our new implicit MC schemes and detail some of its asymptotical properties.

The new Implicit Semi-analog Monte-Carlo (ISMC) solver

In this section, we present the new implicit MC scheme we felt the need to build due to the previously detailed reasons. We here want to build a scheme for which there is no spatial discrepancies at the emission of a radiative MC particle: in other words, we want to cancel, by construction, the error in the sampling of the emitted photons (we insist we want δ x = 0 and not δ x ≈ 0).

Note that, to our knowledge, at least one attempt has been made in order to do so within an IMC and source sampling framework in [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF]. The main idea consists in storing source positions along the flight path of each particles to use the banks to emit photons at the beginning of the next time step. The memory consumption consequently grows with the number of MC particles and the number of banking sites. It also has an additional complexity (how to determine the source sites and how to sample the photons from them) and an overcost [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF]. But overall, the authors still observe their teleportation corrector does not change the competing 7 behaviour between ∆t and δ x . Still, it is interesting noticing that the solver in [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF] has a key ingredient in common with the scheme we develop in the following (namely the SMC scheme of [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF]).

As explained before, we want to cancel the discrepancies at emission. MC schemes for neutronic applications [START_REF] Spanier | Monte Carlo Principles and Neutron Transport Problems[END_REF][START_REF] Lewis | Computational Methods of Neutron Transport[END_REF] generally inject the fission sources into the scattering cross-sections. This allows avoiding having resort to source sampling to model the neutrons emitted during a fission reaction. By construction, the fission neutrons are emitted exactly at the fission banks without the need for storing them in memory. We would like to do the same for the photonic application of this paper. To do so, let us rewrite (1) in a manner allowing to integrate the source term into a scattering part. For this, let us first rewrite (1) with respect to E as

       1 c ∂ t I + ω • ∇I + σ t I = σ a η(E)E + σ s 4π I dω 4π , ∂ t E = 4π cσ a I 4π -η(E)E dω . (11)
In [START_REF] Wollaber | On the stability of the ahrens-larsen or smc equations for thermal radiative transfer[END_REF], we introduced B(E) = η(E)E. We now suggest introducing variable e(x, t, ω) defined by 1 4π 4π e(x, t, ω)dω = E(x, t). Variable e is only an auxiliary variable whose first order moment with 7 See figure 5 of [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF] and the comment "Decreasing the time step [...] increases the amount of teleportation error".

respect to the angular variable ω coincides with E. System (11) can then be rewritten in term of unknowns (I, e) as8 :

     1 c ∂ t I + ω • ∇I = +σ a η(E) 4π e 4π dω -σ t I + σ s 4π I dω 4π , ∂ t e = -cσ a η(E)e + 4π cσ a I 4π dω . ( 12
)
This system is still nonlinear. From now on, let us assume that η(E(x, t)) = η(x, t) for t ∈ [t n , t n + ∆t]. With the previous hypothesis, [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF] becomes

     1 c ∂ t I + ω • ∇I = +σ a η 4π e 4π dω -σ t I + σ s 4π I dω 4π , ∂ t e = -cσ a ηe + 4π cσ a I 4π dω . (13)
Independently of the choice of η, (13) is linear (η 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 + ∆t]. Of course, in practice, a choice must be made for η (explicit, implicit etc.): let us postpone this choice and rather go through few properties of the above linearization (those remains independent of the choice of η).

Few properties of the linearised system (13)

Let us list few properties one can expect asymptotically (i.e. with an important number of MC particles) from an MC resolution of (13) on time step [t n , t n + ∆t].

The multigroup structure of (13) with analytical basis functions

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 in neutronics for example, see [START_REF] Lewis | Computational Methods of Neutron Transport[END_REF][START_REF] Spanier | Monte Carlo Principles and Neutron Transport Problems[END_REF]. But in this case, there are only two groups and the basis functions are analytical, given by φ 0 (v) = δ c (v) and φ 1 (v) = δ 0 (v). To clarify this point, let us put forward the expression of the total and scattering cross-sections of linear system [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF].

To solve [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF], let us build ψ(t, x, ω, v), a new unknown depending on one more dimension and on unknowns (I, e) solutions of [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF]. Variable v is chosen such that ψ(x, t, ω, v) = I(x, t, ω)δ c (v) + e(x, t, ω)δ 0 (v). In the latter expression, δ 0 , δ c are such that

{V } δ c (v)dv = δ V,c and {V } δ 0 (v)dv = δ V,0 ,
where δ V,k is the Kronecker symbol 9 and {V } denotes the singleton V . In fact, v is nothing more than a velocity which can be c for photons or 0 for matter. Let us now build the linear equation satisfied by ψ. Expression (13) can be rewritten (we drop the dependences for conciseness):

∂ t (Iδ c + eδ 0 ) + vω • ∇(Iδ c + eδ 0 ) = -c(σ t Iδ c + σ a ηeδ 0 ) + c σ a δ c η 4π e + σ a δ 0 4π I + σ s δ c 4π I .
It remains to make ψ appear in the collisional part. For the moment the integration is only over the angular distribution. Let us introduce the Kronecker symbols

δ 0,c (v , v) = δ 0 (v)δ c (v ) and δ c,0 (v , v) = δ c (v)δ 0 (v ).
We can then rewrite the above equation as:

∂ t ψ + vω • ∇ψ + c(σ t δ c (v) + σ a ηδ 0 (v))ψ = + c V 4π [σ a δ c,0 (v , v) + σ s δ c,c (v , v) + σ a ηδ 0,c (v , v)] ψ(v , ω )dv dω , (14)
with V = {0, c} a discrete space with only two elements. One can check that performing {c} [START_REF] Mckinley | Comparison of implicit and symbolic implicit monte carlo line transport with frequency weight vector extension[END_REF]dv allows recovering the transport equation of ( 13) and that {0} (14)dv leads to the second equation. We can identify the scattering Σ s (v, v ) and total Σ t (v) cross-sections to rewrite the system above under the general form:

∂ t ψ + vω∇ψ + cΣ t (v)ψ = V 4π cΣ s (v , v)ψ(v , ω )dv dω .
In the above expression, we have:

Σ t (v) = σ t δ c (v) + σ a ηδ 0 (v), and Σ s (v , v) = σ a δ c,0 (v , v) + σ s δ c,c (v , v) + σ a ηδ 0,c (v , v).
Let us rewrite the scattering part as Σ s (v

)P s (v , v) = Σ s (v , v). This implies Σ s (v) = Σ s (v , v)dv , = [σ a δ c,0 (v , v) + σ s δ c,c (v , v) + σ a ηδ 0,c (v , v)] dv , = (σ s + σ a η)δ c (v) + σ a δ 0 (v).
By definition of P s we have:

P s (v , v) = σ s (v , v) σ s (v) , = [σ a δ c,0 (v , v) + σ s δ c,c (v , v) + σ a ηδ 0,c (v , v)] (σ s + σ a η)δ c (v) + σ a δ 0 (v) .
The above expression can be considerably simplified by noticing that

for v = 0, P s (0, v ) = δ c (v ), for v = c, P s (c, v ) = σ s δ c (v ) + σ a ηδ 0 (v ) σ s + σ a η ,
so that P s resumes to

P s (v, v ) = δ 0 (v)δ c (v ) + δ c (v) σ s δ c (v ) + σ a ηδ 0 (v ) σ s + σ a η .
Now, we are interested in a direct resolution of (13) on time step [t n , t n + ∆t]. In other words, we need (cf. [START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF]) to characterise the total and scattering cross-sections of the adjoint form of [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF]. It is given by

-∂ t ψ -vω • ∇ψ + cΣ t (v)ψ = V 4π cΣ S (v , v)ψ(v , ω )dv dω , (15)
where

Σ S (v, v ) = Σ S (v)P S (v, v ) = Σ s (v )P s (v, v ).
Then we have

Σ S (v) = V Σ s (v )P s (v, v )dv , = V ((σ s + σ a η)δ c (v ) + σ a δ 0 (v )) δ 0 (v )δ c (v) + δ c (v ) σ s δ c (v) + σ a ηδ 0 (v) σ s + σ a η dv , = (σ a + σ s )δ c (v) + σ a ηδ 0 (v) = σ t δ c (v) + σ a ηδ 0 (v) = Σ t (v).
Few calculations, similar to the already performed one before to identify P S , show that:

P S (v, v ) = δ 0 (v)δ c (v ) + δ c (v) σ s δ c (v ) + σ a δ 0 (v ) σ s + σ a .
Note that P s and P S only differ from the fact that η is factor of σ a in P s and not in P S .

With the above calculations, we identified the direct and adjoint cross-sections for (13) to be revisited as a 2-group 10 linear transport equation [START_REF] Hilbert | Begrundung der kinetischen gastheorie[END_REF]. The main interest of the above identifications concerns the fact that we here showed the new solver we will build in the next sections can be easily implemented in an already existing MC multigroup simulation code. The relevance of the above calculations will also be highlighted later on when the numerical MC resolution will be tackled.

Conservativity (13) on time step [t n , t n + ∆t]

Independently of the choice of η, the linearized system is conservative as

∂ t V 4π ψdvdω +∇ • V 4π vωψdvdω = 0, ∂ t 4π Idω + 4π edω +∇ • 4π cωIdω = 0.
This implies that if a conservative MC scheme is used to discretise (13), the MC discretisation will be conservative through the successive times steps. Care will be taken to build such MC scheme in the next parts.

4.1.3. The equilibrium diffusion limit on [t n , t n + ∆t] for system [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF] Let us finally show why linearization ( 13) is satisfactory with respect to the equilibrium diffusion limit (4) we aim at capturing accurately. On time step [t n , t n + ∆t], in the limit defined by ( 5) the previously presented linearization resumes to

     δ c ∂ t I + ω • ∇I = + 1 δ σ a η 4π e 4π dω - 1 δ σ t I + σ s 4π I dω 4π , ∂ t e = -1 δ 2 cσ a ηe + 4π 1 δ 2 cσ a I 4π dω . ( 16
)
10 With analytical basis functions.

Let us perform a Hilbert development [START_REF] Hilbert | Begrundung der kinetischen gastheorie[END_REF][START_REF] Densmore | Asymptotic equilibrium diffusion analysis of time-dependent monte carlo methods for grey radiative transfer[END_REF] of the above system with respect to δ by expanding I, e and η as formal power series of δ. System (16) yields for 0 th , 1 st and 2 nd orders:

4π (η 0 e 0 -I 0 )dω = 0, cω • ∇I 0 = cσ a 4π (η 0 e 1 4π + η 1 e 0 4π ) + cσ s 4π I 1 4π dω -cσ t I 1 , ∂ t I 0 + cω • ∇I 1 = cσ a 4π (η 0 e 2 4π + η 1 e 1 4π + η 2 e 0 4π ) + cσ s 4π I 2 4π dω -cσ t I 2 .
It leads to (we used the notations

E i = 4π e i 4π , Φ i r = 4π I i
4π ):

η 0 E 0 = Φ 0 r , c 3 ∇Φ 0 r = -σ t Φ 1 r , ∂ t Φ 0 r cσ -∇ c 3σ t ∇Φ 0 r = σ a (η 0 E 2 + η 1 E 1 + η 2 E 0 -Φ 2 r ).
The same development of the second equation yields:

η 0 e 0 = 4π I 0 4π dω, η 0 e 1 + η 1 e 0 = 4π I 1 4π dω, ∂ t e 0 = -cσ a η 0 e 2 + η 1 e 1 + η 2 e 0 - 4π I 2 4π dω .
Integrating with respect to the angular variable leads to :

η 0 E 0 = Φ 0 r , η 0 E 1 + η 1 E 0 = Φ 1 r , ∂ t E 0 = -cσ a (η 0 E 2 + η 1 E 1 + η 2 E 0 -Φ 2 r
). Some equations are redundant (hence no incompatibility) but we finally obtain the asymptotic limit for the linearized system:

     Φ 0 r = η 0 E 0 + O(δ) = Φ 0 + O(∆t) + O(δ), ∂ t (E 0 + Φ 0 r ) -∇ • c 3σ n t ∇Φ 0 r = O(δ). ( 17
)
We recall the above properties are independent of the choice of η. To sum up, any MC scheme discretising system (13) will capture the equilibrium limit up to a fine enough time discretisation (cf. the O(∆t) term in the first line of ( 17)) and the diffusion limit (second line of ( 17)).

4.2.

The explicit SMC solver of [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF] as a particular case of linearisation [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF] Let us here take few lines to emphasize the fact that the SMC solver described in [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF] is a particular case of linearisation [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF]: it corresponds to the case η = η n and leads to

     1 c ∂ t I + ω • ∇I = +σ a η n 4π e 4π dω -σ t I + σ s 4π I dω 4π , ∂ t e = -cσ a η n e + 4π cσ a I 4π dω . ( 18
)
This means the few properties highlighted in the previous sections also apply to [START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF]. Of course, the multigroup structure and the conservativity were already put forward in [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF] but in some particular cases (namely σ s = 0, cross-sections independent of the temperatures, perfect gas eos). The analysis here is more general and, to our knowledge, original with respect to the equilibrium diffusion limit studied in section (4.1.3). The cross-sections are given by (adjoint form for a direct resolution)

-Σ t (v) = Σ S (v) = σ t δ c (v) + σ a η n δ 0 (v) so that the equivalent of the absorption cross section Σ A (v) = Σ t (v) -Σ S (v) = 0 is zero.
-and

P S (v, v ) = δ 0 (v)δ c (v ) + δ c (v) σsδc(v )+σaδ0(v ) σs+σa .
Having Σ A = 0 corresponds to the particular case where sampling the interaction time from Σ S or Σ t is equivalent. The weight of each MC particle along the flight path remains the same as the weight modification during interval of time [t 0 , t 1 ] is by definition e -vΣ A (v)(t1-t0) = e 0 = 1. With this remark, we only want to point out the fact that the scheme is called 'semi-analog' in [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF] but could have been called 'non-analog' due to the fact both are equivalent for Σ A = 0. Now, the MC scheme to solve [START_REF] Hilbert | Begrundung der kinetischen gastheorie[END_REF], presented also in [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF], implies -letting 'photon' MC particles (v = c) travel at speed of light c until it undergoes a collison. The interaction time is sampled from an exponential distribution of parameter cΣ t (c) = cΣ S (c) = cσ n t .

-A 'matter' MC particle is held motionless (v = 0) at its absorption site until it is re-emitted.

The time to emission is sampled with an exponential distribution of parameter cΣ t (0) = cΣ S (0) = σ n a η n .

-The scattering term P S tells a 'matter' MC particle becomes a photon with probability 1 once the interaction time reached: indeed P S (0, v ) = δ c (v ).

-On another hand, 'photon' MC particles can remain in a 'photon' state (i.e. endure a scattering) with probability P S (c, c) = σs σt or be absorbed (i.e. become a 'matter' MC particle) with probability σa σt . The identification of Σ t , Σ S , P S as in the previous sections allows being in the conditions of theorem 3.2.1 of [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF] and ensures, if the above treatments are performed to the MC particles, the convergence toward system (18) on time step [t n , t n + ∆t] as N M C (the number of MC particles) grows.

With this tracking algorithm, space locations of absorption sites of the particles are never lost (by construction δ x = 0). This remarkable property of this MC solver is to be tempered with the fact it may need severe restriction on the time step due to the poor explicit time discretisation of η as outlined by [START_REF] Wollaber | Four decades of implicit monte carlo[END_REF].

In practice, the SMC scheme is computationally unusable due to the small time steps it needs for stability [START_REF] Wollaber | On the stability of the ahrens-larsen or smc equations for thermal radiative transfer[END_REF]. The solver also has the reputation of not being able to deal with arbitrary equations of state (only perfect gases). The aim of the next section is to build a new scheme combining the advantages of both SMC (δ x = 0) and IMC (implicitness) solvers together with being able to take into account arbitrary equations of state (eos).

Combining the benefits of SMC and IMC: the new ISMC scheme

In this section, we build a new MC scheme combining the properties of both the SMC (with respect to the teleportation error) and IMC (with respect to the time discretisation) ones. Just as Fleck & Cummings suggested replacing B by B n+1 to go from ( 1) to ( 6), we suggest replacing η by η n+1 in [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF]. In the next lines, we are going to go through exactly the same steps as in section 2 but we are going to work on η n+1 instead of B n+1 . The first step yields,

     1 c ∂ t I + ω • ∇I = +σ a η n+1 4π e 4π dω -σ t I + σ s 4π I dω 4π , ∂ t e = -cσ a η n+1 e + 4π cσ a I 4π dω . ( 19
)
Let us identify the equation satisfied by η. Remember, by definition, we have B(E) = η(E)E so that we have

∂ t B(E) = E∂ t η(E) + η(E)∂ t E = (Eη (E) ζ(E) + η(E))∂ t E. ( 20
) Note that ζ(E) = Eη (E) = E dη dE (E) echoes the role of β(E) = dΦ dE (E) (see section 2)
for the IMC scheme. In particular, just as β helps dealing with arbitrary eos for the IMC linearisation, ζ will make sure the solver can be used in the same conditions. Replacing ∂ t E by its expression on both sides of the above equation reads

∂ t η(E) = ζ(E)cσ a 4π I 4π -η(E)E . (21)
Following [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF] (steps recalled in section 2), we propose an estimation of ζ n , η n+1 using an explicitimplicit time discretization:

∂ t η = ζ n cσ n a 1 E 4π I 4π dω -η n+1 . ( 22
)
We now integrate the above equation with respect to time to write

η n+1 = η n + ζ n cσ n a t 0 1 E 4π I 4π -∆tη n+1 , ≈ η n + ζ n cσ n a ∆t 1 E 4π I 4π -η n+1 . ( 23
)
We finally have

η n+1 = η n 1 1 + ζ n cσ n a ∆t + ζ n cσ n a ∆t 1 + ζ n cσ n a ∆t 1 E I 4π , = η n χ n + (1 -χ n ) 1 E 4π I 4π , (24)
with χ n = 1 1+cσ n a ζ n ∆t . Note that χ n echoes the Fleck factor f n = 1 1+cσ n a β n ∆t introduced in section 2. It also plays the same role with respect to the eos via ζ n instead of β n . Remark 4.1. Quantity ζ(E) = Eη (E) plays the same role as β(E) = Φ (E) ≥ 0 in an IMC framework. The positiveness of β is a sufficient condition for the convenient application of an MC scheme to the IMC linearisation (because β n > 0 =⇒ f n ∈]0, 1]). Otherwise, time step restrictions are necessary (in order to have 1 + cσ n a β n ∆t > 0 and ensure the positiveness of the Fleck factor). For the new quantity ζ we introduced, we have

ζ(E) = β(E) -η(E) = β(E) -Φ(E) E , = β(E) -Φ(E)-Φ(0) E-0 , = β(E) -1 E-0 E 0 dΦ dE (e)de, = β(E) -1 E-0 E 0 β(e)de.
The mean value theorem ensures, under smoothness conditions on β, that there exists -for a perfect gas eos, C v (T ) = C v and the positiveness is fulfilled as

Ẽ ∈]0, E[ such that 1 E-0 E 0 β(e)de = β( Ẽ). This implies ζ(E) = β(E) -β( Ẽ). If furthermore β is non-decreasing (i.e. d 2 Φ dE 2 = β (E) ≥ 0), the positiveness of ζ(E) is ensures ∀E ∈ R + as β(E) -β( Ẽ) = ζ(E) ≥
ζ(E) = 3η(E) ≥ 0. -For a cold solid, C v (T ) = C v T θ 3 + K where θ is a Debye temperature and T ∈ [θ, ∞[, the positiveness is fulfilled if K > 0 as β (T ) = 12aT 2 K ρC 2 v (T ) > 0.
The eos is not the main purpose of this paper and the positiveness of ζ may be, as illustrated above, case dependent. All in all, if ζ > 0, this ensures the positiveness of the MC scheme (in the sense it cannot produce negative densities of material and radiative energy). In the following sections, we will call χ n the modified Fleck factor in the following. In the next benchmarks, perfect gas are considered and χ n = 1 1+3cσ n a η n ∆t . To our knowledge, being able to deal with arbitrary eos in the SMC framework is new and original. With this approximation, one can rewrite the transport equation in [START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF] using the new estimation of η on the time step:

1 c ∂ t I + ω • ∇I +σ n t I =χ n σ n a η n 4π e 4π + ((1 -χ n )σ n a + σ n s ) 4π I dω 4π , ∂ t e +χ n cσ n a e =χ n cσ n a 4π
Idω.

(25)

Like in [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF], the implicit time discretization gives rise to a pseudo-scattering term in the transport equation which tends to unstrech the stiff coupling between radiation and matter. Furthermore, with the converging MC tracking algorithm we describe next, space locations of absorption sites of the particles are never lost (by construction δ x = 0). System (25) is closed, linear, explicit (i.e. only quantities at time t n appear) and exactly conservative in total energy. The cross-sections are given by (adjoint form for a direct resolution)

-Σ t (v) = Σ S (v) = σ t δ c (v) + σ a χ n η n δ 0 (v) so that the equivalent of the absorption cross section Σ A (v) = Σ t (v) -Σ S (v) = 0 is zero.
-and

P S (v, v ) = δ 0 (v)δ c (v ) + δ c (v) [σs+(1-χ n )σa]δc(v )+σaχ n δ0(v ) σt .
Having Σ A = 0 corresponds to the particular where sampling the interaction time from Σ S or Σ t is equivalent. This means that once again, the scheme is called 'semi-analog' but could have been called 'non-analog' due to the fact both are equivalent for Σ A = 0. Now, the new MC scheme to solve (15) implies -letting 'photon' MC particles (v = c) travel at speed of light c until it undergoes a collison. The interaction time is sampled from an exponential distribution of parameter cΣ t (c) = cΣ S (c) = cσ n t .

-A 'matter' MC particle is held motionless (v = 0) at its absorption site until it is re-emitted.

The time to emission is sampled with an exponential distribution of parameter cΣ t (0) = cΣ S (0) = σ n a χ n η n .

-The scattering term P S tells a 'matter' MC particle becomes a photon with probability 1 once the interaction time reached: indeed

P S (0, v ) = δ c (v ).
-On another hand, 'photon' MC particles can remain in a 'photon' state (i.e. endure a scattering) with probability P S (c, c) = σs+(1-χ n )σa σt or be absorbed (i.e. become a 'matter' MC particle) with probability σaχ n σt . The identification of Σ t , Σ S , P S as in the previous sections allows being in the conditions of theorem 3.2.1 of [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF]: it ensures, if the above treatments are performed to the MC particles, the convergence with respect to N M C toward system (18) on time step [t n , t n + ∆t]. Note that the new presented implicit solver is easy to implement if one has already access to an SMC implementation: one just has to correct few samplings thanks to the modified Fleck factor χ n (compare the few items of sections 4.2 and 4.3).

In a sense, we could here sum up the original idea of this scheme (and whole paper) is to apply the implicit time discretization introduced by Fleck and Cummings (IMC) in [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF] to the approximate system (SMC) described by Ahrens and Larsen in [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF]. For this reason, we call this new MC scheme the ISMC one (combining IMC [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF] and SMC [START_REF] Ahrens | A semi-analog monte carlo method for grey radiative transfer problems[END_REF]).

We now suggest revisiting the 'fil rouge' problem of this paper together with some benchmarks from the literature with the MC scheme we just described.

Numerical results

In this section, we first come back to our 'fil rouge' application with the new ISMC scheme we presented in the previous section. We then compare our new MC solver with the IMC solver on Marshak wave test-case presented in [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF].

Back to the 'fil rouge' application

One last time, we suggest verifying numerically our analysis on the 'fil rouge' problem of figures (no tilt) 1-2 and (first order tilt) 3-4. Let us revisit the two previous convergence studies (with respect to N x as in figures 1-3 and with respect to ∆t as in figures 2-4) in the same numerical conditions but with the new ISMC scheme. The results are displayed in figures 5 and 6. Figure 5 (to be compared with figures 1-3) presents the convergence study with respect to the spatial discretisation N x . With the new ISMC scheme we suggest in this paper, the equilibrium diffusion limit is captured on a coarse grid N x = 20, see figure 5 (top-left). As the number of cells grows, the noise induced by the MC resolution of (11) also grows: as we keep N M C fixed, the number of particles per cell becomes smaller and smaller. But the new implicit SMC scheme converges faster than the tilted or not IMC ones with respect to N x as, by construction, the spatial discrepancy at emission δ x = 0. As a consequence, a good accuracy can be obtained with few cells but also few particles. Let us revisit the convergence study with respect to ∆t for fixed N x , N M C in exactly the same conditions as in figures 2-4. Figure 6 presents the results obtained with the ISMC scheme. First, for big time steps, the relaxation of the material and radiative temperatures is not fulfilled. This was already the case for the IMC and tilted IMC solvers and is predicted by the numerical analysis [START_REF] Lewis | Computational Methods of Neutron Transport[END_REF]. On another hand, making sure ∆t tends to zero does not degradate the quality of the solution: no teleportation error is produced and the relaxation is better and better captured. Finally, we do not experience a competitive behaviour between the spatial and time discretisation parameters as was the case for the (tilted or not) IMC scheme with source sampling.

In the next section, we revisit one benchmark of the litterature with our new ISMC solver.

Marshak wave test-problem from [3]

In this section, we revisit a benchmark from [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF]. It is, to our knowledge, the only benchmark of the literature putting forward the competing behaviour between the spatial and time discretisation (i.e. with a convergence study with respect to ∆t for a fixed spatial discretisation).

Let us begin by a short description. The test-case corresponds to the study of a 1D Marshak wave [START_REF] Mihalas | Foundations of Radiation Hydrodynamics[END_REF] with dimensionless units. A black body heats the left boundary of the domain x ∈ D = [0, 4] with temperature T (x = 0) = 1. The radiation constant is a = 1 and so is the speed of light c = 1. There is no scattering (i.e. σ s = 0) and σ t (T m ) = σ a (T m ) = 10 T 3 m . Note that this benchmark will demonstrate our new MC solver can be used with temperature dependent opacities. Besides, the test-problem considers a perfect gas eos with ρ = 1 and C v = 7.14. The medium is initially cold as T (x, t = 0) = T 0 (x) = 10 -2 ∀x ∈ D = [0, 4]. We are here interested in the (material and radiative) temperature profiles at t * = 500.

Figure 7 presents a convergence study with respect to the spatial parameter N x = {20, 40, 80, 160, 320, 640, 1280, 2560} for ∆t = 5 × 10 -2 and N M C ≈ 10 7 . Here we have approximately 10 7 MC particles (and not exactly) because we have resort to source sampling (at least for the boundary for the ISMC scheme). Of course, IMC also needs it for the volumic sources. The figure displays a convergence study on the material temperature at t * = 500 for IMC (left) and ISMC (right). Note that a reference solution obtained from a (finely resolved) deterministic solver for the equilibrium diffusion limit (4) is also plotted. Once again, we recover the fact that for fixed ∆t, the IMC solver (see left picture of figure 7) converges toward the reference solution. But the convergence is slow. For coarse discretisations N x = {20, 40, 80, 160}, the error remains very important. The finer the mesh, the noisier the results as we kept N M C almost constant (up to source sampling fluctuations) for this study. On another hand, the ISMC solver (see right picture of figure 7) shows a fast convergence toward the reference solution. For N x = 20, the accuracy of the ISMC solver is equivalent to the one of IMC with N x = 160, hence a gain of a factor 80. Of course, the finer the mesh, the noisier the results, also for ISMC. Figure 8 presents the results obtained on the same quantity of interest but for N x = 10 (i.e. a coarse spatial discretisation as in [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF]) for ∆t ∈ {10 0 , 10 -1 , 10 -2 , 10 -3 }. The reference solution obtained from a deterministic solver (solving the equilibrium diffusion limit (4)) is also displayed. The IMC results are worse and worse as ∆t decreases. For such coarse mesh, ∆x is not small enough with respect to ∆t to compensate the advection term ∝ δx ∆t (see ( 9)): as a consequence, the IMC solutions are way faster than the reference solution. Of course, tilts or the teleportation correction of [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF] tend to mitigate the effect of decreasing the time step: those results are not recalled here but can be found in [START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF] in the same spatial conditions (N x = 10 zones, same scales for the figures) and for a slightly different range of time steps. On another hand, the results obtained with the ISMC solver (left picture of figure 8) tend to show that the new MC scheme is not sensitive to a decrease of ∆t: every curves for ∆t ∈ {10 0 , 10 -1 , 10 -2 , 10 -3 } are indistinguishable. Furthermore, despite the fact that we consider a coarse mesh, the base of the Marshak wave is rather well located. This benchmark confirms the new ISMC scheme has no teleportation error and allows avoiding the competing behaviour of ∆t and δx.

N x = 20 N x = 40 T m ISMC nx=20 T r ISMC nx=20 T m = T r ,

Performance studies for ISMC

Finally, we would like to finish with few performance studies. But before, we insist on the fact that in our opinion, having the above properties (no teleportation error, no competing behaviour between spatial and time discretisation parameters) for an MC scheme are priorities and must be dealt with prior to performance considerations.

Let us revisit the previous test-cases of sections 5.1 and 5.2 but from a performance point of view. Let us begin with the 'fil rouge' problem of section 5.1 and compare the restitution times for fixed parameters ∆t = 10 -12 , N x = 20 and N M C ≈ 10 6 :

-IMC: 228s., -tilted IMC: 336s., -ISMC: 229s.. Of course, the accuracy of the different solver is different for such fixed parameters. To give an idea of their accuracies, the conditions are almost the same as for the top left pictures of figure 1 (for IMC), figure 3 (for the 1 st order tilted IMC) and figure 5 (for ISMC). In term of accuracy, ISMC is better than the 1 st order tilted IMC which is better than IMC. But the overcost of tilted IMC is significative whereas IMC and ISMC have similar restitution times. We recall here that each (tilted or not) IMC computation is carried out with an analog MC scheme. Non-analog MC computations would lead to higher cost without any gains with respect to the teleportation error nor in term of smoothness of the results (indeed, for such coarse grid, the analog MC scheme is converged in term of MC particles).

Performing the same kind of study on the Marshak wave of section 5.2 and choosing ∆t = 5 × 10 -2 , N x = 20 and N M C = 10 6 leads to restitution times given by -IMC: 2854s., -ISMC: 1681s.. Once again, the solvers do not have the same accuracies: to give an idea, the results obtained in the conditions above are presented figure 7 (see the the coarser spatial discretisations). ISMC is both more accurate and faster than IMC.

Let us now compare restitution times for fixed accuracies. Obtaining similar accuracies with three solvers having three discretisation parameters is quite complex. Here is what we suggest. For both problems of sections 5.1 and 5.2, we first select a time discretisation which ensures having T m ≈ T r with a good accuracy (less than 1% error for T m (x, t * ) ≈ T r (x, t * ) with t * the final time). We then choose a spatial discretisation ensuring the chronometry of the radiative wave is converged up to 1% with respect to the reference solution (obtained with a finely resolved deterministic solver for (4)). We then make sure the noise are almost equivalent by choosing the number N M C of MC particles. For the fil rouge problem of section 5.1, we obtain:

-IMC: ∆t = 10 -12 , N x = 1280 and N M C ≈ 10 6 , restitution time 1300s. -ISMC: ∆t = 10 -12 , N x = 80 and N M C ≈ 10 6 , restitution time 305s.

To give an idea of the accuracies of the three above solvers, see figure 9 (left). Note that we should have used more than N M C = 10 6 MC particles for the IMC schemes to reach the accuracy of ISMC.

In other words, ISMC is faster and more accurate. Finally, for the Marshak wave problem of section 5.2, we obtain:

-IMC: ∆t = 5 × 10 -2 , N x = 1280 and N M C ≈ 2 × 10 7 , restitution time 3570s.

-ISMC: ∆t = 5 × 10 -2 , N x = 160 and N M C ≈ 2 × 10 7 , restitution time 1584s.

Figure 9 (right) gives an idea of the respective accuracies reached with the above parameters for IMC and ISMC. Note that once again, we should have used more MC particles for the IMC solver to give equivalent results as ISMC. Once again, ISMC is faster and more accurate.

Conclusion

In this paper, we suggest a new implicit MC scheme for photonics. It combines the advantages of IMC (implicitness) and of SMC (no teleportation error) together with a generalisation for arbitrary equations of state. The scheme is called ISMC. An asymptotic analysis of ISMC scheme yields, to leading order, a valid implicit conservative discretisation of the equilibrium diffusion equation (4). The MC scheme does not introduce 'source sampling' as a resolution strategy and consequently does not suffer the teleportation error (avoiding potential competing discretisation parameters, no cycle-to-cycle error explosion). Dealing with arbitrary eos is as straightforward as for the IMC linearisation. However, the leading order radiation intensity is not given by a planckian at the local end of time-step (see first equation of ( 17)) in general. Thus, the above method only has the diffusion limit, and not the equilibrium diffusion one, or only up to order O(∆t). Still, time steps to ensure equilibrium are affordable (just as with an IMC solver). The scheme also has the following properties:

-It is a convergent MC scheme for the linearized problem on each time step [t n , t n + ∆t].

-It is similar to a multigroup solver except the basis functions within the (2) groups are analytical (δ 0 (v) and δ c (v)). Hence, the solver is easy to plug in an existing multigroup simulation code. It is also similar to an SMC solver and can also be easily implemented in a simulation code implementing SMC (few modification via the modified Fleck factor).

-The solver is conservative: this is mainly due to the choice of the linearized system on time step [t n , t n+1 = t n + ∆t]. In practice, an MC particle switches status between 'matter' and 'radiation' so that the energy is always conserved for the whole system 'matter+photons'.

-The solver, just as the IMC linearization (7), allows capturing the diffusion limit and the equilibrium one up to O(∆t).

-The solver has, by construction, no teleportation error and does not need source tilting: radiative sources are emitted at the emission banks (no need to store them), during the tracking of the MC particles. The operations one must perform to each MC particle are simple (no tilts, not teleportation corrector etc.).

-From an HPC point of view, the solver only needs an 'estimated value' estimator to update (I, E) at each beginning of time step. No need for a track length estimator as for the IMC linearization. This allows gains in term of atomic for shared memory parallelisation, see [START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF]. Note that there is no need anymore for source sampling which is a sequential phase in the whole program (see [START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF]): in other words, the scheme goes in the direction of a better scalability.

-The spatial and time discretisation parameters are not competing (we, at least, experimentally verified it, but this is mainly due to the absence of additional spatial approximations during the MC phase).

-The solver once implemented, does not need any more effort in multidimensional spatial problems and on unstructured grids as one only has to tally the matter and radiative contributions within each cell volume at the end of the time step: no multidimensional gradient reconstructions needed.

-The spectral extension can be carried out naturally, this is a consequence of the structure of the linearized problem [START_REF] Clouët | Asymptotic diffusion limit of the symbolic monte-carlo method for the transport equation[END_REF]. Same remark concerning eventual anisotropic scattering. Taking into account the dependance with respect to the temperature of the opacities is straightforward as testifies both the description of the scheme of section 4.3 and the benchmarks of section 5.

-It can be extended to arbitrary eos under the condition ζ n can be estimated. In fact, ζ = E dη dE echoes the coefficient β = dΦ dE in the IMC framework. Exactly as β for IMC, function ζ may not be known analytically in general. But many numerical approximation strategies can be applied. Depending on the type of gas one needs to consider, the latter approximation may require finer time steps but will definitely not introduce any spatial discrepancies (we will still have δ x = 0 by construction).

-Finally, in a sense, this work also emphasizes the fact that the IMC implicitation of [START_REF] Fleck | An implicit monte-carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF] can be applied to different linearisation: it may benefit other regimes or even other coupled systems needing an MC resolution.

For all those reasons and the numerical results on the different benchmarks, we think we built a promising MC numerical solver. More properties will be studied in future works.

The 'fil rouge' problem considered all along the talk can be described as follows: let us consider a 1D spatial domain such that x ∈ Ω = [0, 1]. The domain is filled with a diffusive media σ t = 2000, with no (physical) scattering, i.e. σ s = 0 and σ t = σ a . Initially, a Heaviside of temperatures at equilibrium is set in the middle [0.4, 0.6] of domain Ω = [0, 1]. In other words, we have at t = 0: T m (x, t = 0) = T r (x, t = 0) = 2.3 × 10 7 1 [0.4,0.6] (x) + 2.3 × 10 4 1 [0,1]\[0.4,0.6] (x).

Note that 1 Ω (x) denotes the indicatrix of domain Ω. The initial condition is displayed in figure A.10 together with the solution of system (4) at final time T = 10 -8 . This reference solution has been obtained solving (4) with a deterministic solver (with a fine mesh).

Note that for time t ∈ [0, T ], the solution does not reach the boundaries. This test-case has been chosen precisely in order to avoid having resort to the sampling of boundary particles so that the only sampled particles for t > 0 are produced by the source term (for the IMC schemes). Benchmarks of the literature often consists of Marshak waves needing boundary particle sampling. Of course, more classical benchmarks are tackled in section 5. Besides, in order to avoid the multiplication of test-cases together with being able to show many different features of the different analysed schemes, we autorised ourselves to change few constants: the radiative constant is set to a = 10 -14 , the speed of light to c = 3 × 10 10 . A perfect gas is considered to that E(T m ) = ρC v T m with ρ = 20, C v = 4 × 10 7 . The configuration may not appear particularly physical but is sitll relevant for real life encountered difficulties.

TTFigure 1 :

 1 Figure 1: Comparisons of the material and radiative temperatures Tm and Tr obtained from (reference solution) a deterministic solver for the equilibrium diffusion limit (4) and the IMC approximation for Nx ∈ {20, 40, 80, 160, 320, 640, 1280, 2560}, ∆t = 10 -12 and N M C ≈ 4 × 10 5 MC particles (analog MC scheme + source sampling).

- 10 T- 12 T- 14 TFigure 2 :

 1012142 Figure2: Comparisons of the material and radiative temperatures Tm and Tr obtained from (reference solution) a deterministic solver for the equilibrium diffusion limit (4) and the IMC approximation for Nx = 20, N M C ≈ 4×10 5 MC particles (analog MC scheme + source sampling) and ∆t ∈ {10 -8 , 5×10 -9 , 10 -9 , 10 -10 , 10 -11 , 10 -12 , 10 -13 , 10 -14 }.

Figure 3 :

 3 Figure 3: Comparisons of the material and radiative temperatures Tm and Tr obtained from (reference solution) a deterministic solver for the equilibrium diffusion limit (4) and the first order tilted IMC approximation for Nx ∈ {20, 40, 80, 160, 320, 640, 1280, 2560}, ∆t = 10 -12 and N M C ≈ 4 × 10 5 MC particles (analog MC scheme + source sampling).

- 10 T- 12 T- 14 TFigure 4 :

 1012144 Figure4: Comparisons of the material and radiative temperatures Tm and Tr obtained from (reference solution) a deterministic solver for the equilibrium diffusion limit (4) and the first order tilted IMC approximation for Nx = 20, N M C ≈ 4 × 10 5 MC particles (analog MC scheme + source sampling) and ∆t = {10 -8 , 5 × 10 -9 , 10 -9 , 10 -10 , 10 -11 , 10 -12 , 10 -13 , 10 -14 }.

  0 and the positiveness of the modified cross-sections is ensured. Otherwise, the positiveness of ζ is closely related to the heat capacity C v (T ) just as the positiveness of β in the IMC framework. Let us consider two examples from [19]:

TTFigure 5 :

 5 Figure 5: Comparisons of the material and radiative temperatures Tm and Tr obtained from (reference solution) a deterministic solver for the equilibrium diffusion limit (4) and the ISMC approximation for Nx ∈ {20, 40, 80, 160, 320, 640, 1280, 2560}, ∆t = 10-12 and N M C ≈ 4 × 10 5 MC particles (analog MC scheme + source sampling).

∆t = 10 - 10 T- 12 T- 14 TFigure 6 :ISMCFigure 7 :

 1010121467 Figure6: Comparisons of the material and radiative temperatures Tm and Tr obtained from (reference solution) a deterministic solver for the equilibrium diffusion limit (4) and the ISMC approximation for Nx = 20, N M C ≈ 4×10 5 MC particles (analog MC scheme + source sampling) and ∆t = {10 -8 , 5×10 -9 , 10 -9 , 10 -10 , 10 -11 , 10 -12 , 10 -13 , 10 -14 }.

Figure 8 :

 8 Figure8: Material temperature profiles T (x, t * = 500) at t * = 500 for the IMC scheme (left) and the ISMC one (right) for Nx = 10 (same conditions as in[START_REF] Irvine | Reducing the spatial discretization error of thermal emission in implicit monte carlo simulations[END_REF]) and ∆t ∈ {10 0 , 10 -1 , 10 -2 , 10 -3 }. The figure also displays a reference equilibrium diffusion solution obtained with a deterministic scheme solving (4).

- 1 Figure 9 :

 19 Figure 9: Material temperature profiles T (x, t * ) for the...

The steeper the derivative of Φ with respect to E at E n (i.e. the greater dΦ dE (E n ) = β n ), the smaller the time step needs to be to ensure equilibrium.

Integrate (12) with respect to 1 4π 4π •dv and recover[START_REF] Wollaber | On the stability of the ahrens-larsen or smc equations for thermal radiative transfer[END_REF].

i.e. is such that δ V,k = 0 if V = k and δ V,k = 1 if V = k.

Appendix A. Details about the 'fil-rouge' problem of this article ( initial and boundary conditions, test-case justifications)

In this paper, we intensively make use of the configuration presented in this appendix. We call it our 'fil rouge' configuration as we solve it with several numerical schemes and progressively highlight the effect of each of them on this toy problem. The initial and boundary conditions together with the problem justifications are provided here for both, the sake of conciseness of the paper and of reproducibility of the results.