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

A gPC-intrusive Monte-Carlo scheme for the resolution of the uncertain linear Boltzmann equation

Keywords: transport, Monte-Carlo, uncertainty quantification, intrusive, non-intrusive, numerical scheme, generalized Polynomial Chaos

In this paper, we are interested in the resolution of the time-dependent problem of particle transport in a media whose composition is uncertain. The most common resolution strategy consists in running, at prescribed points in the uncertain space (at experimental designs points), a simulation device as a black-box and perform the uncertainty propagation resolution (i.e. resolution of SPDE). This kind of strategy is commonly called non-intrusive. The nonintrusive resolution can be carried with Monte-Carlo, quasi Monte-Carlo, generalized Polynomial Chaos etc. The latter is of interest in this document for its fast convergence rate. After going over and illustrating the main drawbacks of the non-intrusive (gPC or not) uncertainty propagation resolution for the linear Boltzmann equation in a simplified configuration, we build a new gPC based MC scheme solving intrusively the uncertain counterpart of the problem. The paper ends with some numerical examples.

Introduction

Polynomial Chaos [START_REF] Wiener | The Homogeneous Chaos[END_REF][START_REF] Ghanem | Stochastic Finite Elements: a Spectral Approach[END_REF][START_REF] Cameron | The Orthogonal Development of Non-Linear Functionals in Series of Fourier-Hermite Functionals[END_REF] and its generalizations (gPC 1 ) [START_REF] Wan | Beyond Wiener-Askey Expansions: Handling Arbitrary PDFs[END_REF][START_REF] Le Maître | Uncertainty Propagation using Wiener-Haar Expansions[END_REF][START_REF] Le Maitre | Multi-Resolution Analysis of Wiener-Type Uncertainty Propagation Schemes[END_REF][START_REF] Vos | Time dependent polynomial chaos[END_REF][START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF][START_REF] Wan | Multi-Element generalized Polynomial Chaos for Arbitrary Probability Measures[END_REF][START_REF] Poëtte | Non Intrusive Iterative Stochastic Spectral Representation with Application to Compressible Gas Dynamics[END_REF][START_REF] Gaël Poëtte | Iterative polynomial approximation adapting to arbitrary probability distribution[END_REF][START_REF] Paulson | Arbitrary polynomial chaos for uncertainty propagation of correlated random variables in dynamic systems[END_REF] have been successfully applied to take into account uncertainties in many physical domains (stochastic elastic materials [START_REF] Ghanem | Stochastic Finite Elements: a Spectral Approach[END_REF], finite deformations [START_REF] Acharjee | Uncertainty Propagation in Finite Deformations -A Spectral Stochastic Lagrangian Approach[END_REF], heat conduction [START_REF] Wan | Stochastic Heat Transfer Enhancement in a Grooved Channel[END_REF], incompressible flows [START_REF] Xiu | Stochastic Modeling of Flow-Structure Interactions[END_REF][START_REF] Najm | A Stochastic Projection Method for Fluid Flow I: Basic Formulation[END_REF][START_REF] Le Maître | A Stochastic Particle-Mesh Scheme for Uncertainty Propagation in Vortical Flows[END_REF], reacting flows and detonation [START_REF] Lucor | Multi-Physics Stochastic Design Optimization: Application to Reacting Flows and Detonation[END_REF], computational fluid dynamics [START_REF] Le Maître | Uncertainty Propagation using Wiener-Haar Expansions[END_REF][START_REF] Wan | Multi-Element generalized Polynomial Chaos for Arbitrary Probability Measures[END_REF][START_REF] Poëtte | Treatment of Uncertain Interfaces in Compressible Flows[END_REF] ...). It is commonly accepted it stands for an efficient alternative to Monte-Carlo (MC) methods in relatively low stochastic dimensions (small number of uncertain parameters). Two observations can be done considering the furnished literature on the subject:

(i) first, if gPC is very efficient in physical applications involving regular/smooth solutions (such as thermal heat [START_REF] Matthies | Galerkin methods for linear and nonlinear elliptic stochastic PDEs[END_REF][START_REF] Wan | Stochastic Heat Transfer Enhancement in a Grooved Channel[END_REF][START_REF] Xiu | A New Stochastic Approach to Transient Heat Conduction Modeling with Uncertainty[END_REF][START_REF] Hien | Stochastic Finite Element Modeling in Linear Transient Heat Transfer[END_REF], structural mechanics and reliability [START_REF] Sudret | Life-Time Reliability based Assessment of Structures submitted to Thermal Fatigue[END_REF][START_REF] Sudret | Stochastic Finite Element Methods and Reliability -A State of the Art Report[END_REF][START_REF] Blatman | Sparse Polynomial Chaos Expansions and Adaptive Stochastic Finite Elements using a Regression Approach[END_REF][START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF], ...), it needs to be wisely adapted/modified for problems involving strong nonlinearities and solutions bearing steep gradient to achieve comparable performances, see [START_REF] Le Maître | Uncertainty Propagation using Wiener-Haar Expansions[END_REF][START_REF] Le Maitre | Multi-Resolution Analysis of Wiener-Type Uncertainty Propagation Schemes[END_REF][START_REF] Wan | Multi-Element generalized Polynomial Chaos for Arbitrary Probability Measures[END_REF][START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF][START_REF] Poëtte | Uncertainty Quantification for Systems of Conservation Laws[END_REF]] amongst many others. In other words, the direct efficiency of gPC is closely related to the structure 2 of the set of partial differential equations (PDEs) modeling the physical phenomenon of interest. In this paper, we tackle uncertainty quantification with gPC applied to the linear Boltzmann equation. The relevance of the latter integrodifferential PDE is not to demonstrate. Amongst the applications (non exhaustive list), one can quote biology [START_REF] Perthame | Transport Equations in Biology[END_REF] with population dynamics, plasma physics (transport of ions and electrons) [START_REF] Duderstadt | Transport Theory[END_REF], photonics [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][START_REF] Lowrie | The coupling of radiation and hydrodynamics[END_REF][START_REF] Hermeline | A discretization of the multigroup pn radiative transfer equation on general meshes[END_REF][START_REF] Kaoua | Solution of the Nonlinear Radiative Transfer Equations by a Fully Implicit Matrix Monte Carlo Method Coupled with the Rosseland Diffusion Equation via Domain Decomposition[END_REF] or neutronics [START_REF] Spanier | Monte Carlo Principles and Neutron Transport Problems[END_REF][START_REF] Lewis | Computational Methods of Neutron Transport[END_REF][START_REF] Dufek | Stochastic approximation for Monte-Carlo Calculation of Steady-State Conditions in Thermal Reactors[END_REF][START_REF] Isotalo | Computational Methods for Burnup Calculations with Monte Carlo Neutronics[END_REF][START_REF] Isotalo | Preventing Xenon Oscillations in Monte Carlo Burnup Calculations by Enforcing Equilibrium Xenon Distribution[END_REF][START_REF] Dufek | The stochastic implicit Euler method: A stable coupling scheme for Monte Carlo burnup calculations[END_REF]. (ii) Second, gPC can be applied intrusively (see for examples [START_REF] Le Maître | Uncertainty Propagation using Wiener-Haar Expansions[END_REF][START_REF] Tryoen | Adaptive Anisotropic Stochastic Discretization Schemes for Uncertain Conservation Laws[END_REF][START_REF] Després | Robust Uncertainty Propagation in Systems of Conservation Laws with the Entropy Closure Method[END_REF]) or non-intrusively (the most common way) and whether one strategy performs better than the other is far from being obvious. An intrusive resolution implies modifying or even rewriting a simulation code. One needs to investigate on the intrications of the uncertainties within the set of PDEs of interest. Considerations about its parallel resolution are case dependent and complex.

A non-intrusive resolution on another hand uses an already existing code as a black-box just as a MC one 3 . It is clearly the most convenient and common strategy. Parallel considerations are immediate: one can launch as many independent runs as available computational devices without need for communications between processors (commonly called embarassingly parallel). The uncertainty analyst does not really need knowing the content of the black-box device, neither the underlying PDE structure it solves, nor the numerical solver embedded. In the following, we build a gPC-intrusive MC (gPC-i-MC ) scheme to solve the uncertain linear Boltzmann equation. This may appear surprising as MC methods seem intrinsically non-intrusive. Care will be taken in the following to highlight why such unconventional and original resolution strategy is relevant regarding this particular set of PDEs. To sum up, we will see that with the new approach, -the uncertainty analysis can be done with only a minor additional cost with respect to only one deterministic MC simulation solving the linear Boltzmann equation, -and can be achieved with minor modifications of an existing MC code solving the deterministic linear Boltzmann equation. In other words, in this context 4 , opening the black-box resolution code can be very efficient.

We are interested in the resolution the uncertain linear Boltzmann equation. It models the time-dependent problem of particle transport in a media whose composition or collisional characteristics are uncertain 5(reaction rates or compositions). We suppose transport to be driven by the linear Boltzmann equation [START_REF] Brunner | An efficient, robust, domain-decomposition algorithm for particle Monte Carlo[END_REF] for particles having position x ∈ D ⊂ R 3 , velocity v ∈ R 3 , at time t ∈ [0, T ] ⊂ R + and where the quantity u(x, t, v) is the density of presence of the particles at (x, t, v):

     ∂ t u(x, t, v) + v • ∇ x u(x, t, v) = -vσ t (η(x, t), v)u(x, t, v) + vσ s (η(x, t), v , v)u(x, t, v ) dv , u(x, 0, v) = u 0 (x, v). (1)
Note that in [START_REF] Brunner | An efficient, robust, domain-decomposition algorithm for particle Monte Carlo[END_REF], we introduced the notation |v| = v to denote the norm of the velocity v and ω = v v for the angular variable of the particle distribution. Equation (1) must come with proper boundary conditions for wellposedness [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF][START_REF] Golse | Transport et Diffusion[END_REF] but we omit them for the sake of conciseness. The left hand side of (1) will be hinted at as the streaming counterpart of (1) whereas its right hand side will be called the collisional one. The above equation is linear and can be used to model the behaviour of particles interacting with a background media defined via both the vector of compositions η(x, t) = (η 1 (x, t), ..., η M (x, t))

t and the microscopic properties of its components denoted by (σ α,m ) α∈{t,s},m∈{1,...,M } . The interaction of particles with matter is described through the macroscopic total interaction probability of particles with media σ t (x, t, v) and a scattering one σ s (x, t, v , v), both characterised by σ t (η(x, t), v) = M m=1 σ t,m (v)η m (x, t), and σ s (η(x, t), v , v) = M m=1 σ s,m (v , v)η m (x, t).

(

) 2
The description of the microscopic collision term (σ s,m ) m∈{1,...,M } can be decomposed into reaction crosssections: for example, it is common to separate the elastic reaction (σ S,m ) m∈{1,...,M } from the multiplicative one (σ f,m ) m∈{1,...,M } , especially in neutronic applications. The macroscopic collision term then reads:

σ s (η(x, t), v, v ) = M m=1 σ S,m (v, v )η m (x, t) + ν f,m (v)σ f,m (v, v )η m (x, t), (3)
where ν f,m (v) designates the multiplicity of reaction f for material m at velocity v. Note that in the following, we will gather every reactions under a unique collision term σ s for the sake of conciseness but the methodology described in this document applies to take into account uncertainties in the different reactions.

We aim at considering uncertainties in the collisional part of the transport equation, i.e. (σ α ) α∈{s,t} are characterised stochastic processes, and we want to quantify its impact on the particle flow. It is common, in uncertainty quantification works [START_REF] Ghanem | Stochastic Finite Elements: a Spectral Approach[END_REF][START_REF] Acharjee | Uncertainty Propagation in Finite Deformations -A Spectral Stochastic Lagrangian Approach[END_REF][START_REF] Wan | Stochastic Heat Transfer Enhancement in a Grooved Channel[END_REF][START_REF] Xiu | Stochastic Modeling of Flow-Structure Interactions[END_REF][START_REF] Najm | A Stochastic Projection Method for Fluid Flow I: Basic Formulation[END_REF][START_REF] Le Maître | A Stochastic Particle-Mesh Scheme for Uncertainty Propagation in Vortical Flows[END_REF][START_REF] Lucor | Multi-Physics Stochastic Design Optimization: Application to Reacting Flows and Detonation[END_REF], emphasizing the fact a quantity is uncertain by introducing explicitly an additional dependence to the unknown of interest with respect to an uncertain parameter (a stochastic process) here denoted by X(x, t, v): in other words, we have σ t (x, t, v, X(x, t, v)) and σ s (x, t, v , v, X(x, t, v)). Note that without loss of generality in the following, we will consider X is a vector X = (X 1 , ..., X Q ) t of Q independent random variables of probability measure dP X = Q i=1 dP Xi rather than a stochastic process: in theory, it is always possible to come back to such framework 6 . The stochasticity can here affect indifferently the media compositions or the microscopic cross-sections (reaction by reaction if wanted), the formulation is general enough. For the sake of conciseness in the following, we may drop the dependences and write

σ t (x, t, v, X(x, t, v)) = σ t (x, t, v, X) and σ s (x, t, v , v, X(x, t, v)) = σ s (x, t, v , v, X).
As a result, the uncertain problem (1) is linear but bears a stochastic process parametered by x, t, v, i.e. u(x, t, v, X), as a solution: solving the uncertain counterpart of (1) consequently resumes to solving the SPDE given by

         ∂ t u(x, t, v, X) + v∇ x u(x, t, v, X) = -vσ t (η(x, t), v, X)u(x, t, v, X) + vσ s (η(x, t), v , v, X)u(x, t, v , X) dv , u(x, 0, v) = u 0 (x, v). ( 4
)
The paper is organized as follow: in section 2, we will recall the most common gPC-based numerical strategy to solve (4), the non-intrusive one. It consists in running a deterministic simulation device to solve (1) at some prescribed (experimental design, see [START_REF] Fedorov | Theory Of Optimal Experiments[END_REF][START_REF] Atkinson | Optimum Experimental Designs, With SAS[END_REF]) points7 denoted by (X i ) i∈{1,...,N } of weights (w i ) i∈{1,...,N } to project the solution u on a gPC basis before performing some post-treatments to estimate the statistical observables of interest (mean, variance, probability of failure, moments etc.). We briefly present some numerical results in a very simple configuration. We even consider several couples of resolution strategies for the black-box code8 solving (1) and for the choice of the design points9 (X i , w i ) i∈{1,...,N } . Care will be taken in this section to illustrate the implications of such decoupled resolution strategies. In section 3, we focus on one of the resolution method for the deterministic system (1) which deserves a particular attention with respect to uncertainties, the Monte-Carlo method. We recall briefly the construction of a MC scheme to solve the deterministic linear Boltzmann equation (1) together with its convergence properties. These can be found in many books but briefly reminding them with UQ-friendly notations eases the description of the gPC-i-MC scheme we present in section 4 to solve [START_REF] Atkinson | Optimum Experimental Designs, With SAS[END_REF]. In particular, it eases the description of the minimal modifications to an existing MC code to do so. It consists in building a new Monte-Carlo scheme allowing the on-the-fly resolution of the uncertain counterpart of the linear Boltzmann equation made possible by -gPC, -the structure of the transport equation, -and the mathematical properties of MC methods. The benefit of the new Monte-Carlo scheme are then illustrated on various examples in section 5 and section 6 is devoted to concluding remarks.

Non-intrusive resolution of the uncertain linear Boltzmann equation

In this section, we very generally describe non-intrusive uncertainty propagation methods [START_REF] Crestaux | Polynômes de Chaos pour la Propagation et la Quantification d'Incertitudes[END_REF][START_REF] Martinez | Modélisation des Incertitudes par Polynômes de Chaos -Étude d'un Écoulement en Milieux Poreux[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF]. The description may, at first glance, look like a recipe but it is representative of its practical use. Recall random variable10 X has probability measure dP X , we aim at estimating a given statistical quantity of interest

I(u) = F (u(x, t, v, ξ)) dP X (ξ) = E [F (u(x, t, v, X))] , (5)
depending on the solution u of ( 4) and on a post-treatment of it, introduced in (5) via F . Post-treatment F can either be a (vectorial) functional or an operator. For example, -if u → F (u) = u then the statistical quantity I corresponds to the mean of u(x, t, v, X).

-If u → F (u) = u 2 it corresponds to its second order moment.

-

If u → F (u) = 1 [U,∞[ (u) then I(u)
becomes the probability of having u(x, t, v, X) beyond threshold U , commonly called a failure probability, see [START_REF] Sudret | Uncertainty Propagation and Sensitivity Analysis in Mechanical Models, Contribution to Structural Reliability and Stochastic Spectral Methods[END_REF].

-If u → F (u) = uφ X k where (φ X k (X)
) k∈N denotes the generalised Polynomial Chaos basis (gPC), see [START_REF] Wan | Long-Term Behavior of Polynomial Chaos in Stochastic Flow Simula-tions[END_REF][START_REF] Wan | Beyond Wiener-Askey Expansions: Handling Arbitrary PDFs[END_REF][START_REF] Wan | Multi-Element generalized Polynomial Chaos for Arbitrary Probability Measures[END_REF][START_REF] Ghanem | Stochastic Finite Elements: a Spectral Approach[END_REF][START_REF] Acharjee | Uncertainty Propagation in Finite Deformations -A Spectral Stochastic Lagrangian Approach[END_REF][START_REF] Wan | Stochastic Heat Transfer Enhancement in a Grooved Channel[END_REF][START_REF] Xiu | Stochastic Modeling of Flow-Structure Interactions[END_REF][START_REF] Najm | A Stochastic Projection Method for Fluid Flow I: Basic Formulation[END_REF][START_REF] Le Maître | A Stochastic Particle-Mesh Scheme for Uncertainty Propagation in Vortical Flows[END_REF][START_REF] Lucor | Multi-Physics Stochastic Design Optimization: Application to Reacting Flows and Detonation[END_REF], associated 11 to measure dP X then I is the k th gPC coefficient of u defined by:

I(u)(x, t, v) = u(x, t, v, X)φ X k (X) dP X = u X k (x, t, v), ∀k ∈ N. (6)
In this case, the P-truncated gPC expansion bears some interesting convergence properties [START_REF] Wiener | The Homogeneous Chaos[END_REF][START_REF] Cameron | The Orthogonal Development of Non-Linear Functionals in Series of Fourier-Hermite Functionals[END_REF] as

u P (x, t, v, X) = P k=0 u X k (x, t, v)φ X k (X) L 2 -→ P →∞ u(x, t, v, X). (7)
We will particularly focus on such developments in the following. -We insist F can also denote an operator: one may be interested not directly in u(x, t, v, X) but in some post-treatments with respect to the physical variables x, t, v. For example, we can have

F (u(•, t, •, X)) = 1 |D| × |F| 1 D (x)1 F (v)u(x, t, v, X) dx dv, (8)
where |D|, |F| denote the volumes of spatial and kinetic spaces D ⊂ R 3 and F ⊂ R 3 . In this case, we especially want to emphasize it is possible applying the material of the previous point to

F P (t, X) = P k=0 F X k (t)φ X k (X) -→ F (u(•, t, •, X)), (9)
to approximate it on a gPC basis. In other words, the next discussion remains very general with respect to the observable of interest, even if we basically focus on u in the following. The non-intrusive methodology then consists in several steps: (i) The first one corresponds to the discretisation of the random variable and its probability measure (X, dP X ) by a numerical integration method with N points:

(X, dP X ) ≈ (X i , w i ) i∈{1,...,N } . (10)
The notation [START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF] for the punctual discretisation of (X, dP X ) is very general and convenient as it can take into account many integration methods: suppose the points (X i ) i∈{1,...,N } are chosen sampled from the probability law of X and (w i = 1 N ) i∈{1,...,N } , then it corresponds to a Monte-Carlo integration [START_REF] Saporta | Probabilités, Analyse de Données et Statistique[END_REF] for the estimation of I. With the same writing, we can conveniently consider Gauss quadrature points, Latin Hypercube Samples, Sparse Grids etc. [START_REF] Saporta | Probabilités, Analyse de Données et Statistique[END_REF][START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF][START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF][START_REF] Halton | Algorithm 247: Radical-inverse quasi-random point sequence[END_REF][START_REF] Van Der Corput | [END_REF][START_REF] Sobol | Uniformly distributed sequences with an additional uniform property[END_REF][START_REF] Gautschi | Orthogonal polynomials: applications and computation[END_REF][START_REF] Trefethen | Is Gauss Quadrature Better than Clenshaw-Curtis?[END_REF]. The latter sets of points in dimension Q differ only by their asymptotic error analysis O(N β(Q) ), i.e. the weak or strong dependence of their convergence rates β(Q) with respect to the number of uncertain parameters Q.

For example β(Q) = β = -1 2 for MC methods: the convergence rate is slow but independent of Q. We have β(Q) = -1 for N Q equidistant points in [-1, 1] Q for example. It is well-known in 1D stochastic dimension, equidistant points (β = -1) outperform a MC resolution (β = -1
2 ). But to obtain a given accuracy, the number N Q of equidistant points grows exponentially fast with the dimension whereas it does not for MC ones. There exists some intermediary solution and alternative, we refer to [START_REF] Saporta | Probabilités, Analyse de Données et Statistique[END_REF][START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF][START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF][START_REF] Halton | Algorithm 247: Radical-inverse quasi-random point sequence[END_REF][START_REF] Van Der Corput | [END_REF][START_REF] Sobol | Uniformly distributed sequences with an additional uniform property[END_REF][START_REF] Gautschi | Orthogonal polynomials: applications and computation[END_REF][START_REF] Trefethen | Is Gauss Quadrature Better than Clenshaw-Curtis?[END_REF] and the reference therein for the interested reader. (ii) The next step consists in running N independent runs of a black-box code (in this paper it refers to the resolution of ( 1)) at the a priori chosen points (X i , w i ) i∈{1,...,N } and gathering a new collection of output points: (u(x, t, v, X i ), w i ) i∈{1,...,N } . It is supposed to bear the main computational effort as (1) must be solved N times. The N runs are independent and their resolutions can consequently be carried on in parallel (it is often called an embarassingly parallel strategy as there are no communication costs between processes, except during the post-treatment), simultaneously if one has access to as many computation devices than runs (N ). Equation ( 1) is solved thanks to a black-box code 12 up to a certain accuracy depending on its numerical solver. Let us denote by ∆ the discretisation parameter of the latter, then for a first order resolution method with respect to ∆ we have

u(x, t, v, X i ) = u BB (x, t, v, X i ) + K BB (x, t, v, X i )∆ + O(∆ 2 ). ( 11
)
Once again, the notation in [START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF] is very general. For example, if (1) is solved by a MC method, ∆ =

1 √ N M C
where N M C is the number of MC particles 13 and K BB is the standard deviation of the process (error estimator) see [START_REF] Lewis | Computational Methods of Neutron Transport[END_REF][START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF]. If (1) is solved with a deterministic scheme, ∆ = max(∆x, ∆t, ∆v). High-order schemes (see [START_REF] Jourdren | Arbitrary High-Order Schemes for the Linear Advection and Wave Equations: Application to Hydrodynamics and Aeroacoustics[END_REF][START_REF] Gautier Dakin | High-order accurate Lagrangeremap hydrodynamic schemes on staggered Cartesian grids[END_REF] for examples) aim at cancelling K BB and the higher coefficients of powers of ∆. Rigorously speaking, [START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF] comes from some numerical analysis and implies choosing a norm, defining a space for the solution etc. but assuming a general form such as [START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF] will ease the later discussions. (iii) Once the N runs obtained, the rest is only postprocessing at the observation points of interest. The estimation of I(u) is mainly made by numerical integration and we have

I(u) = F (u(x, t, v, X)) dP X , = N i=1 F (u(x, t, v, X i ))w i + O(N β ), = N i=1 F (u BB (x, t, v, X i ))w i + O(∆) + O(N β ), = I ∆ N + O(∆) + O(N β ). ( 12
)
At the end of the process, one has access to an approximation I ∆ N . The error between I and I ∆ N can then be decomposed in two main parts :

||I -I ∆ N || = O(N β ) integration error (UQ) + O(∆) numerical error (BB) . (13)
In [START_REF] Brun | Patmos: A prototype monte carlo transport code to test high performance architectures[END_REF], the error of the non-intrusive approximation have (explicitly) two parameters, N for the integration error which is relative to the resolution of the uncertainty propagation counterpart and ∆ relative to the resolution of each deterministic runs. Basically, if the norm in [START_REF] Brun | Patmos: A prototype monte carlo transport code to test high performance architectures[END_REF] is the L 2 one, quantity [START_REF] Brun | Patmos: A prototype monte carlo transport code to test high performance architectures[END_REF] expresses an error on the variance of the observable I: it is clear that if ∆ N β , the estimated variance of I is closer to a numerical error than a variability due to the uncertain parameters. In other words, to perform an uncertainty analysis, one must make sure ∆ N β . The next figure illustrates this fact. Figure 1 presents some convergence curves in a very simple configuration for two different couples of numerical methods:

∆ = ∆t convergence N = N GL convergence N=2 N=4 N=8
-the first line of figure 1 presents convergence studies obtained with a deterministic scheme of parameter14 ∆ = ∆t for the resolution of the black-box code and N Gauss-Legendre (denoted by N GL ) points for the uncertainty propagation. The top left picture presents a convergence study with respect to ∆ = ∆t for fixed values of N GL = 2, 4, 8. The right one displayes a convergence study with respect to N = N GL for fixed values of ∆ = ∆t. -The second line shows the same convergence studies except the deterministic black-box code is solved by a stochastic (MC) scheme 15 of discretisation parameter ∆ =

1 √ N M C
, the uncertain counterpart being also solved with a MC sampling of (X, dP X ) with N = N M C U Q points. The bottom left picture presents a convergence study with respect to ∆ = N M C for fixed values of N = N M C U Q = 100, 1000. The right one displayes a convergence study with respect to

N = N M C U Q for fixed values of ∆ = N M C .
The quantity of interest here is the variance (see (A.3) in appendix A) of the total amount of particles in a closed box x ∈ D, v ∈ R 3 . Monitoring the error on any different statistical quantity would give similar results at the price of more tedious calculations in order to compute the reference solution: the configuration of interest, described in appendix A, ensures we have access to an analytical solution for V[I](t).

Let us comment on figure 1. Independently of the couple of numerical methods to solve the uncertain and the deterministic counterparts, the behaviours are very similar: every curves present, first, a converging behaviour with a slope characteristic of the numerical method applied O(∆t 1 ) figure 1

(top left), O(N β GL ) figure 1 (top right), O(N -1 2 M C ) figure 1 (bottom left) and O(N - 1
2 ) figure 1 (bottom right). Then, the curves present a kink: it corresponds to the point where the general accuracy becomes driven by the second numerical method. After that kink, the general error stagnates because the overall error is driven by the second discretisation parameter and increasing the one relative to the x-axis does not allow any significant gain. In a sense, the locations of the kinks corresponds to optimal parameter choices (∆, N ): increasing the accuracy in one direction without the other induces a loss of computational time. Looking for this optimal set of parameters for efficiency can be complex and is not the purpose of this paper.

Of course, the application of gPC implies an additional parameter P which is the truncation order of the polynomial approximation

u(x, t, v, X) = ∞ k=0 u X k (x, t, v)φ X k (X) ≈ u P (x, t, v, X) = P k=0 u X k (x, t, v)φ X k (X). ( 14
)
This parameter has been intensionally omitted in the previous studies because it will remain common to the new gPC-i-MC scheme we will describe in the next sections. In [START_REF] Cacucci | Handbook of Nuclear Engineering[END_REF], every gPC coefficients are approximated up to an O(∆) + O(N β(Q) ) accuracy as in [START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF]. In this paper, we intend to show that if the black-box code solves (1) with a MC scheme (i.e. ∆ =

1 √ N M C
) then it can be enriched (intrusively) at a relatively low cost to compute on-the-fly the gPC coefficients (of u as in [START_REF] Cacucci | Handbook of Nuclear Engineering[END_REF] or of any general F as in ( 9)) during the MC resolution. Intuitively, it is easy noticing that when using a non-intrusive gPC (ni-gPC ) method on a MC black-box code, basically, N × N M C MC particles are treated for an overall O( 1

√ N M C
) accuracy (as in the bottom left picture of figure 1). Such tensorisation of the MC particles with the experimental design related to the uncertain parameters can be avoided at the cost of minimal modifications to an existing MC solver. In order to identify accurately those modifications, we suggest briefly and formally recalling how a MC solver for (1) is built in the next section 3 and enrich the MC scheme for uncertainties in section 4.

The Monte-Carlo resolution of the linear Boltzmann equation (1)

In this section, we recall the construction of the semi-analog 16 MC scheme to solve (1) and the general structure of a MC code. This MC strategy is commonly used in neutronic applications and is called implicit capture [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF]. The aim is to very concisely present the different steps of the construction and the algorithmic implications. This way, it will be easier, in the next section, to highlight the differences and commons of the resolution schemes once uncertainties/gPC are introduced. We here describe a backward resolution of the transport equation [START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF][START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF] because it saves some calculations and remains very general. We refer to [START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF] for the reader interested in its forward resolution 17 .

The transport equation with deterministic collisional counterpart is given by

∂ t u(x, t, v) + v • ∇u(x, t, v) + vσ t (x, t, v)u(x, t, v) = vσ s (x, t, v) P s (x, t, v , v)u(x, t, v ) dv . (15)
In [START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF], we introduced

σ s (x, t, v) = σ s (x, t, v, v ) dv , P s (x, t, v, v ) = σ s (x, t, v, v ) σ s (x, t, v) .
Let us now rewrite [START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF] in a recursive integral form: to do so, we apply the characteristic method, multiply (15) by exp( s 0 vσ t (x + vα, α, v) dα) and integrate the expression between 0 and t to obtain

u(x, t, v) = u 0 (x -vt, v) exp - t 0 vσ t (x -v(t -α), α, v) dα + t 0 vσ s (x -v(t -s), s, v)u(x -v(t -s), s, v )e - t s vσt(x-vv(t-α)) dα P s (x -v(t -s), s, v , v) dv ds.
( 16) With ( 16), the transport equation is written in a recursive integral form. Using the fact that exp -

t 0 vσ t (x -v(t -s), s, v) ds = exp - t 0 vσ t (x -vα, t -α, v) dα , = ∞ t 1 [0,∞[ (s)vσ t (x -vs, t -s, v) exp - s 0 vσ t (x -vα, t -α, v) dα ds, = ∞ t f τ (x, t, v, s) ds, ( 17
)
where f τ (x, t, v, s) ds is a probability measure18 , (16) becomes

u(x, t, v) =    +1 [t,∞[ (s)u 0 (x -vt, v) +1 [0,t[ (s) u(x -vs, t -s, v )P s (x -vs, t -s, v , v) σ s (x -vα, t -α, v) σ t (x -vα, t -α, v) dv    f τ (x, t, v, s)ds. (18)
Let's now introduce τ, V sampled from the probability measures τ ∼ f τ (x, t, v, s) ds and V ∼ P s (x, t, v, v ) dv and rewrite the integral equation in term of a recursive expectation over these two random variables

u(x, t, v) = E 1 [t,∞[ (τ )u 0 (x -vt, v) + 1 [0,t] (τ )u(x -vτ, t -τ, V) σ s (x -vτ, t -τ, v) σ t (x -vτ, t -τ, v) . ( 19
)
The next step consists in introducing a MC discretization. Formally, the construction of a MC scheme relies on looking for solutions of ( 19) having the particular forms

u p (x, t, v) = w p (t)δ x (x p (t))δ v (v p (t)). ( 20
)
Such solution u p will be commonly called the 'MC particle' p. The MC scheme intensively uses the linearity of equation ( 15): if (u p ) p∈{1,...,N M C } are independent solutions of ( 15) then

N M C
p=1 u p is also solution of [START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF]. Plugging [START_REF] Després | Robust Uncertainty Propagation in Systems of Conservation Laws with the Entropy Closure Method[END_REF] into [START_REF] Gautier Dakin | High-order accurate Lagrangeremap hydrodynamic schemes on staggered Cartesian grids[END_REF] finally leads to the recursive (backward, see [START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF]) treatment

         w p (t) = 1 [t,∞[ (τ )w p (0) +1 [0,t] (τ ) σ s σ t (x p (t -τ ), t -τ, v p (t -τ ))w p (t -τ ), x p (t) = 1 [t,∞[ (τ )(x 0 + vt) +1 [0,t] (τ )(x t-τ + vτ ), v p (t) = 1 [t,∞[ (τ )(v) +1 [0,t] (τ )(V). (21)
In expression [START_REF] Duderstadt | Transport Theory[END_REF], we recognize the very classical operations one must apply to a MC particle to solve (15):

-the 'census event' corresponds to the condition τ ∈ [t, ∞[: the MC particle is transported along the straight line between x(0) and x(t) = x(0) + vt with no change in its attributes (same weight w p (t) = w p (0), same velocity v p (t) = v p (0) = v). -the 'scattering event' corresponds to the condition τ ∈ [0, t[: the MC particle is transported along a straight line between x(tτ ) and x(tτ ) + vτ together with a modification of its weight according to

w p (t) = σ s σ t (x p (t -τ ), t -τ, v p (t -τ ))w p (t -τ ). ( 22
)
It also implies a change of velocity from v to V at the interaction time tτ and position x t-τ . -Of course, if a grid is introduced, a 'cell exit' event is usually introduced. It is not detailed here because easy to deal with (based on the memorylessness of the exponential distribution for the interaction time, see [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF][START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF]), even with the new MC scheme we present later on. By construction of the MC resolution scheme, theorem 3.2.1 of [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF] ensures the narrow convergence of the MC solver toward the solution of ( 15) in the limit N M C → ∞ for the considered time step [0, t]. The central limit theorem ensures its O( 1

√ N M C
) convergence rate. It does not need a mesh or any tesselation of the space if the cross-sections are analytically known. In practice, only the skins of the materials are projected on a grid and cross-sections are constant per cell and time step. In this case, the operations to perform on a MC particle u p have much friendlier expressions:

-the interaction time is sampled from an exponential law of parameter v p σ t (v p ), i.e. we have

τ = - ln(U τ ) v p σ t (v p ) where U τ ∼ U([0, 1]
) and v p is an embedded particle field.

The latter expression has been obtained inversing the cumulative density function of the forementioned exponential law. This is very classical in MC simulations, see [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF][START_REF] Spanier | Monte Carlo Principles and Neutron Transport Problems[END_REF][START_REF] Lewis | Computational Methods of Neutron Transport[END_REF]]. -The outer velocity V is also obtained locally inversing the cumulative density function of P s (x, t, v, v ) dv .

This means we have

U V = V -∞ P s (x p (t -τ ), t -τ, v p (t -τ ), v ) dv where U V ∼ U([0, 1]), (24)
where x p , v p are embedded particle fields and τ is the current sampled interaction time (the change of velocity occurs only at a collision point/time). In [START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF], we intensionally keep the expression implicit as the inversion of the cumulative density function [START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF] strongly depends on the format of the crosssections (multigroup, continuous [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF][START_REF] Spanier | Monte Carlo Principles and Neutron Transport Problems[END_REF][START_REF] Lewis | Computational Methods of Neutron Transport[END_REF]) but the material of this paper remains independent of such considerations. -The weight modification at the interaction point/time on another hand remains punctual, hence given by ( 22). -At cell interfaces, we rely on the memorylessness of the exponential law to stop the particle and resample the interaction time, see [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF][START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF]. Algorithm 1 in appendix B present the general canvas of a MC resolution code. The latter will be useful to highlight the few modifications needed to take into account uncertainties on-the-fly during the MC computations.

The Monte-Carlo resolution of the uncertain linear Boltzmann equation (4)

As explained earlier, applying ni-gPC to a MC code solving the uncertain linear Boltzmann equation leads to a tensorisation of the N experimental design points with the N M C MC particles. Computational ressources are lost in the sense the overall accuracy remains 19 O( 1

√ N M C
) with a cost being O(N × N M C ). On another hand, MC schemes are computationally insensitive to an increase in dimension. The idea of this section is to explain how one can make the most of the MC resolution of the deterministic linear Boltzmann equation to take into account 3(x) + 1(t) + 3(v) + Q(X) = 7 + Q dimensions and treat, on-the-fly during the MC resolution, the effect of uncertainties on the particle flow. Basically, we would like to keep the accuracy O( 1

√ N M C
) with a cost being O(N M C ) as for any MC method. For this, we must build a new MC scheme. We suggest going through the same steps as in the next section and emphasize its subtleties progressively 20 . The transport equation with uncertain collisional counterpart is given by

∂ t u(x, t, v, X) + v • ∇u(x, t, v, X) = -vσ t (x, t, v, X)u(x, t, v, X) +vσ s (x, t, v, X) P s (x, t, v , v, X)u(x, t, v , X) dv , (25)
with

σ s (x, t, v, X) = σ s (x, t, v, v , X) dv , P s (x, t, v, v , X) = σ s (x, t, v, v , X) σ s (x, t, v, X) . ( 26
)
The idea is to go through the same steps as before but having in mind the quantities depend also on X and identify the changes one must perform to the different samplings to take X into account. For example, we must introduce

f τ (x, t, v, s, X) ds = 1 [0,∞[ (s)vσ t (x -vs, t -s, v, X) exp - s 0 vσ t (x -vα, t -α, v, X) dα ds, (27)
which, under some boundedness conditions21 ∀X ∈ Supp(X) where Supp(X) denotes the support of the random variable, remains an exponential probability measure [START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF]. The uncertain counterpart of ( 18) is then given by

u(x, t, v, X) =    +1 [t,∞[ (s)u 0 (x -vt, v, X) +1 [0,t[ (s) u(x -vs, t -s, v , X)P s (x -vs, t -s, v , v, X) σ s (x -vα, t -α, v, X) σ t (x -vα, t -α, v, X) dv    f τ (x, t, v, s, X)ds.
(28) Introduce the set of random variables τ X , V X sampled from the probability measures τ X ∼ f τ (x, t, v, s, X) ds and V X ∼ P s (x, t, v, v , X) dv , then the above integral equation rewritten as a recursive expectation becomes

u(x, t, v, X) = E 1 [t,∞[ (τ X )u 0 (x -vt, v, X) + 1 [0,t] (τ X )u(x -vτ X , t -τ X , V X , X) σ s (x -vτ X , t -τ X , v, X) σ t (x -vτ X , t -τ X , v, X) . ( 29
)
The next step consists in introducing a MC discretization allowing to take into account the uncertain variables. Let us introduce an 'uncertain MC particle' u p defined as

u p (x, t, v, X) = u p (x, t, v)δ X (X p (t)) = w p (t)δ x (x p (t))δ v (v p (t))δ X (X p (t)). ( 30
)
We are now going to identify the operations we must perform to ensure (30) is solution of [START_REF] Fedorov | Theory Of Optimal Experiments[END_REF]. For this, we plug (30) into (29) and make sure u p (x, t, v, X) is a particular solution of [START_REF] Fedorov | Theory Of Optimal Experiments[END_REF]. Plugging u p into (29) leads to the construction of a (compatible) system of equations of unknowns w p (t), x p (t), v p (t), X p (t) given by

                 w p (t) = 1 [t,∞[ (τ X ) w p (0) +1 [0,t] (τ X ) σ s (x p (t -τ X ), t -τ X , v p (t -τ X ), X p (t -τ X )) σ t (x p (t -τ X ), t -τ X , v p (t -τ X ), X p (t -τ X )) w p (t -τ X ), x p (t) = 1 [t,∞[ (τ X ) (x(0) + vt) +1 [0,t] (τ X )(x p (t -τ X ) + vτ X ), v p (t) = 1 [t,∞[ (τ X ) v +1 [0,t] (τ X )(v p (t -τ X ) = V X ), X p (t) = 1 [t,∞[ (τ X ) X p (0) +1 [0,t] (τ X )(X p (t -τ X )). ( 31
)
Let us focus on the last equation: inconditionally with respect to time t, X p (t) is not modified. Indeed, if τ X < t we have X p (t) = X p (tτ X ) until, events after events, the initial condition is reached leading to X p (t) = X p (0) = X p . Remark 4.1 The latter result tells the uncertain variable must be sampled initially for every MC particles and remain unchanged. It also implies a MC particle must transport amongst its attributes the realisation of a random vector of size Q. This has some impact on the memory consumption of the algorithm. Now we know X p (t) = X p , (31) reduces to

           w p (t) = 1 [t,∞[ (τ Xp ) w p (0) +1 [0,t] (τ Xp ) σ s (x p (t -τ Xp ), t -τ Xp , v p (t -τ Xp ), X p ) σ t (x p (t -τ Xp ), t -τ Xp , v p (t -τ Xp ), X p ) w p (t -τ Xp ), x p (t) = 1 [t,∞[ (τ Xp ) (x(0) + vt) +1 [0,t] (τ Xp )(x p (t -τ Xp ) + vτ Xp ), v p (t) = 1 [t,∞[ (τ Xp ) v +1 [0,t] (τ Xp )(v p (t -τ Xp ) = V Xp ), (32)
where, we recall τ Xp , V Xp are sampled from the probability measures τ Xp ∼ f τ (x, t, v, s, X p ) ds and V Xp ∼ P s (x, t, v, v , X p ) dv . System (32) is very similar to system (21) but the samplings depending on Xp may need few more details: assume, for the sake of simplicity, the cross-sections do not depend on x, t locally22 (i.e. within a cell or an element of geometry). Then the probability measure [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF] for the sampling of the interaction time resumes to

f τ (v p , s, X p ) ds = 1 [0,∞[ (s)vσ t (v p , X p )e -vpσt(vp,Xp)s ds. ( 33
)
In practice, this implies sampling τ Xp according to23

τ Xp = - ln(U) vσ t (v p , X p )
where U ∼ U([0, 1]) and X p is an embedded particle field just as v p .

Expression [START_REF] Isotalo | Preventing Xenon Oscillations in Monte Carlo Burnup Calculations by Enforcing Equilibrium Xenon Distribution[END_REF] echoes [START_REF] Dufek | The stochastic implicit Euler method: A stable coupling scheme for Monte Carlo burnup calculations[END_REF]. The same apply to the sampling of the outer24 velocity V Xp and to the weight modification of the uncertain MC particles: the cross-sections at play in ( 26)-( 32) must be used at both the physical (x p (t), v p (t)) and uncertain (X p ) fields of the uncertain MC particle. Now the gPC coefficients can easily be estimated thanks to the uncertain MC particles: the scheme once again intensively uses the linearity of equation ( 25) together with the linearity of the P -truncated gPC approximation defined via (6)- [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]. Indeed, (u p (x, t, v, X)φ X k (X)) p∈{1,...,N M C } , ∀k ∈ {0, .., P } are independent solutions of the projection of the solution of ( 25) onto a Ptruncated gPC basis. This implies the sum over the number of uncertain MC particles verifies ∀k ∈ {0, .., P }

N M C p=1 u p (x, t, v, X)φ X k (X) ≈ u X k (x, t, v).
Applying the operations related to [START_REF] Hermeline | A discretization of the multigroup pn radiative transfer equation on general meshes[END_REF] to any given uncertain MC particles ensures, by construction (see theorem 3.2.1 of [START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF]), the convergence of the MC solver toward the projection of the solution of (25) onto the truncated gPC basis in the limit N M C → ∞. The overall cost remains O(N M C ) together with an O( 1

√ N M C
) accuracy on the gPC coefficients to compute. Note that the computation of the gPC coefficients explicitly appears in the algorithmic presentation 2 of the new MC scheme in appendix C together with an exhaustive description and a discussion on parallel strategies. From the previous description, one must understand the basic idea is to try to avoid a tensorisation between the N experimental design points concerning the random variable X and the N M C samplings related to the physical variables (x, t, v) for the MC particles. This imposes some identified operations to perform a full MC approximation of the gPC coefficient with N M C sampling in the whole space of variables (x, t, v, X). We intensively make use of the insensitiveness of a MC integration with respect to dimension to compute the gPC coefficient of any given output of interest. The new MC scheme is intrusive in the sense one must modify 25-the attributes of the MC particles to take into account a discretisation (X p ) p∈{1,...,N M C } of (X, dP X ), -the call to the cross-sections at those points (X p ) p∈{1,...,N M C } to sample the interaction time, the outer velocity and modify the weight of any uncertain MC particle, -the tallies (to embed the computations of the gPC coefficients and other outputs of interest). The last point may deserve few more details: to approximate any output of interest of the form (5), the post-treatment must also be embedded in the MC resolution. Any other quantity of interest will not directly be available unless every fields of the uncertain MC particles are tracked in some files to be post-treated. We clearly want to avoid such solution because tracking down information with such frequency (many tallies 26 of MC particles per seconds leading to a very important volume of I/O27 ) slows drastically down the computations and can easily make a filesystem collapse. More details will be given in the numerical examples of the next section and in the description of algorithm 2.

At this stage of the discussion, one may also wonder why relying on gPC and consequently remaining sensitive to the dimension Q via the increasing number 28 of coefficient (u X k ) k∈{0,..,P } to be evaluated. To give an element of answer, let us build the PDE satisfied by the moment of order 2 of u, solution of (4). It is defined by

M 2 (x, t, v) = u 2 (x, t, v, X) dP X , = m 2 (x, t, v, X) dP X ,
and certainly corresponds to one of the simplest statistical observable. In this case, quantity m 2 is solution 29of

∂ t m 2 (x, t, v, X) + v • ∇m 2 (x, t, v, X) = -2vσ t (x, t, v, X)m 2 (x, t, v, X) +2u(x, t, v, X) vσ s (x, t, v, v , X)u(x, t, v , X) dv . ( 35
)
The latter equation is nonlinear (see the scattering term). The difficulty to solve [START_REF] Jourdren | Arbitrary High-Order Schemes for the Linear Advection and Wave Equations: Application to Hydrodynamics and Aeroacoustics[END_REF] with a MC method can be compared to the one to solve the quadratic Boltzmann equation [START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF][START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] for example. In other words, to be solved numerically, it may -either need an additional linearisation hypothesis. For example, for a Nanbu-like [START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF] resolution, this implies relying on a time step discretisation and a MC resolution of the explicited equation

∂ t m 2 (x, t, v, X) + v • ∇m 2 (x, t, v, X) = -2vσ t (x, t, v, X)m 2 (x, t, v, X) +2u(x, t n , v, X) vσ s (x, t, v, v , X)u(x, t, v , X) dv , (36)
where u(x, t n , v, X) is the approximated solution at the beginning of the time step ∆t. In other words the convergence depends on N M C but also on the time step ∆t.

-or perform a splitting of operator between the streaming part and the collisional one with an adaptation of Bird's algorithm [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]. This splitting, even if having very good mathematical properties (conservations), also introduces a dependence with respect to a time step ∆t. -or apply an analog MC scheme and keep track of the count rate to take correlations and higher moments into account [START_REF] Cacucci | Handbook of Nuclear Engineering[END_REF]. This is usually done in a file which must be post-treated (binning and linear fit etc. see [START_REF] Cacucci | Handbook of Nuclear Engineering[END_REF]). But analog schemes are known to have a slower convergence rate, to be computationally intensive and inadapted to very multiplicative media (in this case the size of the written file is known to explode). Of course, the above list of alternatives may not be exhaustive. Anyway, the application of gPC does introduce a new parameter (P rather than a time step as in [START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF] or [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]) but the modifications to an existing solver are minor (compare algorithms 1 and 2) and the approximation with respect to P can even be expected to yield spectral convergence for smooth solutions [START_REF] Cameron | The Orthogonal Development of Non-Linear Functionals in Series of Fourier-Hermite Functionals[END_REF]. In the example above, a gPC-i-MC approximation of M 2 is simply given by

M 2 (x, t, v) = ∞ k=0 (u X k (x, t, v)) 2 ≈ P k=0 (u X k,N M C (x, t, v)) 2 ,
where ∀k ∈ {0, .., P } we have

u X k (x, t, v) ≈ u X k,N M C (x, t, v) = N M C p=1 w p (t)δ x (x p (t))δ v (v p (t))φ X k (X p ).
In the following section, we numerically verify and illustrate the previous points and even consider more elaborated statistical outputs of interest (in particular Sobol indices for sensitivity analysis, see [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF][START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF]).

Numerical examples and discussions

In this section, we go through several test-problems emphasizing the strengths and weaknesses of the new MC scheme we suggested in the previous section. We first briefly go back to the simple problem of section 2 which motivated the introduction of the new approach. We then consider three spatially dependent problems for which analytical solutions are not anymore available (to our knowledge). The test-cases are progressive in difficulty and relevance: first we consider a monodimensional uncertain problem (i.e. Q = 1) with simple statistical outputs of interest. Second, we perform a sensitivity analysis with respect to three parameters (Q = 3), implying the approximations of spatial Sobol indices [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF][START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF]. Finally, we tackle a last sensitivity analysis with respect to six parameters (Q = 6) for spatial particle density profiles in an uncertain two-layer medium. For the two first test-problems, we rely on a ni-gPC application, as described in section 2, to produce reference solutions. For the last, the classical ni-gPC method is too costly (this will be justified in section 5.4). To produce the reference solutions (in sections 5.2-5.3), Gauss-Legendre points are used (N = N GL ) for their accuracy 30 and efficiency in relatively low stochastic dimensions (here, up to Q = 3). The test-cases are all monokinetic for the sake of simplicity and ease of the reproducibility of the results. Remark 5.1 The reader may notice the variability of the uncertain input parameters in the benchmarks is always relatively important: accordingly, we make sure the numerical discretisation of the MC simulation code (i.e. the black-box) is not too constraining (i.e. O( 1

√ N M C ) O(N β (Q))
, recall the example of section 2) and that the reference solutions obtained with ni-gPC can be produced in reasonable time. We insist this does not affect the relevance of the discussion. 30 For such integration method, we have, for any continuous functions g g

(x) dP X - N i=1 g(X i )w i = K 2N ! g (2N ) (ξ) = N β ,
in which g (n) denotes the n th derivative of g , see [START_REF] Gautschi | Orthogonal polynomials: applications and computation[END_REF][START_REF] Akhiezer | The Classical Moment Problem[END_REF].

Back to the simple configuration of figure 1

In this section, we briefly go back to the simple problem tackled in section 2. In figure 2, we display the same curves as in figure 1, obtained with a non-intrusive strategy (N = N U Q M C and ∆ = N M C ), together with the one obtained with the gPC-i-MC scheme we described in the previous section. new scheme is such that N U Q M C = N M C : the experimental design is not anymore tensorized with the MC particles and the methodology consequently has one less numerical parameter. Besides, as displayed in figure 2, the approximation obtained with the new MC scheme does not stagnate with the increasing number of samplings. The uncertainty is solved on-the-fly during the MC resolution and the convergence rate for the whole problem remains O( 1

√ N M C
) avoiding the kinks in the curves obtained non-intrusively.

Transport in an uncertain diffusive material

Let us now tackle a new test-problem for which an analytical solution is not available despite the relative simplicity of the configuration. Let us consider x ∈ D = [0, 2], v = 1. Furthermore, we assume the scattering is isotropic (notation dω = 1 S 2 (ω) 1 4π dω = 1) and that the medium is only diffusive (no absorption, i.e. σ t (X) = σ s (X)) even if uncertain. The initial condition is a Dirac mass at x = 1. In this particular case, (4) resumes to

     ∂ t u(x, t, ω, X) + vω∇ x u(x, t, ω, X) = -vσ s (X)u(x, t, ω, X) + vσ s (X)u(x, t, ω , X) dω , u(x, 0, v) = u 0 (x) = δ 1 (x). ( 37
)
Let us consider a monodimensional uncertain parameter (i.e. Q = 1) and assume X ∼ U([-1, 1]) with σ s (X) = σ s + σs X with σ s = 1 and σs = 0.99. The variability is important in this example, see remark 5.1. The fact the medium is only diffusive typically allows focusing on the difficulty tackled in the previous section relative to moment m 2 of u, see [START_REF] Jourdren | Arbitrary High-Order Schemes for the Linear Advection and Wave Equations: Application to Hydrodynamics and Aeroacoustics[END_REF].

Realisations of U (x, t = 0.5, X)

Mean and Variance

u(x, t = 0.5, X1) u(x, t = 0.5, X2) u(x, t = 0.5, X3) u(x, t = 0.5, X4) Fig. 3. Left: four realisations (taken at four non-intrusively obtained Gauss-Legendre points) of U (x, t = 0.5, X). Right: corresponding mean and variance of U (x, t = 0.5, X). The scale for the mean is on the left, for the variance on the right.

0 0.
Let us now comment on the results of figures 3-4: figure 3 (left) presents four non-intrusively obtained realisations of the spatial profile U (x, t = 0.5, X) = u(x, t = 0.5, ω, X) dω. The particles propagate toward the left and right boundaries of the domain, hence the more or less steep fronts of the density U depending on the value of the uncertain parameter X. For X = X 1 , the medium behaves almost as a vacuum whereas for X = X 4 it is very diffusive. From now on we focus on statistical observables such as the mean E[U ] and variance V[U ] of U (x, t = 0.5, X) with respect to x at time t = 0.5. Note that every spatial statistical observables of this section are computed via the approximations of the gPC coefficients

U X k (x, t) = u X k (x, t, ω) dω = u(x, t, ω, X)φ X k (X) dP X dω = U (x, t, X) dP X , (38)
and some of their post-treatments. For example, (P -truncated) approximations of the mean and variance are given by

E[U ](x, t) = U X 0 (x, t), V[U ](x, t) = P k=1 (U X k (x, t)) 2 . ( 39
)
Many other classical statistical quantities can be obtained from post-treatments of the gPC coefficients, see [START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF]. Some examples will be given in the following. 2s. The new gPC-i-MC method recovers the same results as the best compromise solution with only one run of the MC simulation device and very similar computational times (at least for this problem, this will not exactly be the case for the next benchmarks): it testifies of the insensitiveness to the number of dimensions of the tracking32 of the MC particles. The runs were performed on N replication = 32 replicated33 domains for the parallel strategy, in very similar conditions as in [START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF]. The latter remark allows insisting on the fact the parallel strategies applying for a classical MC simulation device also applies with the new gPC-i-MC one. The minor code modifications described in sections 4-C to implement the gPC-i-MC solver do not imply a porting of the parallel counterpart, it is straightforward if already developed. A short discussion on parallel possibilities is provided in appendix C with the description of the new algorithm. The gain in computational time, for this test-case, is approximatively of N GL ≈ 4. In the above example, it remains relatively low: if one has access to 4(N GL ) × 32(N replications ) = 128 processors, which is common nowadays, the computational times are equivalent and the gain is only in term of computational ressources 34 . We here recall that MC simulation codes are known to be computationally intensive and even such low factor (gain of only 4) may be welcome. In the following, we tackle some multidimensional uncertain problems.

E[U ](x, t = 0.5) V[U ](x, t = 0.5) ni-gPC (reference) ni-gPC (best compromise) MC-i-gPC 0 0.5 1

Sensitivity analysis in 3D stochastic dimension

This new example is a 3-dimensional stochastic (i.e. Q = 3) test-problem for which a reference solution with ni-gPC can still be obtained in reasonable times. The set-up is as follow:

-

v = 1, x ∈ D = [0, 1], subdivided into N x = 100 cells ∪ Nx i=1 D i = D.
-Specular boundary condition on left (at x = 0) and vacuum one on the right hand side (at x = 1).

-Initially, the density of particles is homogeneous and deterministic, equal to 1, i.e. u(x, t = 0, ω, X) = u 0 (x, ω, X) = 1 ∀x ∈ D, ∀ω ∈ S 2 . -The medium is pure (i.e. M = 1 and η = η 1 see ( 2)), homogeneous and considered uncertain. It depends on three parameters X = (X 1 , X 2 , X 3 ) affecting the total and scattering cross-sections and the material density as

σ t (x, t, X) = σ t (X 1 ) = σ t + σ t X 1 , ∀x ∈ D, t ∈ R + , σ s (x, t, ω, ω , X) = σ s (X 2 ) = σ s + σ s X 2 , ∀x ∈ D, t ∈ R + , ∀(ω, ω ) ∈ S 2 , η(x, t, X) = η(X 3 ) = η + ηX 3 , ∀x ∈ D, t ∈ R + , (40)
in which (X 1 , X 2 , X 3 ) are independent uniformly distributed random variables on [-1, 1], i.e. ∀i ∈ {1, 2, 3},

X i ∼ U([-1, 1]).
-For the next computations, the mean quantities are set to σ t = 1.0, σ s = 0.9, η = 1.0 and the ones controling the variability to σ t = 0.4, σ s = 0.4, η = 0.4. Note that remark 5.1 also applies here.

-We are interested in the mean E[U ], variance V[U ] and Sobol indices S tot [U ] profiles of U (x, t, X) = u(x, t, ω, X) dω at time t = 1.0. The total and first order Sobol indices [START_REF] Sobol | Uniformly distributed sequences with an additional uniform property[END_REF][START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF] relative to output U are denoted by

S tot [U ] = (S tot 1 [U ], ..., S tot Q [U ]) t and S 1 [U ] = (S 1 1 [U ], ..., S 1 Q [U ]) t .
They represent a powerful but costly (see [START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF]) statistical tools designed to identify, for a given output of interest, which of the uncertain parameters explain the most its variability. They can also be accurately approximated via post-treatments of the gPC coefficients, see [START_REF] Blatman | Sparse Polynomial Chaos Expansions and Adaptive Stochastic Finite Elements using a Regression Approach[END_REF][START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF]. x Fig. 5. Comparison between ni-gPC and gPC-i-MC . Top picture: mean of U (x, t = 1, X) with respect to x. Bottom left: variance of U (x, t = 1, X) with respect to x. Bottom right: total Sobol indices of U (x, t = 1, X) with respect to x.

E[U ](x, t =
Figure 5 compares results obtained with ni-gPC and the gPC-i-MC scheme for the forementioned statistical outputs of interest. Before going through resolution strategy comparisons, let us present briefly the results: figure 5 (top) shows the mean of U (x, t = 1, X). Particles are globally absorbed in the vicinity of x = 0: indeed, we initially have U (x, t = 0, X) = 1, ∀x ∈ D and on figure 5 (top), we have E[U ](x ∼ 0, t = 1) < E[U ](x ∼ 0, t = 0) = 1. This averaged particle absorption occurs despite the probable multiplicative effect (σ s (X) > σ t (X) for some realisations of X) of the medium. Particles are globally lost in the vicinity of x = 1, as

E[U ](x ∼ 1, t = 1) < E[U ](x ∼ 1, t = 0) = 1,
mainly due to the vacuum boundary conditions. Figure 5 (bottom-left) shows the variance of U . The uncertainty is more important in the vicinity of x = 0 and drops of a decade between x = 0 and x = 1. Now we are interested in identifying which of the uncertain parameters explain the most the previous variability: the total Sobol indices for X = (X 1 , X 2 , X 3 ) are displayed figure 5 (bottom-right). Globally, X 3 35 is the lesser important parameter: its total Sobol indice is the lowest at every spatial location x ∈ D of the simulation domain. Parameters X 1 and X 2 , impacting respectively the scattering and total cross-sections, have a globally equivalent influence. The uncertainty on the scattering cross-section (X 2 ) is in particular more important than the one on the total cross-section (X 1 ) in the vicinity of the vacuum boundary condition x = 1. Figure 6 presents a comparison of total (S tot [U ]) and first order (S 1 [U ]) Sobol indices: we recall (see [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]) that the total indice S tot i [U ] relative to X i takes into account the first order indice S 1 i [U ] of X i together with its interaction with every other variables. In particular, in figure 6 (bottom-right), we see that X 3 is important mainly through its interactions 36 with X 1 and X 2 . Hence, even if every of the uncertain parameters are not negligible (relatively important total Sobol indices on figure 6 bottom-right), figure 6 bottom-right attests reducing the uncertainty on X 1 , X 2 will also lead to a reduction of the uncertainty due to X 3 . Going further in interpretation would need the careful study of second order Sobol indices but this is not really the scope of this section. The next one will provide a more pedagocial example. We here wanted to put forward the new gPC-i-MC scheme is able to accurately recover statistical quantities known to be very efficient but costly [START_REF] Bobylev | Theory of collision algorithms for gases and plasmas based on the boltzmann equation and the landau-fokker-planck equation[END_REF].

S tot 1 [U ](x, t = 1) vs. S 1 1 [U ](x, t =
A quick look to figure 5 showed the three resolutions are in agreement, of equivalent accuracies, whatever the statistical observable of interest (mean, variance or Sobol indices). Every resolutions have N M C = 3.2 × 10 7 particles. Let us now focus on their differences: 35 Related to an uncertainty on η. 36 Important difference between S gPC-i-MC with and (P + 1) Q = (2 + 1) 3 = 27 coefficients. First, once again two different levels of discretisation for ni-gPC give equivalent results. The good agreement between the latters testifies the uncertainty quantification counterpart of the decoupled approach is converged. Once again, the fast convergence rate of gPC is put forward: we obtain accurate solutions for low (P = 2) polynomial orders in every directions. The gPC-i-MC scheme gives equivalent results with one run and (P + 1) Q = (2 + 1) 3 = 27 gPC coefficients: it consequently takes advantage of the fast convergence rate of gPC.

Second, let us discuss the average CPU times and costs of the two methods for equivalent accuracies: For the previous test-case (section 5.2), the overall cost of one gPC-i-MC run was very similar to the average CPU time of the N GL ni-gPC ones. For the one of this section, the CPU times present a signicative difference deserving a careful study: one gPC-i-MC run costs about ≈ 1.26 times the average CPU time of the ni-gPC ones. The main difference with the case of section 5.2 is the number of coefficients to be computed: (P + 1) Q = (2 + 1) 3 = 27 here instead of (P + 1) Q = (2 + 1) 1 = 3 for the example of section 5.2. This increase affects two phases of the presented computations:

-the size of the parallel reduction/communication between the N replication replicated domains, -the number of tallies 37 one MC particle must perform. Table 1 compares the CPU times of sequential 38 and parallel runs in comparable conditions (same N M C ). A horizontal reading of table 1 gives an idea of the cost of the reduction phase between the N replication = 8 replicated domains whereas a vertical one gives information on the cost of increasing the number of tallies the uncertain MC particles must perform to apply gPC-i-MC . The main increase in computational time

gPC-i-MC N replication = 1 N replication = 8
t CPU for (P + 1) Q = (0 + 1) comes from the tallying and not the parallel reduction. The tallying phase is more sensitive to an increase of dimension Q. Nonetheless, using N replication = 40 replicated domains with less (N M C = 0.8 × 10 7 ) uncertain MC particles instead of N replication = 32 and N M C = 3.2 × 10 7 ensures recovering the same results as in figure 5 with similar restitution times: in this case, for equivalent accuracies and CPU times, the gain in computational ressources is of a factor

N GL ×32 40 = 64×32 40 = 51.2.
As a concluding remark of this section, we would like to emphasize the numerical accuracy of the gPC-i-MC scheme is still O( 1

√ N M C
) but the cost of the treatment of one uncertain MC particle is more important 37 See algorithm 1 for the definition of the tallying phase. 38 In this case there are no communication cost.

than for a classical one, especially as (P + 1) Q increases. We only briefly tackled domain replication as a (distributed [START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF][START_REF] Palais | Load balancing of parallel Monte Carlo transport calculations[END_REF][START_REF] Brunner | An efficient, robust, domain-decomposition algorithm for particle Monte Carlo[END_REF][START_REF]Monte carlo photon transport on shared memory and distributed memory parallel processors[END_REF][START_REF]Parallel performance study of monte carlo photon transport code on shared-, distributed-, and distributed-shared-memory architectures[END_REF]) parallel strategy but shared memory parallel ones (threads [START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF][START_REF] Palais | Load balancing of parallel Monte Carlo transport calculations[END_REF][START_REF] Brunner | An efficient, robust, domain-decomposition algorithm for particle Monte Carlo[END_REF][START_REF]Monte carlo photon transport on shared memory and distributed memory parallel processors[END_REF][START_REF]Parallel performance study of monte carlo photon transport code on shared-, distributed-, and distributed-shared-memory architectures[END_REF], GP-GPU, vectorisation) could be applied, for example, to accelerate the tallying phase. The idea is in the same vein as the one depicted in [START_REF] Brun | Patmos: A prototype monte carlo transport code to test high performance architectures[END_REF][START_REF] Wang | Optimization of Monte Carlo Neutron Transport Simulations with Emerging Architectures[END_REF] for the on-the-fly Doppler broadening with the SIGMA1 algorithm to differently make the most of new computer architectures. Of course, we do not expect the efficiency for the tallying phase to reach the one of the SIGMA1 algorithm as in [START_REF] Brun | Patmos: A prototype monte carlo transport code to test high performance architectures[END_REF] but it opens to new possibilities.

From now on we will let this aspect (finer parallel strategies) of the discussion for future papers.

A two layered uncertain material

In this section we tackle a last test-problem in 6-stochastic dimensions for which we consider applying ni-gPC is not possible in reasonable times: assuming the same convergence properties as for the other testcases applies (i.e. need for at least N GL = 4 Gauss-Legendre points and P = 2 per directions), we would need N Q GL = 4 6 = 4096 runs a MC black-box code (i.e. N GL × N replication = 4096 × 32 = 131072 processors) to compute (P + 1) Q = 3 6 = 729 gPC coefficients and solve the problem. On another hand, gPC-i-MC can handle it efficiently as presented below.

The general set-up is as follow:

-v = 1, x ∈ D = [0, 1], subdivided into N x = 100 cells ∪ Nx i=1 D i = D.
-Specular boundary condition on the left (at x = 0) and vacuum one on the right hand side (at x = 1). -Homogeneous and deterministic particle density u(x, t = 0, ω, X) = u 0 (x, ω, X) = 1 ∀x ∈ D, ∀ω ∈ S 2 . -The material is composed of two layers of different media, A and B with D A = [0, 1 2 ] and

D B = [ 1 2 , 1] such that D A ∪D B = D = [0, 1]
. Both media are pure (i.e. M = 1 and η = η 1 see ( 2)), homogeneous and considered uncertain. Each depends on three parameters (X i ) i∈{A,B} = (X i 1 , X i 2 , X i 3 ) i∈{A,B} affecting the total and scattering cross-sections and the material density as in the test-problem of section 5.3. We have ∀i ∈ {A, B}

σ t (x, t, X) = i∈{A,B} σ t (X i 1 )1 Di (x) = i∈{A,B} σ i t + σ i t X i 1 1 Di (x), ∀x ∈ D, t ∈ R + , σ s (x, t, ω, ω , X) = i∈{A,B} σ s (X i 2 )1 Di (x) = i∈{A,B} σ i s + σ i s X i 2 1 Di (x), ∀x ∈ D, t ∈ R + , η(x, t, X) = i∈{A,B} η(X i 3 )1 Di (x) = i∈{A,B} η i + η i X i 3 1 Di (x), ∀x ∈ D, t ∈ R + , (41)
in which (X i 1 , X i 2 , X i 3 ) i∈{A,B} are independent uniformly distributed random variables on [-1, 1], i.e. ∀i ∈ {A, B}, j ∈ {1, 2, 3}, X i j ∼ U([-1, 1]). -For the next computations, the mean quantities are set to

σ A t = 1.0, σ A s = 1.3, η A = 1.0, σ B t = 1
.0, σ B s = 0.9, η B = 1.0, and the ones controling the variability to

σ A t = 0.4, σ A s = 0.4, η A = 0.4, σ B t = 0.4, σ B s = 0.4, η B = 0.
4, Note that remark 5.1 also applies here. In the following, we aim at answering several questions: which of the 6 parameters explain the most the variability of the particle density at some prescribed locations? Is it possible to reduce the dimensionality of the problem? By which factor the uncertainty on the main parameters should be reduced to make the remaining ones equally influent?

Figure 7 presents the results obtained with the gPC-i-MC scheme with (P + 1) Q = (2 + 1) 6 = 729 gPC coefficients estimated thanks to N M C = 1.024 × 10 9 uncertain MC particles. The computation has been performed on N replication = 1024 processors and took about 750s. Figure 7 variance profiles of U (x, t = 1, X) = u(x, t = 1, ω, X) dω. The scale for the mean is on the left hand side, for the variance on the right hand side. The profiles present two regimes, on each side of x = 0.5. Material A, on the left hand side, is the most uncertain: the variability drastically drops between x = 0 and x = 1 2 . On another hand, even if smaller in B, the variability remains almost constant on the whole range x ∈ [ 1 2 , 1]. Figure 7 top-right shows the total Sobol indices' profiles at t = 1. The contribution of X 6 to the global variability of the particle density is very small ∀x ∈ D: the dimension of the stochastic problem can consequently be reduced by assuming X 6 ≈ E[X 6 ]. The six other pictures of figure 7 compare the total and first order Sobol indices. They confirm the variability of X 6 can be neglected without impact on the problem together with the fact every other parameters are influent: except for X 3 , they do not depend too much on interactions as the total and first order indices are almost equal. For the next studies, we can consequently keep Q = 5 and choose X 6 = E[X 6 ] = 0. The question arising now is: of how much should we reduce the uncertainty on the main parameters X 1 , X 2 , X 3 (which according to the previous study are the most influent parameters) to have every remaining inputs contributing equally to the global variability? Parameter X 3 being in interaction with other parameters, we can expect its variability to drop by only reducing the one on X 1 , X 2 . We consequently ran several similar calculations with decreasing values of σA t , σA s (controling the variabilities of X 1 , X 2 ) until every total Sobol indices were of the same magnitude (as the ones of X 4 , X 5 ). We performed a dozen calculations and realize having σA t = 0.15, σA s = 0.15 leads to satisfying results. They are displayed in figure 8: it compares the mean, variance and total Sobol indices in the nominal configuration (σ A t = 0.40, σA s = 0.40) and in the newly studied one (σ A t = 0.15, σA s = 0.15). First, figure 8 (left) compares the means and variances in the nominal and new configurations. The mean is slighly impacted but if we focus on the variances per parameters (sobol indices multiplied by the total variance), the study shows that a small decrease (from 0.40 to 0.15) of the variabilities of the identified most influent parameters (X 1 , X 2 ) leads to an important one on the outputs of interest. Besides, on figure 8 (right) as expected, -the magnitude of the total Sobol indices of X 4 , X 5 has not changed (remember their interactions with other variables are negligible). -reducing the uncertainty of X 1 and X 2 has an impact on X 3 : remember we have an important interaction counterpart for this variable, see the big difference between S 1 3 and S tot 3 in figure 7 (bottom left). As a result, its new total Sobol indice is now even below the ones of X 4 , X 5 .

E[U ](x, t = 1) and V[U ](x, t = 1) σA t = 0.
-The total Sobol indices for X 1 , X 2 are now of the same order as the one for X 3 , X 4 .

With this example, we wanted to emphasize the strength of Sobol indices to conduct a sensitivity analysis study. They allow focusing on the most relevant parameters and defining minimal objectives on which the efforts must be carried on to reduce the overall uncertainty. In the previous example, by focusing on X 1 , X 2 and reducing their uncertainty of only a factor 0.15 0.40 = 0.375, it allowed reducing the effect of X 3 of a factor ≈ 0.038 0.01 = 3.80 together with the overall uncertainty on the output of a factor ≈ 0.245 0.025 = 9.80. Working more on X 1 , X 2 would not be efficient as the overall accuracy would remain governed by the uncertainty on X 4 , X 5 . The Sobol indices allowed us identifying the most influent variables, the gPC-i-MC scheme allowed performing the study in reasonable times/with reasonable ressources.

Conclusion

To close this paper, we would like to come back on its main aspects and on future work for the gPC-i-MC solver. This paper introduces a new gPC based intrusive Monte-Carlo method for uncertainty quantification. Combining both intrusiveness for UQ and MC methods is singular and original. The manuscript explains how one can easily modify (hence the intrusiveness) an existing MC simulation device to take into account uncertainties in the linear Boltzmann equation on-the-fly during the MC resolution.

The efficiency of the new gPC based MC scheme is closely related to the insensitiveness of an increase of dimension of the tracking of the MC particles together with the fast convergence rate of gPC. The new method remains sensitive to the curse of dimensionality (as it depends on gPC) mainly during the tallying phase but shows important gains in comparison to non-intrusive gPC. We detailed the minor modifications to apply to the semi-analog/implicit capture MC scheme to take uncertainties into account on-the-fly during the MC resolution. The non-analog, analog and any other MC schemes can be modified the same way without further difficulties applying the material of this paper. In other words, the methodology is not scheme dependent: it also means it can be combined to variance reduction technics without major difficulties.

The next steps consist in considering uncertain nonlinear couplings involving the linear Boltzmann equation: Bateman equations for neutronics [START_REF] Dufek | Stochastic approximation for Monte-Carlo Calculation of Steady-State Conditions in Thermal Reactors[END_REF][START_REF] Isotalo | Computational Methods for Burnup Calculations with Monte Carlo Neutronics[END_REF][START_REF] Isotalo | Preventing Xenon Oscillations in Monte Carlo Burnup Calculations by Enforcing Equilibrium Xenon Distribution[END_REF][START_REF] Dufek | The stochastic implicit Euler method: A stable coupling scheme for Monte Carlo burnup calculations[END_REF][START_REF] Bernede | An Unsplit Monte-Carlo scheme for the resolution of the linear Boltzmann equation coupled to (stiff) Bateman equations[END_REF]], Stefan's law for photonics [START_REF] Lowrie | The coupling of radiation and hydrodynamics[END_REF] or even hydrodynamics with the resolution of the BGK [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF] model, solved with the new gPC-i-MC scheme. The efficiency of the new method is closely related to the unconditional stability of the MC scheme for the linear Boltzmann equation (less frequent parallel reductions). The stakes, for the forementioned couplings/nonlinear models, will reside in the capability for the solvers to take large times steps during the MC phase, see [START_REF] Bernede | An Unsplit Monte-Carlo scheme for the resolution of the linear Boltzmann equation coupled to (stiff) Bateman equations[END_REF] for example. Such applications will certainly be the subject of future papers.

#δt is the time remaining to reach census The gPC-intrusive MC semi-analog scheme described in term of algorithmic operations in order to compute (backward) the gPC coefficients of u(x, t, v, X). The differences with the classical canvas are highlighted in blue: they concern the initial sampling (X p from dP X ), the tallying (estimated the gPC coefficients (u X k (x, t, v)) k∈{0,..,P } intrusively) and the different calls to cross-sections.

Let us briefly tackle parallelism considerations. We insist the operations to perform to one uncertain MC particle remain independent of the others 39 . Consequently, in term of parallel modification, the canvas between algorithm 1 and algorithm 2 does not change. The same parallel strategies apply and the parallel reduction/communications/synchronizations are at the same places within the code. There are still few differences (some are illustrated in the benchmarks of section 5):

-Distributed parallel strategies (such as domain replication or domain decomposition, see [START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF][START_REF] Palais | Load balancing of parallel Monte Carlo transport calculations[END_REF][START_REF] Brunner | An efficient, robust, domain-decomposition algorithm for particle Monte Carlo[END_REF][START_REF]Monte carlo photon transport on shared memory and distributed memory parallel processors[END_REF][START_REF]Parallel performance study of monte carlo photon transport code on shared-, distributed-, and distributed-shared-memory architectures[END_REF]): domain replication has been discussed through the examples of section 5: algorithm 2 does imply a reduction on more arrays (P + 1 quantities to average over the number of replications N replication ). The size of the communication is consequently more important (see table 1). Fortunately, the MC scheme being inconditionally stable, those reduction can occur only at the observation locations/times of interest and are consequently scarse or of reasonable volume during the computation. Domain decomposition has not been tackled in this document. It is not the most common parallel strategy for MC simulation codes. It consists in decomposing the simulation domain D into N decomposition ones and applying the numerical scheme on every one of them. Domain decomposition has not been implemented with gPC-i-MC but should not induce more difficulties than for the classical treatment: particles would still be buffered and communicated between subdomains, see [START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF]. Now, each uncertain particle carries an additional field of size Q, hence a bigger volume to be communicated. -Shared memory strategies have only been hinted at in section 5.3 but surely present an interest, especially vectorisation strategies for the tallying phase which now bears the main computational load (especially as P or Q grow, see section 5.3). Once again, the idea is in the same vein as the one in [START_REF] Brun | Patmos: A prototype monte carlo transport code to test high performance architectures[END_REF] for the vectorisation of on-the-fly Doppler broadening: the gPC basis being orthonormal, there exists a three-term recurrence formulae [START_REF] Szego | Orthogonal Polynomials[END_REF][START_REF] Gautschi | Orthogonal polynomials: applications and computation[END_REF] allowing to express φ X k recursively as a function of φ X k-1 , φ X k-2 . Denote by Φ X P = (φ X 0 , ..., φ X P ) t the vector of P + 1 components of the sequence of orthonormal polynomials associated to X, this ensures there exists a matrix J X P (Jacobi matrix, see [START_REF] Gautschi | Orthogonal polynomials: applications and computation[END_REF]) such that β P φ X P +1 (X)e P = J X P Φ X P (X) -XΦ X P (X), (C.1)

where e P = (0, ..., 0, 1) t of size P + 1 and The latter matrix only depends on the uncertain input distributions and would not change during the whole computation. We here claim no originality, such matrix-vector products are at the basis of Chebyshev (or modified Chebyshev) algorithms [START_REF] Gautschi | Orthogonal polynomials: applications and computation[END_REF] which can take advantage of vectorisation strategies. The gPC basis could then also be computed more efficiently on-the-fly for every uncertain MC particles and the cost of the gPC-i-MC scheme could be reduced. This subject will be the purpose of future studies.

J X P =            α 1 √ β 1 0 0 ... 0 √ β 1 α 2 √ β 2

Fig. 1 .

 1 Fig. 1. Convergence studies with respect to N and ∆ for two couples of numerical methods for the resolution of the uncertain linear Boltzmann equation in a homogeneous configuration.

Fig. 2 .

 2 Fig. 2. Convergence studies with respect to N = N U Q M C and ∆ = N M C as in figure 1 (bottom right) together with a curve obtained with the new gPC-i-MC scheme. For the latter, the uncertain parameters are sampled within the N M C MC particles.

Fig. 4 .

 4 Fig. 4. Comparison between two discretisation of ni-gPC (N GL = 50, P = 6 and N GL = 4, P = 2) and gPC-i-MC (P = 2). Every computations have N M C = 3.2 × 10 8 particles. The results are in perfect agreement. Left: mean of U (x, t = 0.5, X). Right: variance of U (x, t = 0.5, X).

Figure 3

 3 Figure 3 presents the mean and variance profiles of U obtained with ni-gPC . The uncertainty transmitted from the medium to the particles is especially strong on the steep propagation fronts and at x = 1 (i.e. vicinity of the initial condition δ 1 (x)). Now, figure 4 compares the results in term of mean and variance profiles from -(reference) ni-gPC with N M C = 3.2 × 10 8 particles, N GL = 50 and P = 6, -(best compromise) ni-gPC with N M C = 3.2 × 10 8 particles, N GL = 4 and P = 2, -gPC-i-MC with N M C = 3.2 × 10 8 particles and P = 2. The results obtained with the above three options are in very good agreement on both observables. The fact the 'reference' solution (P = 6) and the 'best compromise 31 ' (P = 2) have equivalent accuracies testifies of the fast convergence rate of gPC with respect to P . The fact the gPC-i-MC scheme also does with P = 2 testifies the new method can take advantage of it. Now regarding the cost of the three above options, we have reference : cost = N GL × the averaged CPU time of one run ≈ 50 × 85.0s, -best compromise : cost = N GL × the averaged CPU time of one run ≈ 4 × 85.0s, -gPC-i-MC : cost = 1 × the effective CPU time of one run = 1 × 86.2s. The new gPC-i-MC method recovers the same results as the best compromise solution with only one run of the MC simulation device and very similar computational times (at least for this problem, this will not exactly be the case for the next benchmarks): it testifies of the insensitiveness to the number of dimensions of the tracking 32 of the MC particles. The runs were performed on N replication = 32 replicated 33 domains for the parallel strategy, in very similar conditions as in[START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF]. The latter remark allows insisting on the fact the parallel strategies applying for a classical MC simulation device also applies with the new gPC-i-MC one. The minor code modifications described in sections 4-C to implement the gPC-i-MC solver do not imply a porting of the parallel counterpart, it is straightforward if already developed. A short discussion on parallel possibilities is provided in appendix C with the description of the new algorithm. The gain in computational time, for this test-case, is approximatively of N GL ≈ 4. In the above example, it remains relatively low: if one has access to 4(N GL ) × 32(N replications ) = 128 processors, which is common nowadays, the computational times are equivalent and the gain is only in term of computational ressources34 . We here recall that MC simulation codes are known to be computationally intensive and even such low factor (gain of only 4) may be welcome. In the following, we tackle some multidimensional uncertain problems.

Fig. 6 .

 6 Fig.6. Every results of this figure have been obtained applying the gPC-i-MC scheme. Comparison between total and first order Sobol indices for the three uncertain inputs X 1 (top), X 2 (bottom-left) and X 3 (bottom-right).

1 3 [

 3 U ] and S tot 3 [U ]. -(reference) ni-gPC with N Q GL = 5 3 = 125 and (P + 1) Q = (3 + 1) 3 = 64 coefficients, -(best compromise) ni-gPC with N Q GL = 4 3 = 64 and (P + 1) Q = (2 + 1) 3 = 27 coefficients, -

  GL × the averaged CPU time of one run ≈ 125 × 3 min 50s, -best compromise : cost = N GL × the averaged CPU time of one run ≈ 64 × 3 min 52s, -gPC-i-MC : cost = 1 × the effective CPU time of one run = 1 × 4 min 50s. Each computation were performed on N replication = 32 replicated domains (i.e. processors).

Fig. 7 .

 7 Fig. 7. Every results of this figure have been obtained applying the gPC-i-MC scheme. Top left: mean and variance profiles. Top right: total Sobol indices. Others: comparisons between total and first order Sobol indices for the six uncertain inputs X 1 , X 2 , X 3 ; X 4 , X 5 , X 6 .

Fig. 8 .

 8 Fig. 8. Every results of this figure have been obtained applying the gPC-i-MC scheme. Left: mean and variance profiles. Right: total Sobol indices. Comparisons between the nominal (σ A t = 0.40, σA s = 0.40) and newly (σ A t = 0.15, σA s = 0.15) studied configurations.

13 #move the particle p 14 x 17 # 19 # 21 # 24 else 25 # 26 #move the particle p 27 x 30 #change its weight 31 w 32 # 33 s p ←-s p -τ > 0 34 end 35 end 36 end Algorithm 1 : 16 # 17 #move the particle p 18 x 19 # 20 s p ←-0 21 # 23 #do not change the weight of particle p 24 w 25 #tally the contribution of particle p 26 for 31 # 32 #move the particle p 33 x 36 #change its weight 37 w 38 #

 13141719212425262730313233116171819202123242526313233363738 p (0) = v p δt + x p (s p ), 15 #set the life time of particle p to zero: 16 s p ←-0 do not change the angle or velocity of particle p 18 v p (0) = v p (s p ) = v, do not change the weight of particle p 20 w p (0) = w p (s p ), tally the contribution of particle p 22 u(x, t, v)+ = w p (0) × u 0 (x p (0), v p (0)) 23 end see the recursive treatment in factor of 1 [0,t] (τ ) in (21) p (s pτ ) ←v p τ + x p (s p ), 28 Sample the velocity of particle p from P s(x p , s p , τ, v p , v ) dv 29 v p (s pτ ) = V p (s pτ ) ←-σ s (x p (s pτ ), s pτ, v p (s pτ )) σ t (x p (s pτ ), s pτ, v p (s pτ ))w p (s p ) set the life time of particle p to: The MC semi-analog scheme described in term of algorithmic operations in order to compute (backward) u(x, t, v). It corresponds to the general canvas of MC simulation codes. The while loop is commonly called the tracking of a MC particle. The + = operation is commonly called the tallying. δt = (s pτ ) + δt is the time remaining to reach census p (0) = v p δt + x p (s p ), set the life time of particle p to zero: do not change the angle or velocity of particle p22 v p (0) = v p (s p ) = v, p (0) = w p (s p ), k ∈ {0, ..., P } do 27 u X k (x, t, v)+ = w p (0) × u 0 (x p (0), v p (0), X p (0))φ X k (X p (0)) see the recursive treatment in factor of 1 [0,t] (τ ) in (21) p (s pτ ) ←v p τ + x p (s p ),34Sample the velocity of particle p from P s (x p , s p , τ, v p , v , X p ) dv35 v p (s pτ ) = V p (s pτ ) ←-σ s (x p (s pτ ), s pτ, v p (s pτ ), X p )σ t (x p (s pτ ), s pτ, v p (s pτ ), X p ) w p (s p ) set the life time of particle p to: 39 s p ←s pτ > 0

  Comparison of CPU times for sequential (N replication = 1) and parallel (domain replication with N replication = 8) runs for gPC-i-MC and ni-gPC .

	3 = 01	1 min 44s	2 min 02s
	t CPU for (P + 1) Q = (1 + 1) 3 = 08	2 min 11s	2 min 33s
	t CPU for (P + 1) Q = (2 + 1) 3 = 27	4 min 16s	4 min 33s
	t CPU for (P + 1) Q = (3 + 1) 3 = 64	8 min 29s	8 min 52s
	ni-gPC	N replication = 1 N replication = 8
	t CPU	3 min 28s	3 min 45s
	Table 1

  top-left displays the mean and Mean and Variance of U (x, t = 1, X)

	1.4 1.1 1.2 1.3							MC-i-gPC mean MC-i-gPC variance			0.25 0.2	0.6 0.7 0.8	S tot i [U ](x, t = 1), i ∈ {1, .., 6} Sobol total X3 Sobol total X6 Sobol total X5 Sobol total X4 Sobol total X2 Sobol total X1
	1												0.15	0.5
	0.9													0.4
	0.8												0.1	0.3
	0.7													0.2
	0.6												0.05
	0.5													0.1
	0.4	0	0.1	0.2	0.3	0.4	0.5 x	0.6	0.7	0.8	0.9	1	0	0	0	10	20	30	40	50 x	60	70	80	90	100
			S tot 1 [U ](x, t = 1) vs. S 1 1 [U ](x, t = 1)			S tot 4 [U ](x, t = 1) vs. S 1 4 [U ](x, t = 1)

A description of a non-intrusive applications is detailed section 2 together with some references.

deterministic linear Boltzmann equation solved with a MC scheme.

The uncertainties could impact the initial or boundary conditions, the methodology suggested in this paper would be the same.

At the cost of more or less tedious pretreatments leading to a controled approximation[START_REF] Todor | Karhunen-Loève approximation of random fields by generalized fast multipole methods[END_REF][START_REF] Meyer | Efficient model reduction in non-linear dynamics using the Karhunen-Loève expansion and dual-weighted-residual methods[END_REF][START_REF] Spanos | Stochastic Finite Element Expansion for Random Media[END_REF][START_REF] Mercer | Functions of Positive and Negative Type and their Connection with the Theory of Integral Equations[END_REF][START_REF] Gaël Poëtte | A stochastic surrogate model approach applied to calibration of unstable fluid flow experiments[END_REF] and decorrelation[START_REF] Lebrun | A Generalization of the Nataf Transformation to Distributions with Elliptical Copula[END_REF][START_REF] Lebrun | An Innovating Analysis of the Nataf Transformation from the Copula viewpoint[END_REF].

We insist the notations is general enough: if X is a stochastic process, each X i denotes a realisation of this stochastic process.

deterministic counterpart.

uncertain counterpart.

or vector.

i.e. such that φ X k φ X l dP X = δ kl .

exponant BB appears every time we are considering a numerical approximation obtained from the black-box (BB) code.

Those MC particles will be defined in section 3.

Here, the configuration is such that ∆ = max(∆x, ∆t, ∆v) = ∆t, see appendix A.

described section 3.

The next description can easily be generalized and applied to the analog or the non-analog MC scheme.

It consists in applying the very same steps as above but on the adjoint version of[START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF].

with some proper boundedness properties of the cross-section σt see[START_REF] Papanicolaou | Asymptotic Analysis of Transport Processes[END_REF].

see figure1.

Note that we describe a MC resolution scheme for uncertainty based on the semi-analog (implicit capture) one. The next description can easily be generalized and applied to the analog or the non-analog MC scheme.

Here, the hypothesis ∀X ∈ SuppX is certainly not optimal but sufficient for the property to hold.

This assumption is commonly done.

The next expression is obtained inversing the cumulative density function of an exponential law, this is common in MC computations see[START_REF] Lapeyre | Méthodes de Monte Carlo pour les équations de transport et de diffusion[END_REF].

Inner would be more appropriate as we are here identifying the backward samplings.

This is easier to understand thanks to the algorithmic representations of appendices B and C.

See algorithm 1 for the definition of tallying.

I/O refers to input/output.

The number of coefficients increases with Q.

Multiply (25) by u(x, t, v, X) to obtain[START_REF] Jourdren | Arbitrary High-Order Schemes for the Linear Advection and Wave Equations: Application to Hydrodynamics and Aeroacoustics[END_REF].

Best compromise in term of relative accuracy and restitution times. To find it, we simply ran many tests for different N GL , P, N M C .

See algorithm 1 for the definition of the tracking of one MC particle.

Domain replication[START_REF] Dureau | Hybrid Parallel Programming Models for AMR Neutron Monte-Carlo Transport[END_REF][START_REF] Palais | Load balancing of parallel Monte Carlo transport calculations[END_REF][START_REF] Brunner | An efficient, robust, domain-decomposition algorithm for particle Monte Carlo[END_REF][START_REF]Monte carlo photon transport on shared memory and distributed memory parallel processors[END_REF][START_REF]Parallel performance study of monte carlo photon transport code on shared-, distributed-, and distributed-shared-memory architectures[END_REF] corresponds to the most common parallel strategy for MC simulation codes. It takes advantage of the independence between MC particles, hence populations of MC particles: N replication processors each have a batch of particles and the processors only communicate at the end of the time step to average over the N replication populations.

Note that rigorously speaking, the cost of ni-gPC remains given by the more costly run amongst the N GL but we will keep considering the average CPU time over the N GL calculations as a reference in the following.

This would not be the case for a numerical strategy based on the resolution of (35) for example.

Acknowledgements

The author would like to thank Hervé Jourdren: this work has been initiated within his team. Many thanks also to David Dureau, Adrien Bernede, Xavier Valentin, Boukhmes Mechitoua and Philippe Humbert for valuable discussions about High Performance Computing, Numerical Analysis and physical considerations.

Appendix A. A simple analytical uncertain solution

In this section we build an analytical solution in a simple uncertain configuration. It is used as a reference solution for the convergence studies of section 2. The test-problem is also revisited in section 5.1 with the new gPC-i-MC scheme. In particular, this analytical solution is used in the quantitative results of figures 1 and 2.

The configuration is monokinetic (i.e. v = 1) and homogeneous (i.e. u(x, t, v, X) = u(t, ω, X)). We assume the uncertainty, one-dimensional here for the sake of simplicity, affects the scattering cross-sections σ s = σ s + σs X, where X ∼ U[-1, 1] and σs is closely related to the variance of the uncertain scattering crosssection. Let us introduce U (t, X) = u(x, t, ω, X) dx dω. In the previously described configuration, the uncertain linear Boltzmann equation resumes to the following stochastic ordinary differential equation

satisfied by U . Introduce σ a = σ tσ s , then the solution is given by

The quantity U (t, X) is a random variable indexed by time t, i.e. it is a stochastic process. In this case, mean and variance of the stochastic process (A.2) can be computed analytically and are given by

Of course, higher order moments, probability of failure, complete characterisation of the probability density function of the stochastic process can be calculated but in figures 1-2 we focus on the variance V[U ](t) to perform the convergence studies. Note that in practice, we take v = 1, U 0 = 1, σ t = σ s = 0.1, σs = 0.1. The L 1 -norm of the error is computed at time t = 10. The curves of figure 1 implying a MC scheme are averaged over 128 computations with different seeds.

Appendix B. Algorithmic treatments to solve (1) with a MC scheme

set s p = t #this will be the remaining life time of particle p, it must go down to zero (backward) set s p = t #this will be the remaining life time of particle p, it must go down to zero (backward)

Sample X from the distribution having probability measure dP X .