Johann Guilleminot

Christian Soize

IT Ô SDE-BASED GENERATOR FOR A CLASS OF NON-GAUSSIAN VECTOR-VALUED RANDOM FIELDS IN UNCERTAINTY QUANTIFICATION * †

Keywords: random field, Itô stochastic differential equation, adaptive algorithm, uncertainty quantification

This paper is concerned with the derivation of a generic sampling technique for a class of non-Gaussian vector-valued random fields. Such an issue typically arises in uncertainty quantification for complex systems, where the input coefficients associated with the elliptic operators must be identified by solving statistical inverse problems. Specifically, we consider the case of non-Gaussian random fields with values in some arbitrary bounded or semi-bounded subsets of R n . The approach involves two main features. The first one is the construction of a family of random fields converging, at a user-controlled rate, towards the target random field. Each of these auxialiary random fields can be subsequently simulated by solving a family of Itô stochastic differential equations. The second ingredient is the definition of an adaptive discretization algorithm. The latter allows refining the integration step on-the-fly and prevents the scheme from diverging. The proposed strategy is finally exemplified on three examples, each of which serving as a benchmark, either for the adaptivity procedure or for the convergence of the diffusions.

Introduction.

1.1. Notation. For later reference, the following matrix sets are introduced: (i) M q (R) is the set of all the (q × q) real matrices. (ii) M S q (R) is the set of all the (q × q) symmetric real matrices. (iii) M + q (R) is the set of all the (q × q) symmetric positive-definite real matrices. For x and y in R n , we denote by x, y = n i=1 x i y i and x = x, x 1/2 the Euclidean inner product and the associated norm, respectively. Let {e 1 , . . . , e n } be the canonical basis of R n . The null and identity (q × q) matrices are denoted by [0 q ] and [I q ]. Let a ∧ b = min(a, b) and a ∨ b = max(a, b). Notations c, c x and c x are used to denote normalization constants, the values of which may vary from line to line. Let E denote the mathematical expectation.

Let {Ξ(x), x ∈ R d }, d 1, be a second-order mean-square continuous centered homogeneous R n -valued Gaussian random field, defined on a probability space (Θ, T , P) by a continuous M n (R)-valued normalized correlation function (

x, x ) → [R Ξ (x -x )] such that [R Ξ (x -x )] = diag(R 1 (x -x ), . . . , R n (x -x )) , ∀(x, x ) ∈ R d × R d , (1.1)
wherein each normalized correlation function x → R i (x) is defined from R d into [-1, 1] and is such that R i (0) = 1 for 1 i n (hence, the components of {Ξ(x), x ∈ R d } are mutually independent Gaussian random fields). We denote by x → Ξ(x, θ r ) the r-th independent realization of this Gaussian random field. Accordingly, and for any fixed value of x, Ξ x (θ r ) denotes the r-th realization of the R n -valued Gaussian random variable Ξ x = Ξ(x).

1.2. Preliminaries. This work addresses the random generation of a large class of random fields with values in some subset of R n . Specifically, it is concerned with the simulation of a class of random fields:

• That are non-Gaussian and of second-order.

• For which the family of first-order marginal probability distributions (m.p.d.f.) and spatial correlation lengths must be prescribed. • That contains the particular case of homogeneous random fields (for the translation in R d ). Models belonging to such a class are referred to as Prior Algebraic Stochastic Models (PASMs) in the literature of uncertainty quantification (UQ) and are typically devoted to identification solving underdetermined statistical inverse problems [START_REF] Soize | A computational inverse method for identification of non-Gaussian random fields using the Bayesian approach in very high dimension[END_REF] [START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF]. In order to fix the main ideas and concepts that will be discussed throughout, we first denote by {A(x), x ∈ Ω} the non-Gaussian random field, defined on (Θ, T , P), to be generated. It is assumed that this random field takes its values in some set S ∈ R n and is indexed by an open bounded domain Ω in R d . Let S • , S and ∂S = S\S • denote the interior, closure and boundary of S, respectively. The notation x ↓ ∂S means that x approaches the boundary ∂S from S • . Let {p t (• ; x)} x∈Ω be the prescribed family of first-order m.p.d.f. For most UQ applications, the aforementionned random field corresponds to a vector-valued representation of a random field with values in some set M + q (R) ⊆ M + q (R) (with n q(q + 1)/2), the matrix-valued random field being typically interpreted as the stochastic coefficient of an elliptic operator.

As a matter of illustration, let {[A(x)], x ∈ Ω} be a random field with values in

M + 3 (R) ⊂ M + 3 (R) such that: [A(x)] =   A 1 (x) A 2 (x) 0 A 2 (x) A 3 (x) 0 0 0 A 4 (x)   , ∀x ∈ Ω . (1.2)
Let a = (a 1 , a 2 , a 3 , a 4 ) ∈ R 4 . The random field {A(x), x ∈ Ω} such that A(x) = (A 1 (x), A 2 (x), A 3 (x), A 4 (x)) for all x in Ω thus takes its values in

S = {a ∈ R 4 | a 1 > 0 , a 3 > 0 , a 1 a 3 -a 2 2 > 0 , a 4 > 0} , (1.3)
and S = {a ∈ R 4 | a 1 0 , a 3 0 , a 1 a 3 -a 2 2 0 , a 4 0} . (1.4) This basic example shows that even for a simple case, the boundary ∂S can be rather complex, and it is worth mentioning that this complexity can fairly be expected to grow significantly as the dimension n increases.

1.3. Overview of the approach. The generation of non-Gaussian stochastic processes or random fields has been extensively discussed in the literature, at least for the case S = R n . A first class of algorithms is based on the construction of nonlinear memoryless transformations acting on Gaussian models (see [3] [42] for a recent survey, among others). Such transformations can be defined, for instance, through Hermite series expansions (see [18] [22] [41] and the references therein). The parameters of the underlying Gaussian model(s), as well as the coefficients of the expansions, are then numerically calibrated so that the image by the nonlinear mapping exhibits some target probabilistic properties -in most cases, a target marginal distribution and a specified correlation structure. In the present case where only the system of first-order marginal distributions is concerned, and since A(x) is a second-order vector-valued random variable, a convenient way to define such a mapping is to have recourse to a pointwise polynomial Chaos expansion (PCE) [START_REF] Ghanem | Stochastic finite elements: A spectral approach[END_REF] (see also [START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF] and [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF]):

A(x) = +∞ α,|α|=0 a α (x)h α (Ξ x ) , ∀x ∈ Ω , (1.5)
in which α = (α 1 , . . . , α n ) ∈ N n is a multi-index with modulus |α| = i α i and {h α } α are the multidimensional Hermite polynomials. In effect, the mean-square convergent expansion defined by Eq. (1.5) is truncated after polynomials of degree N pce , that is:

A(x) Npce α,|α|=0 a α (x)h α (Ξ x ) , ∀x ∈ Ω . (1.6)
From a random generation point of view, the use of such an expansion is certainly very efficient, as it can easily be implemented using fast matrix-based operations. In addition, the generation of a realization of {A(x), x ∈ Ω} requires the simulation of a single realization of the Gaussian random field {Ξ(x), x ∈ R d } only -as opposed to the second class of algorithms detailed below. However, one should note that the use of PCE rises some practical numerical difficulties which may degrade the sampling quality, among which the identification of the deterministic fields {x → a α (x)} α (even for moderate values of n) and numerical convergence issues to name a few.

A second general class of techniques involves transformations with memory and essentially relies on the construction of stochastic differential equations. The algorithm which is proposed and numerically studied in this paper belongs to this second class of numerical strategies and is devoted to the case where S is any bounded or semi-bounded part of R n . Broadly speaking, the approach consists in prescribing the family of target first-order marginal probability distributions as the family of invariant measures associated with a family of diffusion processes, the latter being defined as the stationary solutions of a family of Itô stochastic differential equations (ISDEs). Conceptually, such an approach amounts to alternatively define the random field {A(x), x ∈ Ω} through the following mapping:

A(x) = H x ({W (r, x), r ≥ 0}) , ∀x ∈ Ω , (1.7)
with H x a nonlinear operator and {W (r, x), r ≥ 0} a R n -valued normalized Wiener process defined on (Θ, T , P). The family {{W (r, x), r ≥ 0}} x∈Ω of Wiener processes thus introduced will be defined more precisely in § 2.1.

1.4. Outline. The outline of the paper, along with the intended contributions of this work, are as follows:

• In § 2, the theoretical background used for constructing each ISDE in the so-called free case (that is, for S = R n ) is first recalled for completeness. Specifically, we propose a particular definition of a family of Wiener processes that allows to generate spatial dependendies while solving the family of ISDEs. The definition of the latter is further addressed, and general assumptions related to the existence and uniqueness of invariant measures are recalled.

• Section 3 is concerned with the constrainted situation where S is any part of R n . In order to handle this case, an approximating family of non-Gaussian random fields, the definition of which is provided in § 3.2, is introduced. The strategy involves the construction of a family of approximating potential functions through a regularization procedure defined in § 3.3. • In § 4, we propose a novel adaptive integration algorithm for the family of ISDEs. The adaptation is performed on-the-fly, depending on the current values of the drift coefficients. It prevents from using a too small integration step -while ensuring a reasonably fast convergence towards the stationary solutions. • The proposed numerical strategy is finally exemplified through three applications in § 5. The first one deals with the restriction of a Gaussian bivariate random variable to a compact support. The latter, as well as the imposed variances, are calibrated such that the particle reaches the boundary of admissible domain S with a quite large probability, hence providing a benchmark for the regularization and adaptivity procedures. The second application addresses a regularization on a matrix manifold. In particular, a reference generator is used to illustrate the convergence in probability distribution for the proposed algorithms. The last exemple exemplifies the simulation of a non-Gaussian random field with values in the positive-definite cone.

2. Theoretical background for the free case.

2.1. Definition of a family of Wiener processes. Let W = {W (r, x) = (W 1 (r, x), . . . , W n (r, x)), r ∈ R + , x ∈ Ω} be a R n -valued centered second-order Gaussian random field such that:

(i) For all x in Ω, W (0, x) = 0 a.s.

(ii) The generalized r-derivative D r W of W is the cylindrical normalized Gaussian white noise N (see e.g. [START_REF] Kree | Mathematics of Random Phenomena[END_REF]). The covariance generalized function [C N ] of N is written as

∀(x, x ) ∈ Ω × Ω , ∀τ ∈ R, [C N (x, x , t + τ, t)] = δ 0 (τ )[R N (x, x )] , (2.1)
where δ 0 is the Dirac generalized function at the origin of R. In Eq. (2.1), [R N ] denotes the continuous M n (R)-valued function defined on Ω × Ω as

[R N (x, x )] = [R Ξ (x -x )] , (2.2)
where [R Ξ ] is the correlation function of the homogeneous Gaussian random field {Ξ(x), x ∈ R d } defined in § 1.2. It follows that for all x fixed in Ω, {W (r, x), r ≥ 0} is a normalized R n -valued Wiener process, that is, {W (r, x), r ≥ 0} satisfies the following properties: (i) The n real-valued stochastic processes {W 1 (r, x), r 0}, . . ., {W n (r, x), r 0} are mutually independent. (ii) W (0, x) = 0 a.s. (iii) The process {W (r, x), r 0} has independent increments. (iv) For all 0 s < t < +∞, the increment ∆W (s, t, x) = W (t, x) -W (s, x) is a R n -valued second-order random variable which is Gaussian, centered and whose covariance matrix reads as (t -s)[I n ]. Note that {W (r, x), x ∈ R d } is then a homogeneous colored Gaussian stochastic field for any fixed value of r ∈ R + .

2.2. Construction of a family of ISDEs. Throughout this section, we set S = R n . It is assumed that for all x in Ω, each element p t (• ; x) of the family of prescribed m.p.d.f. takes the following algebraic form:

p t (u; x) = c x exp{-Φ(u; x)} , ∀u ∈ R n , (2.3)
in which {Φ(• ; x)} x∈Ω is a family of continuous real-valued functions defined on R n . In what follows, those functions are classically referred to as the potential functions. Let F ISDE be a family of Itô stochastic differential equations indexed by x ∈ Ω, each element of which writes for r > 0:

dU (r, x) = V (r, x) dr , dV (r, x) = -∇ u Φ(U (r, x); x) - f x 2 V (r, x) dr + f x dW (r, x) , (2.4)
where {f x } x∈Ω is a family of free R + -valued parameters, U (0, x) = U 0 x and V (0, x) = V 0

x almost surely (a.s.) and {W (r, x), x ∈ Ω, r ∈ R + } is the Gaussian random field defined in § 2.1. The probability distribution P U 0 x ,V 0 x (du, dv; x) of (U 0

x , V 0 x ) is assumed to be known. Note that a deterministic initialization vector (u 0

x , v 0 x ) ∈ R n ×R n can be used, in which case

P U 0 x , V 0 x (du, dv; x) = δ 0 (u-u 0 x )⊗δ 0 (v -v 0 x )
, with δ 0 the Dirac measure at the origin of R n . We denote by {b

x : R n × R n → R 2n } x∈Ω and {[σ x ] ∈ M S
2n (R)} x∈Ω the families of drift vectors and diffusion matrices associated with the family of solution diffusion processes, defined for any x fixed in Ω as

b x (u, v) = v -∇ u Φ(u; x) - f x 2 v (2.5)
and

[σ x ] = [a x ][a x ] T , [a x ] = [0 n ] [0 n ] [0 n ] √ f x [I n ] (2.6)
for all (u, v) ∈ R n ×R n . We first assume that each ISDE of the family F ISDE admits a unique solution which is defined almost surely for all r ≥ 0 and that for r → +∞, there is an asymptotic stationary solution which is denoted as {(U s (r, x), V s (r, x)), r 0}. We denote as {P s (•, • ; x)} x∈Ω the family of invariant measures associated with the family of stationary solutions {{(U s (r, x), V s (r, x)), r 0}} x∈Ω and introduce the family {ρ s (•, • ; x)} x∈Ω of probability density functions such that for x fixed in Ω:

P s (du, dv; x) = ρ s (u, v; x) du dv . (2.7)
The construction of a unique asymptotic stationary second-order solution for any x ∈ Ω requires the family {Φ(• ; x)} x∈Ω to satisfy a few properties. Specifically, each element of the the family {Φ(• ; x)} x∈Ω of potential functions is assumed to satisfy the following properties, for any x fixed in Ω:

(H 1 ) u → ∇ u Φ(u; x) is a locally bounded function from R n into R n . (H 2 ) inf u R n >R Φ(u; x) → +∞ as R → +∞. (H 3 ) inf u∈R n Φ(u; x) = Φ min x , Φ min x ∈ R. (H 4 ) R n ∇ u Φ(u; x) R n p t (u; x)du < +∞.
Under the above assumptions, each element ρ s (• ; x) of the family {ρ s (• ; x)} x∈Ω is the unique solution of a steady-state Fokker-Planck equation [START_REF] Soize | The Fokker-Planck Equation for Stochastic Dynamical Systems and its Explicit Steady State Solutions[END_REF], hence implying the existence and uniqueness of the associated invariant measure P s (•, • ; x) (see Eq. (2.7)). It can further be shown1 [START_REF] Soize | Construction of probability distributions in high dimension usign the maximum entropy principle: Applications to stochastic processes, random fields and random matrices[END_REF] that for all x in Ω:

ρ s (u, v; x) = c x exp{- 1 2 v 2 -Φ(u; x)} , ∀(u, v) ∈ R n × R n . (2.8)
Therefore, one has

p t (u; x) = R n ρ s (u, v; x) dv , ∀u ∈ R n (2.9)
in view of Eq. (2.3). It follows that for any x fixed in Ω,

lim r→+∞ U (r, x) = A(x) (2.10)
in probability distribution. It can be deduced that solving the family of ISDEs allows sampling a family {U s (r * , x} x∈Ω (with r * sufficiently large) of random variables which defines a non-Gaussian random field admitting {p t (• ; x)} x∈Ω as its family of first-order m.p.d.f. and whose correlation structure results from the nonlinear transformation of the n correlation functions {R 1 , . . ., R n }.

3. Theoretical background for the constrained case. Let us now consider the general situation where S is any part, bounded or semi-bounded, of R n . In this case, the potentiel functions of the family {Ψ(• ; x)} x∈Ω do not satisfy assumption (H 1 ). We first provide in § 3.1 a typical example where such a situation is met in practice. The proposed strategy is then outlined in § 3.2. In effect, it relies on the definition of regularized potential functions, the construction of which is addressed in § 3.3.

3.1. An illustrative example in stochastic linear elasticity. In the framework of uncertainty quantification, the development of PASMs has been pioneered in [START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF] for the class of elliptic operators, opening up the way for extensions and generalizations to various situations (such as the modeling of permeability or elasticity tensor random fields with symmetry constraints) in [START_REF] Guilleminot | Non-Gaussian positive-definite matrix-valued random fields with constrained eigenvalues: application to random elasticity tensors with uncertain material symmetries[END_REF] [10] [START_REF] Guilleminot | Generalized Stochastic Approach for Constitutive Equation in Linear Elasticity: A Random Matrix Model[END_REF] [12] [START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF]. Such models are intended to ensure that the modeled random quantity does satisfy all the required fundamental mathematical properties (e.g. positive-definiteness) and are tailored to allow for an inverse identification through a limited parametrization. In this context, one is typically concerned with the construction of the family {p A(x) } x∈Ω of first-order m.p.d.f. associated with the S-valued random field {A(x), x ∈ Ω}, with2 S ⊂ R n . Following [START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF], a convenient way to achieve such a construction is to invoke the maximum entropy (MaxEnt) principle [START_REF] Jaynes | Information Theory and statistical mechanics[END_REF] [START_REF] Jaynes | Information Theory and statistical mechanics. II[END_REF]. The available information on random variable A(x) is then assumed to be synthesized by the following constraint (which does not include the normalization condition):

E{g (A(x))} = h x , ∀x ∈ Ω , (3.1)
in which a → g(a) is a measurable mapping defined from a given set V ⊇ S into R m , m 1, and {h x } x∈Ω is a family of vectors in R m . It can be shown that the MaxEnt-based p.d.f. a → p A(x) (a) takes the form:

p A(x) (a) = 1 S (a) c x exp {-λ x , g(a) } , ∀x ∈ Ω , (3.2)
wherein a → 1 S (a) is the indicator function for the support S of p A(x) and {λ x } x∈Ω is a family of vector-valued parameters defined on a subspace H λ of R m , the definition of which ensures the integrability of the p.d.f. given by Eq. (3.2). It is seen that the obtained p.d.f., defined by Eq. (3.2), takes a form that is similar (except for the presence of the indicator function) to the one stated by Eq. (2.3).

3.2. Definition of an approximating family for the non-Gaussian random field {A(x), x ∈ Ω} using the free space construction. Following the above discussion, let us consider the situation where each element p t (• ; x) of the family of prescribed marginal p.d.f. takes the form

p t (u; x) = c x exp{-Ψ(u; x)} , ∀u ∈ S , (3.3)
in which S denotes a semibounded (or bounded) subset of R n . We assume the following algebraic form (which is algebraically similar to the one derived from a MaxEnt approach; see Eq. (3.2)) for the potential function Ψ(• ; x):

Ψ(u; x) = -log(1 S (u)) + λ x , g(u) , ∀u ∈ S , (3.4)
with λ x ∈ H λ a given vector and g : V → R m a measurable mapping. The target p.d.f. thus reads as:

p t (u; x) = 1 S (u) c x exp{-λ x , g(u) } , ∀u ∈ S . (3.5)
Because of the indicator function, the potential function defined by Eq. (3.4) cannot be readily handled by using the derivations introduced in § 2 and necessitates the development of a specific methodology. For this purpose, we propose to proceed in two steps. First of all, an additional family {Ψ (• ; x)} x∈Ω of potential functions, with 0 < 1, is introduced. Specifically, this family is such that for all x in Ω, each function Ψ (• ; x):

(i) Is defined by

Ψ (u; x) = -log (1 S (u)) + λ x , g(u) , ∀u ∈ V , (3.6)
where 1 S denotes a positive regularized indicator function (defined in § 3.3) such that lim ↓0 1 S = 1 S and V is defined more precisely in the assumption below.

(ii) Satisfies hypotheses (H 1 -H 4 ). (iii) Is such that lim ↓0 Ψ (• ; x) = Ψ(• ; x) in S • .
Throughout the paper, it is assumed that the family {Ψ(• ; x)} x∈Ω satisfies the following assumption. Assumption 1. Let x be fixed in Ω. The potential function Ψ(• ; x) defined by Eq. (3.4) is assumed to be:

• Either repulsive on ∂S, in which case it is further assumed that the mapping g is defined on V = S • .

• Or non-repulsive on ∂S, in which case it is assumed that g is defined on V = R n and such that:

R n exp{-Ψ (u; x)}du < +∞ , ∀ > 0 . (3.7)
Let A (x) denote the random variable defined by the p.d.f.

p t (u; x) = c x exp{-Ψ (u; x)} , ∀u ∈ V , (3.8)
such that:

lim ↓0 p t (• ; x) = p t (• ; x) , ∀x ∈ Ω . (3.9)
For any fixed value of , we then consider solving a family F ISDE of ISDEs whose elements read as dU (r, x) = V (r, x)dr ,

dV (r, x) = -∇ u Ψ (U (r, x); x) - f x 2 V (r, x) dr + f x dW (r, x) , (3.10)
for r > 0 and x fixed in Ω. Eq. (3.10) is supplemented with suitable initial conditions. All the parameters (initial conditions, and so on) associated with the above family F ISDE are denoted using notations similar to those introduced in § 2, except for the subscript . For any fixed value of and following the derivations in § 2, one has:

lim r→+∞ U (r, x) = A (x) , ∀x ∈ Ω . (3.11)
Thus, the family {U s (r * , x} x∈Ω defines, in the stationary regime (i.e. for a sufficiently large value of r * ), a non-Gaussian random field {A (x), x ∈ Ω} that exhibits {p t (• ; x)} x∈Ω as its family of first-order m.p.d.f. The nonlinear transformation (see e.g. Eq. (1.7) for the free case) and the correlation function of the centered Gaussian random field {Ξ(x), x ∈ R d } completely define the system of marginal distributions of random field {A (x), x ∈ Ω}. When the first-order m.p.d.f. is invariant under translation in R d , this mapping is independent of x and random field {A (x), x ∈ Ω} is homogeneous for the translation in R d . Finally, for all x fixed in Ω, if Eq. (3.9) holds, then the sequence {A (x)} >0 converges (in probability distribution) towards A(x) as ↓ 0.

In a second step detailed in § 4, the use of an adaptivity procedure for the discretization of each ISDE is investigated. To this aim, we build on the algorithm proposed by Lemaire in [START_REF] Lemaire | An adaptive scheme for the approximation of dissipative systems[END_REF] for the approximation of diffusion processes with nonglobally Lipschitz drift vector fields. A key feature of the algorithm (which is detailed in § 4) is that the stochastic time step is automatically refined whenever the current state goes closer to the boundary of the admissible space, hence avoiding the use of a too small integration step that may slow down the convergence towards the stationary solution (and thus, increase the global CPU-time). At this stage, it should be noted that this paper is not concerned with the general study of some optimal adaptation mechanism or scheme for the discretization of stochastic differential equations (the optimality being understood in the sense of some error metric, e.g. a L 2 norm). The interested reader is referred to the large literature on this topic; see for instance [START_REF] Hofmann | Optimal approximation of stochastic differential equations by adaptive step-size control[END_REF] [27] [START_REF] Müller-Gronbach | The optimal uniform approximation of systems of stochastic differential equations[END_REF] for the definition of algorithms in the case of scalar stochastic differential equations with an additive noise; see e.g. [START_REF] Lamba | An adaptive Euler-Maruyama scheme for SDEs: convergence and stability[END_REF] and the references therein for a discussion and the derivation of an alternative adaptation approach controlled by the drift coefficient.

Remark 1. The first case in Assumption 1 is especially relevant to physically sound UQ analysis, in which one is usually concerned with the modeling of random variables with values in M + q (R). In this case, it is usually required that the vector representations A(x) and A(x) -1 are both second-order random variables, that is,

E{ A(x) 2 } < +∞ , E{ A(x) -1 2 } < +∞ , ∀x ∈ Ω .
(3.12)

It can then be shown that a sufficiently fast decrease of the probability measure near the origin 0 ∈ ∂S of R n is required for Eq. (3.12) to hold [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF], meaning that the gradient of the potential function significantly increases in the neighbourhood of ∂S. This important feature, combined to the additional repulsion generated by the regularized indicator function, makes the associated potential function strongly repulsive near ∂S, so that U is trapped inside the admissible domain S and never reaches the boundary in practice (see § 5.2). For V = R n , the point U is allowed to go beyond ∂S (in a way that depends on the selected value of ) and is then pulled back into S in finite time by the repulsion generated by the sole regularized indicator function (see § 5.1), hence avoiding the use of a rejection-like generation algorithm. This repulsive effect allows, through a simple particle tracking, to accomodate the case where the generator is subsequently plugged into a formulation which requires that U ∈ S • or U ∈ S almost surely at any point x in Ω. The tracking procedure may for instance consists in selecting, for each simulation run, a time instant r * r * such that U ( r * , x) ∈ S • (or U ( r * , x) ∈ S) for all x in Ω -such a strategy may, however, increase the computational time for simulations performed over large lattices.

Remark 2. Below, we briefly discuss an alternative approach following a strategy similar to that introduced above, except that the formulation is based on a classical Metropolis-Hastings algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]. Assume first that the random field {A(x), x ∈ Ω} under consideration has been discretized and must be accordingly sampled over a grid defined by a set {x 1 , . . . , x N } of N non-overlapping points. Upon using the same generation structure as the one followed for the discretization of each ISDE (see § 4), a Metropolis-Hastings-like generation scheme for the random field can be readily obtained as

A(x i ) = lim k→+∞ U k (x i ) , ∀i ∈ {1, . . . , N } , (3.13)
where

U k+1 (x i ) = U k (x i ) + σ n    N j=1 1 [0;min(1,q (k+1) j )] α (k+1) j    Ξ x (i) (θ k+1 ) , (3.14)
and U 0 (x i ) = u 0 x i , with u 0 x i an arbitrary initialization vector in S • . In Eq. (3.14), σ n > 0 denotes the standard deviation for the Gaussian transition kernel (see [START_REF] Roberts | Weak convergence and optimal scaling of random walk Metropolis algorithms[END_REF] for a discussion regarding the tuning of this parameter), {α

(k+1) j } N
j=1 is a set of N independent realizations of a random variable that is uniformly distributed over [0, 1] and

q (k+1) j = p t x U k (x j ) + σ n Ξ x (j) (θ k+1 ) p t x (U k (x j ))
, ∀j ∈ {1, . . . , N } .

(3.15)

Besides its simplicity (for implementation in particular), the main advantage of such a strategy is that it can readily handle any type of target m.p.d.f. without requiring the convergence of a sequence of approximating solutions to be studied. Nevertheless, it should be pointed out that such an algorithm necessitates the acceptance of the states at all the points x 1 , . . ., x N (i.e. that at least, U k (x j ) + σ n Ξ x (j) (θ k+1 ) ∈ S for all j in {1, . . . , N }) for a new sample to be generated (see Eq. (3.14)), hence resulting in a cumbersome computational effort for typical values of N . In practice, such a scheme may thus be restricted to the case where S occupies a large region of R n (so that the probability of reaching ∂S is reasonably small), as well as to moderate values of N .

Construction of the approximating family {Ψ (• ;

x)} x∈Ω of potential functions. Following § 3.2, let {Ψ(• ; x)} x∈Ω be the family of potential functions such that for all x in Ω:

Ψ(u; x) = -log (1 S (u)) + λ x , g(u) , ∀u ∈ S .
(3.16)

Let 1 S be the regularized indicator function defined as

1 S = K * 1 S , (3.17)
in which z → K (z) is a regularizing kernel on R n , parametrized by 0 < 1, and * denotes the convolution operator:

1 S (u) = R n K (z)1 S (u -z)dz , ∀u ∈ R n . (3.18)
Here, and without loss of generality, a regularization by a centered Gaussian kernel

K (z) = 1 n (2π) n/2 exp - 1 2 z 2 2 , ∀z ∈ R n , (3.19)
is used. It can be deduced that 1 S equivalently reads as

1 S (u) = E{1 S (u -Z )} , ∀u ∈ R n , (3.20)
where Z is a R n -valued centered Gaussian random variable whose p.d.f. is given by the right-hand side of Eq. (3.19). The approximating family {Ψ (• ; x)} x∈Ω of potential functions is then defined by:

Ψ (u; x) = -log (1 S (u)) + λ x , g(u) , ∀u ∈ V . (3.21)
Remark 3. Several approaches can be pursued in order to compute the indicator function 1 S . A brute-force approach consists in evaluating the regularization given by Eq. (3.17), regardless of the current value of u. From a computational standpoint, it is clear that such a strategy is certainly not optimal, unless the neighborhood of boundary ∂S is expected to be reached with a quite large probability. One may alternatively choose to compute the regularized function only when the state u is "sufficiently" close to ∂S. However, one should note that the numerical expense induced by such a strategy, which requires the distance between u and its projection onto ∂S to be computed, is conditionned by the nature of the boundary.

4. Definition of a discretization scheme with a stochastic adaptive step sequence. For k 0, we let

U k (x) = U (r k , x), V k (x) = V (r k , x) and W k (x) = W (r k , x)
, where {r 0 , . . . , r M } constitutes some stochastic discretization (see § 4.2) of time interval [0; T ], with T = r M . For x fixed in Ω and k 0, we denote as L k (x) the R n -valued random variable defined as:

∀j ∈ {1, . . . , n} , (L k (x)) j = - ∂Ψ(u; x) ∂u j u=U k (x) . (4.1)
Similarly, (L k (x)) j denotes the j-th component of the gradient of Ψ (• ; x) evaluated at u = U k (x). When a comparison between either numerical schemes or several parametrizations is performed (at a given point x), we use the same initial conditions and the same sample path for the increments of the normalized Wiener processes.

4.1. Numerical comparison between the explicit Euler and Störmer-Verlet schemes for the free case construction. Before proceeding to the construction of the general algorithm, we numerically compare below the convergence of the explicit Euler and Störmer-Verlet (SV) schemes, the former serving as a basis for Lemaire's derivations [START_REF] Lemaire | An adaptive scheme for the approximation of dissipative systems[END_REF]. Details about the SV algorithm for the integration of ordinary differential equations can be found in [START_REF] Vogelaere | Methods of integration which preserve the contact transformation property of the Hamiltonian equations[END_REF] -see also [START_REF] Hairer | Geometric numerical integration illustrated by the Störmer/Verlet method[END_REF] [16] [START_REF] Milstein | Stochastic Numerics for Mathematical Physics[END_REF] for surveys about geometric integration; see [START_REF] Milstein | Quasi-symplectic methods for Langevin-type equations[END_REF] and e.g. [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF] for geometric integration applied to Langevin-type equations. For illustration purposes, we consider a potential function derived in [START_REF] Soize | Construction of probability distributions in high dimension usign the maximum entropy principle: Applications to stochastic processes, random fields and random matrices[END_REF] (with S = R 7 ) for the MaxEnt-based probabilistic modeling of a band M + 4 (R)-valued random matrix [A].

Definition of the target probability distribution.

Let A be the R 7valued random variable corresponding to a vector representation (see below) of the nonzero entries of [A]. Let p A denote its p.d.f., obtained by invoking the MaxEnt principle. Let u → e(u) and u → [G(u)] be the two unique non-linear deterministic mappings defined from R 7 into R 7 and from R 7 into M + 7 (R) respectively, such that ∀u ∈ R 7 , e(u) = (u 4 1 , u

and [G(u)] = [L(u)] T [L(u)], with [L(u)] =     u 2 1 u 2 0 0 0 u 2 3 u 4 0 0 0 u 2 5 u 6 0 0 0 u 2 7     . (4.3)
It is assumed that the available information on random vector A takes the form:

E{e(A)} = e , E{ [G(A)] 2 F } = 7.8575 , E{ [G(A)] -1 2 F } = 5.6 , (4.4)
with e = (1, 0, 1, 0, 1, 0, 1). The associated potential function Φ : R 7 → R such that p A (a) = c exp{-Φ(a)} for all a in R 7 is then given by: Φ(u) = λ (1) , e(u) + λ (2)

[G(u)] 2 F + λ (3) [G(u)] -1 2 F , ∀u ∈ R 7 . (4.5)
In Eq. (4.5), the Lagrange multipliers λ (1) ∈ R 7 , λ (2) ∈ R + * and λ (3) ∈ R + * are such that the constraints given by Eq. (4.4) are fullfiled and are found to be (see [START_REF] Soize | Construction of probability distributions in high dimension usign the maximum entropy principle: Applications to stochastic processes, random fields and random matrices[END_REF]): λ (1) = (0.7381, 4.1697, 1.2465, -0.9248, 0.8998, 4.2584, 0.8714), λ (2) = 2.2293 and λ (3) = 1.7749.

Discretization schemes and numerical results.

Let U 1 = u 0 and V 1 = v 0 , where u 0 and v 0 are two given deterministic vectors in R 7 . For k = 1, . . . , M -1, the explicit Euler scheme reads as

U k+1 = U k + ∆r V k , (4.6) V k+1 = 1 - f ∆r 2 V k + ∆r L k + f ∆W k+1 , (4.7)
where the increment ∆W k+1 = W k+1 -W k of the normalized Wiener process W is a second-order Gaussian centered R n -valued random variable with covariance matrix ∆r[I n ].

For k = 1, . . . , M -1, the SV scheme is defined as

U k+1/2 = U k + ∆r 2 V k , (4.8
)

V k+1 = 1 -β 1 + β V k + ∆r 1 + β L k+1/2 + √ f 1 + β ∆W k+1 , (4.9
)

U k+1 = U k+1/2 + ∆r 2 V k+1 , (4.10)
with β = f ∆r/4. In order to provide a convergence analysis of the two above schemes, we introduce the function Conv : N → R + defined as:

Conv(M c ) = 1 M c Mc k=1 U k 2 , ∀M c 1 . (4.11)
This function measures the convergence of the ergodic estimator for the second-order moment of U k . We then consider the following error function ErrConv :

N → R + : ErrConv(M c ) = |Conv(M c ) -m ref 2 |/m ref 2 ∀M c 1 , (4.12)
in which m ref 2 denotes the reference value for the second-order moment of U . In practice, such a value is obtained by using the explicit Euler scheme (with a sufficiently small time step, e.g. ∆r = 10 -5 ) and an ergodic estimator over a large integration interval (typically M c = 10 9 ). The precision and rate of convergence (towards the stationary solution) are then illustrated below by studying the convergence of mapping M c → ErrConv(M c ). This convergence is displayed in Fig. 4.1 for several values of the integration step ∆r and for the two different schemes defined by Eqs. (4.6) and (4.8). As expected, it is seen that for large values of ∆r, the Störmer-Verlet scheme yields a better approximation than the explicit Euler one, whereas the approximations coincide for a sufficiently small time step. These results, together with other numerical evidences which are not reported here for brevity, motivate the combination of Lemaire's adaptivity strategy with the SV scheme for the random field case.

Construction of an adaptive

Störmer-Verlet scheme. The aim of this section is to propose an adaptive discretization algorithm for solving the family F ISDE of ISDEs (the parameter being fixed) introduced in § 3.2. This scheme involves (i) a specific definition of the family of Wiener increments (in accordance with § 2.1) and (ii) a stochastic step sequence which essentially depends on two different ingredients, namely: -a dominating deterministic nonincreasing step sequence, which is used whenever the current state is "far enough" from ∂S. -a stochastic counterpart that prevents the scheme from diverging and controls the adaptivity procedure as the particle reaches the neighbourhood of the boundary ∂S. Following [START_REF] Lemaire | An adaptive scheme for the approximation of dissipative systems[END_REF], let (γ k ) k≥1 be a positive nonincreasing sequence such that:

γ k = γ 0 k -1/τ , ∀k ≥ 1 , (4.13)
where γ 0 > 0 is a parameter corresponding to the initial time step and τ ∈ N * is a free parameter controlling the rate of decay for the deterministic step sequence (see the remark below). From Eq. (4.13), it is readily seen that lim k→+∞ γ k = 0 (4.14)

and k≥1 γ k = +∞ . (4.15)
Furthermore, we introduce a family of stochastic processes {χ k x , k ≥ 1} x∈Ω defined for all x as

χ k x = 2 V x (U k (x), V k (x)) b x (U k (x), V k (x) 2 ∨ 1 , ∀k ≥ 1 . (4.16)
In Eq. (4.16), V x is the Lyapunov function defined at point x as [START_REF] Mattingly | Ergodicity for SDEs and approximations : locally Lipschitz vector fields and degenerate noise[END_REF] V

x (u, v) = H x (u, v) + f x 4 u, v + f x 2 16 u 2 + 1 , ∀(u, v) ∈ V × R n , (4.17)
in which H x is the following Hamiltonian:

H x (u, v) = v 2 /2 + Ψ (u; x) , ∀(u, v) ∈ V × R n . (4.18)
The stochastic step sequence ( γ k+1 ) k≥0 is then defined as

γ k+1 = min x∈Ω γ k+1 ∧ χ k x , ∀k ≥ 1 , (4.19)
with γ 0 = γ 0 . For any x in Ω and k = 1, . . . , M -1, the adaptive explicit Störmerverlet (ASV) scheme then reads as:

U k+1/2 (x) = U k (x) + γ k+1 2 V k (x) , (4.20)
V k+1 (x) = 1 -β x 1 + β x V k (x) + γ k+1 1 + β x L k+1/2 (x) + √ f x 1 + β x ∆W k+1 (x) , (4.21)
U k+1 (x) = U k+1/2 (x) + γ k+1 2 V k+1 (x) , (4.22)
with β x = f x γ k+1 /4. For any x fixed in Ω and following § 2.1, the increment of the normalized Wiener process {W (r, x), r 0} between time instances r k and r k+1 is finally defined as:

∆W k+1 (x) = W k+1 (x) -W k (x) = ( γ k+1 ) 1/2 Ξ x (θ k+1 ) . (4.23)
It is worth noticing that such a definition allows spatial dependencies to be introduced while solving the family of ISDEs and is therefore well-adapted to the generation of the class of random field prior models.

Remark 4. From Eq. (4.13), it is seen that setting a large value of τ amounts to consider a very slowly decreasing deterministic counterpart in the step sequence. In effect, this allows for a fast exploration of the admissible space S, but might slow down the convergence towards the stationary regime when the neighborhood of ∂S is reached with moderate to large probabilities. Conversely, retaining a small value of τ decreases the number of stochastic refinements (in the transient regime in particular) and yields a better (asymptotic) approximation of the invariant measure [START_REF] Lamberton | Computation of the Invariant Distribution of a Diffusion[END_REF].

Numerical experiments.

In this section, we illustrate the proposed approach through three applications. The first two applications do not address the random field case (the spatial indexation is then temporarily dropped in § 5.1 and § 5.2) and specifically investigate the behavior of the algorithm for two different kinds of supports. An application to the random field case is finally provided in § 5.3, where a homogeneous system of first-order m.p.d.f. is considered. Assuming that the random field is sampled at N non-overlapping points x (1) , . . ., x (N ) , let Conv be the following function:

Conv(M c , x (i) ; ) = 1 Mc k=M0 γ k Mc k=M0 γ k U k-1 (x (i) ) 2 , ∀M c 1 , (5.1)
where M 0 1 may be selected to speed up the convergence by removing the transient part in the ergodic estimator. For an homogeneous system of first-order m.p.d.f. and when a reference value m ref 2 for the second-order moment of U is known, the convergence is characterized by using the following error function M c → ErrConv(M c ; ):

ErrConv(M c ; ) = max 1 i N |Conv(M c , x (i) ; ) -m ref 2 |/m ref 2 ∀M c 1 , (5.2)
with Conv given by Eq. (5.1). All routines are programmed in C, making use of the FFTW library [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF] for FFT computations. The underlying random generator is the 64-bit Mersenne Twister generator from Matsumoto and Nishimura [START_REF] Matsumoto | Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator[END_REF]. of the admissible space. Specifically, the adaptive step is refined 36,410 times in average, out of one million iterations (the average being estimated from a set of 1, 000 independent realizations). Note that the number of refinements does not change that much from one realization to another. A sample path for k → (U k 1 , U k 2 , γ k ) for = 0.04 and = 0.01 is displayed in Fig. 5.6 for a large integration range, hence showing the incidence of this parameter. The graphs of the p.d.f. p U1 of random variable U 1 (in the stationary regime) are finally shown in Figs. 5.7 and 5.8 for different values of parameter , hence illustrating the convergence of the sequence as ↓ 0.

Application 2.

5.2.1. Problem setting. This second illustration is concerned with a regularization on the matrix set M + q (R). Specifically, we consider sampling from the random matrix ensemble SG + constructed in [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF] (for which a reference generator is available). Let [G] ∈ SG + and n = q(q + 1)/2. The target p.d.f. p of [G] is then defined as (see the reference above)

p([G]) = c exp {-ψ([G])} ∀[G] ∈ M + q (R) , (5.6)
where the potential function ψ :

M + q (R) → R is given by ψ([G]) = -log 1 M + q (R) ([G]) -(λ -1) log (det([G])) + q -1 + 2λ 2 tr([G]) (5.7)
for all [G] in M + q (R). In Eq. (5.7), the scalar parameter λ > 0 can be shown to control the level of statistical fluctuations of the random matrix [G]. Let Mat : R n ⊃ S → M + q (R) be the so-called vector-matrix mapping such that for all

[G] in M + q (R), ∃ u ∈ S such that [G] = [Mat(u)] with [Mat(u)] =       u 1 u 2 . . . u q(q-1)/2+1 u 3 . . . sym. . . . u n       . ( 5.8)
The admissible space S ⊂ R n is thus defined as:

S = u ∈ R n | [Mat(u)] ∈ M + q (R) .
(5.9)

It follows that the potential function Ψ : S → R can be written as

Ψ(u) = -log (1 S (u)) + g(u) , (5.10)
with g : S → R such that

g(u) = -(λ -1) log (det([Mat(u)])) + q -1 + 2λ 2 q k=1 u k(k+1)/2 .
(5.11)

The approximation Ψ of Ψ, parametrized by 0 < 1, is then given by: Ψ (u) = -log (1 S (u)) + g(u) , ∀u ∈ S .

(5.12)

For numerical purposes, we let q = 2 (n = 3), λ = 3.67 and f = 10. This value of λ corresponds to a large level of statistical fluctuations, that is δ = 0.6.

Numerical results.

The graph of the error function n → ErrConv(n) is first displayed in Fig. 5.9 for = 0.01, γ 0 = 2 -6 and τ = 10 6 . For this sample path, 223 refinements are performed, out of 80, 000 iterations. It is seen that the convergence towards the stationary solution is reached for M c 60, 000, with a relative error that is less than 6%. For any k 1, we further introduce the R + * -valued random variable ρ k such that

ρ k = min S [Mat(U k )] ,
(5.13) with S([A]) the spectrum of [A] ∈ M S q (R). At the k-th iteration, ρ k is used to quantitatively characterize in part the distance to the boundary ∂S. Sample paths of k → ρ k and k → γ k are displayed in Fig. 5.10. As for the first application, it is readily observed that the step is properly refined whenever the current state gets closer to the boundary of the admissable space S. Further, it is seen that the probability of reaching the boundary ∂S turns out to be small, with only 1.6% of the steps being refined for = 0.01: this property stems from the repulsion of the potential function in the neighbourhood of ∂S, as discussed in § 3.2. The convergence of the proposed algorithm is illustrated in Figs. 5.11, 5.12 and 5.13, where the firstand second-order marginal p.d.f. estimated with both the proposed and reference generators are displayed. Whereas an overall good match is observed, regardless of the random variables involved, some small discrepancies still remain for the second order marginal distributions. This fact may be due to an undersampling of the joint distributions, which are estimated from a set of 60, 000 independent realizations.

Application 3.

In this section, we finally exemplify the algorithm on the M + q (R)-valued random field whose first-order marginal p.d.f. is defined through the potential function given by Eq. (5.12). For simplicity, a two-dimensional domain Ω = (]0, 10[) 2 is considered and the correlation function [R Ξ ] of the stationary Gaussian for all 1 i n, ij c being interpreted as the correlation length of the Gaussian random field {Ξ i (x), x ∈ Ω} along the direction defined by the unit vector e j . The choice of such a Matérn-type (squared exponential) correlation function follows from its wide use in applications (especially in geophysics; see [START_REF] Apanasovich | A Valid Matérn Class of Cross-Covariance Functions for Multivariate Random Fields With Any Number of Components[END_REF] for instance), as well as from the properties (e.g. infinite differentiability) inherited by the associated Gaussian field. Here, the generation of the stationary R-valued Gaussian random fields is performed by using the circulant embedding method introduced in [START_REF] Dietrich | A Fast and Exact Simulation of Stationary Gaussian Processes Through Circulant Embedding of the Covariance Matrix[END_REF] and [START_REF] Wood | Simulation of Stationary Gaussian Processes in [0; 1] d[END_REF] (see also [START_REF] Chan | Simulation of stationary Gaussian vector fields[END_REF]; see [START_REF] Gneiting | Fast and exact simulation of large Gaussian lattice systems in R 2 : Exploring the limits[END_REF] [49] for extensions). For illustration purposes, the random field is sampled over a regular grid with 8 points along each direction and the simulation is performed with the following parameters: ij c = 2 for 1 i 3 and 1 j 2; f = 10; τ = 10 6 ; γ 0 = 2 -8 . A convergence analysis (similar to those performed in the preceding sections; see e.g. § 5.2.2) shows that an independent realization is obtained after 60, 000 iterations. The CPU time required to generate one realization of the random field is 15.98 seconds on a 2.8 GHz single core processor -note that this CPU time mainly depends on the algorithms that are used for simulating a R-valued centered homogeneous Gaussian random field (which has to be generated at least n × r * times; see e.g. [START_REF] Gneiting | Fast and exact simulation of large Gaussian lattice systems in R 2 : Exploring the limits[END_REF] and the references therein for a discussion about computational complexity). Since numerical evidences show that the algebraic form of the input correlation functions is indeed preserved through the nonlinear transformation [START_REF] Guilleminot | Prior modeling and generator for non-Gaussian tensor-valued random fields with symmetry properties. Application to elastic random media[END_REF], the estimated correlation function (along a given direction) can be fitted to a squared exponential one by a least-square algorithm, hence providing estimates for the correlation lengths associated with the generated fields. The comparison between the estimated correlation function and its least-square fit in the direction specified by e 1 is shown in Fig. 5.14 (left panel) for the random field {U 1 (x), x ∈ Ω} (the statistical estimator being computed from a set of 200 realizations). The graphs of the fitted correlation functions (along the same direction) for the three generated random fields {U 1 (x), x ∈ Ω}, {U 2 (x), x ∈ Ω} and {U 3 (x), x ∈ Ω} can be found in Fig. 5.14 (right panel) as well. The correlation length along the considered direction for the random field {U 1 (x), x ∈ Ω} (resp. {U 2 (x), x ∈ Ω} and {U 3 (x), x ∈ Ω}) is found to be 1.84 (resp. 1.88 and 1.94).

Conclusion.

In this work, we addressed the construction of a generator for a large class of non-Gaussian random fields with values in some subset of R n . The numerical strategy involves two main ingredients. The first one is the definition, through a regularization procedure, of a sequence of auxiliary random fields converging towards the target random field. The second ingredient is the derivation of a discretization algorithm which is defined by combining an extension of an adaptivity procedure recently proposed in the literature with an explicit Verlet scheme. The scheme thus obtained allows refining the integration step on-the-fly, hence preventing the scheme from diverging as the current state goes closer to the boundary of the admissible set. The proposed strategy was finally exemplified on three examples that illustrate both the adaptivity procedure and the convergence of the algorithm.

Fig. 4 . 1 .

 41 Fig. 4.1. Convergence analysis: graph of the mapping Mc → ErrConv(Mc) for ∆r = 1.5×10 -3 (left panel) and ∆r = 0.5 × 10 -3 (right panel). Solid line: explicit Euler scheme. Dashed line: Störmer-Verlet scheme.

Fig. 5 . 2 .

 52 Fig. 5.2. Convergence analysis: graph of the mapping Mc → ErrConv(Mc) for several values of parameter , γ 0 = 2 -8 and τ = 10. Black line: = 0.01; red line: = 0.04; blue line: = 0.07; magenta line: = 0.1.

Fig. 5 . 3 .

 53 Fig. 5.3. Illustration: sample path of k → (U k 1 , U k 2 ) for k ∈ 5 × 10 5 ; 6 × 10 5 .

Fig. 5 . 4 .

 54 Fig. 5.4. Illustration of the adaptivity procedure: sample path of k → γ k for k ∈ 5×10 5 ; 6×10 5 (left panel) and k ∈ 5.013 × 10 5 ; 5.017 × 10 5 (right panel).

Fig. 5 . 5 .

 55 Fig. 5.5. Illustration of the adaptivity procedure: sample path of k → (U k 1 , U k 2 , γ k ) for k ∈ 5 × 10 5 ; 6 × 10 5 .

Fig. 5 . 6 .

 56 Fig. 5.6. Illustration: sample path of k → (U k 1 , U k 2 ) for k ∈ 5 × 10 5 ; 1 × 10 6 for = 0.04 (left panel) and = 0.01 (right panel).

Fig. 5 . 7 .

 57 Fig. 5.7. Convergence analysis: graph of the p.d.f. U → p U 1 (U ) for = 0.1 (left panel) and = 0.07 (right panel): theoretical p.d.f. (thick line) and ergodic estimator (thin line).

Fig. 5 . 8 .

 58 Fig. 5.8. Convergence analysis: graph of the p.d.f. U → p U 1 (U ) for = 0.04 (left panel) and = 0.01 (right panel): theoretical p.d.f. (thick line) and ergodic estimator (thin line).

Fig. 5 . 9 .

 59 Fig. 5.9. Convergence analysis: graph of the mapping Mc → ErrConv(Mc) for = 0.01, γ 0 = 2 -6 and τ = 10 6 .

Fig. 5 . 10 .

 510 Fig. 5.10. Illustration of the adaptivity procedure: sample paths of k → ρ k (red circles) -left vertical axis -and k → γ k (black solid line) -right vertical axis -for = 0.01, γ 0 = 2 -6 and τ = 10 6 .

Fig. 5 . 11 .

 511 Fig. 5.11. Graph of p.d.f. u → p U 1 (u) (left panel) and u → p U 2 (u) (right panel): reference algorithm (solid line) and proposed approach (dashed line).

Fig. 5 . 12 .

 512 Fig. 5.12. Graph of p.d.f. u → p U 3 (u) (left panel) and contour plot of (u, v) → p U 1 U 2 (u, v) (right panel): reference algorithm (solid/thick line(s)) and proposed approach (dashed/thin line(s)).

Fig. 5 . 13 .

 513 Fig. 5.13. Contour plots of (u, v) → p U 1 U 3 (u, v) (left panel) and (u, v) → p U 2 U 3 (u, v) (right panel): reference algorithm (thick lines) and proposed approach (thin lines).

Fig. 5 . 14 .

 514 Fig. 5.14. Plot of the correlation function along the direction defined by e 1 . Left panel: correlation function associated with {U 1 (x), x ∈ Ω} (circle marker: estimated values; solid line: leastsquare fit). Right panel: fitted correlation functions associated with {U 1 (x), x ∈ Ω}, {U 2 (x), x ∈ Ω} and {U 3 (x), x ∈ Ω}.

In order to avoid any possible confusion, it is recalled here that cx denotes a normalization constant whose value may vary from line to line.

For some particular cases, one may recover a free case construction by using an alternative algebraic construction based on the exponential map; see[START_REF] Guilleminot | Prior modeling and generator for non-Gaussian tensor-valued random fields with symmetry properties. Application to elastic random media[END_REF] for instance.

Acknowledgments. This research was funded by the French National Research Agency (ANR) under TYCHE (ANR-2010-BLAN-0904) and MOSAIC (ANR-12-JS09-0001-01) contracts.

This work was supported by the French National Research Agency (ANR) under TYCHE and MOSAIC contracts.

5.1. Application 1.

Problem setting.

In this first application, we consider a two-dimensional case, n = 2, and define the target p.d.f. p t as

in which the potential function is given by

where the regularized indicator function 1 S is defined by Eq. (3.17). The graph of u → 1 S (u) is displayed in Fig. 5.1 for = 0.01 (left panel) and = 0.04 (right panel).


A Pasko

F Mazaleyrat

L K Varga

P Stamenov

J M D Coey

High-field magnetization behavior of Mn-Al-C alloys

Keywords: Magnetic anisotropy, manganese-aluminum alloys, permanent magnets, spark plasma sintering

Dense Mn-Al-C permanent magnets have been produced using melt-spinning and spark plasma sintering techniques. The samples have been magnetically characterized from 4 K to 300 K in a field of up to 14 T. The ahysteretic parts of magnetization loops were analyzed within a model for the approach to saturation which allows for possible crystallographic texture. The saturation magnetization and magnetocrystalline anisotropy energy have been determined taking into account the estimated fraction of ferromagnetic phase. The obtained magnetic parameters have been compared with available theoretical and experimental data.

I. INTRODUCTION

Recent years have seen a renewed interest in the Mn-Al alloys, first reported by Kono [START_REF] Kono | On the ferromagnetic phase in manganese-aluminum system[END_REF] and Koch et al. [START_REF] Koch | New material for permanent magnets on a base of Mn and Al[END_REF], as promising candidates for permanent magnet applications. These materials belong to the large family of ferromagnetic phases with L1 0 ordered structure and possess a high uniaxial magnetocrystalline anisotropy arising from a breakdown of the cubic crystal symmetry to a tetragonal one. The escalating prices of rare earths present an industrial opportunity for new rare-earth-free magnets such as the Mn-Al-based alloys having good technological properties (light, strong, machinable) and low cost. With proper optimization, they would be capable of filling the gap between cheap hexagonal ferrites and highperformance Nd-Fe-B magnets [START_REF] Coey | Permanent magnets: Plugging the gap[END_REF]- [START_REF] Kramer | Prospects for non-rare earth permanent magnets for traction motors and generators[END_REF]. The ferromagnetic tetragonal τ-phase is metastable and can be obtained in the Mn 0.5+x Al 0.5-x alloy (x = 0-0.1) alloy by quenching the hightemperature hexagonal ε-phase with subsequent annealing at 673-973 K. In order to stabilize the ferromagnetic phase, carbon has been introduced into the crystal lattice, which also improves coercivity due to fine Mn 3 AlC precipitates [START_REF] Ohtani | Magnetic properties of Mn-Al-C permanent magnet alloys[END_REF]. The mechanism of phase formation comprises either an atomic ordering followed by a martensitic (diffusionless) phase transition or a massive (thermally activated compositioninvariant) structural transformation [START_REF] Wiezorek | Grain boundary mediated displacive-diffusional formation of τ-phase MnAl[END_REF]. The magnetic hysteresis of Mn-Al-C alloys is therefore sensitive to the microstructure and lattice defects developed during the τphase formation and is greatly influenced by the manufacturing route.

Textured Mn-Al-C samples with optimized functional properties were produced in Japan by hot extrusion [START_REF] Ohtani | Magnetic properties of Mn-Al-C permanent magnet alloys[END_REF]. Various methods have been proposed to make isotropic magnets: melt-spinning [START_REF] Sakka | Rapid quenching and properties of hard magnetic materials in MnAl-X (X = Ti, Cu, Ni, C, B) systems[END_REF]- [START_REF] Chaturvedi | A comparison of τ-MnAl particulates produced via different routes[END_REF], gas atomization [START_REF] Chaturvedi | A comparison of τ-MnAl particulates produced via different routes[END_REF]- [START_REF] Lee | Effect of mechanical milling and heat treatment on the structure and magnetic properties of gas atomized Mn-Al alloy powders[END_REF], mechanical alloying [START_REF] Obi | Magnetic and structural properties of heattreated high-moment mechanically alloyed MnAlC powders[END_REF]. In addition, a ball milling step was introduced after quenching [START_REF] Zeng | Structural and magnetic properties of nanostructured Mn-Al-C magnetic materials[END_REF] or rapid solidification [START_REF] Chaturvedi | A comparison of τ-MnAl particulates produced via different routes[END_REF], [START_REF] Lee | Effect of mechanical milling and heat treatment on the structure and magnetic properties of gas atomized Mn-Al alloy powders[END_REF] for experiments with nanostructured materials. Despite relatively high magnetocrystalline anisotropy, it has not been easy to develop coercivity in the Mn-Al-C alloys. In this paper we report on the high-field magnetic characterization of the Mn-Al-C bulk samples produced by spark plasma sintering (SPS). The idea behind SPS is simultaneous application of heat and pressure for more efficient consolidation of powders. Unlike hot pressing, SPS uses a resistive heating by a pulsed electric current passing directly through the sample [START_REF] Munir | The effect of electric field and pressure on the synthesis and consolidation of materials: A review of the spark plasma sintering method[END_REF]. When applied to the synthesis of dense Mn-Al-C magnets, this technique can be advantageous, by combining in one step the heat treatment necessary for the τ-phase formation and the powder sintering itself, while avoiding excessive grain coarsening and decomposition of the metastable phases. Magnetic measurements in a high-field region enable us to determine the saturation magnetization and magnetocrystalline anisotropy energy, within a model for the approach to saturation.

II. EXPERIMENTAL

Ternary Mn 54 Al 44 C 2 and Mn 49 Al 49 C 2 alloys were prepared from pure elements by induction melting in a water-cooled copper crucible. The ingots were then cast into ~ 3 mm wide and ~ 25 µm thick ribbons using the melt-spinning technique under protective atmosphere at a wheel speed of 25 m•s -1 [START_REF] Fazakas | Preparation of nanocrystalline Mn-Al-C magnets by melt spinning and subsequent heat treatments[END_REF]. The ribbons appeared brittle and were subsequently ground by hand with pestle and mortar. The obtained powders were compacted in a graphite die with internal diameter of 8 mm and sintered in a Sumitomo Dr. Sinter Lab SPS machine under vacuum with an applied uniaxial pressure of 70 MPa. Two thermal regimes were chosen with a maximum temperature of 923 K or 1223 K, heating/cooling rate of 6000 K•s -1 and a dwell time of 300 s. At lower temperatures sintering is accompanied by the ε→τ phase transition; at higher temperatures the material is sintered in the stable ε-phase state, and transformation to the τ-phase occurs on cooling [START_REF] Pasko | Spark plasma sintering of Mn-Al-C hard magnets[END_REF]. Crystal structures were determined by a PANalytical X'Pert Pro X-ray diffractometer (XRD) equipped with an X'Celerator detector. Magnetic measurements were carried out using a Quantum Design PPMS instrument in a temperature range of 4-300 K and in an applied field of up to 14 T.

III. RESULTS

According to the phase diagram of the Mn-Al system [START_REF] Shukla | Thermodynamic assessment of the Al-Mn and Mg-Al-Mn Systems[END_REF], two equilibrium paramagnetic phases can be formed in the Manuscript received March 7, 2014. Corresponding author: A. Pasko (email: pasko@satie.ens-cachan.fr).

Digital Object Identifier inserted by IEEE compositional range concerned: the Mn-rich cubic β-phase and the Al-rich trigonal γ 2 -phase. In addition, two metastable phases may appear in non-equilibrium conditions: disordered hexagonal ε-phase and ordered tetragonal τ-phase. Stable at higher temperatures, the paramagnetic ε-phase can be preserved on quenching and converted into the ferromagnetic τ-phase by a composition-invariant structural transformation. The resulting phase constitution and related magnetic properties are therefore strongly dependent on the heat treatment regime.

FIG. 1 HERE

Magnetic hysteresis loops of the Mn 49 Al 49 C 2 and Mn 54 Al 44 C 2 alloys, sintered at 923 K or 1223 K, at room temperature in a field of up to 5 T are compared in Fig. 1. In the Mn-Al system an excess of Mn is usually required to stabilize the ferromagnetic τ-phase. Hence, the Mn 54 Al 44 C 2 samples show higher magnetization, in spite of partial occupation of Al sites by Mn atoms with antiferromagnetic Mn-Mn coupling. However, the Mn 49 Al 49 C 2 samples have higher coercivity due to the presence of a parasitic phase. The increase of sintering temperature leads to denser material with higher content of the τ-phase, but adversely affects its coercivity. The XRD phase analysis confirms the correlation between the τ-phase mass fraction and the saturation magnetization value.

FIG. 2 HERE

In view of its better combination of magnetic properties, the Mn 54 Al 44 C 2 alloy, sintered at 923 K, was selected for further detailed high-field characterization. Fig. 2 shows the magnetization as a function of temperature, measured in a constant field of 14 T, together with a quadratic in T fitting curve. The specific moment extrapolated to 0 K is σ s = 114 A•m 2 •kg -1 . According to the comparative analysis of structural and magnetization data, this sample contains approximately 85 % of the ferromagnetic phase [START_REF] Pasko | Spark plasma sintering of Mn-Al-C hard magnets[END_REF]. Taking into account the τ-phase density 5086 kg•m -3 , we can estimate the equivalent single-phase saturation magnetization as M s = 0.68 MA•m -1 (μ 0 M s = 0.86 T).

FIG. 3 HERE

Magnetic hysteresis loops of this alloy in a field of up to 14 T at 4 K and 300 K are presented in Fig. 3. The large reversible segments of these curves have been analyzed using a model describing the approach to saturation in textured materials. The total magnetic moment as a function of applied field is calculated by averaging the numerical solution of the torque equation over possible orientations of anisotropy axes. The orientation distribution is described by a simple expression with minimum number of free parameters. The resulting model function is used to fit directly, using nonlinear least-squares regression, the experimental data.

FIG. 4 HERE

Fitting of the ahysteretic part of magnetization loops for the Mn 54 Al 44 C 2 alloy, sintered at 923 K, is shown in Fig. 4. The approach to saturation at 300 K corresponds well to isotropic material. However, the magnetization curve at 4 K can be better fitted with an assumption of anisotropic texture. The specific moment is σ s = 114 A•m 2 •kg -1 at 4 K and σ s = 98 A•m 2 •kg -1 at 300 K; the anisotropy field is H a = 4.0 MA•m -1 (μ 0 H a = 5.0 T) at 4 K and H a = 3.1 MA•m -1 (μ 0 H a = 3.9 T) at 300 K. Using again the estimated mass fraction of the τ-phase, we found the corrected values of saturation magnetization M s and magnetocrystalline anisotropy constant K 1 presented in Table 1. The presence of light crystallographic texture can be attributed to the uniaxial pressure applied to the sample during SPS.

TABLE 1 HERE

IV. DISCUSSION Electronic structure and magnetic properties of the equiatomic Mn-Al compound were computed using density functional techniques (LMTO, FLAPW) in a number of works. The magnetocrystalline anisotropy energy of the τphase depends on its c ∕ a axial ratio so that a tetragonal distortion from a cubic structure stabilizes the system in terms of the Jahn-Teller effect [START_REF] Sakuma | Electronic structure and magnetocrystalline anisotropy energy of MnAl[END_REF]. The ferromagnetic L1 0 phase may be interpreted to originate from a nonmagnetic B2 phase through a two-step process: an electronic transition in the cubic lattice due to spin polarization, followed by a structural deformation into the tetragonal system [START_REF] Kurtulus | Electronic structure, chemical bonding, and spin polarization in ferromagnetic MnAl[END_REF]. Different first principles calculations give similar values for the magnetocrystalline anisotropy constant: 1.5 MJ•m -3 [START_REF] Sakuma | Electronic structure and magnetocrystalline anisotropy energy of MnAl[END_REF], 1.525 MJ•m -3 [START_REF] Park | Saturation magnetization and crystalline anisotropy calculations for MnAl permanent magnet[END_REF]. The saturation magnetization (assuming a density of 5008 kg•m -3 ) is obtained as 0.806 MA•m -1 [START_REF] Park | Saturation magnetization and crystalline anisotropy calculations for MnAl permanent magnet[END_REF].

Magnetic properties of polycrystalline Mn-Al and Mn-Al-C alloys within the τ-phase existence domain were experimentally studied using the singular point detection technique by Pareti et al. [START_REF] Pareti | Magnetic anisotropy of MnAl and MnAlC permanent magnet materials[END_REF]. It was confirmed that uniaxial magnetocrystalline anisotropy of this material is well described with only a second-order parameter K 1 . The anisotropy field H a at room temperature increases with increasing Mn content and decreases with increasing C content, the variation being from 3.1 MA•m -1 to 4.2 MA•m -1 . It is remarkable that the saturation magnetization M s exhibits an inverse dependence on the chemical composition, so that the anisotropy constant K 1 = μ 0 M s H a ∕ 2 remains indeed almost constant in the considered compositional range. From figures in [START_REF] Pareti | Magnetic anisotropy of MnAl and MnAlC permanent magnet materials[END_REF] these parameters for Mn 54 Al 44 C 2 alloy can be estimated as follows: M s = 0.6 MA•m -1 , K 1 = 1.0 MJ•m -3 , in good agreement with our data.

V. CONCLUSION

Dense Mn-Al-C permanent magnets have been produced from rapidly solidified alloys by spark plasma sintering without additional annealing. The samples have been magnetically characterized in a field sufficient to saturate this material. The ahysteretic parts of magnetization loops were analyzed within a model for the approach to saturation which takes into account possible crystallographic texture. The measured magnetocrystalline anisotropy energy is close to the available predictions of the density functional theory (in T → 0 K limit) and consistent with other experimental studies. For many materials the observed coercivity amounts to ~ 10 % of the anisotropy field. Isotropic Mn-Al-C magnets have not yet attained this level.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Magnetic hysteresis loops of the Mn54Al44C2 alloy, sintered at 923 K, in a field of up to 14 T at different temperatures.

Fig. 1 .

 1 Fig. 1. Magnetic hysteresis loops of the Mn49Al49C2 and Mn54Al44C2 alloys sintered at different temperatures. The samples with higher coercivity have lower specific moment because of smaller content of the τ-phase.

Fig. 2 .

 2 Fig. 2. Specific moment as a function of temperature for the Mn54Al44C2 alloy sintered at 923 K. The experimental data follow well a quadratic (without linear term) law.


A E Gheribi

D Corradini

L Dewan

P Chartrand

C Simon

P A Madden

M Salanne

Prediction of the thermophysical properties of molten salt fast reactor fuel from first-principles

Molten fluorides are known to show favorable thermophysical properties which make them good candidate coolants for nuclear fission reactors. Here we investigate the special case of mixtures of lithium fluoride and thorium fluoride, which act both as coolant and fuel in the molten salt fast reactor concept. By using ab initio parameterized polarizable force fields, we show that it is possible to calculate the whole set of properties (density, thermal expansion, heat capacity, viscosity and thermal conductivity) which are necessary for assessing the heat transfer performance of the melt over the whole range of compositions and temperatures. We then deduce from our calculations several figures of merit which are important in helping the optimization of the design of molten salt fast reactors.

I. INTRODUCTION

The molten salt fast reactor (MSFR) is currently the most advanced concept of nuclear fission reactor involving the use of molten salts both as coolant and as fuel [START_REF] Locatelli | Generation IV nuclear reactors: Current status and future prospects[END_REF] . A great deal of research has recently been devoted to the optimization of its neutronic [START_REF] Mathieu | The thorium molten salt reactor: moving on from the MSBR[END_REF][START_REF] Aufiero | Calculating the effective delayed neutron fraction in the molten salt fast reactor: Analytical, deterministic and monte carlo approaches[END_REF] , safety [START_REF] Brovchenko | Design-related studies for the preliminary safety assessment of the molten salt fast reactor[END_REF] , recycling [START_REF] Delpech | Reactor physics and reprocessing scheme for innovative molten salt reactor system[END_REF][START_REF] Hamel | Reduction process of uranium (IV) and uranium (III) in molten fluorides[END_REF][START_REF] Gibilaro | Electrochemical extraction of europium from molten fluoride media[END_REF][START_REF] Taxil | Lanthanides extraction processes in molten fluoride media: Application to nuclear spent fuel reprocessing[END_REF][START_REF] Gibilaro | Direct electrochemical reduction of solid uranium oxide in molten fluoride salts[END_REF] and thermodynamic [START_REF] Benes | Thermodynamic properties and phase diagrams of fluoride salts for nuclear applications[END_REF][START_REF] Benes | Thermodynamic assessment of the LiF-CeF3-ThF4 system: Prediction of PuF3 concentration in a molten salt reactor fuel[END_REF][START_REF] Capelli | Thermodynamic investigation of the LiF-ThF4 system[END_REF] characteristics. Its advantages are an intrinsic safety of the core, due to very negative feedback coefficients, a high working temperature which ensures a good efficiency, and the possibility to perform online refuelling and on-site treatment. But of particular interest is its high compositional flexibility. In particular, our precognized liquid fuel is a mixture of lithium fluoride and thorium fluoride (LiF-ThF 4 ). Being able to develop fission reactors based on the thorium cycle ensures a good availability of the resources for many countries. From the point of view of basic research, a consequence of the MSFR compositional flexibility is the need to generate substantial databases providing the variation of the properties of the melt with composition and temperature.

Due to the high temperature and to the corrosive character of molten fluorides, measuring their thermophysical properties is not a straightforward task, and the experimental results are generaly limited to pure salts or binary mixtures at the eutectic composition [START_REF] Benes | Molten salt reactor fuel and coolant[END_REF][START_REF] Benes | A comprehensive study of the heat capacity of CsF from T = 5 K to T = 1400 K[END_REF][START_REF] Kohklov | Heat capacity and thermalconductivity of molten ternary lithium, sodium, potassium, and zirconium fluorides mixtures[END_REF] . An efficient alternative is to use computer simulations. In particular, molecular dynamics (MD) simulations are able to provide all the necessary properties, such as the density, the heat capacity, the viscosity and the thermal conductivity. The necessary ingredients are a good model for the ionic interactions, and its associated parameters. In recent years, we have developed from first-principles, i.e. without involving any empirical information during the parameterization procedure, such interaction potentials for molten fluorides, in the framework of the polarizable ion model [START_REF] Madden | From first-principles to material properties[END_REF][START_REF] Salanne | Polarization effects in ionic solids and melts[END_REF][START_REF] Salanne | Including many-body effects in models for ionic liquids[END_REF] . They were tested against many physical properties, as well as on diffraction and spectroscopy experiments. An excellent agreement was obtained for the full set of data [START_REF] Heaton | Turq A first-principles description of liquid BeF2 and its mixtures with LiF: 1. Potential development and pure BeF2[END_REF][START_REF] Pauvert | In situ experimental evidence for a non-monotonous structural evolution with composition in the molten LiF-ZrF4 system[END_REF][START_REF] Levesque | Structure and dynamics in yttrium-based molten rare earth alkali fluorides[END_REF] . In particular, for LiF-ThF 4 eutectic mixture, which corresponds to a ThF 4 mole fraction x = 0.22, the largest deviation from experimental data observed was 5 % for the density [START_REF] Dewan | Molecular dynamics simulation of the thermodynamic and transport properties of the molten salt fast reactor fuel LiF-ThF4[END_REF] . Such a good predictive power from MD simulations has not been obtained yet for water, for example. It is thus possible to perform simulations under a large variety of thermodynamic conditions in order to feed the databases. In this work, we focus on the LiF-ThF 4 binary system, for which we calculate all the relevant physical properties over the whole range of compositions and temperatures. This allows us to map the variations of the thermophysical figures of merit which are necessary for the engineering of MSFRs.

II. MODEL AND METHODS

The pair additive form of the interionic potential includes charge-charge, charge-dipole, dipole-charge, dipole-dipole, repulsion, and two dispersion interaction terms, as follows :

U ij (r ij ) = q i q j r ij + q i r ij • µ j r 3 ij f ij 4 (r ij ) - µ i • r ij q j r 3 ij f ji 4 (r ij ) + µ i • µ j r 3 ij - 3(r ij • µ i )(r ij • µ j ) r 5 ij +B ij exp (-α ij r ij ) - C ij 6 r 6 ij f ij 6 (r ij ) - C ij 8 r 8 ij f ij 8 (r ij ) (1)
where q i and µ i are the charge and dipole moment of particle i, respectively, and f ij n the damping function for short-range correction of interactions between charge and dipole and dispersion interactions [START_REF] Tang | An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients[END_REF] :

f ij n (r ij ) = 1 -c ij n e -b ij n ×r ij n k=0 (b ij n × r ij ) k k! . ( 2
)
{B ij , α ij , C ij 6 , C ij 8 , b ij n , c ij n
} are a set of parameters which were determined from first-principles calculations in our previous works [START_REF] Dewan | Molecular dynamics simulation of the thermodynamic and transport properties of the molten salt fast reactor fuel LiF-ThF4[END_REF] . The dipole moment induced on each ion is a function of the polarizability and electric field on the ion caused by the charges and the dipole moments of all the other ions. The instantaneous dipole moment is determined self-consistently at every time step by minimization of the total energy using the conjugate gradient method. The charge-charge, charge-dipole, and dipoledipole contributions to the potential energy and forces on each ion are evaluated under periodic boundary conditions by using the Ewald summation technique [START_REF] Aguado | Ewald summation of electrostatic multipole interactions up to the quadrupolar level[END_REF] .

A first series of MD simulations was performed in the N P T ensemble, where N is the particle number included in the simulation cell, P is the pressure which was always set to 1 atm and T is the temperature. The barostat and thermostat relaxation times were both set to 10 ps. This allowed us to determine the equilibrium density at a given thermodynamic point. Then a second series of simulations was performed at fixed density, i.e. in the N V T ensemble where V is the volume. The same relaxation time was used for the thermostat. The simulation timestep was 0.5 fs, and the total simulation time was 0.25 ns and 7.5 ns for the N P T and N V T simulations, respectively. All the calculation conditions are summarized in the supplementary information; in total the system was simulated at 62 thermodynamic points with this procedure.

III. RESULTS AND DISCUSSION

Our ab initio-parameterized interaction potentials have been validated against many experimental data for similar systems. More precisely, they are able to reproduce with very high accuracy the measured EXAFS spectra [START_REF] Pauvert | In situ experimental evidence for a non-monotonous structural evolution with composition in the molten LiF-ZrF4 system[END_REF] (LiF-ZrF 4 , NaF-ZrF 4 and KF-ZrF 4 ), the Raman spectra 26 (LiF-BeF 2 ), the X-ray diffraction 27 (LiF-BeF 2 ), the electrical conductivity 21 (LiF-NaF-ZrF 4 and LiF-YF 3 ), the diffusion coefficients [START_REF] Levesque | Structure and dynamics in yttrium-based molten rare earth alkali fluorides[END_REF][START_REF] Sarou-Kanian | Diffusion coefficients and local structure in basic molten fluorides: in situ NMR measurements and molecular dynamics simulations[END_REF] (LiF-KF and LiF-YF 3 ), the heat capacity [START_REF] Benes | A comprehensive study of the heat capacity of CsF from T = 5 K to T = 1400 K[END_REF] (CsF) and the viscosity [START_REF] Salanne | A first-principles description of liquid BeF2 and its mixtures with LiF: 2. network formation in LiF-BeF2[END_REF] (LiF-BeF 2 ). For the thermal conductivity, this quantity is very difficult to measure experimentally and our simulations have shown an excellent agreement in the case of molten chlorides [START_REF] Ohtori | Calculations of the thermal conductivities of ionic materials by simulation with polarizable interaction potentials[END_REF] , for which accurate measurements could be made by using a method based on forced Rayleigh scattering. In the case of LiF-ThF 4 , in our previous work [START_REF] Dewan | Molecular dynamics simulation of the thermodynamic and transport properties of the molten salt fast reactor fuel LiF-ThF4[END_REF] we have evaluated the density, the viscosity and the electrical conductivity at the eutectic composition. Our MD results showed a very good agreement with the corresponding experimental data [START_REF] Benes | Molten salt reactor fuel and coolant[END_REF] . In the present work, we can therefore consider our simulations to provide us accurate predictions.

In order to compare the performances of a series of coolants, it is useful to define generalized heat-transfer metrics. Several figure of merits (FOM) have been proposed by Bonilla to evaluate the properties of molten fluorides as coolants for nuclear reactors 25 . They take the following general form:

FOM = η s β t ρ u C v p,m λ w r ( 3
)
where η is the viscosity, β the thermal expansion, ρ the density, C p,m the massic heat capacity at fixed pressure and λ the thermal conductivity of a given melt. The set of exponents r, s, t, u, v and w vary depending on the conditions. The corresponding values are summarized in Table I. Our objective is to map the variations of all the properties which are necessary for determining the FOMs over the whole range of temperatures and compositions. We have therefore tried to fit all the simulation data with simple analytic functions, by finding the parameters which allow to minimize

χ 2 (a) = 1 N sim Nsim i=1 a i (x, T ) -a MD i a MD i 2 (4)
where a MD i is the value extracted from the MD simulation for the property a at a given thermodynamic point (N sim is the total number of simulated points) and a i (x, T ) is the corresponding value for the analytical function (x is the mole fraction of ThF 4 in the melt and T the temperature), which contains the parameters which need to be fitted.

First we start with the density. For this quantity, a small temperature dependence has to be included since this quantity usually varies linearly with temperature at fixed composition [START_REF] Benes | Molten salt reactor fuel and coolant[END_REF][START_REF] Dewan | Molecular dynamics simulation of the thermodynamic and transport properties of the molten salt fast reactor fuel LiF-ThF4[END_REF] . An excellent fit of the simulation data was obtained by using the following polynomial expression:

ρ(x, T ) = 2.2247 -0.0004329 × T +(15.412 -0.0015748 × T ) × x +(-8.1106 -0.0088374 × T ) × x 2 +(-11.380 + 0.020484 × T ) × x 3 +(9.2358 -0.010673 × T ) × x 4 (5)
The average error is χ(ρ) = 1.5 %. The ρ(x, T ) map is plotted on Figure 1; unsuprisingly the largest variations are due to the composition changes, and larger densities are obtained for high ThF 4 concentrations due to the large difference of molar weight between the two cations. The thermal expansion, which corresponds to the material's volume change as a function of temperature, was derived from the simulation results using the following equation:

β = - 1 ρ ∂ρ ∂T P (6)
Its variation with temperature at fixed composition is negligible, and the following polynomial relation was fitted for the composition changes:

β(x, T ) = β(x) = 0.0001×(3.0829-1.3187×x-0.43862×x 2 ) (7) yielding a value of 0.8 % for χ(β).

Similarly, the heat capacity of molten salts is usually constant with respect to temperature variations at fixed composition. This is not the case for the corresponding solid phases 14 . Our simulations, in which the molar heat capacity C p is calculated by differentiation of the molar enthalpy H (extracted from the N P T ensemble simulations) with respect to the temperature,

C p = ∂H ∂T P , (8)
confirm this behavior. Unlike other molten fluoride mixtures [START_REF] Beilmann | Excess heat capacity in liquid binary alkali-fluoride mixtures[END_REF] , for which a deviation from ideality was measured, in the LiF-ThF 4 melts C p varies almost linearly with composition. However, a fourth degree polynomial analytic function was fitted, The viscosity η was determined by the Green-Kubo method, i.e. by integrating the correlation function of the stress tensor, as follows : (11) in which σ αβ is one of the components of the stress tensor [START_REF] Hansen | Theory of simple liquids[END_REF] . In order to improve the statistics, the stress autocorrelation function is averaged over the five independent components (σ xy , σ xz , σ yz , σ xx-yy and σ 2zz-xx-yy ). Much longer simulation times are needed to compute the viscosity, because enough statistics needs to be accumulated in order to yield a plateau for the integral defined in Equation 11.

η = V k B T ∞ 0 σ αβ (0)σ αβ (t) dt
In molten salts, the viscosity typically follows an Arrhenius law for the variation with temperature. We have therefore used a more complex analytic form than for the other properties,

η(x, T ) = 0.21412 × exp 2383.3 T +(0.022784 × exp 7000.1 T ) × x +(94.189 × exp -4807.0 T ) × x 2 +(-3.0284 × exp 2281.7 T ) × x 3 +(0.76363 × exp 4016.4 T ) × x 4 (12)
Although a larger overall error is observed, χ(η)=9 %, it remains within the error bars of the MD simulation values. The variation of η with respect to x and T is shown on Figure 2. Unlike other properties, the composition variation appears to be small. This behavior is different from other melts such as LiF-BeF 2 mixtures, for which viscosities larger by several orders of magnitudes are obtained on the BeF 2 -rich side [START_REF] Salanne | Simulation of the liquid-vapor interface of molten LiBeF3[END_REF] . This is due to the strong network-forming ability of BeF 2 , which is not observed for ThF 4 . The largest viscosities are thus obtained for the lower temperatures, around the eutectic composition. Finally, the thermal conductivity is the most complex quantity to calculate. Like the viscosity, it is obtained using the Green-Kubo method [START_REF] Ohtori | Calculations of the thermal conductivities of ionic materials by simulation with polarizable interaction potentials[END_REF] , but in that case several phenomelogical coefficients need to be calculated from the integration of the corresponding microscopic fluxes auto-correlation functions. Nevertheless, it was shown by Ohtori et al. that its temperature variation is very small, and is a consequence of the change of the density of the melt only [START_REF] Ohtori | Thermal conductivity of molten alkali halides: temperature and density dependence[END_REF] . We used the same analytic function as for the density, and the following expression was obtained. with χ(λ) = 5.9 %, which is again well within the simulation error bars. The λ(x, T ) map is plotted on Figure 3, interestingly we observe a minimum for a composition of x = 0.6. Unlike the viscosity, the thermal conductivity of liquids does not show substantial variations when the structure of the melt changes (upon formation of a network, for example). Some efforts are currently made in order to improve our understanding of this quantity from the physical chemistry point of view in molten salts [START_REF] Ohtori | Thermal conductivity of molten alkali halides: temperature and density dependence[END_REF] as well as in molecular solvents such as water [START_REF] Roemer | Heat conduction and thermomolecular orientation in diatomic fluids: a non-equilibrium molecular dynamics study[END_REF][START_REF] Roemer | Nonequilibrium molecular dynamics simulations of the thermal conductivity of water: A systematic investigation of the spc/e and tip4p/2005 models[END_REF] . It is then possible to calculate the FOMs as defined by Equation 3 and the sets of exponents given in Table I at any thermodynamic point. The example of FOM1 in the eutectic region is shown on Figure 4; similar plots are provided for FOM2, 3 and 4 in the supplementary information. The smaller the FOM is, the better heat-transfer performance. We note that the FOMs of LiF-ThF 4 are relatively similar to those of LiF-NaF-KF and NaF-ZrF 4 eutectic mixtures [START_REF] Salanne | Heattransport properties of molten fluorides: determination from first-principles[END_REF] , which are among the candidates for being used in another molten salt-based reactor, the advanced high temperature reactor [START_REF] Forsberg | The advanced high-temperature reactor: high-temperature fuel, liquid salt coolant, liquid-metalreactor plant[END_REF] . The figure reveals that a LiF-rich melt is favorable for the thermal properties. Since the thorium is the fertile material, from which the fissile 233 U is formed during the reactor operation, a compromise would have to be found for choosing an optimal composition for the reactor performance. Furthermore, one of the most important technological difficulties which needs to be overcome is to have container materials with a very high resistance to corrosion. A consequence of this is that the temperature has to be kept as low as possible, leading to composition choices close to the eutectic. The spatial variation of the temperature may be relatively large inside the reactor core. Here we observe that the FOMs tend to decrease with increasing temperature, which should ensure a good extraction of the heat at all the conditions.

IV. CONCLUSION

The development of reliable interaction potentials, with ab initio accuracy, and the availability of large com-puter resources now allow us to build reliable thermophysical databases for molten salts. Such an approach is complementary to experimental measurements, yielding a large number of data at a very reasonable cost. In the present work we have illustrated this by calculating the density, the thermal expansion, the heat capacity, the viscosity and the thermal conductivity of LiF-ThF 4 mixtures for a wide range of compositions and temperatures. By fitting the calculated data with analytic functions, we have been able to map the variations of these properties for any thermodynamic point, and to determine the figures of merit which are necessary for designing molten salts reactor.

Although LiF and ThF 4 are the main components of the MSFR, it is sure that small amounts of other fluorides (UF 3 , UF 4 and PuF 3 , mainly) will also be present. In addition a large number of species will be formed due to the fission and the neutron capture reactions, and it remains impossible to extend the present work for all the elements. Therefore only the most important ones can be considered. However this difficulty can be overcomed by using thermodynamic models, which have up to now been fitted on experimental data only. The improvement of the thermodynamic models by the inclusion of the large amount of data provided by MD simulations should allow for the extension of the range of their applicability in the future.

Figure 1 :

 1 Figure 1: Variation of the density of the LiF-ThF4 melts with composition and temperature. The white region corresponds to temperatures below the liquidus.

CFigure 2 :

 2 Figure 2: Variation of the viscosity of the LiF-ThF4 melts with composition and temperature. The white region corresponds to temperatures below the liquidus.

Figure 3 :

 3 Figure 3: Variation of the thermal conductivity of the LiF-ThF4 melts with composition and temperature. The white region corresponds to temperatures below the liquidus.

  λ(x, T ) = 3.6014 -0.0014497 × T +(-24.140 + 0.017325 × T ) × x +(87.337 -0.072664 × T ) × x 2 +(-126.47 + 0.10924 × T ) × x 3 +(61.298 -0.053182 × T ) × x 4 ,(13)

Figure 4 :

 4 Figure of merit 1!

Table I :

 I Coefficients for the figures of merit used to evaluate the heat-transfer properties of molten fluorides25 . FOM4 is used for the heat-exchanger area.

	FOM Convection Régime	r	s t u v w
	1	Forced Turbulent 1.00 0.2 0.0 2.0 2.8 0.0
	2	Natural Turbulent 0.36 0.2 1.0 2.0 1.8 0.0
	3	Natural	Laminar 0.50 1.0 1.0 2.0 1.0 0.0
	4	-	-	1.00 0.2 1.0 0.3 0.6 0.6

Acknowledgements

This work was supported by the EVOL project in the 7 th Framework Programme of the European Commission (Grant agreement No. 249696) and by funding from the Natural Sciences and Engineering Research Council of Canada (NSERC).Computational resources for this work were provided by the Réseau Québécois de Calcul de Haute Performance (RQCHP).