J Wouters
email: jeroen.wouters@uni-hamburg.de

F Bouchet
email: freddy.bouchet@ens-lyon.fr

Rare event computation in deterministic chaotic systems using genealogical particle analysis

In this paper we address the use of rare event computation techniques to estimate small over-threshold probabilities of observables in deterministic dynamical systems. We demonstrate that genealogical particle analysis algorithms can be successfully applied to a toy model of atmospheric dynamics, the Lorenz '96 model. We furthermore use the Ornstein-Uhlenbeck system to illustrate a number of implementation issues. We also show how a time-dependent objective function based on the fluctuation path to a high threshold can greatly improve the performance of the estimator compared to a fixed-in-time objective function.

Introduction

Rare events may have a large impact on the dynamics of geophysical turbulent flows and the climate. In bistability situations, a rare transition can drastically change the structure of the flow, like for instance the bistability of the Kuroshio current [1,2] or a change of polarity of the Earth's magnetic field due to the turbulent dynamics of the Earth metal core [3]. Rare events can also be extremely important because of their impact on society, ecosystems or the economy. There are many such examples in climate dynamics, for example extreme droughts, heat waves, rainfalls and storms [START_REF]Intergovernmental Panel on Climate Change 2012 Managing the risks of extreme events and disasters to advance climate change adaption: special report of the Intergovernmental Panel on Climate Change[END_REF]. The probability and the impact of these events is likely to change in the future due to a changing climate [START_REF]Intergovernmental Panel on Climate Change 2012 Managing the risks of extreme events and disasters to advance climate change adaption: special report of the Intergovernmental Panel on Climate Change[END_REF]. The magnitudes of possible changes are however still uncertain [START_REF] Barriopedro | [END_REF].

On the one hand, for climate dynamics, there is a lack of sufficient reliable empirical data [6]. How could one assess faithfully the probability of events with recurrence times longer than one decade with only one or two century long reliable data? In the last decade, many methods have been developed to extract the most information possible from this too short time-series. For instance extreme value statistics [START_REF] De Haan | Extreme value theory: an introduction Springer series in operations research[END_REF][START_REF] Leadbetter | Extremes and related properties of random sequences and processes Springer series in statistics[END_REF] has interestingly allowed to extrapolate from the information available from empirical observation [START_REF] Kharin | [END_REF]10].

Another approach would be to critically study and understand the dynamics of rare events produced by complex climate models. This second approach seems to be the only available one for events with a recurrence time longer than decades or centuries. However, the current scientific state of the art does not yet allow to obtain many results by following this route. The first critical issue is a sampling problem. Indeed, if one wants to study events with century or millennial recurrence time and assess the reliability of the model dynamics to produce those events, a direct numerical simulation would require to have model runs of at least hundreds of thousands of years long in order to get reliable statistics on both the probability of and the dynamics leading to these events. As it is not always reasonable to trade computational length with model complexity, especially for the turbulent part which is responsible for most of the fluctuations, it is clear that we are facing an extremely difficult scientific challenge.

Is there a way to produce reliable statistics of specific rare events of a given model, without having to rely on prohibitively long direct numerical simulations? The same issue has been faced in many other scientific fields and has led to the development of some interesting approaches. Indeed, many of the complex systems studied in different branches of science feature events that are very rare but nevertheless very relevant due to their high impact. Take for example buffer overflows in digital communication networks, the insolvency of an insurer or bank, collisions in planetary systems, the dynamics of phase transitions in condensed matter, the long time dynamics of complex molecules in chemistry or biology, to name but a few. In recent years a number of promising algorithms have been developed to tackle these problems [START_REF]Rare event simulation using Monte Carlo methods[END_REF][START_REF] Bucklew | An introduction to rare event simulation[END_REF][START_REF] Moral | Mean field simulation for Monte Carlo integration[END_REF][START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF]. These rare event simulation algorithms can drastically reduce the error made on the estimation of small probabilities.

Generally speaking, the objective of the algorithms is to make rare events less rare, either by altering the dynamics (importance sampling) or by targeted killing and cloning of an ensemble of realizations (genealogical particle analysis or interacting particle algorithms). Upon estimation the intervention of the algorithm is then taken into account to obtain an estimate for quantities of the original system.

Within the class of genealogical particle analysis algorithms, a number of different strategies exist. A first crucial difference is the type of quantity one aims at estimating. One can be interested in the distribution of first entrance time to a set or to sample transition paths [START_REF] Cérou | [END_REF][START_REF] Rolland | Computing transition rates for the 1-D stochastic Ginzburg-Landau-Allen-Cahn equation for finite-amplitude noise with a rare event algorithm[END_REF][START_REF] Rolland | [END_REF]18], the probability of a rare event [19,20,21] or expectation values of long time averages such as the scaled cumulant generation function [22,23,24,25,26,27] as it appears in large deviation theory. For these different aims again different algorithms exist, for example geneological algorithms with fixed [19,20] versus variable [21] particle numbers, minimum action algorithms [18], or milestoning [28]. These algorithms have already been applied to a wide range of systems, for example percolation problems [29], in complex chemistry [30], polymer and biomolecule dynamics [31,32,33,34], magnetism [35,36], Burger turbulence [37,38].

The aim of this work is to make a first step in the application of those approaches to climate dynamics problems. Climate dynamics has specificities that make past approaches not directly adaptable. First, the climate is clearly out of equilibrium (without time reversal symmetry or detailed balance), therefore only non-equilibrium approaches can be considered. Second, the phenomenology of geophysical turbulent flows is dominated by large scale synoptic scales and is rather different from other complex dynamics, for instance molecular dynamics. And third, most climate models are deterministic models, or sometimes include a stochasticity that does not affect directly the synoptic scales.

The aim of this paper is to consider the latter specificity of many climate model. We address the following question: can rare event algorithms based on genealogical particle analysis be used effectively and efficiently for deterministic dynamics? Most algorithms rely on a Markov assumption, which is verified for deterministic models. However at the cloning stage, a new trajectory is branched from another one in order to produce a new ensemble member. For a strictly deterministic system, the offspring trajectory will not be different from its parent. To ensure separation of the two trajectories, one has therefore to add either a very small noise on the overall dynamics or a small change on the initial condition to the offspring trajectory and rely on the dynamics chaoticity.

A key issue is then to verify a posteriori that the noise is small enough so as not to distort the measured statistics and probabilities. A test using different decreasing noise strengths and checking for stability should therefore be used.

In order to perform the first study of the effectiveness of these approaches for chaotic deterministic dynamics, we have chosen to study a simple chaotic system with many degrees of freedom, and of relevance for climate dynamics: the Lorenz '96 model [START_REF] Lorenz | Predictability: A problem partly solved[END_REF]40]. We have also chosen the conceptually simplest and most robust genealogical algorithm that allows to sample invariant measure or transition probabilities: the genealogical particle analysis algorithm. We give a detailed heuristic presentation of the algorithm and a benchmark on the Ornstein-Uhlenbeck process in section 4.

The genealogical particle analysis algorithms explore the statistics of solutions of the dynamical systems by running an ensemble of realizations, interrupting the ensemble simulation at given times and killing ensemble members that do not perform well as measured by a weight or objective function and cloning the ones with a high weight. This selective procedure explains the terminology genealogical particle analysis. The individual realizations are also sometimes referred to as particles. The design of a good objective function is then arguably the main design issue one faces when using genealogical particle analysis algorithms. Other choices that have to be made are the number and timing of interactions and the number of particles to use. We will address these practical issues in a detailed study of the genealogical particle analysis algorithms on the Ornstein-Uhlenbeck process. This process is easy to simulate numerically and allows for analytic expressions to be derived; it is therefore well suited for the purpose of illustration and testing.

Another aim of this paper is to propose a systematic approach and procedure to get reliable results and error estimates. We propose to build the tail of the cumulant distribution funciton of interest by gluing together pieces of results obtained for different cloning parameter by a systematic study of the most reliable one, through an empirical estimate of the algorithm variance. Moreover, we propose a procedure to test empirically this class of algorithms against the real dynamics. Indeed, for a model like the Lorenz '96 model, we have no theoretical results that can serve as a benchmark.

The paper is organized as follows. In Section 2 we discuss how the need for rare event simulation techniques arises, what the objective of such algorithms is (making rare events typical) and how this goal can be achieved for stochastic processes by implementing a genealogical particle analysis simulation. In Section 3 we present a brief discussion of the theory of large deviations and what it can say about the way in which rare events are reached by a process. This theory can be used to implement a more efficient rare event sampling method. In Section 4 we proceed by implementing the genealogical particle analysis simulation to the Ornstein-Uhlenbeck system. We discuss in depth the selection of the parameters in the algorithm. In Section 5 we then present the implementation of the genealogical particle analysis simulation on a chaotic deterministic dynamical system. Finally, we present our conclusions in Section 6.

Rare event computation for Markov dynamics

In sections 2.1 and 2.2 we present a classical discussion of the inefficiency of brute force Monte Carlo simulation for estimating small probabilities. This motivates the need for rare event computation techniques. We introduce the genealogical particle analysis algorithm and the related theory in Section 2.3.

Motivation

The goal of rare event computation techniques is to make the numerical estimation of small probabilities more efficient. The necessity of using such techniques is demonstrated by the sampling of the tail of a distribution P using independent samples identically distributed according to the distribution P . Say one wants to estimate a small probability γ A = P (X ∈ A)

1 by means of a brute force Monte Carlo estimate

γA = 1 N N i=1 1 A (X i ) (1)
where 1 A is the indicator function on the set A. The estimator γA is an unbiased estimator of γ A since the expectation value of γA is clearly γ A . When the number of samples N is large enough for γA to follow a central limit theorem, the statistical error of the estimator can be quantified by its variance V ar(γ A ) = V ar(1 A (X))/N . Furthermore

V ar(1 A (X)) = E((1 A (X) -γ A ) 2 ) = E(1 A (X)) -γ 2 A = γ A -γ 2 A ≈ γ A (2)
when γ A is small. The relative error of the estimator RE being proportional to the standard deviation divided by the estimated quantity, we have RE ∼ 1 √ γ A N . The relative error quickly becomes large as γ A goes to zero for fixed sample size N . Fortunately there exist methods for estimating small probabilities more efficiently.

Importance sampling

The main ingredient of rare event computation techniques is a sampling from a modified distribution together with an adapted estimator to counteract this change of measure. This method to lower the estimator variance of a rare event probability is termed importance sampling. Again the example of the sampling from independent identically distributed random variables provides valuable insights.

Say we want to estimate

γ A = ˆdXρ(X)1 A (X) 1
where ρ is the density for our random variable X. Instead of doing a straightforward sampling of X as in (1), assume we can sample from a modified measure ρ for which ρ(X) > 0 whenever X ∈ A and ρ(X) > 0. In such a case, the probability we want to estimate can be rewritten as

γ A = ˆdX ρ(X) ρ(X) ρ(X) 1 A (X) = Ẽ(L(X)1 A (X)) L(X) := ρ(X) ρ(X) whenever 1 A (X)ρ(X) > 0
and we can therefore estimate γ A using the estimator

γA = 1 N N i=1 L( Xi )1 A ( Xi ) (3)
on samples Xi distributed according to ρ.

The variance for such an estimator is

Ṽ ar(L(X)1 A (X)) = Ẽ(L 2 (X)1 A (X)) -γ 2 A = E(L(X)1 A (X)) -γ 2 A
If we could take ρ as the conditional measure with ρ(X) = ρ(X)/γ A for X ∈ A and zero elsewhere, such that L(X) = 1 A (X)γ A , this would result in a zero variance estimator. This estimator is however not practically implementable, since for this we would need to know the value of γ A , which is the value we seek to calculate. This calculation demonstrates some important points however. First of all, it shows that a change of measure can indeed reduce the variance of the estimator. Although the ideal change of measure is not feasible in practice, a change of measure that is in a sense close to it should also give a substantial variance reduction. This modified measure should therefore have most of its weight on the set of interest A. On the other hand, this also implies that one needs to have some understanding of the shape of the set A and the distribution on it to construct an efficient importance sampling.

Skewing a normal distribution

To illustrate how a change of measure can provide significant variance reductions, even if the modified measure is not the ideal conditional measure, we discuss an example for normally distributed random variables. This example will also be useful to illustrate and validate our rare event algorithm for dynamical systems.

Say we want to estimate the probability of the rare event A = {x > a} for a normally distributed random variable x ∼ N 0,1 with zero average and standard deviation equal to one. Assume that we can skew the distribution with an exponential function ρ

(X) = ρ(X) exp(CX)/E(exp(CX)) = 1 √ 2π exp -(X-C) 2 2
which constitutes of a shift of the average by C. Since L(X) = E e CX / exp(CX) = exp (-CX + C 2 /2), the variance of the terms in the importance sampling estimator is now

Ṽ ar(L(X)1 A (X)) = P -c,1 (x > a)e C 2 -γ 2 A ( 4
)
where P µ,σ denotes probabilities under a normal distribution with mean µ and variance σ. The standard deviation, the square root of the variance (4), is plotted for a = 2 in Figure (1). The standard deviation has a single minimum, which is obtained for a value of C which is close to C = 2, for which the mean of the tilted value coincides with the threshold. This basic example illustrates how importance sampling can lower dramatically the estimator variance.

As Fig. (1) shows, the relative error, which is proportional to the plotted quantity Std(γ A )/γ A , can be reduced by a factor of more than 4 for the case where A = [2, +∞). Since the error decreases as 1/ √ N , this means a at least 16-fold longer brute force simulation would be necessary to obtain a similarly accurate result. For graphical purposes a relatively low threshold 2 was chosen here. For higher thresholds, the performance gains increase drastically, with a reduction of computational effort by a factor 6 × 10 5 when A = [5, +∞).

Genealogical particle analysis algorithm

The motivation for rare event simulation and the discussion of importance sampling have shown that it is necessary to make rare events less rare. This concept can be applied to stochastic processes such as the paths followed by either stochastic dynamics or chaotic deterministic dynamics. In those cases the objective is to alter the probability of certain paths that are connected to the rare event one wants to study.

Two different strategies are employed to alter the path sampling in stochastic dynamical systems. The first one is to alter the dynamical equations of the system by introducing a forcing term [41]. By tuning a parameter of the added forcing term, Figure 1: The ratio of the standard deviation to estimated probability of an exponentially tilted gaussian importance sampling estimator for a threshold a = 2 with N = 1 one can then attempt to decrease the variance of the rare event estimator. A second strategy consists in calculating an ensemble of realizations of the stochastic system in parallel and manipulating the ensemble members by performing selections at a finite number of selection times so as to bias the population.

Here we will use the second strategy, by employing a variant of the so called genealogical particle analysis algorithms. The selections applied to the ensemble consist of dynamical trajectories, called particles, being copied and killed depending on weight factor asigned to every ensemble member. This strategy has the advantage of not altering the dynamical trajectories themselves, such that their dynamics can be studied a posteriori.

Extensive analysis of the convergence of genealogical particle analysis algorithms can be found for example in [START_REF] Moral | Mean field simulation for Monte Carlo integration[END_REF][START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF]. In the following sections, we perform a simpler calculation, assuming a mean field approximation, to demonstrate the evolution of the expected particle distributions in a genealogical particle analysis. The validity of this mean field approximation for large particle number is the subject of the complete proofs given in [START_REF] Moral | Mean field simulation for Monte Carlo integration[END_REF][START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF]. Before going to a truly interacting genealogical particle analysis in Section 2.3.3, we first get some insight by looking at an algorithm where particles are reweighted, but by a factor depending only on the evolution of the particle itself, in Section 2.3.1.

A non-interacting genealogical particle analysis

We calculate rare events of a continuous time Markov chain. P (2) (y|x, ∆t) denotes the transition probabilities from configuration x to y over a time interval ∆t. We are interested in the probability of Algorithm 1 Non-interacting weighted particle system (i) Initiate M particles in configuration x 0 : ξ i,0 = x 0 for 1 i N t 0 = M (ii) For every time step k ∈ {1, . . . , n} (a) Propagate ξ i,t k-1 under the dynamics, resulting in ζ i,t k distributed according to

P (2) (ζ i,t k |ξ i,t k-1 , ∆t k ) with ∆t k = t k -t k-1 (b) Calculate weights W i,k for particle i: W i,k = W (ζ i,t k , ξ i,t k-1 ) := exp(V (ζ i,t k ) -V (ξ i,t k-1 )) for a suitably chosen weight function V (c) Generate a new particle distribution ξ j,t k consisting of N i,k copies of particle with configuration ζ i,t k where N i,k is chosen at random such that E(N i,k ) = W i,k (note that N t k = i N i,k ) (iii) Finally, for any F , calculate F = 1 M Σ Nt n i=1 F (ξ i,tn ) to estimate E 0,tn (F e V )e -V (x 0 ) (to estimate γ A take F (x) = F A (x) := 1 A (x) exp(V (x 0 ) -V (x)))
being in a set of configurations A at a time t = t n , given that the process started in configuration x 0 at time t = 0. N t denotes the number of particles at time t, whereas {ξ i,t } 1 i Nt denote the particle configuration at time t. E 0,t denotes expectation values under the original Markov dynamics P (2) at time t. The algorithm to generate the particles is described in the box Algorithm 1.

The algorithm can be summarized as follows: after initialization (step i) the ensemble members are evolved forward in time (step (ii)a) and weight values are calculated from the previous and current configurations of the particle, ζ i,t k and ξ i,t k-1 respectively (step (ii)b). Based on these weight values, ensemble members are killed or cloned (step (ii)c). Repeating this procedure results in a reweighted sample of paths, from which expectation values of the unweighted path distribution can be estimated (step iii). The rare event probability for a set A can be obtained by taking as observable

F (x) = 1 A e -V (x) e V (x 0 ) such that E 0,tn (F e V )e -V (x 0 ) = E 0,tn (1 A ).
Note that the random number N i,k generated in step (ii)c can be zero, such that particles can be killed as well as cloned (when N i,k > 1). A way to generate the random number described in step (ii)c is to take N i,k = W i,k + u where u is uniformly distributed on [0, 1] and x is the floor of x (the largest integer smaller than x).

Unbiased estimator

We first show that Algorithm 1 provides an unbiased estimator for the quantity E 0,tn (F e V )e -V (x 0 ) , i.e. the algorithm results in a random estimate whose expectation value equals the quantity to be estimated:

E 1 1 M Σ Nt n i=1 F (ξ i,tn ) = E 0,tn (F e V )e -V (x 0 )
where E 1 is the expectation over the random variables in the algorithm.

Write N (x, t) the particle number at configuration x, i.e. N (x, t)dx is the number of particles with x x i,t < x + dx:

N (x, t k-1 ) = Σ Nt k-1 i=1 δ(x -ξ i,t k-1 )
According the algorithm 1, if a particle sits at ξ i,t k-1 at time step k -1, N i,k copies are created of ζ i,t k at the next time step. Hence, the particle number at the next time step will be

N (x, t k ) = Σ Nt k-1 i=1 N i,k δ(x -ζ i,t k ) (5)
One step in the algorithm involves the generation of two sets of random variables, the updated particle configurations ζ i,t k , which is conditioned on ξ i,t k-1 , and the number of particle copies N i,t k , which depends on both ζ i,t k and ξ i,t k-1 . The expectation value of functions depending on the particle configurations ξ i,t k at step k can therefore be expressed as the expectation value

E ξ i,t k (•) = E ξ i,t k-1 (E ζ i,t k |ξ i,t k-1 (E N i,k |ζ i,t k ,ξ i,t k-1 (•)))
Applying this expression to Eq. 5 and using the probabilities for the updated particle configurations

P (ζ i,t k |ξ i,t k-1 ) = P (2) (ζ i,t k |ξ i,t k-1 , ∆t k )
and that the number of particle copies

E N i,k |ζ i,t k ,ξ i,t k-1 (N i,k ) = W i,k (ζ i,t k , ξ i,t k-1 ), we have E 1 (N (x, t k )) = E ξ i,t k-1 (E ζ i,t k |ξ i,t k (Σ Nt k-1 i=1 W i,k (ζ i,t k , ξ i,t k-1 )δ(x -ζ i,t k ))) = E ξ i,t k-1 Σ Nt k-1 i=1 ˆdyP (2) (y|ξ i,t k-1 , ∆t k )W i,k (y, ξ i,t k-1 )δ(x -y) = E ξ i,t k-1 Σ Nt k-1 i=1 P (2) (x|ξ i,t k-1 , ∆t k )W i,k (x, ξ i,t k-1 ) = E ξ i,t k-1 ˆdzP (2) (x|z, ∆t k )W (x, z)Σ Nt k-1 i=1 δ(z -ξ i,t k-1 ) = ˆdzP (2) (x|z, ∆t k )W (x, z)E 1 (N (z, t k-1 )) = ˆdzP (2) (x|z, ∆t k )e V (x)-V (z) E 1 (N (z, t k-1 ))
This equation relates the expected particle density at step k to the density at step k -1. By iteration we can relate the density at step k to the density at the start of the algorithm, which is M δ(x -x 0 ):

E 1 (N (x, t k )) = ˆdx k-1 . . . dx 1 dzP (2) (x|x k-1 , ∆t k ) . . . P (2) (x 1 |z, ∆t k ) × e V (x)-V (z) M δ(z -x 0 ) = M P (2) (x|x 0 , t k -t 0 )e V (x)-V (x 0 )
The expectation value of the quantity calculated at the end of the algorithm in step iii is therefore

E 1 1 M Σ Nt n i=1 F (ξ i,tn ) = 1 M E 1 ˆdxΣ Nt n i=1 δ(x -ξ i,tn )F (x)
Algorithm 2 Interacting genealogical particle analysis (i) Initiate M particles in configuration x 0 , ξ i 0 = x 0 for 1 i N 0 = M (ii) For every time step k ∈ {1, . . . , n} (a) Propagate ξ i,t k-1 under the dynamics, resulting in ζ i,t k distributed according to

P (2) (ζ i,t k |ξ i,t k-1 , ∆t k ) with ∆t k = t k -t k-1 (b) Calculate weights for particle i: Wi,k = W i,k (ζ i,t k , ξ i,t k-1 ) Z k Z k = 1 N t k Σ i W i,k (ζ i,t k , ξ i,t k-1 ) (c) Store the value of the normalizing factor Z k (d) Generate a new particle distribution ξ j,t k consisting of N i,k copies of particle with configuration ζ i,t k where N i,k is chosen at random such that E(N i,k ) = Wi,k (iii) Finally calculate 1 M Σ Nt n i=1 F (ξ i,tn ) n k=1 Z k to estimate E 0 (F e V )e -V (x 0 ) (for γ A take F (x) = 1 A (x) exp(V (x 0 ) -V (x))) = 1 M ˆdxE 1 (N (x, t n ))F (x) = ˆdxF (x)e V (x)-V (x 0 ) P (2) (x|x 0 , t n -t 0 ) = E 0,tn (F e V )e -V (x 0 )
Note that the expected total particle number

E 1 (N (t k )) = ´dxE 1 (N (x, t k )) = M ´dxP (2) (x|x 0 , t k -t 0 )e V (x)-V (x 0
) is in general not preserved over time. The particle number can strongly increase, which entails a large numerical cost. The solution to this problem is to renormalize the weights calculated in step (ii)b of the algorithm, hence introducing an interaction between the ensemble members. We discuss this new algorithm in the next section. As we will see, the interaction complicates the algorithm analysis.

Interacting particles

We now add interaction to the weights of the particle system, so as to control the particle number. A similar analysis as in the previous section can still be carried out, if one assumes that the number of particles used in algorithm is large enough, such that averages over particle configurations can be replaced by an expectation value under the law of large numbers (mean field approximation). The corresponding algorithm is described in the box Algorithm 2. By applying the algorithm to a function F A (x) = 1 A (x) exp(V (x 0 ) -V (x)) with 1 A the indicator function of the set A, estimates γA of the probability γ A can be obtained.

We again perform an analysis of the evolution of the expected particle distribution in Algorithm 2. For simplicity of the derivation we assume that the number of particles

Z k = 1 N t k Σ Nt k i W (ζ i,t k , ξ i,t k-1 ) ≈ E 1 (W (ζ i,t k , ξ i,t k-1 )) (6) = ˆdx dy E 1 (N (x, t k-1 )) E 1 (N t k-1 ) P (2) (y|x, ∆t k )W (y, x) (7)
Using this estimate and the same reasoning as for the non-interacting particle algorithm, we have for the expected particle distribution that

E 1 (N (x, t k )) ≈ ˆdzP (2) (x|z, ∆t k ) W (x, z) E 1 (W ) E 1 (N (z, t k-1 )) (8)
The expectation value of W in the denominator can be substituted using Eq. 6. The particle number is now constant:

E 1 (N t k )) = ˆdxE 1 (N (x, t k )) = E 1 (N t k-1 ) = . . . = E 1 (N t 0 ) = M
We therefore have that

E 1 (W (ζ i,t k , ξ i,t k-1 )) ≈ ˆdxdyP (2) (x|y, ∆t k )W (x, y) E(N (x, t k-1 )) M Inserting this into 8, we have E(N (x, t k )) ≈ M ´dyP (2) (x|y, ∆t k )W (x, y)E(N (y, t k-1 )) ´dxdyP (2) (x|y, ∆t k )W (x, y)E(N (y, t k-1 ))
Therefore by iteration

1 M Σ N i=1 F (x i,tn ) ≈ E 0 (F (x n )W (x n , x n-1 ) . . . W (x 1 , x 0 )) E 0 (W (x n , x n-1 ) . . . W (x 1 , x 0 )) = E 0 (F (x n )W (x n , x n-1 ) . . . W (x 1 , x 0 ))
Z n where Z n = E 0 (W (x n , x n-1 ) . . . W (x 1 , x 0 )). From Eqs. 6 and 8 we see that

Z k ≈ Z k Z k-1 and therefore 1 M Σ N i=1 F (x i,tn ) n k=1 Z k ≈ E 0 (F (x n )W (x n , x n-1 ) . . . W (x 1 , x 0 )) = E 0 (F e V )e -V (x 0 )
The above reasoning can also be extended to show that path dependent quantities (such as E[x(τ )|x(T ) > a] for τ < T ) can be estimated from the ancestral paths of the particle system.

Time-dependent weighting

The weighting function W (x, y) = exp(V (x) -V (y)) results in a particle distribution tilted by exp(V (x)) at all selection times t k . More flexibility can be obtained by using time-dependent weighting, for example with a weighting function of the form

W (t k , x, y) = exp(V t k (x) -V t k-1 (y)) (9)
This way the telescoping canceling of V t k is preserved in products of weights that appear in the calculation of the tilted measure. For example,

W (t k , x, y)W (t k-1 , y, z) = e Vt k (x)-Vt k-1 (y) e Vt k-1 (y)-Vt k-2 (z) = exp(V t k (x) -V t k-2 (z))
The result is again a particle distribution tilted by exp(V t k (x)) at time t k , as with the time-independent weight function. However, paths up to the final time will have different weights, which can make a large difference in the algorithm performance, as we will demonstrate in Section 4.2.4.

Fluctuation paths and the weighting function

The ideal change of measure discussed in Section 2.2 suggests to make the rare event that is the least rare the most probable one under the reweighted dynamics. This rationale extends not only to the distribution of the system at the final time, but also to the entire path up to the final time. This means that variance can be reduced if the least unlikely path leading to a high threshold is made more likely under the particle system dynamics.

For stochastic differential equations in the weak noise limit, the least unlikely path from an attractor can be calculated from Freidlin-Wentzell type large deviation theory and is called a fluctuation path (also sometimes an instanton). The particle system dynamics can be made to more closely follow the fluctuation path by using the time-dependent weighting discussed in Section 2.3.4. Even if the particle distribution at the final time is the same as the particle distribution obtained with a constant weighting, there is still a variance reduction since less particles are killed, increasing the independence of the particle and thus the effective particle number.

Fluctuation paths

The probability of a given path in a stochastic differential equation with small noise,

dX = b(X )dt + √ dW,
where W is a Brownian motion, can be estimated using the Freidlin-Wentzell large deviation theory. The theory determines the probability of seeing a path that is close to a specified continuous function in the limit of going to zero. It roughly states that lim

→0 log P [X ∈ F ] = -inf ω∈F I(ω)
where F is any closed subset of the set of continuous trajectories and the rate functional I is called the action. The action is given by

I(ω) = 1 2 ˆT 0 dt( ω(t) -b(ω(t))) 2 = ˆT 0 dtL[ω, ∂ω ∂t ] (10)
L[ω, ∂ω ∂t ] := 1 2 ( ω(t) -b(ω(t))) 2 (11)
The distribution of paths leading to rare fluctuations then concentrates around action minima as decreases, with given constraints. If the set of paths F contains the evolution along the deterministic dynamics ẋ = b(x), this path will obviously minimize the above action, hence the need for constraints to obtain more interesting results. For example, in the simple case were the deterministic dynamics ẋ = b(x) has a single attractor x 0 , the distribution of the paths conditioned on X(0) = X 0 , X 0 = x 0 concentrate close to the minima of the action ´0 -∞ dtL[ω, ∂ω ∂t ] with the boundary conditions X(-∞) = x 0 and X(0) = X 0 . Such a path is called a fluctuation path leading to X 0 (it is also sometimes called an instanton, but instanton usually rather refers to those fluctuation paths that connect attractors to saddle points).

Rare event simulation for a stochastic process: the Ornstein-Uhlenbeck process

We now illustrate some of the practical issues arising when implementing a genealogical particle analysis algorithm for rare event estimation. We start off with a stochastic process for which we can calculate explicitly all of the probabilities that we want to estimate, for pedagogical reasons, and so that we can compare the numerical results to the analytic expressions.

Description of the Ornstein-Uhlenbeck process

We consider the Ornstein-Uhlenbeck process

dx = -λxdt + σdW (12)
As the transition probabilities P (x(t)|x(0)) are Gaussian, the Ornstein-Uhlenbeck process preserves Gaussianity. Using the Itô formula, one can derive that the mean m(t) = E(x(t)) and the variance v(t) = E((x(t) -m(t)) 2 ) evolve according to the equations

ṁ = -λm v = σ 2 -2λv (13)
which can be easily solved explicitly. The probability that x exceeds a certain threshold a at a time t, given that the process started at x(0) = 0 at time zero can be calculated explicitly as

P (x(t) > a|x(0) = 0) = ˆ∞ a dxN m(t),v(t) (x), (14)
where m and v solve [START_REF] Moral | Mean field simulation for Monte Carlo integration[END_REF] and N is the probability density function of the normal distribution with mean m and variance v. We consider in the following the estimation of this probability through a genealogical particle analysis algorithm. Below we will use the parameter values λ = 1 and σ = 1.

Algorithm implementation

Let us assume we seek to estimate a small probability γ A , for instance γ A = P (x(t) ∈ A|x(0) = 0). We denote M the number of particles for each realization of the algorithm.

Then each independent realization i of the algorithm, with M particle each, will give an estimate γA,i . According to a theorem discussed in [19], asymptotically for large M , the random number γA,i is distributed according to a Gaussian distribution with standard deviation σ A (M ) = σ A / √ M and a corresponding relative error RE(M ) = σ A (M )/γ A . The value of the estimator relative error RE(M ) is essential as it quantifies the relative error one should expect for each realization of the algorithm, and thus the quality of the result. How the estimator relative error RE(M ) depends on the number of selections, on their timing, and the type observables are critical questions that we analyze in this section.

Number of particles

The result in [19] proves the existence of the central limit theorem, but does not give a value for the estimator variance RE(M ). In order to get an estimate of RE(M ), we compute it empirically by performing K independent algorithm realization and using the estimator

RE(M ) 1 K K i=1 (γ A,i -γ A ) 2 /γ A . (15)
In this formula the value of γ A will be either the theoretical value when it is available, for instance for the Ornstein-Uhlenbeck process, or the estimated value of the probability by averaging γA,i over K realizations. In the following, by an abuse of notation, RE denotes either the theoretical estimator variance of the estimator variance evaluated from [START_REF] Cérou | [END_REF], which one should be clear from the context.

We first study the estimator variance RE for the OrnsteinUhlenbeck case. We first test whether or not the regime of the central limit theorem has been reached by changing the number of particles M , and verifying whether RE (15) reduces by the corresponding √ M factor. Figure 3a shows the expected decrease in relative error as the number of particles is increased for a range of thresholds a (see eq. ( 14)). The inverse square root behavior of the error with increasing number of particles is demonstrated for a fixed threshold in Figure 3b. The parameters are specified in the figure captions.

We study how the estimator variance RE depends on the other numerical parameters in the following sections.

Number of selections and their timing

Since little theoretical analysis has been performed on the optimal number of selection steps, this is the most heuristic choice to be made. Some numerical analysis of this issue has been performed in [42]. For the problem they investigate, changing the number of selections, and using equidistant in time selections, the estimator variance clearly shows a minimum for a certain number of selections.

This result can be interpreted as follows. Selections shouldn't be performed too frequently, as cloning increases correlations between the particles and therefore reduces the effective number of independent particles, increasing the estimator variance. If not performing selections frequently enough however, the particle distribution relaxes to the unbiased particle measure, leading to the poor brute force Monte-Carlo variance. This can be seen in Figure 4: for low thresholds a, importance sampling is useless and estimations with a small number of selections have the lowest estimator relative error RE. Due to the large time between selections, the particles have relaxed to the particle measure of brute force Monte Carlo simulations and therefore have a similar estimator variance. For higher thresholds, for instance for a > 1.7, it becomes advantageous to kill a larger number of particles to obtain a more skewed final particle distribution, in order to lower the variance. For the threshold value a = 2 the optimal number of interactions among the values in the figure is N = 16. For higher thresholds there is a small reduction in error by increasing the number of selections, although increasing the number of selections much further beyond N = 64 results in an overall increase of error.

N i =2 N i =8 N i =16 N i =64 (
N i =16 N i =256 N i =1024 N i =4096 (b)
Figure 4 also illustrates the large estimator variance improvement for the genealogical particle analysis algorithm compared with Monte-Carlo sampling, as soon as a ≥ 2.

Besides the number of selections, there also seems to be little theoretical understanding of the optimal timing for selections. One strategy to selection timing is to calculate on-the-fly a criterion on the distribution of particle weights (such as the squared coefficient of variation or entropy) and only perform selection if a fixed threshold is exceeded. The convergence of such adaptive selection strategies is discussed in [43].

Estimating a range of over-threshold probabilities

In the following the weight function W (x, y) = exp(C(x -y)) is used. From the point of view of the estimator variance RE [START_REF] Cérou | [END_REF], to each value of the threshold a corresponds an optimal value of C, denoted C * (a), or equivalently for each value of C the estimator variance has a minimum for a given value of a, denoted a * (C). For instance Figure 4 shows that the value C = 4 is optimal for a 2.5 = a * . In simple cases, we expect C * (a) to increase monotonically with a.

There is an optimal value of a for each value of C, however the estimate is good for a range of thresholds around this optimum. When instead of a particular over-threshold probability one is interested in the tail of the complete distribution probability, one can perform a number of genealogical particle analysis simulations each with different value of C, and select for each threshold the value corresponding to the lowest estimator variance RE. Figure 5 illustrates how the tail of P (a) = P (x(t) ≥ a|x(0) = 0) can be estimated this way, for x the Ornstein-Uhlenbeck process. For large values of the threshold (above a ≈ 4.6) all estimates have a high error and the highest value of C is chosen by default. As can be seen on this figure there is very good agreement with the theoretical value up to probabilities as low as 10 -10 . Using this strategy, we can accurately estimate the tail of the over-threshold distribution down to probabilities as small as 10 -10 , with relative error lower than one.

Selections along the fluctuation paths

We have discussed in Section 2.3.4 that a time-dependent weighting function can be used. In this way the particle distribution can be weighted with different exponential factors C(t k ) all along the path, but still lead to the same exponentially tilted final particle distribution. Furthermore, in Section 3.1 we have discussed how for small noises, most of the paths leading to a rare event will concentrate around fluctuation paths that minimize the action functional. The aim of this section is to demonstrate the interest of using fluctuation paths to construct time-dependent weighting functions in order to increase the efficiency of the genealogical particle analysis algorithm.

Since the Ornstein-Uhlenbeck process is linear, taking limits of higher thresholds is equivalent to taking a weak-noise limit through a rescaling of the x coordinate. Hence, for fixed noise intensity σ, paths starting at x 0 = 0 conditioned on reaching the final threshold a will concentrate around the fluctuation paths in the limit a → ∞. The action (10) for the Ornstein-Uhlenbeck process ( 12) is given by I [X] = ´T 0 dτ ( Ẋ +λX) 2 . Taking as boundary conditions X(0) = 0 and X(T ) = a the fluctuation paths are easily computed to be X f (t) = a sinh λt sinh λT . By using the potential function W (t, x, y) = exp(C(t k )x -C(t k-1 )y) with a weight parameter C(t k ) dependent on the selection time t k , we can control μ(t k ), the mean particle position at t k , by fixing C(t k ). The expected particle distribution for the Ornstein-Uhlenbeck process tilted with this weighting function after the selection at

t k is exp (C(t k )x) N 0,v(t k ) (x)/ ˆdx exp (C(t k )x) N 0,v(t k ) (x)
as discussed in Section 2.3.4. The corresponding expected mean particle position is therefore μ(t k ) = C(t k )v(t k ) where v(t) = (1 -exp(-2t))/2 is the variance of the Ornstein-Uhlenbeck process at time t (the solution of Eq. 13 with v(0) = 0). Choosing C(t k ) = X f (t k )/v(t k ), μ(t) follows X f (t) and the algorithm particle distribution closely follows the fluctuation path leading to the threshold a.

Figure 6 shows the effect of using a weight function based on a fluctuation path versus an exponential weight function. The bottom two plots show how using the fluctuation path significantly decreases the fraction of particles that are killed during the selection steps (N (-) k ) to the number of particles at that time step (N k ). This is also illustrated in the plots in the top and middle rows. The top plots show the paths from the initial state for all surviving particles at the final time. Paths that have been killed during the process are not shown. We call these paths the ancestral paths. As can be seen on the top left plot, only few trajectories from the initial stage of the algorithm are ancestors of the final positions. This is not the case for the top right plot. The algorithm using a weight based on a fluctuation path has a much larger number of ancestors. This richer ancestral tree results in a decreased estimator error for the over-threshold probabilities, as is demonstrated in Figure 7.

Note that for both the exponential weighting function and the weighting based on the fluctuation path, the paths reaching the threshold follow the fluctuation path. Other paths reaching the threshold are so rare that few of them are generated, even if they are more likely to survive selection in the case of exponential weighting. Note that the killed paths, partially shown in the middle row of Fig. 6, tend to have a negative change in position before being killed. The higher target path in the exponential makes for a higher average dissipative force -x on the particles, leading to a large discrepancy to the total number of particles N k (bottom) for genealogical particle analysis algorithms with either exponential weighting with C = 6.0 (left) or weighting based on the fluctuation path (right) for the fluctuation path ending at a = 3.0 at the final time T = 2. The dashed black lines in the top plots show the fluctuation path. The dash-dotted line in the top left plot shows the mean of the target particle distributions after selection (equals Cv(t)). The average number of particles for both simulations is M = 10 4 and the number of selections steps is 32. For graphical purposes a randomly selected sample of 2% of the ancestral and killed paths are shown in the first two rows between the actual particle distribution and the target distribution at selection times.

Genealogical particle analysis algorithm for a deterministic dynamical system: the Lorenz '96 model

The Lorenz '96 model is a deterministic dynamical system that is often used as a toy model in the meteorology community. It was proposed by Lorenz as part of a study on Figure 7: The over-threshold probability P (x > a) as estimated by the genealogical particle analysis algorithm either with an exponential weight (blue long-dashed line) or a weight based on fluctuation paths (red medium-dashed line). The two shortdashed lines, at equal distance from the estimated averages, correspond to a 2 standard deviation interval of the estimator. The full line is the analytic result. Both implementations use the same number of particles N = 1e4 and 32 selections steps and both have roughly the same computational cost. error growth and predictability for chaotic dynamical systems [START_REF] Lorenz | Predictability: A problem partly solved[END_REF]40].

A crucial difference between the famous Lorenz '63 model and the less well-known Lorenz '96 model is that the latter has a large number of degrees of freedom. Indeed macroscopic variables of deterministic systems with a large number of degrees of freedom often behave qualitatively similar to solutions of stochastic differential equations with much less degrees of freedom. Such results can be proven for some specific types of models (feauturing separation of time scale, independence, ...).

It is believed however that similar results remain true for a wide range of models and observables even though mathematical proofs are out of reach. If this conjecture is correct, then the sampling of rare events through genealogical particle analysis algorithms should be applicable to macroscopic variables of deterministic systems with a large number of degrees of freedom. In this section, we demonstrate empirically, through numerical simulation, that genealogical particle analysis algorithms can indeed efficiently sample the tail of the energy distribution for the Lorenz '96 model.

Description of the model

The Lorenz '96 model consists of L variables x i on a ring i ∈ {0, .., L-1}, with dynamics ẋi = x i-1 (x i+1 -x i-2 ) + R -x i where indices i are in Z L , i.e. the index i is identified with i mod L if i / ∈ {0, .., L -1}. The non-linear part of the dynamics x i-1 (x i+1 -x i-2 ) conserves the energy

E(x) = 1 2L L i=1 x 2
i , while R is a forcing and -x i a linear dissipation. The dynamics is chaotic for R ≥ 8 [START_REF] Lorenz | Predictability: A problem partly solved[END_REF]40]. We will estimate the probability of reaching a certain energy threshold after a time t, starting from the zero vector x 0,i = 0 ∀i. A small perturbation N 0,1 is added to the particle configuration to make the trajectories diverge. For large enough times T , the system will therefore relax to its physical invariant measure and it makes sense to determine probabilities of exceedances of macroscopic observables γ Et := P (E(x(T )) > E t ). Throughout the article we will use a number of variables L = 32 and a forcing R = 2 8 = 256.

Figure 8 shows a plot of the over-threshold probabilities of the energy of the Lorenz '96 system, estimated through a brute force Monte-Carlo simulation from randomly perturbed initial conditions. Given that we have finite computer resources at our disposal, assume we can generate at most M = 10 5 independent measurements of the energy. If the maximal relative error that we are willing to tolerate is for example 0.5 then since RE = γ Et -γ 2

Et /( √ M γ Et ) ≈ 1/( M γ Et ) the lowest probability that we can estimate is approximately γ Et = 1/M (RE) 2 = 4.10 -5 . From Figure 8 we can deduce that the corresponding highest energy threshold obtainable lies around an energy threshold E t = 1785. Beyond this threshold the use of rare event algorithms becomes necessary.

Algorithm implementation

We use the following settings for the genealogical particle analysis simulation. The initial condition is set to x i = N 0,1 . The total integration time per realization is T = 1, 27. This time interval corresponds to roughly 5 times the decorrelation time of the energy observable. The standard deviation of the estimator

1 K K i=1 (γ A,i -γ A ) 2
is estimated from K = 10 independent runs of the algorithm and the truth γ A is taken from a long brute force Monte Carlo simulation. The number of interaction is set to 64.

Weight function

For simplicity, we have employed an exponential weight function W = exp(C∆E) where ∆E is the change in energy between two interactions. This choice doesn't require any a priori knowledge of the dynamics and is easy to implement. This weight function is not optimal, but, as we will show, it already gives good results. For the value of the forcing parameter R = 2 8 the distribution of the energy values is roughly Gaussian. One can therefore estimate the mean µ E and the variance σ 2 E from a brute force Monte Carlo simulation and use these values along with the reasoning of Section 2.2.1 to determine an appropriate value of the exponential weighting factor C in the weighting function W = exp(C∆E). One can then choose a value C = ∆µ E /σ 2

E

where ∆µ E is the desired change of the mean energy of the final particle distribution.

The values of C in the weighting function W = exp(C∆E) for the calculations presented in this section are taken as C r = r/(2σ E ) with r ∈ {1, 2, 3, , 4} and σ E being the standard deviation of the energy so as to increase the mean energy ∆µ E by steps of size σ E /2.

Noise perturbation

For deterministic dynamical systems, in order for two trajectories to have different dynamics after selection, a small perturbation can be added. This can be achieved by adding for example a weak Brownian perturbation at all times, or by adding a small instantaneous perturbation to offspring at the selection times. The former approach provides a simpler mathematical framework. Indeed the study of the noise effect would amount to the study of the stochastic differential equation properties in the weak noise limit, independent of the genealogical particle analysis algorithm. By contrast the latter approach intertwines the random perturbation with the genealogical particle analysis algorithm effects and is therefore more complicated to analyze. The latter approach, however, has the practical advantage of being computationally simpler. In this study, as we will proceed purely empirically, we have opted for the latter approach. The clones are perturbed by N 0,1 where N 0,1 is a standard L-dimensional Gaussian random variable, i.e. the noise acts independently on all of the variables.

The small noise perturbation invariably adds an error to the estimates of the tail probabilities. To obtain a rough upper bound on the strength of the perturbation that can be added without significantly perturbing the tail, we first perform a brute force simulation with the added noise for different noise strengths and verify that the tail probabilities do not change significantly compared to the sampling error of the brute force calculation. A set of independent realizations is performed like in the brute force Monte Carlo approach, the only difference being that at the selection times t k the same noise perturbations is added as in the genealogical particle analysis simulation. No selection is performed however in these runs. This way we can estimate the effect of the noise on the final time particle distribution. Figures 9a and9b show that below a perturbation strength of = 0.87 and for thresholds higher than a = 1600, the noise does not have a significant effect on the over-threshold probabilities. More complex schemes of noise perturbation could be implemented to assure that the perturbed trajectory remains close to the attractor, for example by storing a configuration at a time point before t k , adding a small perturbation to it and evolve it up to t k to have the perturbation relax towards the attractor.

Furthermore, after performing the genealogical particle analysis algorithm, we check that the perturbing noise intensity is small enough by decreasing and checking that the estimates of the over-threshold tail statistics are consistent. Figure 9b shows that for = 0.1 the results remain stable upon halving the noise intensity.

Estimating the tail of the energy distribution

We use for the Lorenz '96 model the procedure described in Section 4.2.3: we increase the values of the weight parameter C and use for each threshold value the best estimate from the point of view of the empirical estimator variance. The result is shown in Figure 10. As there is no analytic expression for the energy distribution tail of the Lorenz '96 system, we use a long brute force Monte Carlo estimation as comparison. The estimator variance markedly decreases when using the genealogical particle analysis algorithm. When constructing the overthreshold probability, we see that the tail can be reliably reproduced when compared to the longer brute force calculation.

The improvements in efficiency from using a rare event simulation scheme can be quantitatively estimated from Figure 10. The plot of the empirical relative error shows how for a threshold around a = 1800 a brute force Monte Carlo calculation yields a relative error of 0.5, whereas the genealogical particle analysis simulation yields a relative error of approximately 0.05. A reduction in relative error by a factor 10 is achieved. Since the brute force Monte Carlo error scales as 1/ √ M , a similar reduction by a raw increase of processing power would require M to increase by a factor of 100. For higher thresholds and with more fine tuning of the selection process, a much larger reduction is likely to be achievable. system with perturbations of varying strengths at times t k , without performing killing and cloning, with 10 000 independent realisations. The dotted line shows the estimated 2σ interval for over-threshold probability of the energy of the Lorenz '96 system, as estimated from two realizations of the genealogical particle analysis algorithm with different perturbing noise strength upon cloning.

Conclusion

In this paper we have addressed the use of rare event computation techniques to estimate small over-threshold probabilities of observables in deterministic dynamical systems. We have demonstrated that the genealogical particle analysis algorithms can be successfully applied to a toy model of atmospheric dynamics, the Lorenz '96 model as presented in Section 5.1. We have furthermore used the Ornstein-Uhlenbeck system to illustrate a number of implementation issues. The example of the Ornstein-Uhlenbeck has illustrated the importance of the choice of the objective function for the performance of the genealogical particle analysis algorithm estimator. We have shown how a time-dependent objective function based on the fluctuation path to a high threshold can greatly improve the performance of the estimator compared to a fixed-in-time objective function. Furthermore we have discussed how the number of particles and the number of selection steps influence the performance of the estimator.

For the deterministic chaotic system a complication arises in that a stochastic perturbation needs to be added to the system to make identical clones of one parent diverge and explore the system's path space. We have demonstrated in this example how the estimates of the rare event simulation are stable for small perturbations and agree with results from brute force Monte Carlo estimations. We therefore can have The over-threshold probability tail as estimated from the genealogical particle analysis algorithm compared to a brute force reconstruction. The number of particles used is 1, 000 for the genealogical particle analysis simulation and 10, 000 for the brute force simulation. confidence in the correctness of these estimates.

For the example of a deterministic chaotic system that we have studied we have not yet used the fluctuation path approach, since this would require information on the dynamics to the rare event that we a priori do not possess. This lack of knowledge can be improved by iterating the estimation procedure, where one uses estimates of an initial brute force simulation to estimate the fluctuation path, after which an genealogical particle analysis simulation based on this path can be used to estimate a higher fluctuation path, which can be used for a next iteration of the algorithm. However, the results of the straightforward implementation of the rare event simulation already shows significant improvements compared to brute force estimation.

Figure 2 :

 2 Figure 2: Schematic representation of the genealogical particle analysis

Figure 3 :

 3 Figure 3: (a) The estimator variance RE =

Figure 4 :

 4 Figure 4: The estimator variance RE =

Figure 5 :

 5 Figure 5: (a) The estimator variance RE for different weight factors C for a range of thresholds. The brute force error is computed as γ A -γ 2 A /( √ M γ A ) where A = [a, +∞). (b) The estimated over-threshold probability P (a) compared to the analytic result. For each value of the threshold a, the estimate corresponding to the value of C with lowest estimator variance is chosen.

Figure 6 :

 6 Figure 6: Ancestral paths (top), the final portion of killed paths, plotted only between t k-1 and t k if the path is killed at t k (middle) and the fraction of the number of particles killed N (-) k

Figure 8 :

 8 Figure 8: Over-threshold probability γ Et = P (E(t) > E t ) estimated from a brute force simulation and estimator variance RE = γ Et -γ 2 Et /( √ M γ Et ) of the energy E of the

Figure 9 :

 9 Figure9: (a) The over-threshold probability P (E > a) of the energy E of the Lorenz '96 system with perturbations of varying strengths at times t k , without performing killing and cloning, with 10 000 independent realisations. The dotted line shows the estimated 2σ interval for over-threshold probability of the energy of the Lorenz '96 system, as estimated from two realizations of the genealogical particle analysis algorithm with different perturbing noise strength upon cloning.

Figure 10 :

 10 Figure 10: (a) The empirical relative error RE for different weight factors C, for a range of energy thresholds a, for the Lorenz '96 model. The number of particles is M = 1, 000. The brute force Monte Carlo relative error (in blue) is estimated with the same number of realizations M as γa -γ2 a /(γ a √ M ) (b)The over-threshold probability tail as estimated from the genealogical particle analysis algorithm compared to a brute force reconstruction. The number of particles used is 1, 000 for the genealogical particle analysis simulation and 10, 000 for the brute force simulation.

Acknowledgements

The authors would like to thank Eric Vanden-Eijnden, Eric Simonnet, Joran Rolland, Pascal Yiou, Bérengère Dubrulle, Francesco Ragone and Takahiro Nemoto for fruitful discussions.

JW acknowledges the support of the AXA Research Fund.

The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement


T Naohara

H Aono

T Maehara

S Matsutomo

H Hirazawa

Y Watanabe

Simulation of Heat Generation Behavior in Ti-coated Ablation Needle for Cancer Therapy using High-frequency Induction Technique

Keywords: high-frequency induction technique, cancer therapy, ablation treatment, heat transfer simulation

published or not. The documents may come

Simulation of Heat Generation Behavior in Ti-coated Ablation Needle for Cancer Therapy using High-frequency Induction Technique

T Naohara, H Aono, T

Introduction

As is well known, the high-frequency induction technique finds its greatest application in the metal-processing industries. The primary uses are classified into the major categories of heating prior to metalworking, heat treatment, welding and metal melting. This technique offers a number of advantages over the conventional furnace process such as energy savings, rapid heating, and resultant high production rate [START_REF] Zinn | Elements of Induction Heating: Design, Control, and Applications[END_REF]. In addition to these industrial fields, another application of this technique has been sought and developed as a novel thermotherapy of solid cancer in the medical fields [START_REF] Watanabe | [END_REF]3]. Thus, we have attempted to evaluate the applicability of the novel ablation therapy for oral cancer utilizing a highfrequency induction technique [4]. As reported in a previous paper [5], we have fabricated a prototype ferromagnetic ablation needle, which is totally covered with a Ti layer, as demonstrated in Fig. 1 (a) and(b). According to experiments performed by changing the inclination angles (θ= 0 o , θ= 45 o and θ= 90 o ) relative to the magnetic flux direction, the newly designed ablation needle had a complete non-oriented heating property relative to the magnetic flux direction, suggesting that precise control of the ablation temperature is possible in clinical use without a lead-wire connected to a fiber-optic thermometer. The purpose of the present study was to simulate the heat generation property of a prototype ablation needle by a finite element method. In addition, the heat transfer simulation was also carried out to visualize the temperature distribution in the tissue-mimicking phantom subjected to the insertion of the Ti-coated ablation needle.

Magnetic field and heat transfer simulation analyses

The heating property of the Ti-coated ablation needle was investigated using the electromagnetic field analysis software, JMAG Studio, ver. 10.0 (The JRI Solution, Ltd.) [6]. Fig. 1(c) shows the appearance of the specimen used in the present simulation study. The two-dimensional schematic view of the simulation model is shown on the left-side of Fig. 2. The simulation model was completed by considering the physical parameters of the outer Ti component, inserted carbon steel rod and tissue-mimicking phantom, as well as the frequency, current condition, magnetic field strength, turn number of the coil and the initial temperature. These simulation conditions are summarized on the right-side of Fig. 2. The simulation was carried out for the three different inclination angles of θ= 0 o , 45 o and 90 o relative to the magnetic flux direction. The magnetic permeability value of the ferromagnetic carbon steel rod used in the simulation was 2000, while this value was 1 for the non-magnetic outer Ti component.

Heat generation property of the Ti-coated ablation needle

Fig. 3 shows the simulation images of the magnetic flux density for the prototype Ti-coated ablation needle subjected to the insertion of the ferromagnetic carbon steel rod under the θ= 0 o condition relative to the magnetic flux direction. In the color display bar presented on the right portion of the figure, the color transition (red→ green→ blue) reflects the continuous decreases in the value of the magnetic flux density, referring to the numerical value on the its left-side. As indicated by the red and dark-blue colors in Fig. 3(a), the magnetic flux density is the highest at the surface and in the vicinity of the inserted carbon steel rod. However, we note in Fig. 3(b) that no significant magnetic flux density is observed in the outer Ti component with the non-magnetic property.

In the high-frequency induction field [START_REF] Bozorth | Ferromagnetism[END_REF][START_REF] Saito | [END_REF], the magnetic total loss (P total ) is composed of the hysteresis loss (P h ), eddy c current loss (P e ) and residual loss (P r ) in the magnetic materials as given by Eq. [START_REF] Zinn | Elements of Induction Heating: Design, Control, and Applications[END_REF].

P total = P h + P e + P r [W/m 3 ] (1)
Moreover, the P h and P e values, denoted by Eqs. ( 2) and (3), respectively, significantly affect the heating property of the ferromagnetic metals and alloys. Since the magnetic flux concentrates only on the surface and in the vicinity of the inserted carbon steel rod, the effect of P h plays an important role in its heat generation property.

P h = k h • f • B m 1.6 [W/m 3 ] (2)
where k h : a constant value, f: frequency (Hz), B m : maximum magnetic flux density

P e = k e (t • f • B m ) 2 / ρ [W/m 3 ] (3)
where k e : a constant value, t: thickness (m), f: frequency (Hz), B m : maximum magnetic flux density (T), ρ: electrical resistivity (Ωm) Fig. 4 shows the simulation images of the eddy current density at 350kHz for the prototype Ti-coated ablation needle subjected to the insertion of the ferromagnetic carbon steel rod under the θ= 0 o condition relative to the magnetic flux direction. As indicated by the red and dark-blue colors in Fig. 4(a), the eddy current density is the highest at the surface and in the vicinity of the ferromagnetic carbon steel rod, but drastically decreases with the depth due to the skin effect, showing a similar tendency with the results in Fig. 3(a). It is appropriate to consider that the high eddy current flowing at the outer surface of the carbon steel rod causes a high heat generation ability through the effect of P e given by Eq. ( 3). Additionally, the eddy current of approximately 2.5×10 6 A/m 2 flows on the non-magnetic outer Ti component, as indicated by the light-green color in Fig. 4(b). The region colored red is observed in the upper portion having a greater diameter, since the high eddy current of 5.0×10 6 A/m 2 flows on this portion, where a notch is made to fit the claw of Kocher's forceps. Fig. 5 shows the simulation images of the Joule loss density at 350kHz for the prototype Ti-coated ablation needle subjected to the insertion of the ferromagnetic carbon steel rod under the θ= 0 o condition relative to the magnetic flux direction. In comparison to Fig. 4, it seems likely that the behavior of the Joule loss density is similar to that of the eddy current density. The high value of 10.0×10 6 W/m 3 is obtained at the surface and in the vicinity of the inserted carbon steel rod, as presented in Fig. 5(a), suggesting that the surface eddy current results in the considerable heat generation through the effect of P e denoted in Eq. ( 3). Moreover, it is obvious from Fig. 5(b) that the red color region is seen in the upper part of the Ti component possessing a greater diameter to some extent, corresponding to the gripping part of the prototype Ti-coated ablation needle. Fig. 6 shows the simulation images of the current density at different inclination angles (θ= 0 o , θ=45 o and θ=90 o ) relative to the magnetic flux direction for the prototype Ti-coated ablation needle. It is important to note that a fairly enlarged area perpendicular to the magnetic flux direction is achieved for the Ti-coated ablation needle with the increasing inclination angles up to θ= 90 o . The complete encapsulation provides an uninterrupted closed circuit for the outer Ti component surrounding the enlarged area of the ferromagnetic carbon steel rod. As denoted by the red color in the figure, a high eddy current flows near the surface of the outer Ti component, suggesting that P e plays a major role in the heating property for the θ= 90 o specimen. These results allowed us to conclude that the magnetic field analysis is an effective method to visually clarify the origin of the heat generation behavior in the Ti-coated ablation needle for oral cancer treatment.

Heat transfer property of the prototype Ti-coated ablation needle Fig. 7 shows the heat transfer simulation images from the prototype Ti-coated ablation needle to the tissue-mimicking phantom taken for the longitudinal section. The numerical value, which corresponds to the color transition, is presented as a color display bar on the right-side of the figure. As seen in Fig. 7(a), the temperature distribution is colored darkblue at the initial temperature of 20 o C which corresponds to no increase in temperature (ΔT). It is obvious from Fig. 7(b) that the generated heat transfers in the lateral direction during the induction time of 300s. Referring to the color display bar on the right-side, the heat transfer region is not very wide as the light-green color, and the ΔT value is at most 4.0~6.0 o C. Taking note of the longitudinal section, no significant heat transfer is achieved to the underside of the specimen, which is situated in the longitudinal direction. Thus, it is considered that no ablation effect is expected to the underside of prototype Ti-coated ablation needle, even though it is inserted into the tumors during clinical use. Fig. 8 shows the simulation images of the heat transfer from the prototype Ti-coated ablation needle to the tissue-mimicking phantom at the induction times of 600s, 900s and 1200s. When compared to Fig. 7, it is obvious that the heat transfer region gradually spreads with the increasing induction time. In Fig. 8(a), the specimen is colored red at the induction time of 600s, suggesting that the ΔT value reaches approximately 10.0 o C, as demonstrated by the color display bar. Moreover, the region colored light-green significantly spreads in the tissue-mimicking phantom, as compared to the image obtained at the

Conclusion

With the aim of studying the applicability for the novel ablation therapy of oral cancer, the heat generation and transfer behaviors of a prototype Ti-coated ablation needle subjected to the insertion of a carbon steel rod was simulated in a high-frequency induction field at 350kHz. According to the magnetic field analysis, both the magnetic flux density and eddy current density were the highest at the surface and in the vicinity of the inserted ferromagnetic carbon steel rod.

In the heat transfer simulation images, the temperature distribution was visualized in the longitudinal section of the tissue-mimicking phantom. The heat transfer concentrically occurred in the lateral direction of the tissue-mimicking phantom. However, no significant heat transfer was achieved in the underside of the specimen, which is situated in the longitudinal direction.

Fig. 1 :

 1 Fig.1: Appearance of the inserted ferromagnetic carbon steel rod (a), the newly designed prototype Ti-coated ablation needle (b) and the specimen (c) used in the present simulation study.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Simulation images of magnetic flux density at 350kHz from the inserted carbon steel rod (a) and prototype Ti-coated ablation needle (b) under the θ= 0 o condition. Fig. 4: Simulation images of eddy current density at 350kHz from the inserted carbon steel rod (a) and prototype Ti-coated ablation needle (b) under the θ= 0 o condition.

Fig. 2 :

 2 Fig. 2: Simulation model of heat generation and transfer properties used in the present study.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Simulation images of Joule loss density at 350kHz from the inserted carbon steel rod (a) and prototype Ti-coated ablation needle (b) under the θ= 0 o condition. Fig. 6: Simulation images of current density at 350kHz from the prototype Ti-coated ablation needle with different three inclination angles relative to the magnetic flux direction.

Fig. 7 :Fig. 8 :

 78 Fig. 7: Heat transfer simulation images showing the temperature distribution at 0s (a) and 300s (b) in the longitudinal section of the tissue-mimicking phantom.Fig. 8: Heat transfer simulation images showing the temperature distribution at 600s (a), 900s (b) and 1200s (c) in the longitudinal section of the tissue-mimicking phantom. the induction time of 300s given in Fig. 7(b). It is noted in Fig. 8(b) that the region colored red spreads from the specimen to the tissue-mimicking phantom, accompanying the enlargement of the right-green region at the induction time of 900s. After the induction time of 1200s, the generated heat significantly transfers over a wide region of the tissue-mimicking phantom, as seen in Fig. 8(c).

Acknowledgment

The present study was supported by a Grant-in-Aid from The Ministry of Education, Science, Sport and Culture of Japan (No.23500559: T. Naohara).


D Zhang

L Nastac
email: lnastac@eng.ua.edu

Numerical Modeling of the Nanodispersion during Ultrasonic and Electromagnetic Stirring Processing of Al-Based Nanocomposites

Keywords: Ultrasonic and electromagnetic stirring, Numerical modeling of dispersion of ceramic nanoparticles, 6061 nanocomposites, SiC nanoparticles, Molten metal processing, Unidirectional solidification

The metal-matrix-nano-composites (MMNCs) in this study consist of a 6061 alloy matrix reinforced with 1.0 wt.% SiC 50 nm diameter nanoparticles that are dispersed within the molten alloy matrix using ultrasonic cavitation and induction melting technologies. The required ultrasonic parameters to achieve the required cavitation for adequate degassing and refining of the aluminium alloy as well as the fluid flow characteristics for uniform dispersion of the nanoparticles into the 6061 matrix are being investigated in this study by using an inhouse developed magneto-hydro-dynamics model. The magneto-hydro-dynamics model accounts for turbulent fluid flow, heat transfer and solidification, electromagnetic field as well as the complex interactions between the solidifying alloy and nanoparticles by using ANSYS Maxwell and ANSYS Fluent Dense Discrete Phase Model (DDPM) and a particle engulfment and pushing (PEP) model. The PEP model accounts for the Brownian motion. Scanning electron microscope (SEM) analysis was performed on the as-cast MMNC coupons processed via ultrasonic cavitation dispersion technique (UCDS) and confirmed the distribution of the nanoparticles predicted by current model.

Introduction

Recently, to overcome the counterpart that the ductility of the aluminium-based metal matrix composites (MMCs) deteriorates with high ceramic particle concentration [1][2][3][4][5], more attention has been drawn to metal matrix nanocomposites (MMNCs), since the properties of metallic alloys reinforced by ceramic nanoparticles (with dimensions less than 100 nm) would be enhanced considerably while the ductility of the matrix is retained [6][7][8][9][10][11][START_REF] Ferguson | Light Metals TMS[END_REF]. Currently, several fabrication technologies including high-energy ball milling [9,[START_REF] Ferguson | Light Metals TMS[END_REF], in-situ synthesis [8], electroplating [START_REF] Sautter | [END_REF], and ultrasonic technology (UST) [5,6] are most commonly used, among which UST is supposed to be more reliable and cost effective.

During the induction heating process, the dominant electromagnetic field will generate heat (affecting the temperature field) and momentum (controlling the fluid flow field). Meanwhile, all these three fields may strongly influence one another, which may significantly complicate the problem [START_REF] Mach | International Scientific Colloquium[END_REF]. An improved understanding of the heat and electromagnetically driven flow mostly requires mathematical modeling, as experimental studies, especially for high temperature materials processing, may be restricted because of the physical and chemical properties of the melts [START_REF] Schwarze | [END_REF].

During the solidification process, the particles will either be pushed or engulfed by the solidification front. Traditional PEP models [16][17][18][19][20] cannot explain the evidence in MMNCs that nanoparticles can indeed be engulfed and distributed throughout the material. As proposed by Ferguson [START_REF] Ferguson | Light Metals TMS[END_REF], for sufficiently small particles, "Brownian Motion" can partially or completely counteract forces such as viscous drag, gravity and thermal/concentration gradients, thus leading to engulfment rather than pushing.

In this study, ANSYS Fluent Dense Discrete Phase Model (DDPM) [START_REF]Fluent 6.3: User's Guide Manual Fluent Inc. and Ansys's Fluent[END_REF] was adapted. The DDPM accounts for turbulent fluid flow, heat transfer, electromagnetic forces, and complex interactions between the molten alloy and nanoparticles during the melting and unidirectional solidification processes. Based on the theory proposed by Ferguson, all of the nanoparticles are assumed to be engulfed by the solidification front and no entrapment will occur. The dispersion of SiC nanoparticles with and without induction stirring has been investigated in detail.

Model Description

The geometries of the electromagnetic field and CFD model are shown in Due to harmonic nature of EM field and induced eddy currents, the Lorentz force F Lor can be decomposed into a steady and harmonic part that oscillates with double frequency. Because of much greater inertia times of melt in comparison to alternate EM field timescale, only the steady part of the Lorentz force is taken into account [START_REF] Spitans | [END_REF].

Assuming the fluid flow does not affect the electromagnetic field, the Lorentz force can be obtained according to Ampere's law, Faraday's law, and magnetic flux continuity law [START_REF] Pal | COMSOL Users Conference[END_REF].

The fluid flow model, solidification model, particle tracking model, and solution procedure are presented in detail in [START_REF] Zhang | [END_REF][START_REF] Zang | Presented at Modeling of Casting, Welding, and Advanced Solidification Processes (MCWASP XIV)[END_REF].

Simulation Results and Discussion

Fig. 2 shows the fluid flow and particle distributions after 1.0s, 3.0s, 10s and 30s, respectively, when the induction stirring is ignored [START_REF] Zang | Presented at Modeling of Casting, Welding, and Advanced Solidification Processes (MCWASP XIV)[END_REF]. It can be seen that the flow is much stronger at the center of the furnace. It demonstrates that the particles have little effect on the fluid flow because of the one-way coupling, and the gravity direction is the same as the direction of the stronger flow in the center, so it will help to develop the fluid flow. Therefore the fluid flow after 3.0s is almost the same as that after 1.0s because both of them are fully developed (see Figs. 2(a) and 2(b)). The fluid flow is damped in the mushy zone, where the particles will move relatively slow and stop moving in the solid (see Figs. 2(c) and 2(d)). At 3.0s, the particles have been dispersed pretty well from the bottom to the top, but a little more particles tend to stay near the wall, and fewer nanoparticles at the center where the flow is stronger, which indicates the nanoparticles couldn't disperse well in strong flows (see Fig. 2(f)). Additionally, the particle distribution stays almost the same henceforth which means the distribution becomes stable. Then, the solidification process starts at 3.0s, and in order to model the unidirectional solidification process, all the walls except the bottom one are treated as adiabatic. After 10s, we can see that there are some agglomerations in the solid (see Fig. 2(g)) because the nanoparticles approaching the solid/liquid interface will be engulfed based on the assumption that the Brownian motion will dominate during the particle engulfment/pushing process. When the induction stirring is taken into consideration, the fluid flow and particle distributions after 1.0s, 3.0s, 10s, and 30s are shown in Fig. 3, respectively. Compared to Fig. 2, the fluid flow pattern is basically the same, because the induced flow by induction stirring is relatively weak compared to that induced by ultrasonic stirring. But there are still some differences especially at the center of the furnace that the flow at the bottom becomes weaker than that at the top because of the induction stirring. It can also be seen that there are considerable number of particles in the center of the furnace, which indicates that the particle distribution becomes even better than before. It can be explained that induction stirring modifies the fluid flow, making it more uniform. When solidification starts, the fluid flow is damped in the mushy zone as well, and particles approaching the solidification front are captured. As it can be seen in Fig. 3(h) that, after 30s, the nanoparticles are dispersed quite well throughout the metal matrix, which is also validated by SEM analysis [START_REF] Liu | [END_REF]27].

Concluding Remarks

A magneto-hydro-dynamics model was developed in ANSYS Maxwell and ANSYS Fluent. The DDPM model coupled with the -turbulence model was used to investigate the distribution of SiC nanoparticles into a 6061 matrix under ultrasonic and induction stirring and molten metal and unidirectional solidification processing conditions. It was determined that a strong flow can disperse the nanoparticles relatively well. Also, it was demonstrated that induction stirring can help improve the fluid flow characteristics, making the nanoparticle distribution even more uniform.

Fig. 1 .Fig. 1 :

 11 Fig. 1: Geometry model.The ultrasonic probe has a diameter of 40 mm. The liquid aluminium alloy is 6061. It has a density of 2685 kg/m 3 . The SiC nanoparticles with an average particle size of 55 nm and density of 3216 kg/m 3 are treated as inert-particles. The mass flow rate of the SiC nanoparticles is 0.014 kg/s. Thus, 1.0 wt. % of SiC nanoparticles can be injected at about 15 mm beneath the ultrasonic probe for 1 sec. The multiphase CFD model accounts for turbulent fluid flow, heat transfer, and the complex interaction between the molten alloy and nanoparticles by using the ANSYS Fluent DDPM and -turbulence model[START_REF]Fluent 6.3: User's Guide Manual Fluent Inc. and Ansys's Fluent[END_REF]. The CFD model is described in detail below.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Fluid Flow and particle distributions (without induction stirring): (a) and (e) after 1.0s; (b) and (f) after 3.0s; (c) and (g) after 10s; and (d) and (h) after 30s.


B Paya

P Teixeira

Electromagnetic characterization of magnetic steel alloys with respect to the temperature

Keywords: Magnetic properties, electrical properties, temperature dependence, numerical modeling, measurement

Because of a lack of data, numerical models of induction heating use approximations to describe the temperature dependence of the magnetization curve. Thanks to our measuring equipment acquired in 2013, we are now able to measure the true behavior of ferrous alloys from room temperature to 1,200 °C, sometimes far away from the assumptions commonly used in many software. Key physical quantities and their evolution with temperature are obtained: saturated magnetization, maximal susceptibility and resistivity. They are compared to analytical approximations proposed by commercial software which are not in agreement with the measurements. New analytical approximations are proposed. We also describe the strange behavior of some special alloys designed for specific mechanical properties. The electromagnetic properties are not reversible with regard to the temperature, describing a kind of hysteretic loop during heating and cooling phases.

Introduction: lack of physical data for numerical simulation of induction heating

Numerical modeling of induction heating device often suffers from a lack of physical data, especially the temperature dependence of the magnetization curves or the resistivity. To fill this gap, EDF acquired a unique equipment supplied by the German company Magnet-Physik Steingroever GmbH, able to measure the magnetization curve and the resistivity during the same thermal cycle up to 1,200 °C. After a short description of the measuring device, we present first the results of the characterization of 38MnSiV5 magnetic steel. The measures are compared to analytical approximations commonly used in commercial software. We also show unexpected behavior of special alloys designed for specific mechanical properties.

Description of the measurement equipment

Fig. 1: schematic overview of the measurement unit This equipment presented in detail in [START_REF] Paya | Measurement of Electrical and Magnetic Properties of Steels at Elevated Temperature[END_REF] is consisting of 4 main functions (Fig. 1): -Permeameter: device providing the measurement of magnetic polarization curves. A bar of the specimen is placed inside two concentric coils: the exciting coil providing the magnetic field thanks to a controlled current source, and the measuring coil collecting the magnetic flux to the flux meter. The measuring coil is specially designed in a socalled "J-compensated configuration" that compensates the air magnetization and gives directly the magnetization of the sample [START_REF] Steingroever | Magnetic measuring techniques[END_REF]. -Resistivimeter: device providing the measurement of electrical resistivity using a 4-points configuration. Type K thermocouples are spot-welded on both ends of the specimen and connected to the current source and the microvoltmeter. -Furnace: device providing the heating from room temperature up to 1,200 °C. The heating element is a bifilar platinum resistor configured to minimize magnetic field generation when heating. A water jacket prevents the permeameter's coils from heating. Furnace temperature is controlled by a PID regulation during heating and cooling phases. -Driving and supervising device: computer providing manual or automatic driving of the other equipments, the follow-up of the running measurement campaign and the transfer of measures to other computers. Dedicated software is developed for this purpose. The different temperature steps are defined before the campaign starts. At each step, the resistivity measurement is done first and the magnetization curve, then. Measurement results are available in text files, Excel® files or by magnetization curves drawing.

Example of a measurement campaign: characterization of 38MnSiV5 steel

The specimen to be characterized had a parallelepiped shape 4x4x150 mm 3 . Two thermocouples used for the 4 points resistivity measurement are spot-welded on both ends. Another thermocouple is also spot-welded in the middle for the sample temperature measurement. The heating cycle starts from room temperature up to 1,000 °C, then down to room temperature: the whole campaign lasts about 6 hours. Two samples of the same alloy are characterized to make sure of the reproducibility of the measures. Magnetization curves and electrical resistance are recorded at each temperature step. From these measures, physical quantities are calculated: saturated magnetization, maximal susceptibility and resistivity. These values are plotted in the following figures. The shape of the dot corresponds to the characterized sample; the red and blue colors correspond respectively to the heating and cooling phase.

Analytical approximation of the physical properties

Commercial software, such as FLUX®, proposes analytical models to describe the physical properties of material [START_REF] Cedrat | User guide Flux® 12. Volume 2 Physical description, solving & postprocessing[END_REF]. So, the magnetization J can be described using a univocal relationship with the magnetic field H and the temperature T:

𝐽 𝐻, 𝑇 = 2 𝐽 𝑠𝑎𝑡 𝜋 • tan -1 𝜒 𝑚0 𝜋 2 𝐽 𝑠𝑎𝑡 • 𝜇 0 𝐻 • 𝐹 𝑇 (1)
and 𝐹 𝑇 = 1 -𝑒 -𝑇-𝑇 𝑐 𝐶 for 𝑇 ≤ 𝑇 𝐶 (Curie temperature)

In this case, both the saturated magnetization 𝐽 𝑠𝑎𝑡 𝑇 and the maximal susceptibility 𝜒 𝑚 𝑇 follow the same evolution with the temperature corresponding to the equation [START_REF] Steingroever | Magnetic measuring techniques[END_REF]. The corresponding curves are drawn in orange in Fig. 2 and Fig. 3; as you can see, they do not fit with the experimental data. To better describe the temperature behavior, we propose the following relationships:

𝐽 𝐻, 𝑇 = 2 𝐽 𝑠𝑎𝑡 𝑇 𝜋 • tan -1 𝜒 𝑚 0 𝑇 𝜋 2 𝐽 𝑠𝑎𝑡 𝑇 • 𝜇 0 𝐻 (3)
with

𝐽 𝑠𝑎𝑡 𝑇 = 𝐽 𝑠𝑎𝑡 0 • 1 -𝑎 𝑇 • 1 -𝑒 - 𝑇-𝑇 𝑐 𝐶 𝐽 (4)
and

𝜒 𝑚 𝑇 = 𝜒 𝑚 0 -𝜒 𝑚 𝑚𝑖𝑛 • 𝑇 𝑇 𝑚𝑖𝑛 -1 2 + 𝜒 𝑚 𝑚𝑖𝑛 • 1 -𝑒 - 𝑇-𝑇 𝑐 𝐶 𝜒 (5)
The corresponding curves are drawn in green in Fig. 2 and Fig. 3. The saturated magnetization curve fits very well with the experiments. The susceptibility curve fits better than the previous one especially at temperature close to Curie temperature. One must keep in mind that the measure of the susceptibility with our equipment is not very accurate.

Commercial software do not propose analytical model of the resistivity which present an inflection point; the closest model available has an exponential shape.

𝜌 𝑇 = 𝜌 𝑎 + 𝜌 𝑏 • 1 -𝑒 -𝑇-𝑇 0 𝜏 (6)
When using this equation ( 6), it is possible to describe the resistivity, either in the low (Fig. 4, purple curve) or high (Fig. 4, orange curve) temperature range. Unfortunately, the continuity of the function and its derivate is not respected, especially close to the transition temperature T trans . We propose another way to write this relationship:

𝜌 𝑇 = 𝜌 0 • 𝛼 0 𝑒 𝑇 𝜏 0 + 1 -𝛼 0 for 𝑇 ≤ 𝑇 𝑡𝑟𝑎𝑛𝑠 𝜌 𝑇 = 𝜌 ∞ • 𝛼 ∞ 𝑒 - 𝑇 𝜏 ∞ + 1 for 𝑇 ≥ 𝑇 𝑡𝑟𝑎𝑛𝑠 (7)
The continuity of the function and its derivate is obtained by verifying the two following equations giving relationship between the low and high temperature parameters. The corresponding curve drawn in green in Fig. 4 better agrees with the experiments.

𝜌 0 𝛼 0 𝜏 0 𝑒 𝑇 𝑡𝑟𝑎𝑛𝑠 𝜏 0 = - 𝜌 ∞ 𝛼 ∞ 𝜏 ∞ 𝑒 - 𝑇 𝑡𝑟𝑎𝑛𝑠 𝜏 ∞ 𝜌 ∞ = 𝜌 0 𝛼 0 𝑒 𝑇 𝑡𝑟𝑎𝑛𝑠 𝜏 0 • 1 + 𝜏 ∞ 𝜏 0 + 1 -𝛼 0 (8)
The temperature hysteresis: an unexpected behavior Some special steel alloys are designed for specific mechanical properties. When characterizing them for numerical modeling of induction heating, we observed an unexpected behavior. The electric and magnetic properties were sensitive to the history of the heating. During the heating phase, the specimen loses its magnetic property at a high Curie temperature but it recovers it more slowly during cooling phase. Resistivity evolution has also a hysteretic behavior almost in the same temperature range. To show an example of this phenomenon, characterization of 16MND5 steel is presented in Fig. 5 to Fig. 7. In this case, the resistivity presents a hysteresis between Curie temperature (750 °C) down to 300 °C approximately. At room temperature, the saturated magnetization is almost the same before and after the thermal cycle but the maximal susceptibility is much lower: the heat treatment has shifted down the magnetization curve.

To measure relevant properties for numerical simulation, we should so take care of the preparation of the specimen. It should be extracted from the piece before the heating process. We should avoid heating during the cutting process: electro-erosion is a relevant way to cut the piece without heating it. Otherwise, the measured characteristics may differ from reality and give wrong results in numerical simulation.

Conclusion:

Thanks to our measuring equipment acquired in 2013, we are now able to measure the true behavior of ferrous alloys from room temperature to 1,200 °C. Key physical quantities and their evolution with temperature are obtained: saturated magnetization, maximal susceptibility and resistivity. Analytical approximations are proposed for numerical models. We also show that some special alloys designed for specific mechanical properties present a kind of hysteretic loop respect to the temperature during heating and cooling phases. Consequently, caution should be exercised for the material characterization.

Fig. 2 :

 2 Fig. 2: saturated magnetization of 38MnSiV5 Fig. 3: Maximal susceptibility of 38MnSiV5

Fig. 5 :Fig. 6 :Fig. 7 :

 567 Fig. 5: saturated magnetization of 16MND5

Acknowledgment

Characterization of MnSiV5 steel takes place in Optipro-Indux Project (ANR-2010-RMNP-011), with financial support of the French Research National Agency (ANR).