Loris Ibarrart
email: loris.ibarrart@mines-albi.fr

Stéphane Blanco

Cyril Caliot

Jérémi Dauchet

Simon Eibner

Mouna El-Hafi

Olivier Farges

Vincent Forest

Richard Fournier

Mouna El Hafi

Jacques Gautrais

Raj Konduru

Léa Penazzi

Jean-Marc Tregan

Thomas Vourc'h E ,

Daniel Yaacoub

Advection, diffusion and linear transport in a single path-sampling Monte-Carlo algorithm : getting insensitive to geometrical refinement

come

Introduction

Applicative contexts are numerous where diffusion, advection and linear transport physics interact within systems that require quite refined geometric descriptions. Engi-neers attempting to closely understand and optimize the transfer of heat by conduction, convection and radiation in electronics and electric engines of increasing powers and decreasing sizes is among the most highlighted examples, but complexifying the geometric description is also a major requirement in domains such as housing insulation, energy reception in concentrated-solar plants, process optimization in steel industry, etc. Outside heat transfer, similar needs are associated to many types of transport physics, e.g. neutron transport for design and operation of nuclear reactors, charge carriers transport in solar cells, etc. In all such contexts numerical simulation imposes constraints in terms of geometrical description that engineers translate into trade-offs between accuracy and computer time or computer power. There is therefore a strong demand for numerical approaches that would face no limit as far as geometrical refinement is concerned.

When only linear transport is considered it can be stated that the problem is at least partly solved: thanks to two decades of research on acceleration techniques for pathtracing Monte Carlo algorithms, the film industry now renders highly complex scenes with such an ease that they free the artists of all previous constraints. The scenes can be conceived without any care about how the rendering algorithms will behave: it is acted that there is very little impact of geometry on accuracy and computation times. Physicists have successfully adapted and even upgraded these tools for the handling of other types of linear transport questions, e.g. for solar and infrared radiative heat transfer in planetary atmospheres where the simulation of 3D multiple scattering in cloudy scenes is now strictly insensitive to the level of detailed description of both the cloud volume and the ground surface [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF].

We address the very same need for one-speed linear transport coupled with diffusion and advection. This implies that 1. new paths are defined to account for the coupling of these three physics; 2. these paths preserve the features that allow path-tracing acceleration in computer graphics.

For pure linear transport, path-sampling Monte Carlo algorithms are designed via the available path-integral formulation of the general solution of the linear Boltzmann equation, which is possible whatever the geometry. It is indeed quite straightforward to switch from this specific partial differential equation (PDE) to an integral over particlepaths that are successions of continuous lines (straight lines for linear transport in the absence of external force, or for radiation in the absence of refraction effect) interrupted at collision events (absorption and scattering within the volume, absorption and reflection at the boundary). These paths will here be called linear-transport paths.

Advecto-reacto-diffusive PDEs can be handled similarly using the Feynmann-Kac formula [START_REF] Kac | On some connections between probability theory and differential and integral equations[END_REF], which gives the exact solution as an expectation of an Itô process driven by a stochastic differential equation (SDE) [START_REF] Graham | Stochastic simulation and Monte Carlo methods: mathematical foundations of stochastic simulation[END_REF]. The Itô process being Markovian, Monte-Carlo techniques can be used to solve these PDEs by simulating random paths of the stochastic process, the same way linear-transport paths are sampled to solve the linear Boltzmann equation. Feynmann-Kac formula already handles a first level of coupling: the physics of advection, reaction and diffusion within a unique geometrical domain are handled using paths that will here be called advecto-diffusive paths (or shifted-brownian paths). The name of these paths does not include reaction because reaction occurs either at all location along the path, or at the beginning and the end of each path (depending on the formulation): the paths themselves translate advection and diffusion only. In our case, if the Boltzmann equation solution were known at all location within the volume and at the boundary, then linear transport would enter the description via such a reactive term and the advecto-diffusive paths would transport the corresponding sources throughout the domain as expected along Feynamnn's vision of propagative processes, in strict correspondence with Green formalism.

In the applicative contexts listed above, the coupling goes further in three directions:

• in the advecto-reacto-diffusive equation, sources that obey Boltzmann equation are unknown,

• these sources depend on the solution of the advecto-reacto-diffusive equation itself, since Boltzmann equation includes source terms obeying the advecto-reactodiffusive equation;

• the domain is typically divided into subdomains where there is no advection and no contribution of linear transport physics (solids in the following) and subdomains with all three physics interacting (fluids in the following).

We therefore need a formalism in which advecto-diffusive paths interact with lineartransport paths, and also advecto-diffusive paths interact with pure diffusive paths at the interface between subdomains (depending on the retained physics, typically ensuring continuity of the flux density when combining solid-diffusion on one side, fluid-diffusion on the other, and linear transport particles emission/absorption at the interface itself). This question was already addressed by several authors with the objective of simulating the transfer of heat in 3D porous media. [START_REF] Gerard | A hybrid random walk method for the simulation of coupled conduction and linearized radiation transfer at local scale in porous media with opaque solid phases[END_REF], [START_REF] Caliot | Combined conductive-radiative heat transfer analysis in complex geometry using the monte carlo method[END_REF] and [START_REF] Sans | Solving transient coupled conductive and radiative transfers in porous media with a monte carlo method: Characterization of thermal conductivity of foams using a numerical flash method[END_REF] dealt with thermal diffusion inside opaque solids interacting via radiation, with also convective heat transfer along the solid surfaces for a fluid at a known temperature: advection and diffusion within the fluid is not part of the simulation. [START_REF] Penazzi | Stardis: Propagator evaluation for coupled heat transfer in large geometric models[END_REF] adds the simulation of convection within fluid cells that are assumed perfectly mixed, which means that all three heat transfer modes are coupled but with a simple fluid modeling in which no advection is required. Among these researches, [START_REF] Fournier | Radiative, conductive and convective heat-transfers in a single monte carlo algorithm[END_REF] started to formulate the corresponding path integrals, using a doublerandomization approach, introducing paths that are recursively defined as successions of sub-path corresponding to each of the three heat transfer modes, still with a convection model requiring no advection, and only at the stationary limit. The theoretical foundations that make this coupling possible are fully exposed in [START_REF] Marc Tregan | Coupling radiative, conductive and convective heat-transfers in a single monte carlo algorithm: a general theoretical framework for linear situations[END_REF]. We will here make use of this very same formalism for problems involving pure diffusion in some parts of the geometry (solids), advecto-diffusion in the remaining part (fluid) and one-speed linear transport within the fluid and between the solid-fluid interfaces (exchanges between fluid volumes, between a fluid volume and the solid surfaces, and between the solid surfaces through the absorbing and scattering fluid). This will allow us to test the corresponding advection-diffusion/linear-transport path-sampling Monte Carlo algorithms on realistic configurations inspired from solid-fluid porous heat exchangers.

As far as preserving the features that make path-tracing acceleration efficient in computer graphics, one-speed linear transport can be handled exactly as radiative transfer in computer graphics and the recursive switching from one sub-path to the other is independent of the detail of the geometrical description. The question reduces therefore to accelerating advecto-diffusive paths inside a geometrical domain defined using a very large number of geometric primitives. A very extensive amount of work is reported concerning the fast sampling of diffusive or advecto-diffusive paths for Monte Carlo algorithms based on Feynmann-Kac formula. They are essentially gathered around the walk on sphere denomination [START_REF] Kk Sabelfeld | Integral formulation of the boundary value problems and the method of random walk on spheres[END_REF] (and recently walk on rectangular parallelepipeds [START_REF] Deaconu | A random walk on rectangles algorithm[END_REF]). We put a great deal of efforts into the implementation and testing of these solutions together with the computer-graphics acceleration techniques, but at the present stage we could not reach convincing performances and were led to make other choices 1 . The main reason is that walk on sphere requires the knowledge of the closest distance to the boundary, which is typically a missing information when making use of the data structuring and access choices behind today advanced acceleration strategies (hierarchical grids added on top of unstructured sets of geometrical primitives). This remains a very open field of research [START_REF] Wald | Embree: a kernel framework for efficient cpu ray tracing[END_REF] and may lead to significant further acceleration potentials, but this point will not be further addressed here: we will construct an advecto-diffusive path sampling procedure as successions of straight-line jumps that may not be the best choice for long-term perspectives but are perfectly suited to available accelerators and will allow us to illustrate how such path-sampling Monte Carlo strategies have the same features as for pure linear transport: computation time is insensitive to the number of geometrical primitives.

In order to prepare this illustration, Sec. 1 introduces the physical model and Sec. 2 shows how the coupling of two physics is translated into the recursive sampling of subpaths, each dedicated to a single physics. An approximate brownian walk is then introduced in Sec. 3 that recovers all the features of linear-transport paths as far as tracking acceleration is concerned. Finally Sec. 5 provides simulation examples for a configuration emblematic of the heat transfer literature: a porous heat exchanger composed of an increasing number of kelvin cells.

Model

We consider a motionless solid domain Ω S of boundary ∂Ω S (only diffusion, no advection, no linear transport) and an adjacent absorbing and scattering fluid domain Ω F of boundary ∂Ω F (all three processes active at all locations). The physical properties are uniform within each part, the velocity v ≡ v(x) is known at each location x within the fluid, and we act on the geometrical refinement of the geometry by increasing the characteristic scale ratios of ∂Ω S and ∂Ω F (smallest physically significant scale divided by system scale), typically by increasing the number of pores in a porous geometry.

For the pure diffusion part and the advection-diffusion part of the model, the physical observable is noted η S ≡ η S (x) for each location x within the solid (x ∈ Ω S ) and η F ≡ η F (x) for each location within the fluid (x ∈ Ω F ). For the one-speed linear transport part, the physical observable is noted f ≡ f (x, u) at each location x ∈ Ω F and each direction u of the unit sphere (noted "4π"). In pure theoretical terms η and f can be interpreted as the density and the distribution function of two interacting species. The interaction between the two species is compatible with an equilibrium state where f is isotropic: • the emission direction is isotropically distributed; the corresponding emission rate is 1 4π times the diffusing species absorption (or reaction) rate ν a η F , • the absorption (or reaction) frequency ν a is equal for both the diffusing species and the linear transport species.

An example translation of this model for heat-transfer applications is provided in Appendix F. The boundary of the system is noted B, i.e. B = (∂Ω S ∪∂Ω F )-(∂Ω S ∩∂Ω F ). At each location y on B the incoming distribution function f i is known, i.e. f (y, u) = f i (y, u) for all u.n(y) > 0 where n(y) is the inward normal to the boundary at y. For the diffusive species at the boundary we note η B (y) the value of η at each boundary location and B is split in two parts, B = B D ∪ B R : on B D we use a Dirichlet boundary condition (η B is known), on B R the diffusion flux is null (∇η.n = 0). At the interface ∂Ω S ∩ ∂Ω F between the solid and the fluid, we write the flux continuity for the two species: the diffusion flux on the solid side, -D S ∇η S .n, equals the diffusion flux on the fluid side, -D F ∇η F .n, plus the linear-transport flux, j T .n where j T = 4π f cu du. There is no advective flux because velocity is assumed null at the solid-fluid interface and density is assumed continuous through the interface.

Altogether the model is

         ∆η S = 0, x ∈ Ω S D F ∆η F -∇.(η F v) -∇.j T = 0, x ∈ Ω F cu.∇f = -ν a f - η F 4π -ν s f - 4π p s (u|u )f du , x ∈ Ω F ; u ∈ 4π (1) (2) (3)
where D S and D F are the solid and fluid diffusion coefficients, c is the particles speed in the linear transport model, ν a is the absorption frequency, ν s is the single-scattering frequency, p s is the single-scattering phase function and f ≡ f (x, u ). Still within the system, at the solid-fluid interface,

           -D S ∇η S (y).n(y) = -D F ∇η F (y).n(y) + j T (y).n(y), y ∈ ∂Ω S ∩ ∂Ω F η S (y) = η F (y), y ∈ ∂Ω S ∩ ∂Ω F f (y, u) = α η S (y) 4π + u .n(y)<0 (1 -α)f p r (u|u )du , y ∈ ∂Ω S ∩ ∂Ω F ; u.n(y) > 0 (4) (5) (6)
where α is the surface absorptivity (i.e. one minus reflectivity) and p r the probability density of the direction after surface reflection. At the boundary of the system,

               f (y, u) = f i (y, u), y ∈ ∂Ω F ∩ B; u.n(y) > 0 η S (y) = η B (y), y ∈ ∂Ω S ∩ B D η F (y) = η B (y), y ∈ ∂Ω F ∩ B D ∇η S (y).n(y) = 0, y ∈ ∂Ω S ∩ B R ∇η F (y).n(y) = 0, y ∈ ∂Ω F ∩ B R (7) (8) (9) (10) (11)
Note that Robin boundary conditions would make no difference as far as the present numerical discussion is concerned.

Interacting random walks

Monte Carlo algorithms solving Eqs. 1, 2 and 3, each separately, are well established. They all start with the definition of a random process, the solution of the addressed model being shown to be the expectation of this process. For Eq. 1 (pure diffusion) the process is brownian motion. For Eq. 2 (diffusion-advection-reaction2 ), in accordance with Feynmann-Kac theory, the process is an exponentially interrupted brownian motion with drift. For Eq. 3 (linear Boltzmann equation), the process is a backward multiplescattering and reflection walk. The remaining question is therefore coupling.

When coupling is linear, the main idea behind the combined use of two Monte Carlo algorithms is double randomization [START_REF] Sylvain | Stochastic finite differences for elliptic diffusion equations in stratified domains[END_REF][START_REF] Bossy | Monte carlo methods for linear and non-linear poisson-boltzmann equation[END_REF]. If

• a Monte Carlo algorithm is available for solving problem 1 with N samples of process 1,

• sampling process 1 requires the solution of problem 2,

• a Monte Carlo algorithm is available for solving problem 2 with N samples of process 2, then there is no need to sample N times process 2 each time it is needed for process 1, which would imply N × N samples: it suffices to use only one sample of process 2. If for instance process 1 is an interrupted brownian motion that requires the solution of process 2 at the interruption location, and if process 2 is a backward multiple-scattering walk, then the overall process, representing the coupling of the two underlying physics, is simply the following: a brownian motion that, when interrupted, switches to a multiplescattering path. This remains meaningful when the coupling is both ways, i.e. sampling process 2 requires the solution of problem 1. This is the case when addressing the coupling of Eq. 2 and Eq. 3: when interrupted, the brownian motion switches to a multiplescattering path (because the reaction term involves f that is unknown), and when this multiple-scattering path is interrupted at an absorption/emission location, it switches back to brownian motion (because the emission term involves η F that is unknown). Figure 2 illustrates the practice of such a combination of brownian paths and multiplescattering linear-transport paths in the simplified case of an infinite uniform static fluid. The resulting Monte Carlo algorithm evaluates η F (x obs ) at a given location x obs assuming that η F is known at a given distance from x obs . It starts with the brownian walk by setting the initial location x B to x obs and the initial time t B to t obs . The first step is then the sampling of an exponentially distributed duration δt B for the brownian walk before interruption by emission. The expectation of this exponential distribution is 1 νa . Then the location of the brownian walker at t B +δt B is sampled according to a three-dimension gaussian distribution centered on x B of standard deviation σ = √ 6Dδt B . At this new location x T a linear-transport path is initiated, which starts with the isotropic sampling of a direction u. Then the algorithm is a standard multiple-scattering algorithm:

• a free-flight duration δt T is sampled according to an exponential distribution of expectation 1

νa+νs

• the multiple-scattering path is started with a straight line from x T to x T + cδt T u

• at this collision location a test is made to decide between absorption and scattering

• if scattering is retained a new direction u s is sampled according to p s and the linear-transport path is continued in this new direction.

The linear-transport path ends when absorption is retained, and at this absorption location a new brownian walk is started. This creates a succession of a brownian walk, a multiple-scattering path, a brownian walk, etc, until one of the multiple-scattering path is interrupted at a location x B,i where η F is known. The Monte Carlo weight w is w = η F (x B,i ) and the Monte Carlo estimate is the average of N such weights. In this simple configuration, both the brownian and linear-transport walks can be sampled exactly and the Monte Carlo estimate is strictly unbiased. The fact that the estimate is strictly unbiased in this example will not be preserved in the rest of this article. This is due to the interface and boundary conditions that will lead to brownian walks in confined domains: we will handle these confined walks via a numerical approximation that will introduce a bias. But appart from this approximation (that reduces to zero when lowering the numerical parameter to zero), the exact same approach will be used hereafter for the coupling of linear transport with diffusion and advection. In particular, the computation times indicated in Table A.1 illustrate a well known feature of Monte Carlo simulations that we will also observe in all the forthcoming simulation examples: in the low Knudsen number regime (or high optical thickness regime), i.e. when scattering widely dominates absorption and the mean free path is very small compared to the dimension of the system, then a large number of scattering events are required before either the linear-transport walk is interrupted to initiate a brownian walk, or it reaches the limit of the domain. This feature will not vanish with the following developments: when we will illustrate an insensitivity to the complexification of geometry, this will only mean that the refinement of the geometrical description has no impact on the computation time requirements, but these requirements will remain direct functions of the single-scattering albedo and the Knudsen number.

Brownian walks in confined domains

The work of Feynmann and Kac in the 50's [START_REF] Kac | On some connections between probability theory and differential and integral equations[END_REF] plays a key role in establishing the links between parabolic partial differential equations and random processes. Feynmann-Kac's formula settles indeed that the solutions of advection-diffusion-reaction equations (Eqs. 1 and 2 in our context) can be rewritten as expectation of stochastic processes. Together with the usual translation of linear transport theory in statistical terms (for Eq. 3), this is all the theoretical basis required to justify Monte Carlo algorithms such as the one illustrated in the preceding section. However, this statement leaves aside the huge amount of theoretical and numerical work that was required to deal with boundary conditions, before such algorithms became practical.

The family of Random Walk on Sphere (RWS) algorithms dominates the corresponding literature. The starting point is the observation that exact sampling of confined brownian walks is only possible for simple geometries and then the question becomes the following: how is it possible to cover any geometry using simple patterns such as spheres, sample exact brownian paths inside these patterns in a computationally efficient way, make numerical approximations in the vicinity of the boundaries (any boundary cannot be adjusted exactly with a sphere for instance), and ensuring that the parameters of these approximations can be tuned to reach the required accuracy (ensuring convergence with computational costs remaining satisfactory)?

The historic RWS method proposed by Brown in [START_REF] Brown | Modern mathematics for engineers[END_REF], and rigorously justified by Muler in [START_REF] Muller | Some continuous monte carlo methods for the dirichlet problem[END_REF], makes use of the distribution of first passages at the sphere, for a brownian motion starting from the center of the sphere: this distribution is uniform. Using the same logic as in the preceding section, this allows to write the solution of the stationary diffusion equation as the expectation of the density at random locations on the sphere (which is also the harmonicity property of the laplacian operator). In the present context, this solves Eq. 1 at the center of any sphere as soon as η S is known at all locations on the sphere:

η S (x) = 4π 1 4π du η S (x + δu), x ∈ Ω S , δ max δ ({S(x, δ) ∩ ∂Ω S } = ∅) (12)
where S(x, δ) is the sphere of radius δ centered at x. If η S is not known at the location x + δu sampled on the sphere, then double randomization is used with a new sphere centered on this new location, etc. The system boundary is then handled the following way:

• the sphere is chosen as the largest sphere entirely within the domain (the radius is the closest distance to the boundary, i.e. δ = max δ ({S(x, δ) ∩ ∂Ω S } = ∅)) and when the location sampled on the sphere is at a distance lower than a numerical parameter ε, then the random walk is considered to have reached the boundary;

• the corresponding approximation was very deeply investigated and convergence is satisfied, the number of successive spheres required to reach the boundary increasing only as the logarithm of 1 when reducing .

The corresponding algorithm used to evaluate the density at any probe position x obs is showcased in figure 3.a. This proposition has been extended in very numerous directions, notably with walks on rectangular parallelepipeds [START_REF] Deaconu | A random walk on rectangles algorithm[END_REF]. Considering the present context we can briefly highlight the fact that drifted brownian walks can be handled very similarly. For instance, Sabelfeld [START_REF] Sabelfeld | Random walk on spheres method for solving drift-diffusion problems[END_REF] used first passage statistics of drifted brownian motion in a sphere to extend the RWS method to coupled advection and diffusion for cathodoluminescence and electron beam induced current imaging. At the stationary state, for a homogeneous velocity field, it was proven that the distribution of first passage positions in a sphere corresponds to the von Mises-Fisher distribution. Equation 13shows how this distribution can be applied to equation 2. It can be noted that in the case of a zero velocity field, (a) this expression leads to equation 12, and to a Dirac distribution centered in v = u|v| at the limit of an infinite local Peclet number. For Eq. 2 without coupling with linear transport, where P e = δ|v| D F . When refining the geometrical description, the main point in Eqs. 12 and 13 is the fact that each step requires the computation of the local closest distance to the boundary (to ensure δ = max δ ({S(x, δ) ∩ ∂Ω S } = ∅) and accessing this information can be computationally very expensive. At least, although very significant research efforts are devoted to this question of accessing the shortest sphere-surface intersection [START_REF] Wald | Embree: a kernel framework for efficient cpu ray tracing[END_REF][START_REF] Shellshear | Fast distance queries for triangles, lines, and points using sse instructions[END_REF], the costs of today's available algorithms remain much higher than those of the algorithms accessing the shortest intersection with the geometry for a straight line in a given direction (linesurface intersection). So we can state that, at the present stage of these researches, the insensitivity to geometrical refinement illustrated in [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF] for radiative transfer, is directly connected to the fact that accessing the geometry is reduced to finding line-surface intersections for which computer graphics research has produced acceleration techniques based on precomputed recursive grids that are now very practically available for common usage by numerical scientists [START_REF] Pharr | Physically based rendering: From theory to implementation[END_REF][START_REF] James | The rendering equation[END_REF].

+ x + x1 + x2 + + x 3 (b) + x + x1 + x2 +x3 + x4 + x5 + x6 + x7 + x 8 δS
η F (x) = 4π p M F (u|v)du η F (x + δu), x ∈ Ω F , δ max δ ({S(x, δ) ∩ ∂Ω F } = ∅) (13
In brief, before sphere-surface intersection algorithms reach a similar maturity, RWS algorithms must be adapted to match this requirement of only accessing information about the geometry via line-surface intersections.

An example proposition compatible with ray-tracing

We make here a proposition for an algorithm inspired of the RWS literature, but with the constraint of using no sphere-surface intersection, only line-surface intersections. Then the acceleration techniques described for instance in [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF] will be directly available for a numerical implementation handling highly refined geometries. This is only one example proposition. Many alternatives could be formulated, aiming for instance at increasing the convergence order or reducing the computation time associated to the computation of line-surface intersections in the vicinity of the boundary, but they would all lead to the same conclusion as far as the objective of the present article is concerned: insensitivity to geometrical refinement can still be achieved when introducing brownian motion in a confined domain. Each of the following algorithms are detailed in Appendix D. They are validated against available analytical solutions in Appendix E and against numerical solutions in Appendix F.3 and Appendix F.4.

A finite difference scheme for pure diffusion in the solid

For the handling of geometrical confinement with only line-surface intersections, we made the choice of designing approximate brownian walks as direct statistical translations of standard finite difference or finite volume approaches to advection-diffusion-reaction partial differential equations. Regarding pure diffusion, the formulation of equation 12 is kept unchanged, but with a fixed step δ S (instead of the closest distance to the boundary). Therefore, at each step, given a sampled direction u on the unit sphere (noted 4π), the walk should shift to the position x + δ S u. However, in order to account for boundaries, the rays defined by (x, u) and (x, -u) are traced; let δ u and δ -u be the distances between x and the boundary in the directions u and -u respectively. Then, the walk is only shifted by the minimum of δ S , δ u , and δ -u , which leads to a position in Ω S or directly on its boundary ∂Ω S . This algorithm corresponds to equation 14 (see figure 3.b and Algorithm 2 in Appendix D).

η(x) = 4π 1 4π du η(x + δ(x, u) u), x ∈ Ω S , δ = min(δ S , δ u , δ -u ) (14)
Under the assumption that the density field is locally trilinear within S(x, δ S ), this equation exactly solves equation 1 (see a demonstration in Appendix B).

A Patankar's scheme for coupling diffusion with advection and linear transport in the fluid

Regarding coupled advection and diffusion, Patankar's work [START_REF] Suhas | Numerical Heat Transfer and Fluid Flow[END_REF] suggests a numerical scheme based on the analytical solution of equation 2 in a 1D homogeneous media. By interpreting the coefficients in Patankar's scheme as probabilities, we propose in this section an approximate random walk with drift for advection-diffusion. As in Sec. 2, coupling with linear transport will be handled thanks to an interruption time at which brownian motion switches to a linear-transport walk. In the present context, this interruption time must preserve all the numerical features inspired by Patankar work, ensuring that for a uniform velocity field the walk recovers the available exact solution of an exponential density profile in the velocity direction and a linear profile in the orthogonal plane. This is achieved by introducing a probability p T to switch to linear transport at each Patankar step. We give all the details of the formal developments leading to these probabilities in Appendix C. In a local orthonormal basis (O, e 1 , e 2 , e 3 ) chosen such that v = e 1 |v|, it leads to equation 15 where δ is a fixed step and u i = (-1) i e i+1 2 ( being the floor function, so that u 1 = -e 1 = -v/|v| and u 2 = e 1 = -u 1 ).

η(x) = p T η T (x) + (1 -p T ) 6 i=1 p i η(x + δu i ), x ∈ Ω F (15) with                                p T = δ 2 νa D F 4 + Pe(1+e Pe )
e Pe -1

+ δ 2 νa D F p 1 =
P e e P e P e (e P e + 1) + 4 (e P e -1)

p 2 =
P e P e (e P e + 1) + 4 (e P e -1)

p 3 = p 4 = p 5 = p 6 =
(e P e -1) P e (e P e + 1) + 4 (e P e -1) [START_REF] Suhas | Numerical Heat Transfer and Fluid Flow[END_REF] where η T = f du is the density of the linear-transport species described by equation 3 and P e = δ v D F . This equation can be directly interpreted as a random walk, which has the probability 1 -p T to move in 6 directions according to their 6 respective probabilities, and the probability p T to switch to a linear-transport walk (see Algorithm 3 in Appendix D). The above expression is valid for any orientation of the local orthonormal basis around e 1 = v/|v|. Therefore, equation 15 can be averaged over the direction e 2 (and subsequently u 3,4,5,6 ), leading to equation 17. In the algorithmic translation of this equation, at each step of the walk a diffusion direction u is uniformly sampled on the unit circle C(v) in the plane perpendicular to v. This way, when dealing with the boundaries, the diffusion directions can be treated exactly as in section 4.1. Following the same guidelines as previously, each step now requires 4 ray tracing to evaluate the length of the next step : 2 in the diffusion direction to compute δ u and δ -u , and 2 in the velocity direction to compute the distances δ v and δ -v between x and the first surface in the directions v and -v respectively (see Fig. 4). Those 4 distances δ v , δ -v , δ u , δ -u are then compared to the fixed step δ F (that is equivalent to δ S in the solid).

η(x) =p T η T (x) + (1 -p T ) p 1 η x -δ v(x) ||v(x)|| + p 2 η x + δ v(x) ||v(x)|| +4 p 3 C(v) 1 2π du η(x + δu) , x ∈ Ω F , δ = min(δ F , δ v , δ -v , δ u , δ -u ) (17)
where p T , p 1 , p 2 , and p 3 are defined as mentioned for equation 16 (but with an updated value of δ = min(δ F , δ v , δ -v , δ u , δ -u )).

𝛿 -𝑣

𝛿 𝑣 𝒗 𝒖 𝛿 𝑢 𝛿 -𝑢

𝒞(𝒗)

Figure 4: Example of ray tracing to evaluate the length of the next step in the approximate random walk based on Patankar scheme: starting from the current location, two rays are drawn in the velocity direction v and its opposite, and two others in a random diffusion direction u (and its opposite) in the plane orthogonal to v. These four rays enable to compute δ = min(δ F , δv, δ -v , δu, δ -u ).

Boundary and interface conditions

Using the same approach as in [START_REF] Lejay | Monte carlo simulations in media with interfaces[END_REF] for the solid-fluid interface, the density gradients in Eq. 4 can be approximated by finite differences using the same small steps δ S and δ F as in the previous sections. The source term j T (y) • n(y) in Eq. 4 can be expressed as the difference between the incident and emitted flux densities in the hemisphere around the surface normal n(y):

j T (y) • n(y) = α 4 c η(y) - α 4 c 2π u • n(y) π 4πf (y, -u) du (18)
which leads to

η(y) = p S η S (y-δ S n)+p F η F (y+δ F n)+p Ti 2π u • n(y) π 4πf (y, -u) du, y ∈ ∂Ω S ∩∂Ω F (19)
where

                   p S = δ F D S δ F D S + δ S D F + δ S δ F α 4 c p F = δ S D F δ F D S + δ S D F + δ S δ F α 4 c p Ti = δ S δ F α 4 c δ F D S + δ S D F + δ S δ F α 4 c (20)
Again, in Monte Carlo terms, these probabilities are to be understood using the double-randomization concept leading to successive sub-paths, as in Sec. 2. When a random walk hits a solid-fluid interface Eq. 19 makes the walk switch to a diffusive walk starting in the solid at (y -δ S n) with a probability p S (see Sec. 4.1), an advecto-diffusive walk starting in the fluid at (y + δ F n) with a probability p F (see Sec. 4.2), or a lineartransport walk starting at the point y in a direction -u toward the fluid, distributed according to the Lambertian law of surface emission (or according to an other law if surface emissivity depends on direction).

The same approach is used when dealing with boundary conditions involving a diffusion flux, i.e. Eq. 10 and 11, which leads to a normal reinjection:

η(y) = η S (y -δ S n(y)), y ∈ ∂Ω S ∩ B R η(y) = η F (y + δ F n(y)), y ∈ ∂Ω F ∩ B R (21)
In practice, even if the lengths δ S and δ F are supposed small enough, another surface might be present between y and its point of re-injection in the solid or fluid (e.g. in sharp corners). To avoid this problem, two rays are traced from the point y in directions n and -n to redefine the lengths δ S and δ F for this one step only : if a surface is detected between y and its point of re-injection, the length is reduced so that the re-injection point ends up strictly in the midpoint between the initial surface and the detected one, ensuring the reinjection occurs in the expected domain. Those corrected lengths can be formally written as δ Sc and δ F c :

         δ Sc = min δ S , δ -n 2 δ F c = min δ F , δ n 2
where δ n and δ -n are the distances between y and the next surface in direction n and -n respectively. These values of δ Sc and δ F c are finally used instead of δ S and δ F to compute the probabilities p S , p F and p Ti with Eq. 20 (see Algorithm 5 in Appendix D).

Insensitivity to geometrical refinement

Before discussing the computation times, we need to describe how we make use of computer graphic tools for accelerating the line-surface intersections when geometrical description involves large numbers of geometrical primitives (triangles as far as we are concerned). We will be very brief because the very same approach is implemented, with the same libraries, in [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF]. All required references and technical details can therefore be found in this reference. In short, before the simulation is launched the geometry is processed (all the triangles defining the system boundary and solid-fluid interface) and a recursive grid is constructed that will be used each time a line-surface intersection is computed. This grid has no impact on the final result: it is only meant to improve the computation times. The main ideas are that the grid is locally refined so that the final elements include only a limited number of geometrical primitives (only a limited number of intersections need to be computed when entering these final elements), and the structure of the grid is designed to optimize data access in memory and to allow efficient out-of-core procedures when needed (when facing large geometrical data sets). The open-source libraries implementing these acceleration tools are star-engine [START_REF]Star-engine[END_REF] and embree [1].

A brief summary at this stage:

• thanks to double randomization, coupling is handled with simple successions of sub-paths dedicated to each physics (Sec. 2),

• brownian walk algorithms are available that deal with confined domains with linesurface intersections only (Sec. 3),

• and acceleration tool are available so that the computation of line-surface intersection will take approximately the same time for small or large numbers of geometrical primitives [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF].

These are the three ideas that lead to the insensitivity that the present article aims to illustrate. We make this illustration using the porous structure of Figures 5 and6. An elementary structure is set as a stack of four Kelvin cells in direction e 1 and is deployed in the two orthogonal directions e 2 and e 3 to form a solid foam that will be crossed by a flow. At the inlet face the flow is uniform, parallel to e 1 . The three dimensional velocity field inside the foam is obtained numerically using a standard fluid mechanics solver. A symmetry condition is used for the lateral faces of the foam which allows that the velocity field be computed only once using the elementary structure (see Figure 5): the symmetry imposes that the flow remains the same in each stack of four Kelvin cells whatever the number of time this structure was deployed. We therefore define a porous geometry that can be enlarged by multiplying the number of pores as required, without the need to recompute the internal flow for each new configuration. As far as solid/fluid densities and the distribution function f are concerned, the same symmetry assumption is made at the four lateral faces of the deployed foam: null diffusive and advecto-diffusive fluxes, as well as specular reflection for linear transport. The inlet and outlet faces are shifted away from the foam (once the strand thickness at the inlet, twice the Kelvincell thickness at the outlet). At these two faces the required boundary conditions are therefore only for the density in the fluid and the distribution function:

• the density η F in the fluid is uniform at a known value η B = 0 along the inlet face

∂Ω F ∩ B D ,
• the diffusion flux in the fluid is null at the outlet and lateral faces ∂Ω F ∩B R (e.g. at the outlet face, the gradient ∇η F .e 1 for the density η F is null in the e 1 direction),

• the distribution function f for incoming directions at the inlet face is isotropic at

f i = 1 4π exp -y 2 2σ 2
where y is the 2D center position on inlet face and σ is five times the size of a single pore (see Figure 6),

• incoming distribution function at the outlet is uniform and isotropic at f i = 0,

• incoming distribution function on the lateral faces is provided by specular reflection boundary conditions,

• the diffusion flux in the solid is null at the lateral faces ∂Ω S ∩ B R (there is no solid at the inlet and outlet faces, viz. ∂Ω S ∩ B D = ∅).

Appendix F.4 provides an additional validation for the numerical schemes of Sec. 3 using one single stack of four Kelvin cells and comparing the Monte Carlo estimates to a full simulation made with a standard deterministic solver. In Fig. 7 we display Monte Carlo estimates that start from this four-cells configuration and increase the deployment up to more than ten thousand cells, the problem being truly three-dimensional because of the gaussian-shape incident distribution function f i at the inlet. Typical solving times for a one percent accuracy are 10s using a standard laptop, but as announced in introduction, the main result here is the fact that this computation time remains the same whatever the number of cells. Of course, as announced at the end of Sec. 2, the computation time remains a decreasing function of the Knudsen number (or equivalently, an increasing function of the optical thickness).

Conclusion

To be achievable, this insensitivity fundamentally requires each of the three items at the begin of the above section. First, the coupled problem is addressed in the same Monte Carlo terms as for a single physics : the solution of the problem is the expectancy of a random variable attached to a single random path. Combining several physics (linearly) preserves the main idea of starting from the observation location (starting from the camera in computer graphics, starting from the outlet in previous example) and following a random path that visits the system until it reaches the sources, constructing therefore in a backward manner (with an adjoint approach) how the sources propagate throughout the system to contribute to the observation. Of course, the path is more elaborate, with sub-paths initiated the one after the other. We insisted on the computational costs of multiple scattering in the low Knudsen number regime, but we could have also illustrated an increase of the computational cost when increasing the number of Kelvin cells along 7: Evolution of the averaged outlet density η and computation time for an increasing number of kelvin cells. Loading time corresponds to the time required to generate and copy the geometry in the RAM memory as well as to build the acceleration grid, whilst the solving time is the time required to sample 10 4 coupled random paths in the system. Diffusion coefficients are D F = D S = 5 10 -5 , particles speed is c = 0.05 m/s, surface absorptivity is α = 1, single-scattering and absorption frequencies are νs = νa = 12.5 s -1 , the single-scattering phase function is isotropic, viz. ps = 1 4π . Loading and computation times are obtained with an Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz CPU. the foam thickness. This effect is moderated by the backward following of advection, that tends to bring the paths back to the inlet, but we would still have shown that dealing with 5 cells instead of 4 would have increased the computation times. There is something like a "coupling-Knudsen number" that plays a role similar to the Knudsen number and that is intrinsic to today's Monte Carlo simulations, but we would still have observed an insensitivity to the transversal deployment of the foam.

This insensitivity means that if we leave aside the modification of the physics itself (i.e. the paths have similar structures in each simulation) then the computation time is the same whatever the number of triangles. We insisted on the second item: this is because we managed to design path-sampling algorithms for brownian walkers in confined spaces using only the geometrical computation of line-triangle intersections. This could very well be revisited if the ongoing computer graphics research on sphere-triangle intersections reaches the same maturity as for line-triangle intersections. At present, we needed this step of slightly revisiting brownian walks before we could benefit of the third essential point: the use of acceleration structures that are such that once the geometry has been loaded in memory, the computation times associated to the sampling of the path is strictly identical with four or several thousand cells. Altogether, we reach the very same conclusions as those of [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF] concerning the handling of heavily detailed volumedescriptions, only now the physics is not restricted to radiative transfer. In what we have here presented, there is however a strong restriction associated to the linearity of each of the considered phenomenologies and the linearity of their coupling: this is the condition under which double randomization allows the definition of single paths dealing with the coupled problem. We are exploring the possibilities of similar insensitivities with nonlinear physics along the proposition of [START_REF] Dauchet | Addressing nonlinearities in monte carlo[END_REF], but this will involve branching path statistics.

Appendix A. Simulation examples involving diffusion/linear-transport paths

An infinite uniform static fluid is considered (no drift) for which the fluid density η F is known outside a rectangular cuboid. The addressed quantity is η F at a given location x obs inside the cuboid. The algorithm is the one from Sec. 2 (illustrated in Fig. 2): a brownian walk starts at x obs and is interrupted exponentially in time; at this interruption location a linear-transport path starts until absorption; at the absorption location a new brownian walk is started, etc, until one of the brownian walks is interrupted outside the rectangular cuboid where η F is known; the Monte Carlo weight is the value of η F at this final location and η F (x obs ) is estimated as the average of N such weights. Table A.1 provides simulation examples for the following conditions:

• The faces of the cuboid are perpendicular to each of the unit vectors of an orthonormal basis (e 1 , e 2 , e 3 ), the two extreme summits are (0, 0, 0) and (L, 2L, 3L), and the observation location is x obs = (L/2, L/2, L/2).

• For x outside the cuboid η • The scattering Knudsen number c νsL is varied from 0.01 to 100.

F (x) = a • x + b with a = ( 2 3L 4 , 2 3L 4 , 2 3L
• The dimensionless diffusion coefficient D F cL is varied from 0.01 to 100. • The number of Monte Carlo samples is N = 10 5 .

In this configuration, the exact solution is η F,exact (x obs ) = 2 L 3 , independently of Knudsen number c νsL and dimensionless diffusion coefficient D F cL . On the other hand, both the number of required samples for a 1% accuracy and the computation time depend on these parameters. The computation time is notably sensitive to the scattering Knudsen number (each multiple-scattering path involve more scattering events when decreasing c νsL ) and to the dimensionless diffusion coefficient (when lowering D F cL each diffusive/linear-transport path involves more diffusion/linear-transport switches before exiting the cuboid). Monte Carlo samples and σ is its associated statistical uncertainty as provided by the standard error (take 3σ for a 0.9995 confidence). T is the computation time recorded when performing the N = 10 5 Monte Carlo samples, N 1% is the required number of Monte Carlo samples in order to reach σ = 0.01η F,exact (x obs ) and T 1% is the associated computation time. Computations were made using a x86i 64 Intel(R) Core(TM) i9-9880H CPU 2.30GHz.

Appendix B. Trilinear profile

Let the trilinear profile η(x ) = a • x + b be the local expression of the density inside a sphere S(x, δ S ) of radius δ S around x and check that Equation 14gives exactly the density at x when using the adaptative step δ(x, u) defined in section 4.1 (thanks to symmetry):

4π 1 4π du η(x + δ(x, u)u) = 2π 1 2π du 1 2 η(x + δ(x, u)u) + 1 2 η(x -δ(x, -u)u) = 2π 1 2π du 1 2 (a.(x + δ(x, u)u) + b + a.(x -δ(x, -u)u) + b) = 2π 1 2π du 1 2 (2a.x + 2b + a.u(δ(x, u) -δ(x, -u) =0 )) = 2π 1 2π du η(x) = η(x)
where the unit sphere is noted 4π and 2π is any hemisphere of that sphere.

These results indicate good convergence even for step lengths as big as half of the domain's size (i.e. δ S = 0.5). Such a good convergence is granted by the simplicity of the geometry of the domain Ω along with relatively smooth Dirichlet conditions on ∂Ω. In general, a step length of about 1/20th of the domain's size should be retained in order to ensure a minimum convergence (i.e. δ S = 0.05). anymore (see result in red). δ F plays therefore the role of standard numerical parameter for which convergence must be specifically tested (the value of 1/20th of the domain's size suggested above for pure diffusion can not be retained as such).

δ S x η(x) η(x) σ 0.5 0

Appendix F. Application to heat transfer

The model studied in the body of the article is quite general and can be applied to different applicative fields, e.g. neutron transport, charge carrier transport in semiconductors etc., as discussed in introduction. Among these potential applications, here we detail how the present work can be directly used to address heat transfer problems.

The example in section 5 is directly inspired of an academic porous exchanger from the heat transfer literature [START_REF] Kumar | Simultaneous determination of intrinsic solid phase conductivity and effective thermal conductivity of kelvin like foams[END_REF]. The main assumption in such standard heat transfer applications is that radiative transfer can be linearized as function of temperature, which imposes that the relative temperature differences remain limited. This assumption is detailled in Appendix F.1. Appendix F.2 provides the resulting correspondences between the model addressed in the body of the article and heat transfer. Appendix F.3 and Appendix F.4 provide numerical solutions, obtained from standard heat transfer solvers, that are used for validation.

Appendix F.1. Radiance temperatures

We start from the stationary radiative transfer equation formulated in terms of monochromatic specific intensity I ν ≡ I ν ( x, u, t) at position x, in direction u at time t and frequency ν :

u. ∇I ν = -k ν a I ν + k ν a I eq ν -k ν s I ν + k ν s 4π p s,ν ( u| u )I ν du (F.1)
where k ν a is the absorption coefficient, k ν s the scattering coefficient, p s,ν the scattering phase function, I ν ≡ I ν ( x, u , t) and I eq ν ≡ I eq ν (θ( x, t)) the specific equilibrium intensity at temperature θ( x, t) within the medium. Although the radiative transfer is stationary, I ν depends on t due to the evolution of the solid and fluid temperature.

Assuming that for all times and positions, the solid and fluid temperature remains close to a reference temperature θ ref , the temperature dependence of the specific equilibrium intensity can be linearized:

I eq ν (θ) ≈ I eq ν (θ ref ) + ∂ θ I eq ν (θ ref )(θ -θ ref ) (F.2)
Moreover, equilibrium properties allow to write :

0 = -k ν a I eq ν + k ν a I eq ν -k ν s I eq ν + k ν s 4π p s,ν ( u| u )I eq ν du (F.3)
Introducing the notation Ĩν = I ν -I eq ν (θ ref ) for the perturbations and subtracting equations (F.1) and (F.3), the radiative transfer equation under the assumption (F.2) can be written as follows:

u. ∇ Ĩν ≈ -k ν a Ĩν + k ν a ∂ θ I eq ν (θ ref )(θ -θ ref ) -k ν s Ĩν + k ν s 4π p s,ν ( u| u ) Ĩ ν du (F.4)
We choose to rewrite this equation using the radiance temperature θ ν R, u in the direction u. This radiance temperature is a spectral and directional quantity defined as the temperature for which the equilibrium specific intensity is equal to the specific intensity: For didactic reasons we here make the assumption of grey absorbing and scattering fluid (there would be no difficulty associated to the full preserving of spectral dependancies), which writes :

I eq ν (θ ν R, u ( x, t)) = I ν ( x,
u. ∇θ R, u ≈ -k a θ R, u + k a θ -k s θ R, u + k s 4π p s ( u| u )θ R, u du (F.9)
with θ R, u ≡ θ ν R, u at all frequencies ν.

Similarly at the boundary : where α is the surface absorptivity (or emissivity), θ the surface temperature and p r (u|u ) the distribution for reflection directions. Equation F.9 and F.10 define a closed system for θ R, u from wich the volumic radiative power ψ R and the surfacic radiative power ϕ R write ψ R = -ζ θ - , where λ is the thermal conductivity. Equation 1Steady state heat equation for opaque solids. Equation 2Steady state advection-diffusion heat equation for semitransparent fluids. Equation 3Radiative transfer equation F.9 Appendix F.4. Kelvin-cell porous heat exchanger Algorithm 5 along with algorithms 2, 3 and 4 are used to solve the 3D configuration presented in figure 5: a solid stack of 4 kelvin cells is heated by a radiative temperature θ hot set on the inlet face and cooled by the fluid at θ cold which enters by the same face. On the outlet face, the radiative temperature is set at θ cold to represent the environment. As in section 5, symmetry boundary conditions are used on the lateral faces. Radiative transfer is here reduced to only direct exchange between surfaces (the fluid is perfectly transparent, i.e. 5. Results for locations on the center line going down the flow; the inlet face is at adimensionned position 0 and the outlet at adimensionned position 1. In order to obtain a curve showing high variations of temperature in the fluid, the solid's and fluid's thermal diffusivities are respectively set at 80.10 -6 m 2 .s -1 and 130.10 -6 m 2 .s -1 . The inlet velocity field is uniform at 1 m.s -1 .

Figure 1 :

 1 Figure 1: Illustration of the advection, diffusion and linear transport configuration studied in the article. The solid domain Ω S is shown in gray, the fluid domain Ω F is shown in light blue. The model is: pure diffusion in the solid; advection-diffusion in the fluid; linear transport in the absorbing and scattering fluid. The boundary B = B D ∪ B R of the system is the union of B D , shown in grey, where Dirichlet boundary condition are formulated (the density η is prescribed), and B R , shown in red, where the diffusion flux is null. The incident distribution function f obeying linear transport is prescribed over the whole boundary B. See a practical implementation of this general configuration in Figures 5 and 6.

Figure 2 :

 2 Figure 2: A diffusive/linear-transport path in an infinite uniform static fluid for which the density η F is known outside a rectangular cuboid. N such paths are sampled to estimate η F (x obs ). Each path produces a Monte Carlo weight equal to the value η F at the location of the first brownian-walk interruption occurring outside the parallelepiped. η F (x obs ) is estimated as the average of these N weights. Blue sub-paths are brownian walks. Red sub-paths are multiple-scattering linear-transport paths. All algorithmic details are provided in the text of Sec. 2. Simulation examples are provided in Appendix A and Table A.1.

Figure 3 :

 3 Figure 3: Sampling of confined brownian walks to estimate the density η at point x. (a) Using the RWS method.In order to end the algorithm, the studied domain's surface is thickened by a small value in which the position is projected on the border. (b) Using the proposition made here, that is compatible with ray tracing. All the spheres have the same radius δ S and are distorted when getting close to the boundary.

  ) with the probability density of von Mises-Fisher distribution p M F (u|v)

Figure 5 :

 5 Figure 5: Stack of 4 kelvin cells and surface representation of the fluid flow for an homogeneous inlet fluid flow on the bottom right face and an outlet at the top left face. Fluid velocity at inlet face is 0.35 m/s; maximum velocity is 1 m.s -1 . Pore size is 4 mm and strand thickness is 0.7 mm. This configuration is an emblematic example to study heat transfer (see Appendix F for the correspondences with such problems).

Figure 6 :

 6 Figure 6: a) Representation of the deployed porous media submitted to an incident isotropic distribution function f i with gaussian spatial profile at the inlet face. b) Representation of a typical path in the 3D geometry, which is constituted by three type of subpaths: advective subpaths are in green, conductive subpaths in blue and linear transport subpaths in red.

Figure

  Figure7: Evolution of the averaged outlet density η and computation time for an increasing number of kelvin cells. Loading time corresponds to the time required to generate and copy the geometry in the RAM memory as well as to build the acceleration grid, whilst the solving time is the time required to sample 10 4 coupled random paths in the system. Diffusion coefficients are D F = D S = 5 10 -5 , particles speed is c = 0.05 m/s, surface absorptivity is α = 1, single-scattering and absorption frequencies are νs = νa = 12.5 s -1 , the single-scattering phase function is isotropic, viz. ps = 1 4π . Loading and computation times are obtained with an Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz CPU.

  4 ) and b = 1 L 3 . • The absorption Knudsen number c νaL is 1. • The single scattering phase function is Heiney-Greenstein phase function with asymmetry parameter g = 0.5.

3 :

 3 Results of algorithms 1, 3, and 5 for N = 10000 random samples on an analytic boundary value problem η(x) in a unit cube for various values of position x.e 1 = x.e 2 = x.e 3 = x = ||x|| √ 3 and step length δ F . Monte Carlo results are shown in column η along with their corresponding standard error σ.

θ

  R, u = αθ + (1 -α) u .n(y)<0 p r (u|u ) θ R, u d u (F.10)

θ

  -h R θu .n(y)<0 | u. n| π θ R, u du (F.12) with ζ = 16k a σθ 3 ref , h R = 4ασθ 3 ref , αis the surface emissivity and σ the Stefan-Boltzmann constant. Model in section 1 Heat transfer problems η ρC θ (temperature); We choose ρ F C F = ρ S C S = ρC for the sake of simplicity (this is not a restriction of the approach) 4πf ρC θ R, u (radiance temperature) ref (reference temperature for radiation linearisation); this correspondence is of application for the definition of the particles speed c. D a = λ ρC (thermal diffusivity)

(Figure F. 8 :Comsol -y = φl 2 MC -y = φl 2 Figure F. 9 :

 8229 Figure F.8: 2D-slice of an isolated duct structure used as a heat exchanger between a thermal source at θ hot and an incoming cold Poiseuille circulation at θ cold . Retained values for the simulation are: L l = 10, φ = 0.8, where φ is the ratio of inner duct radius over outer radius, λ S λ F = 10, where λ S and

  ν a = ν s = 0). Results are compared with an ANSYS Fluent simulation in figure F.10

Figure F. 10 :

 10 Figure F.10: Comparison of the temperature obtained with ANSYS Fluent and the Monte Algorithm for the stack of 4 Kelvin cells presented in Figure5. Results for locations on the center line going down the flow; the inlet face is at adimensionned position 0 and the outlet at adimensionned position 1. In order to obtain a curve showing high variations of temperature in the fluid, the solid's and fluid's thermal diffusivities are respectively set at 80.10 -6 m 2 .s -1 and 130.10 -6 m 2 .s -1 . The inlet velocity field is uniform at 1 m.s -1 .

  TableA.1: Simulation results for the diffusion/linear-transport algorithm of Sec. 2 with the parameters of Appendix A. L 3 ηF (x obs ) is the estimation of L 3 η F,exact (x obs ) = 2 that is provided by N = 10 5

	c νsL	D F cL	L 3 ηF (x obs )	L 3 σ	T (s)	N 1%	T 1% (s)
	100 0.01	2.001	0.004201 0.6247	4411	0.02756
	100	1.	1.990	0.006390 0.6897 10207	0.07040
	100 100.	1.927	0.05176 0.6919 669758	4.634
	1.	0.01	1.999	0.003586 0.7919	3214	0.002546
	1.	1.	1.993	0.006148 0.6196	9450	0.05855
	1.	100.	1.927	0.05165 0.6028 666982	4.021
	0.01 0.01	2.000	0.001850 99.04	855	0.8476
	0.01	1.	2.004	0.005765 16.37	8309	1.360
	0.01 100.	2.020	0.05145	13.16 661672	87.10

Table E .

 E

					δ F x	η(x) η(x) σ
	.1	1.3	1.297 0.0026	0.5 0.1 1.12 1.123 0.0029
	0.5 0.3	1.9	1.889 0.0077	0.5 0.3 1.48 1.474 0.0077
	0.5 0.5	2.5	2.479 0.0109	0.5 0.5	2	2.005 0.0094
	0.5 0.7	3.1	3.099 0.0104	0.5 0.7 2.68 2.682 0.0082
	0.5 0.9	3.7	3.706 0.0053	0.5 0.9 3.52 3.533 0.0040
	0.1 0.1	1.3	1.302 0.0025	0.1 0.1 1.12 1.117 0.0028
	0.1 0.3	1.9	1.895 0.0075	0.1 0.3 1.48 1.476 0.0078
	0.1 0.5	2.5	2.499 0.0105	0.1 0.5	2	2.000 0.0098
	0.1 0.7	3.1	3.097 0.0101	0.1 0.7 2.68 2.675 0.0085
	0.1 0.9	3.7	3.705 0.0051	0.1 0.9 3.52 3.525 0.0040
	0.05 0.1	1.3	1.301 0.0025	0.05 0.1 1.12 1.116 0.0027
	0.05 0.3	1.9	1.907 0.0075	0.05 0.3 1.48 1.490 0.0078
	0.05 0.5	2.5	2.503 0.0106	0.05 0.5	2	2.005 0.0098
	0.05 0.7	3.1	3.097 0.0101	0.05 0.7 2.68 2.671 0.0085
	0.05 0.9	3.7	3.702 0.0051	0.05 0.9 3.52 3.522 0.0041
	Table E.2: Results of algorithms 1, 2, and 5
	for N = 10000 random samples on an analytic
	boundary value problem η(x) in a unit cube
	for various values of position x.e 1 = x.e 2 =
	x.e 3 = x =	||x|| √ 3	and step length δ S . Monte
	Carlo results are shown in column η along with
	their corresponding standard error σ.

Table F .

 F 4: Correspondences between the model addressed in the body of the article and heat transfer problems.

Appendix F.2. Correspondences of physical quantities

Table F.

By leaving aside walk on sphere algorithms, we lose there remarkable convergence features when reducing the boundary thickness to zero.

From Equation

3, -∇.j T = -νa η F -f du which has the form of a standard reaction term provided that Eq. 2 is decoupled, i.e. provided that f is known.

Appendix C. A Patankar scheme including sources due to linear transport In Sec. 4.2 we proposed an approximate random walk with drift inspired of Patankar work for advection diffusion in confined domains, with a probability to switch to linear transport at each Patankar step. This appendix provides the developments leading to this set of probability. The resulting path-sampling algorithm is presented in the appendix Appendix D (see algorithm 3).

We address the following advection-diffusion equation with a source due to the coupling with linear transport (see equation 2 where -∇.j T = -ν a η -f du according to equation 3; we note η T = f du as in section 4.2):

where the flux density vector is

Using a finite volume approach, the balance on a cubic volume element of edge δ centered at x writes

where the densities η and η T in the source term of Eq C.1 have been taken uniform and equal to the densities η(x) and η T (x) at the center of the volume element. The fluxes J k=1,2,3,4,5,6 through the 6 faces of the cube are approximated according to Patankar approach, that uses the analytic solution of one dimensional advectiondiffusion equations -∇.j = 0. 6 independent one-dimensional equations are constructed by projecting the flux density vector on the outward-pointing normals u k=1,2,3,4,5,6 on each face of the cube:

with Dirichlet boundary conditions at x and x + δu k (η is known). The solution of this equation is

with the Peclet numbers Pe k = v.u k δ D F . The projections j • u k of the flux density vector (see Eq C.2) are obtained by differentiating Eq C.5 with respect to σ to write the diffusion term -D F ∇η • u k = -D F ∂η ∂σ and by substituting Eq C.5 into the advection term v • u k η:

This expression is used in the balance Eq C.3, assuming uniform flux densities on each face of the cubic volume element equal to the flux density at the center of the face (J k = δ 2 j • u k ):

Finaly, we choose to orient the cubic volume element such that u 1 = -v/|v| and u 2 = v/|v|, leading to : where Pe = ||v|| δ D F . Doing so (and using lim Pe→0 Pe exp (Pe)-1 = 1) we obtain the following expression for η(x):

Working on the ratios of exponential functions, those probabilities can be reformulated as in Eq. 16.

Appendix D. Algorithms

This appendix gathers all the propositions of Sec. 4 under an algorithmic form. These algorithms are those used in Sec. 5.

Algorithm 1: Monte Carlo algorithm used to evaluate η(x)

Call algorithm 3 else // on the boundaries or interfaces Call algorithm 5 end Store weight Uniform sampling of a random direction u on the unit circle C(v) orthogonal to the direction of v(x) Evaluate lengths δ v , δ -v , δ u , and δ -u with ray tracing (see Fig. 4) for which equation E.2 is solution everywhere in the field.

Results in table E.2 are obtained with algorithm 1, along with algorithms 2 and 5 (where we always have y ∈ B D here), in the case where Ω is the unit cube defined by ∀x ∈ Ω, 0 ≤ x.e i ≤ 1 where (O, e 1 , e 2 , e 3 ) is an orthonormal basis.

Appendix E.2. Validation of the random path compatible with ray tracing against an analytic solution for advection-diffusion Let us consider the following boundary value problem :

x ∈ Ω (E.4)

for which equation E.5 is solution everywhere in the field.

Results in table E.3 are obtained with algorithm 1, along with algorithms 3 and 5 (where we always have y ∈ B D here), in the case where Ω is the unit cube defined by ∀x ∈ Ω, 0 ≤ x.e i ≤ 1 where (O, e 1 , e 2 , e 3 ) is an orthonormal basis.

By comparison with the pure diffusion example of section Appendix E.1, present results indicate good convergence but not for step lengths as big as half of the domain's size where we used -∇.j T = ψ R (see equation F.11); θ S and θ F are the temperatures within the solid and the fluid respectively. At the solid-fluid interface,

where we used j T (y).n(y) = ϕ R (see equation F.11). At the boundary of the system,

Note that using the definition of the particles speed c as a function of θ ref provided in Table F.4, the probability set for path-switching at the boundaries in equation 20 becomes


Emmanuel Richaud
email: emmanuel.richaud@ensam.eu

Leticia Pires

Sébastien Roland

Laurianne Truffault

Lucien Laiarinandrasana

Alice Corfa

Julien Pinaud

Sylvain Caillol

Vincent Ladmiral

Christophe Beauvais

Sebastien Roland

Thermal Degradation of Fluorosilicone Elastomers

Keywords: ageing, elastomers, crosslinking

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

Introduction

Fluorosilicone rubbers (FVMQs) are known for their high performances in harsh environment, because of their high resistance to non-polar solvents, mineral oils, fuels compared to non-fluorinated silicones [1,2] and good level of gas permeability. However, they present a lower resistance to moderately polar solvents (esters, ketones…) and a lower thermal stability compared to silicones [3] . Due to their relatively high cost compared to other rubbers [ 4 ], they are usually employed in industrial domains such as aeronautics or defense where high and long lasting performances are required in fluid sealing, compression resistance, or stability to chemicals. Scientific literature mainly deals with their synthesis and the (relatively) wide range of structures and properties which can be improved by using comonomers [ 5 , 6 ]. Their mechanisms of degradation have been much less studied than that of other elastomers and assessment of their lifetime from predictive approaches such as Arrhenius extrapolations [7] is poorly reliable. The aim of the present work is thus to address a multiscale study of FVMQs thermal degradation as a preliminary step in order to identify some trackers of ageing and, in a close future, to establish the structure relationships involved in their ageing induced failure.

Experimental

FVMQ (Figure 1) was supplied as a crude rubber, which was cured using a dicumyl peroxide in a Gibitre press (170°C -4 min) so as to get 150 µm thin films. Vinyl group content was found close to 1.8% by 1 H NMR, corresponding to a double bond concentration about 0.12 mol kg -1 . DCP concentration was about half the concentration of double bond. Since each DCP decomposes into 2 radicals, it can be assumed that each double bond is reacted after curing. Number and weight average molar mass (in PS equivalent) were found close to Mn = 300 and Mw = 740 kg mol -1 . The samples with various ageing times were characterized by tensile tests to plot the stress strain curves and to derive the evolution of the initial Young's modulus with respect to the ageing time. Samples were punched as H4 dogbones samples prior to ageing. Tensile tests were conducted in standard conditions (20 °C, 50% HR) at a 20 mm min -1 elongation rate using an Instron 5966 tensile machine. Elastic modulus was determined in the 0-5% deformation region

Si CH 3 C H O CH 2 Si CH 3 C H 2 O CH 2 CF 3 n m << n
In situ isothermal ageing experiments at 280 °C were performed using TGA (Q500 -TA Instruments) under 100% O 2 or 100% N 2 supplied by a 50 ml min -1 gas flow (samples being previously heated from room temperature to 280 °C under nitrogen at a 20 °C min -1 ).

Characterization at molecular scale was performed by FTIR in transmission mode: spectra were obtained by averaging 16 spectra obtained at a 4 cm -1 minimal resolution using a Frontier 100 apparatus (Perkin Elmer). Elemental analyses were subcontracted at Institut des Sciences Analytiques (Lyon, France) and Crealins (Lyon, France).

Results

The engineering stress-strain curves of thin films after several ageing durations is given in Figure 2, where F is the load, S 0 is the initial cross section area, L is the applied displacement and L 0 is the gauge length. Nonlinear profile typically encountered for elastomers is observed (Figure 2a.) The initial Young's moduli (E 0 ) were then determined and plotted with respect to the ageing duration (Figure 2b). Young's modulus is increasing with the duration of ageing time. The ageing continuously stiffens the FVMQ material. It seemed to be interesting to investigate the mechanisms responsible for those architectural changes. For that purpose, TGA experiments were performed either under 100% O 2 or 100% N 2 atmosphere, in order to enhance the specific effects of each associated mechanism (either thermal ageing with possibly unzipping, or pure oxidation) (Figure 3):

-The existence of a mass loss even under nitrogen suggests the existence of an unzipping process.

-Under oxygen, kinetics of mass loss are relatively different: in the earlier times, residual mass under oxygen is higher than under nitrogen, suggesting that oxygenated species accumulate in the FVMQ matrix and display a certain stability even at 280 °C. At longer time, the sudden drop of mass suggests that those species decompose and generate volatile byproducts. Their nature remains to be determined. FTIR spectroscopy was used to highlight the structural changes. Spectra were recorded after several ageing durations under air (Figure 4). Contrarily to the case of many other polymers, there is no evidence of carbonyl species remaining in the polymer matrix after ageing, possibly because they evaporate as seen above. Interestingly, the decrease of C-H absorbance is observed and occurs simultaneously with mass loss. To better understand the chemical mechanisms involved in the thermal degradation of FVMQs, elemental analyses were carried out for samples aged either under air or under nitrogen (Table 1).

Table 1. Elemental analysis of aged FVMQs (percentage values in mass).

Those results call for the following comments:

-at first, for unaged samples, H and Si contents are in good agreement with the theoretical structure given in Figure 1. C content is lower than expected (ca 30-31%).

-under air, the content in Silicon increases whereas it stays closer to initial value for ageing under nitrogen.

-reversely, C and H content decrease under air whereas the value would remain closer to the value observed for unaged sample. Despite the incertitude on carbon content (see above), it can be noted that the ratio between C and H content remains about 6 for all ageing under air, suggesting that the mechanisms responsible for their loss must be strongly linked.

Discussion

According to experimental results, it seems that: -under nitrogen, an "unzipping" (depolymerization) mechanism would be plausible. The analysis of volatiles evolved under N 2 must be helpful to check its existence.

This preliminary mechanism remains to be fully validated basing for example on the analysis of evolved molecules by TGA-GC-MS or TGA-FTIR having in mind that the nature of some gaseous compounds (trifluoroethanal or trifluoropropanal) generated from the decomposition of alkoxy radicals by a β-scission process would allow to identify the "parent" hydroperoxide.

The analysis of residual aged polymer in a wider range of temperatures would also be helpful in particular to investigate the possible presence of carbonyl products.

Conclusions

This preliminary work reports a study on the thermal degradation of FVMQs. According to a multiscale study, ageing under nitrogen induces a non-negligible mass loss meanwhile elemental analysis suggests the overall composition remains the same. Reversely, under air, a This results in an increase of the Young's modulus of the material. A tentative mechanism is proposed and will be used to derive a kinetic model in a future step.

Figure 1 .

 1 Figure 1. Possible structure of FVMQ under study.

Figure 2 .

 2 Experimental results for FVMQs aged at 280°C under air: Engineering stress-strain curves (a) and evolution of the initial Young's modulus E 0 with respect to the ageing time (b).

Figure 3 .

 3 Figure 3. In situ ageing experiments monitored by TGA under N 2 or O 2 .

N 2 O 2 Figure 4 .

 224 Figure 4. Changes in the C-H stretching region of FVMQs after ageing at 280 °C under air.

  is accompanied by the release of hydrocarbon segments and crosslinking.

  under air, samples undergo a degradation mechanism leading to volatile release together with crosslinking. According to Table1, the volatiles would be constituted of hydrocarbonated groups, consistently with the observed correlation between mass loss and depletion of absorbance relative to C-H groups. A proposal of degradation mechanism is given below, but it remains now to validate this one by further experimental approach, in particular by evaluating the stability of peroxide species (formed in  or  position of Si atom) and the nature of crosslink bridges (Si-Si, Si-CH 2 -CH 2 -Si or Si-CH 2 -Si).

			0 h
			112 h
			285 h
			490 h
	absorbance	1.5 2	705 h 1015 h 1630 h 1870 h
		1
		3200	3100			2800
				wavenumber (cm -1 )
			virgin	100 h air	200 h air	300 h air	100 h N 2
		% Si	23.57	28.02	30.24	29.47	26.39
		% C	24	15.85	13.52	12.99	23.74
		% H	4	2.38	2.25	2.01	3.11

-

Acknowledgements

ITC Elastomères is gratefully acknowledged for having supplied the FVMQs samples.

Agence Nationale de la Recherche is gratefully acknowledged for having granted the FLuorosilicones Ageing Mechanisms and Synthesis (FLAMS ANR-20-CE06-0024) project.


S G Bobkov
email: bobkov@math.umn.edu

N Gozlan

C Roberto
email: croberto@math.cnrs.fr

P.-M Samson

POLAR ISOPERIMETRY. I. THE CASE OF THE PLANE

Keywords: Isoperimetry, Sobolev-type inequalities, Rényi Divergence Power, Relative Fisher Information

This is the first part of the notes with preliminary remarks on the plane isoperimetric inequality and its applications to the Poincaré and Sobolev-type inequalities in dimension one. Links with informational quantities of Rényi and Fisher are briefly discussed.

Isoperimetry on the plane and the upper half-plane

The paper by Diaz et al [D-H-H-T] contains the following interesting Sobolev-type inequality in dimension one.

Proposition 1.1. For any smooth real-valued function f on [0, 1],

(1.1)

1 0 f (x) 2 + 1 π 2 f (x) 2 dx ≥ 1 0 f (x) 2 dx 1/2 .
More precisely, this paper mentions without proof that (1.1) is a consequence of the isoperimetric inequality on the plane R 2 . Let us give an argument, which is actually based on the isoperimetric inequality (1.2) µ + (A) ≥ √ 2π (µ(A)) 1/2 , A ⊂ R 2 + (A is Borel), in the upper half-plane R 2 + = {(x 1 , x 2 ) ∈ R 2 : x 2 ≥ 0}. Here, µ denotes the Lebesgue measure restricted to this half-plane, which generates the corresponding notion of the perimeter µ + (A) = lim inf ε→0 µ(A + εB 2 ) -µ(A) ε (cf. e.g. [B-Z]).

Inequality (1.2) follows from the Brunn-Minkowski inequality in R 2 µ(A + B) 1/2 ≥ µ(A) 1/2 + µ(B) 1/2 along the same arguments as in the case of its application to the usual isoperimetric inequality (see [B-Z]). Indeed, applying it with a Borel set A ⊂ R 2 + and B = εB 2 (ε > 0), we get

µ(A + εB 2 ) ≥ µ(A) 1/2 + µ(εB 2 ) 1/2 2 = µ(A) 1/2 + π 2 1/2 ε 2 = µ(A) + √ 2π (µ(A)) 1/2 ε + O(ε 2 ),
and therefore (1.2) from the definition of the perimeter.

The relation (1.2) is sharp and is attained for the upper semi-discs

A ρ = {(x 1 , x 2 ) ∈ R 2 : x 2 1 + x 2 2 ≤ ρ 2 , x 2 ≥ 0}, ρ > 0.
In this case, µ(A ρ ) = 1 2 πρ 2 is the area size between the upper part of the circle x 2 1 + x 2 2 = ρ 2 and the x 1 -axis x 2 = 0, while the µ-perimeter is just the length of the half-circle µ + (A ρ ) = πρ.

To derive (1.1), one may assume that the function f is non-negative and is not identically zero on [0, 1]. Then we associate with it the set in R 2 + described in polar coordinates as

A = {(x 1 , x 2 ) : 0 ≤ r ≤ f (t), 0 ≤ t ≤ 1}
with x 1 = r cos(πt), x 2 = r sin(πt). Integration in polar coordinates indicates that, for any non-negative Borel function

u on R 2 , (1.3) R 2 u(x 1 , x 2 ) dx 1 dx 2 = π 1 -1 ∞ 0 u r cos(πt), r sin(πt) rdr dt.
Applying it to the indicator function u = 1 A , we get

µ(A) = π 2 1 0 f (t) 2 dt.
On the other hand, µ + (A) represents the length of the curve

C = {(x 1 (t), x 2 (t)) : 0 ≤ t ≤ 1} parameterized by x 1 (t) = f (t) cos(πt), x 2 (t) = f (t) sin(πt). Since x 1 (t) 2 + x 2 (t) 2 = f (t) 2 + π 2 f (t) 2 ,
we find that

µ + (A) = 1 0 x 1 (t) 2 + x 2 (t) 2 dt = 1 0 f (t) 2 + π 2 f (t) 2 dt.
As a result, the isoperimetric inequality (1.2) takes the form

1 0 f (t) 2 + π 2 f (t) 2 dt ≥ √ 2π π 2 1 0 f (t) 2 dt 1/2
. which is the same as (1.1). Note that the condition f ≥ 0 may easily be removed in the resulting inequality.

One can reverse the argument and obtain the isoperimetric inequality (1.2) on the basis of (1.1) for the class of star-shaped sets in the upper half-plane.

The same argument may be used on the basis of the classical isoperimetric inequality

(1.4) µ + (A) ≥ √ 4π (µ(A)) 1/2 (A is Borel)
in the whole plane R 2 with respect to the Lebesgue measure µ. It is attained for the discs

A ρ = {(x 1 , x 2 ) ∈ R 2 : x 2 1 + x 2 2 ≤ ρ 2 }, ρ > 0, in which case µ(A ρ ) = πρ 2 and µ + (A ρ ) = 2πρ.
Starting from a smooth non-negative function f on [-1, 1] such that f (-1) = f (1), one may consider the star-shaped region

A = {(x 1 , x 2 ) : 0 ≤ r ≤ f (t), -1 ≤ t ≤ 1}, x 1 = r cos(πt), x 2 = r sin(πt), enclosed by the curve C = {(x 1 (t), x 2 (t)) : -1 ≤ t ≤ 1} with the same functions x 1 (t) = f (t) cos(πt), x 2 (t) = f (t) sin(πt).
Integration in polar coordinates (1.3) then yields a similar formula as before,

µ(A) = π 2 1 -1 f (t) 2 dt,
and also the perimeter µ + (A) represents the length of C, i.e.,

µ + (A) = 1 -1 x 1 (t) 2 + x 2 (t) 2 dt = 1 -1 f (t) 2 + π 2 f (t) 2 dt.
As a result, the isoperimetric inequality (1.4) takes the form

1 -1 f (t) 2 + π 2 f (t) 2 dt ≥ √ 4π π 2 1 -1 f (t) 2 dt 1/2 , or equivalently, (1.5) 1 2 1 -1 1 π 2 f (t) 2 + f (t) 2 dt ≥ 1 2 1 -1 f (t) 2 dt 1/2 .
To compare with (1.1), let us restate (1.5) on the unit interval [0, 1] by making the substitution f (t) = u( 1+t2 ). Then it becomes 1 2

1 -1 1 4π 2 u 1 + t 2 2 + u 1 + t 2 2 dt ≥ 1 2 1 -1 u 1 + t 2 2 dt 1/2 . Changing x = 1+t
2 , replacing u again with f , and removing the unnecessary condition f ≥ 0, we arrive at:

Proposition 1.2. For any smooth real-valued function f on [0, 1] such that f (0) = f (1), (1.6) 1 0 f (x) 2 + 1 4π 2 f (x) 2 dx ≥ 1 0 f (x) 2 dx 1/2 .
As we can see, an additional condition f (0) = f (1) allows one to improve the coefficient in front of the derivative, in comparison with (1.1). It should also be clear that (1.6) represents an equivalent form of the isoperimetric inequality (1.4) for the class of star-shaped regions.

Relationship with Poincaré-type inequalities

It would be interesting to compare Propositions 1.1-1.2 with other popular Sobolev-type inequalities such as the Poincaré-type and logarithmic Sobolev inequalities. Starting from (1.1) and (1.6), a simple variational argument yields:

Corollary 2.1. For any smooth real-valued function f on [0, 1], (2.1) Var µ (f ) ≤ 1 π 2 1 0 f (x) 2 dx,
where the variance is understood with respect to the uniform probability measure dµ(x) = dx on the unit segment. Moreover, if f (0) = f (1), then

(2.2) Var µ (f ) ≤ 1 4π 2 1 0 f (x) 2 dx.
The constants 1 π 2 and 1 4π 2 in (2.1)-(2.2) are optimal and are respectively attained for the functions f (x) = cos(πx) and f (x) = sin(2πx) (cf. also [B-G-L]).

For the proof, let us note that an analytic inequality of the form (2.3)

1 0 f (x) 2 + cf (x) 2 dx ≥ 1 0 f (x) 2 dx 1/2
with a constant c > 0 becomes equality for f = 1. So, one may apply it to f ε = 1 + εf , and letting ε → 0, one may compare the coefficients in front of the powers of ε on both sides. First,

1 0 f ε (x) 2 dx = 1 + 2ε 1 0 f (x) dx + ε 2 1 0 f (x) 2 dx,
so, by Taylor's expansion, as ε → 0,

1 0 f ε (x) 2 dx 1/2 = 1 + ε 1 0 f (x) dx + ε 2 2 1 0 f (x) 2 dx - 1 8 2ε 1 0 f (x) dx + ε 2 1 0 f (x) 2 dx 2 + O(ε 3 ) = 1 + ε 1 0 f (x) dx + ε 2 2 1 0 f (x) 2 dx - ε 2 2 1 0 f (x) dx 2 + O(ε 3 ).
On the other hand, since

f ε (x) 2 + cf ε (x) 2 = 1 + 2εf (x) + ε 2 f (x) 2 + cf (x) 2 ,
we have

f ε (x) 2 + cf ε (x) 2 1/2 = 1 + εf (x) + ε 2 2 f (x) 2 + cf (x) 2 - 1 8 2εf (x) + ε 2 f (x) 2 + cf (x) 2 2 + O(ε 3 ) = 1 + εf (x) + cε 2 2 f (x) 2 + O(ε 3 ). Hence 1 0 f ε (x) 2 + cf ε (x) 2 1/2 dx = 1 + ε 1 0 f (x) dx + cε 2 2 f (x) 2 dx + O(ε 3 ).
Inserting both expansions in (2.3), we see that the linear coefficients coincide, while comparing the quadratic terms leads to the Poincaré-type inequality

c f (x) 2 dx ≥ 1 0 f (x) 2 dx - 1 0 f (x) dx 2 .
Thus, the isoperimetric inequality on the upper half-plane implies the Poincaré-type inequality (2.1) on [0, 1], while the isoperimetric inequality on the whole plane implies the restricted Poincaré-type inequality (2.2), with optimal constants in both cases.

Sobolev inequalities

If f is non-negative, then f (x) = 0 ⇒ f (x) = 0 and thus f (x) 2 + cf (x) 2 = 0. Hence, applying Cauchy's inequality, from (2.3) we get

1 0 f (x) 2 dµ(x) ≤ 1 0 f (x) f (x) + c f (x) 2 f (x) 1 {f (x)>0} dx 2 ≤ 1 0 f (x) dx 1 0 f (x) dx + c 1 0 f (x) 2 f (x) 1 {f (x)>0} dx .
Therefore, Propositions 1.1-1.2 also yield:

Proposition 3.1. For any non-negative smooth function f on [0, 1] with

1 0 f (x) dx = 1, (3.1) Var µ (f ) ≤ 1 π 2 1 0 f (x) 2 f (x) 1 {f (x)>0} dx,
where the variance is with respect to the uniform probability measure µ on the unit segment. Moreover, if f (0) = f (1), then

(3.2) Var µ (f ) ≤ 1 4π 2 1 0 f (x) 2 f (x) 1 {f (x)>0} dx.
Recall that there is a general relation between the entropy functional

Ent µ (f ) = f log f dµ -f dµ log f dµ (f ≥ 0)
and the variance, namely

(3.3) Ent µ (f ) f dµ ≤ Var µ (f ).
It is rather elementary; assume by homogeneity that f dµ = 1. Since log t ≤ t -1 and therefore t log t ≤ t(t -1) for all t ≥ 0, we have

f (x) log f (x) ≤ f (x) 2 -f (x).
After integration it yields (3.3).

Using the latter in (3.1)-(3.2), we arrive at the logarithmic Sobolev inequalities.

Corollary 3.2. For any non-negative smooth function f on [0, 1], with respect to the uniform probability measure µ on the unit segment we have

(3.4) Ent µ (f ) ≤ 1 π 2 1 0 f (x) 2 f (x) 1 {f (x)>0} dx. Moreover, if f (0) = f (1), then (3.5) Ent µ (f ) ≤ 1 4π 2 1 0 f (x) 2 f (x) 1 {f (x)>0} dx.
Replacing here f by (1+εf ) 2 and letting ε → 0, we return to the Poincaré-type inequalities (2.1) and ( 2.2) with an extra factor of 2. The best constant in (3.4) is however 1 2π 2 and in (3.5) is 1 8π 2 [B-G-L, Proposition 5.7.5]. On the other hand, the inequalities (3.1)-(3.2) are much stronger than (3.4)-(3.5).

Informational quantities and distances

The inequalities (3.4)-(3.5) may be stated equivalently in terms of informational distances to the uniform measure µ on the unit segment. Let us recall that, for random elements X and Z in an abstract measurable space Ω with distributions ν and µ respectively, the Rényi divergence power or the Tsallis distance from ν to µ of order α > 0 is defined by

T α (X||Z) = T α (ν||µ) = 1 α -1 p q α p dλ -1 = 1 α -1 f α dµ -1 ,
where p and q are densities of ν and µ with respect to some (any) σ-finite dominating measure λ on Ω, with f = p/q being the density of ν with respect to µ (the definition does not depend on the choice of λ). If α = 1, we arrive at the Kullback-Leibler distance or an informational divergence

T 1 (X||Z) = D(X||Z) = p log p q dλ = f log f dµ,
which is the same as Ent µ (f ). For α = 2 the Tsallis T 2 -distance is the same as the χ 2 -distance.

If α ≥ 1, necessarily T α (X||Z) = ∞ as long as ν is not absolutely continuous with respect to µ. In any case, the function α → T α is non-decreasing; we refer an interested reader to the survey [E-H].

In the case of the real line Ω = R, and when the densities p and q are absolutely continuous, the relative Fisher information or the Fisher information distance from ν to µ is defined by

I(X||Z) = I(ν||µ) = ∞ -∞ p p - q q 2 p dλ = ∞ -∞ f 2 f dµ,
still assuming that the probability measure ν is absolutely continuous with respect to µ and has density f = p/q. This definition is commonly used when q is supported and is positive on an interval ∆ ⊂ R, finite or not, with the above integration restricted to ∆. With these notations, Proposition 3.1 corresponds to the order α = 2 and therefore takes the form (4.1) T 2 (X||Z) ≤ 1 π 2 I(X||Z), T 2 (X||Z) ≤ 1 4π 2 I(X||Z), holding true for an arbitrary random variable X with values in [0, 1]. Here the random variable Z has a uniform distribution µ on [0, 1], and we use an additional constraint f (0) = f (1) in the second relation.

There is also another non-distance formulation of (3.6) in terms of classical informational quantites such as the Rényi entropy power and the Fisher information

N α (X) = ∞ -∞ p(x) α dx -2 α-1 , I(X) = ∞ -∞ p (x) 2 p(x) dx.
Here the case α = 2 defines the quadratic Rényi entropy power N 2 (X). If µ is supported and has an absolutely continuous positive density q on the interval ∆ ⊂ R, one may also define the restricted Fisher information

I 0 (X) = ∆ p (x) 2 p(x) dx.
For example, if Z is uniformly distributed in the unit interval, so that q(x) = 1 for 0 < x < 1, we have I(Z) = ∞, while I 0 (Z) = 0. In this case, if X has values in [0, 1], we have

T 2 (X||Z) = 1 0 p(x) 2 dx -1 = N 2 (X) -1/2 -1, I(X||Z) = I 0 (X).
Hence, the first inequality in (3.6) may be written as the following.

Corollary 4.1. For any random variable X with values in [0, 1], having there an absolutely continuous density, we have

(4.2) N 2 (X) 1 + 1 π 2 I 0 (X) 2 ≥ 1.
This relation is analogous to the well-known isoperimetric inequality for entropies,

N (X) I(X) ≥ 2πe,
where N (X) = N 1 (X) = e 2h(X) is the entropy power, corresponding to the Shannon differential entropy

h(X) = - ∞ -∞ p(x) log p(x) dx.
The functional I 0 (X) may be replaced with I(X) in (4.2) (since I 0 ≤ I), and then one may remove the assumption on the values of X. Moreover, with the functional I(X), this inequality may be considerably strengthened. Indeed, the relation N 2 (X)(1 + 1 π 2 I(X)) 2 ≥ 1 is not 0-homogeneous with respect to X, and therefore it admits a self-refinement when applying it to the random variables λX, λ > 0. Optimizing over this parameter, we will obtain an equivalent 0-homogeneous relation that is, N 2 (X)I(X) ≥ 1.

Observe that another inequality involving the quadratic Rényi entropy power N 2 (X) and some generalisation of Fisher information can be extracted from [L-Y-Z], namely for all 1 ≤ q < ∞, N 2 (X) q |p | q p ≥ C q for an optimal constant C q . However it's unclear how to related this inequality to (4.3).

  (4.3) N 2 (X)I(X) ≥ c,with c = π/4. But, it is obviously true that with c = 1. To see this, first note that, by the Cauchy inequality, for all x ∈ R,

Research was partially supported by the NSF grant DMS-1612961 and by the Bézout Labex, funded by ANR, reference ANR-10-LABX-58 and the Labex MME-DII funded by ANR, reference ANR-11-LBX-0023-01.


Marie Riebel
email: marie.riebel@unistra.fr

Odile Rohmer

Luisa Weiner

• People with autism spectrum disorder (ASD) are frequently exposed to public stigma 1 which can lead to interiorized stigma or self-stigma 2

• Self-stigma leads to numerous negative health, occupational and social outcomes 3

• Although adults with ASD present high levels of self-stigma 3 , there are no interventions targeting self-stigma in ASD.

• By targeting shame and self-blame 4 compassion can decrease self-stigma. Compassion Focused Therapy (CFT) could thus be useful in the reduction of the self-stigma associated with autism.

Aim: Evaluating the feasibility and efficacy of CFT for reducing the self-stigma of autism.

Single case experimental design (SCED)

Patient : male with two children, self-

WEEKLY SELF COMPASSION SCORES (SCS)

A Phase B Phase

single case experimental design of Compassion Focused Therapy to target the self-stigma associated with Autism Spectrum Disorder Methods Preliminary Results

												Randomized
												starting point
												A : No intervention	B : CFT sessions
		3
		2,8
		2,6
		2,4
	OUTCOME MEASURE	1,8 2 2,2
		1,6
		1,4
		1,2
		1
		1	3	5	7	9	11	13	15	17	19	21
		30
		WEEKLY INTERNALISED STIGMA SCORES (ISMI)
		29										S l e
												f -
		28										s t i g
												m
		26 27 OUTCOME MEASURE										a
		25
		24
									A Phase	B Phase
		23
		1	3	5	7	9	11	13	15	17	19	21

employed • Age : 56 yo • Age of ASD diagnosis: 53 yo • Previous therapy: DBT Repeated measures : • Self-stigma (ISMI) 5 (weekly) • Self-compassion (SCS) 6 (weekly) • Shame (daily) Introduction A Discussion 1. Overton & Medina, 2008; 2. Corrigan & Watson, 2002; 3. Dubreucq et al., 2021; 4. Gilbert, 2009, Kirby et al., 2017; 5. Hammer & Holand, 2017; 6. Neff, 2003 9 weeks 15 weekly one-hour sessions