Rémi Helleboid
email: remi.helleboid@gmail.com

Denis Rideau

Isobel Nicholson

Jeremy Grebot

Bastien Mamdy

Gabriel Mugny

Marie Basset

Megan Agnew

Dominique Golanski

Sara Pellegrini

Jérôme Saint-Martin

Marco Pala

Philippe Dollfus

A Fokker-Planck-based Monte Carlo method for electronic transport and avalanche simulation in single-photon avalanche diodes

Keywords: Advection-Diffusion Monte Carlo (ADMC), avalanche breakdown probability, electronic avalanche, electronic transport, drift-diffusion, jitter, Monte Carlo simulation, singlephoton avalanche diode (SPAD), technology computed aided design (TCAD), Fokker-Planck

We present an efficient simulation method for electronic transport and avalanche in single-photon avalanche diodes (SPAD). Carrier transport is simulated in the real space using a particle Monte Carlo approach based on the Fokker-Planck point of view on an advection-diffusion equation, that enables us to reproduce mobility models, including electric fields and doping dependencies. The avalanche process is computed thanks to impact ionization rates implemented using a modified Random Path Length algorithm. Both transport and impact ionization mechanisms are computed concurrently from a statistical point of view, which allows us to achieve a full multi-particle simulation. This method provides accurate simulation of transport and avalanche process suitable for realistic three-dimensional SPADs, including all relevant stochastic aspects of these devices, together with a huge reduction of the computational time required, compared to standard Monte Carlo methods for charge carrier transport. The efficiency of our method empowers the possibility to precisely evaluate SPADs figures of merit and to explore new features that were untrackable by conventional methods. An extensive series of comparisons with experimental data on state-of-the art SPADs shows a very good accuracy of the proposed approach.

I. INTRODUCTION

Single-photon avalanche diodes are avalanche photodiodes, operated above breakdown voltage, in which a single photon absorption can lead to a self-sustained electron-hole avalanche in the junction. These devices are used in many applications such as Fluorescence Lifetime Imaging Microscopy [START_REF] Perenzoni | A 160 x 120 Pixel Analog-Counting Single-Photon Imager With Time-Gating and Self-Referenced Column-Parallel A/D Conversion for Fluorescence Lifetime Imaging[END_REF], Timeof-Flight (TOF) and Light Imaging Detection and Ranging (LiDAR) [START_REF] Bronzi | 100 000 Frames/s 64 × 32 Single-Photon Detector Array for 2-D Imaging and 3-D Ranging[END_REF], [START_REF] Lindner | A High-PDE, Backside-Illuminated SPAD in 65/40-nm 3D IC CMOS Pixel With Cascoded Passive Quenching and Active Recharge[END_REF]. Predictive numerical simulation of SPAD behavior is a key step to improve the performance of these devices in terms of output parameters such as Photon Detection Efficiency (PDE), timing jitter, afterpulsing and dark-countrate [START_REF] Bronzi | SPAD Figures of Merit for Photon-Counting, Photon-Timing, and Imaging Applications: A Review[END_REF]. A great deal of effort has already been dedicated to finding accurate methods to simulate the SPAD operation and to extract reliable electrical characteristics of a given device architecture without the need of expensive and timeconsuming manufacturing and characterization. The state-ofthe-art methodologies rely on solving the Boltzmann transport equation (BTE) by means of particle Monte Carlo simulation in the phase-space, either under frozen field to investigate the details of avalanche process [START_REF] Petticrew | Simple Monte Carlo Simulator for modelling Linear Mode and Geiger Mode Avalanche Photodiodes in C++[END_REF], [START_REF] Dolgos | Full-band Monte Carlo simulation of high-energy carrier transport in single photon avalanche diodes: Computation of breakdown probability, time to avalanche breakdown, and jitter[END_REF] or self-consistently coupled to Poisson's equation to study all stages of SPAD operation [START_REF] Cazimajou | Quenching Statistics of Silicon Single Photon Avalanche Diodes[END_REF]. This stochastic approach of transport can reproduce the statistical behavior of the SPAD, but it is known to be computationally very intensive, which makes it unsuitable for SPADs with sizes of several micrometers. To address this issue, many empirical models have been developed to extract SPAD's characteristics [START_REF] Oldham | Triggering phenomena in avalanche diodes[END_REF]- [START_REF] Gao | 3D Electro-optical Simulations for Improving the Photon Detection Probability of SPAD Implemented in FD-SOI CMOS Technology[END_REF]. While these models have proven themselves to be accurate for simple architectures, they are not able to reproduce the intrinsic statistical behavior of SPAD devices. Recently, some attempts to address a statistical calculation with a reduced simulation time have been made [START_REF] Rosset | A model for the jitter of avalanche photodiodes with separate absorption and multiplication regions[END_REF], [START_REF] Pilotto | An improved random path length algorithm for p-i-n and staircase avalanche photodiodes[END_REF], but at the cost of strong simplifications such as the need to artificially split the device into an absorption and a multiplication region. These models assume that the electric field is constant in the multiplication region and treat the avalanche as a one-dimensional process. Hence, such methods cannot succeed in simulating state-of-the-art SPAD devices in which the electrostatics is highly non-uniform and bent in order to increase the active absorption volume [START_REF] Van Sieleghem | A Near-Infrared Enhanced Silicon Single-Photon Avalanche Diode With a Spherically Uniform Electric Field Peak[END_REF]. This work aims at gathering efficient empirical models within a stochastic Monte Carlo method, able to simulate accurately and efficiently the statistical figures of SPADs of any architecture, without any simplification of the device geometry. Our approach relies on interpreting a transport advectiondiffusion equation as a Fokker-Planck equation and to solve it for each particle with a Monte Carlo method operating only in the real space. Concurrently, the impact ionization process is considered through a probabilistic Random Path Length algorithm, which randomly creates new particles, according to an ionization coefficient model. The coupling of these two stochastic numerical processes leads to a multi-particle simulation that allows us to extract figures of interest for the SPAD operation such as avalanche breakdown probability and timing jitter.

The strength and accuracy of our method is underlined and confirmed by comparing numerical results with characterization results for state-of-the-art SPADs. In what follows, for clarity purposes, our Monte-Carlo approach will be referred to as the Advection-Diffusion Monte Carlo (ADMC).

II. MODEL

A. Transport model

Carrier transport in semiconductor devices has been widely studied through the Boltzmann transport equation, which describes the evolution of the density of carriers in real and phase spaces:

∂f ∂t (r, k, t) + v • ∇ r f + 1 F • ∇ k f = Q coll f (1)
where f (r, k, t)dr dk dt is the probability of finding a particle, at time t, position r, with a wave vector k within a small volume dr × dk × dt. v is the group velocity at position r and time t, F is the total external force applied to the carriers and Q coll is the collision integral operator. Although this equation was very successful to model all regimes of semi-classical transport, its integration by the means of BTE Monte Carlo methods is known to be computationally heavy, which may prevent its use in many-particle simulations of large devices. The complexity of this phase and real space equation can be lowered by integrating over all the possible k vectors, in the reciprocal space. This leads to the balance equation, which can, at the cost of some simplifications be turned into the well-known Jacoboni drift-diffusion equation [START_REF] Jacoboni | Theory of Electron Transport in Semiconductors[END_REF]:

∂f ∂t (x, t) = -v 0 ∂f ∂x (x, t) + D 0 ∂ 2 f ∂x 2 (x, t) (2)
where v 0 is the carrier drift velocity, D 0 the carrier diffusion coefficient and f (x, t) the carrier density probability, said differently, if the studied system is made of N e electrons, then the local charge density can be extracted as:

n(x, t) = N e f (x, t)
This equation ( 2) is valid for a bulk semiconductor with constant doping level and uniform electric field, in a onedimensional simulation. Yet, it can be generalized to a threedimensional device simulation with varying-in-space doping level and electric field, and then becomes:

∂f ∂t (r, t) = -div(v(r, t)f (r, t)) + div(D(r, t)∇f (r, t)) (3)
with v(r, t) and D(r, t), the space and time dependent driftvelocity and diffusion coefficient at, respectively, position r and time t. The function f in (2) and ( 3) is basically the same as in (1) but the information on the wave vector k is lost, because f was integrated over all the reciprocal space.

To be complete, one must add initial and boundary conditions to this equation. The initial condition can be any density probability function, which represents the distribution of charges at the beginning of the simulation. The boundary conditions describe the behavior of charges at specific boundaries of the simulation domain. They can be either reflecting boundary conditions, for example at interfaces with insulators, or absorbing/emitting boundary conditions at metallic contact interfaces. Usually, the drift-diffusion equation ( 3) is solved as a standard partial differential equation, by the means of deterministic numerical methods such as finite difference, finite element or box methods [START_REF]The Stationary Semiconductor Device Equations[END_REF]. Those methods were explored in-depth to achieve good accuracy while keeping computational time within a reasonable range. Yet, in SPAD simulations, analysis of the statistical behavior is a crucial aspect of device operation, which cannot be achieved using deterministic methods.

Our model relies on the reformulation of the generalized equation (3) as a stochastic differential equation (SDE). Indeed, (3) can be interpreted as a particular form of a Fokker-Planck equation, which describes the spatio-temporal evolution of the probability density function of a continuous random process where both diffusive and dragging forces are present.

A time-continuous random process is a stochastic process, i.e. a sequence of random numbers, indexed by a continuous set of values, it will be written, for example, (Z t ) t≥0 . One can define the probability density of such a process at a given time t by a function g(•, t), for which g(r, t)dr is the probability for the process to occur in a small volume dr at position r at time t. Said differently, if we define the probability law P Zt of the position of the stochastic process of (Z t ) at time t, as:

P Zt : A → P(Z t ∈ A) (4)
with A a subset of R 3 , the probability density function of the random variable Z t , named g(•, t), associated to P Zt is a function such as:

P Zt (A) = A g(ν, t)dν (5)
We now expose the link between (3) and a time-continuous stochastic process. The so-called generalized drift-diffusion can be interpreted as a Fokker-Planck equation, for which solutions are temporal probability density function of the stochastic process (X t ), defined as follows [START_REF] Pavliotis | Stochastic Processes and Applications : Diffusion Processes, the Fokker-Planck and Langevin Equations[END_REF]:

dX t = v(X t , t)dt + σ(X t , t)dB t (6)
The process X t can be directly interpreted as the position of a particle in real space. We thus have X t the random variable representing the position of the particle at time t in the real space, v(X t , t) the drift velocity of the particle at its position and at time t, σ(X t , t) the diffusion matrix of the dynamics at position X t and time t and B t is the three-dimensional standard Brownian motion. The diffusion factor of ( 6) is linked to the electronic diffusion coefficient through the following relation:

σ(X t , t) = 2D(X t , t).
The differential dynamics (6) can be reformulated under its integral form:

X t = t 0 v(X s , s)ds + t 0 σ(X s , s)dB s (7)
The Brownian motion (B t ) is a time-continuous random process describing a succession of infinitesimal random steps that follow a centered Gaussian law with variance dt. Each step is independent of the previous ones, which results in a continuous random walk with no bias and a variance: Var(B t ) = t. It is possible to construct a Brownian motion of any dimension by stacking independent one-dimensional Brownian motions in a vector of the desired dimension.

Our model relies on the following proposition: a process defined by the SDE ( 6) or [START_REF] Cazimajou | Quenching Statistics of Silicon Single Photon Avalanche Diodes[END_REF] has a probability density function, as defined in [START_REF] Petticrew | Simple Monte Carlo Simulator for modelling Linear Mode and Geiger Mode Avalanche Photodiodes in C++[END_REF], that is a solution of the PDE defining the generalized drift-diffusion [START_REF] Lindner | A High-PDE, Backside-Illuminated SPAD in 65/40-nm 3D IC CMOS Pixel With Cascoded Passive Quenching and Active Recharge[END_REF]. Proof of this statement and of existence of such density function can be found in [START_REF] Gallager | Stochastic Processes: Theory for Applications[END_REF]. The particular form of the Fokker-Planck equation we are using can also be found under the name "Kolmogorov forward equation" which is often treated as a so-called "Kolmogorov backward equation with reverse time". For this proposition to be true, the boundary and initial conditions of the PDE and of the SDE must match. The most common initial condition for the stochastic process is to set deterministically its value at time t = 0, for example Z t=0 = r 0 , which amounts setting the initial condition of (3) to the Dirac delta function :

f : (r, t = 0) -→ δ r0 (r) (8)
The process (6) can hence be seen as a random walk, with standard deviation σ(X t , t), which is linked to the diffusion of carriers, and biased by the drift velocity v(X t , t) due to external forces. This view enables us to see the random process (X t ) t≥0 as the position of a carrier over time.

We have thus built a path between the deterministic macroscopic description of the density of carriers (with the driftdiffusion equation) and a stochastic description of individual carriers submitted to drift and diffusion. This is of first importance in the context of SPAD simulation, where the modeling of the statistical behavior of the device, resulting from the stochastic behavior of carriers, is mandatory to estimate all their relevant figures of merit.

In the present study, we have assumed that the transport is isotropic, and thus that the velocity and the diffusion coefficient can be reduced to single scalar functions that depend on the electric field strength and the doping concentration, through carrier mobility. We show in section III that this choice is reasonable for device simulation of Silicon-based SPAD devices at room temperature. Yet, for other applications or simulation parameters, such as cryogenic temperatures, where the transport might be highly anisotropic, one can keep the diffusion coefficient under a full 3x3 matrix-form and adjust the velocity with respect to the crystallographic orientation of the electric field [START_REF] Jacoboni | The Monte Carlo method for the solution of charge transport in semiconductors with applications to covalent materials[END_REF].

The drift-velocity for electrons and holes is computed by using the field, temperature and doping dependent model of Arora [START_REF] Arora | Electron and hole mobilities in silicon as a function of concentration and temperature[END_REF], together with a high-field saturation regime model from Canali's work [START_REF] Canali | Electron and hole drift velocity measurements in silicon and their empirical relation to electric field and temperature[END_REF]. Furthermore, the Einstein relation is assumed to hold:

D e,h ( F , N total ) = µ e,h ( F , N total )k B T q (9)
with D e,h the diffusion coefficient, µ e,h the charge mobility, for electron (e) and hole (h), k B the Boltzmann constant, T the temperature and q the elementary charge. The mobility depends on the electric field strength F and the total doping concentration N total . The diffusion coefficient thus depends on the local electric field and doping concentration, through the mobility. If necessary, one can straightforwardly replace the Einstein relation by a more sophisticated model for diffusion in order to take into account, for example, high-field effects [START_REF] Cheung | The diffusion of electrons in semiconductors in high electric fields[END_REF], [START_REF] Arora | Drift diffusion and Einstein relation for electrons in silicon subjected to a high electric field[END_REF].

The stochastic formulation of the drift-diffusion equation enables us to use stochastic numerical methods to simulate the transport of charges within the semiconductor material and thus gathering results over many simulations to extract the figures of interest. In this case the Monte-Carlo Euler-Maruyama method [START_REF] Pavliotis | Stochastic Processes and Applications : Diffusion Processes, the Fokker-Planck and Langevin Equations[END_REF] is particularly well-suited, the algorithm is straightforward and goes as follows. We choose a constant time step dt, and we start from a given position x 0 . We call X n the numerical approximation of X t at time n • dt.

The Euler-Maruyama Monte-Carlo scheme is then given by:

X 0 = x 0 X n+1 = X n + v(X n , n • dt)dt + σ(X n , n • dt)W n (10)
where (W n ) n≥1 are random variables that follow a centered Gaussian distribution with variance dt. The transport properties of the carrier does not depend on the chosen time step dt, as the Brownian motion construction guarantee to have a standard deviation of √ t where t is the simulated time. Yet, because the Euler-Maruyama deterministic part is a numerical scheme of order one in time, the time step must stay relatively small, we found that a value of 10 fs was a good compromise between accuracy and computational time. Our transport algorithm has been validated by checking that the random process results in the desired electron and hole velocities as displayed in figure 1, which demonstrates the accuracy of our method.

B. Impact ionization model

The peculiarity of our model is that the impact ionization modeling is fully coupled with the transport model. We have implemented a modified version of the "Random Path Length Algorithm" (RPLA) [START_REF] Ong | A simple model to determine multiplication and noise in avalanche photodiodes[END_REF] [START_REF] Pilotto | An improved random path length algorithm for p-i-n and staircase avalanche photodiodes[END_REF].

In our model, the path integral [START_REF] Sun | A Simple Analytic Modeling Method for SPAD Timing Jitter Prediction[END_REF] is computed on the fly, concurrently with the transport simulation, at each time step of the Monte-Carlo process. Thus, our model does not require the assumption that the device is divided into an absorption region and an avalanche region as in other previous works [START_REF] Rosset | A model for the jitter of avalanche photodiodes with separate absorption and multiplication regions[END_REF]. It allows us to simulate the full electron-hole avalanche process three-dimension space without the need for severe approximations on the device electrostatics in the so-called avalanche region. Carrier multiplication is modeled through field-dependent rates of impact ionization [START_REF] Sze | Imapct ionization physics of Semiconductor Devices[END_REF], which is the number of impact ionization events a carrier experienced per distance unit. The probability that an electron will not create an impact ionization event between two positions z and z is given by :

P ii (z |z) = 1 -exp - z z α(x) dx ( 11
)
where α(x) is the impact ionization coefficient at position x [START_REF] Mcintyre | A new look at impact ionization-Part I: A theory of gain, noise, breakdown probability, and frequency response[END_REF].

With the RPLA the occurrence of an impact ionization event is modeled by a random variable that follows the probability P ii . A way of simulating such a probability law is to generate a random number variable r uniformly distributed between 0 and 1, then the integral expression ( 11) is computed along the carrier's trajectory, to get the position z for which P ii (z |z) = r. If such a z exists, then we consider that an impact ionization occurred at this location. The integral P ii is then reset to zero for the particle and a new random number r is drawn. The created electron and hole will follow the same process, and might, eventually also experienced impact ionization.

In [START_REF] Rosset | A model for the jitter of avalanche photodiodes with separate absorption and multiplication regions[END_REF], the authors back-traced the timing of the impactionization event by assuming that the carrier's path in a socalled "avalanche-region" is a one-dimensional straight line at constant velocity. In our algorithm, we compute the integral P ii (z |z) on the fly, consistently with the trajectory, and whenever P ii (z |z) reaches or exceeds the random value r, an electron-hole pair is created at the current position of the parent carrier. The capacity of the multiplication model to be used consistently with the transport simulation is mandatory to accurately model the operation of state-of-the-art devices with highly bent electric fields, for which considering straight lines for carrier trajectories is very inaccurate [START_REF] Helleboid | Comprehensive modeling and characterization of photon detection efficiency and jitter tail in advanced SPAD devices[END_REF]. The extraction of impact ionization coefficients in Silicon was performed by several authors either from measurements or from BTE Monte Carlo simulations. In this work we have used the experimental coefficients given by Van Overstraeten and De Man [START_REF] Van Overstraeten | Measurement of the ionization rates in diffused silicon p-n junctions[END_REF], but one can straightforwardly replace them by any other data.

Our multiplication modeling was tested by extracting the impact ionization coefficients back from simulation, at different electric fields. This is done by counting the number of impact ionization events on a given traveled distance and comparing the obtained coefficient with the values of the Van Overstraeten data used as an input, both for electrons and holes. The results are reported in figure 2 and show that our stochastic algorithm perfectly reproduces the input model of multiplication. It was noted by several authors that modeling particle multiplication with impact ionization rates leads to the possibility Impact ionization coefficient (cm

-1 )
Model electron Model hole ADMC electron ADMC hole Fig. 2. Comparison of the carriers' impact ionization rates between input model [START_REF] Van Overstraeten | Measurement of the ionization rates in diffused silicon p-n junctions[END_REF] and the ADMC results. Error bars represent the 99% confidence interval. ADMC coupled with a RPLA accurately reproduces the input impact ionization model.

for a cold electron injected in a high electric field to create an impact ionization immediately after being injected. This is in contradiction with the fact that an impact ionization can occur only if the electron has gained enough energy from the electric field to be able to knock another electron out of the valence band [START_REF] Okuto | Ionization coefficients in semiconductors: A nonlocalized property[END_REF], [START_REF]Threshold energy effect on avalanche breakdown voltage in semiconductor junctions[END_REF]. To accumulate this minimal threshold energy, the electron is required to travel a minimal mandatory distance in a high electric field. This length is called "dead-space" as no ionization can occur within this distance. Those considerations led to the creation of different models that take into account this energy threshold, among which the dead space model [START_REF] Hayat | Effect of dead space on gain and noise of double-carrier-multiplication avalanche photodiodes[END_REF] and the effective field model [START_REF] Mcintyre | A new look at impact ionization-Part I: A theory of gain, noise, breakdown probability, and frequency response[END_REF] are the most used. The minimal energy required for a carrier to trigger an impact ionization was computed using theoretical considerations on the band structure and conservation of mass and momentum [START_REF] Hauser | Threshold Energy for Avalanche Multiplication in Semiconductors[END_REF], [START_REF] Anderson | Threshold Energies for Electron-Hole Pair Production by Impact Ionization in Semiconductors[END_REF]. It is also common to use effective values for this threshold energy, computed through conventional Monte-Carlo simulations and fitting procedures [START_REF] Spinelli | Dead space approximation for impact ionization in silicon[END_REF]. Those threshold energy values are often used as a fitting parameter, within an acceptable range of values. In this study we have used the following values for Silicon, the threshold energy for an electron to impact ionize is E e th = 1.8eV and E h th = 2.4eV for a hole, following [START_REF] Tan | Avalanche noise measurement in thin Si p+-i-n+ diodes[END_REF]. For a particle injected with no energy in the system, we compute its energy gain from the electric field by extracting the electrostatic potential it traveled through, which is equivalent to integrating the electric field F over the particle path:

E gained = P ath --→ F(z) dz = n k=0 ---→ F(X k ) v(X k ) • dt (12)
While E gained is lower than the threshold energy E e th the impact ionization coefficient α in (11) is zero, and no impact ionization can occur. This non-local impact ionization model is called the "hard dead space model" because of the sharp transition between zero and non-zero impact ionization probability. The need to avoid such a sharp transition, that is not observed in BTE Monte Carlo simulation results, has been discussed by several authors and is still under debate [START_REF] Sano | Soft and hard ionization thresholds in Si and GaAs[END_REF]. One might also want to take into account the residual energy of the carrier that engender the impact ionization, as well as the initial energy of the created electron and hole [START_REF] Dolgos | Full-band Monte Carlo simulation of single photon avalanche diodes[END_REF]. We did not focus on those topics as they would not significantly affect BrP or jitter result in SPADs. Yet we show that within our model, it is straightforward to implement models that would cope with those small-scale physical effects. Indeed, by randomly selecting the initial energy of the carriers on a given statistical distribution, one can retrieve a so-called "soft dead-space" [39] model. As an example, we have used a Gaussian distribution whose mean and standard deviation were set to roughly match results from BTE Monte Carlo. When using a dead space impact ionization model, one must tweak the commonly used function for impact ionization rates for electrons α e and holes α h . Indeed, those values, as they come from measurements, intrinsically include the presence of the dead space. Following [START_REF] Spinelli | Mean gain of avalanche photodiodes in a dead space model[END_REF], the impact ionization coefficients for dead space models are transformed as such:

α dead-space e,h = 1 1 α e,h + d e,h (13)
where d e and d h are the dead space for electrons and holes, respectively.

The dead space can be extracted by simulating multiplication at different fields and looking at the distance traveled by an electron (d e ) or a hole (d h ) to reach its minimal threshold energy. We have checked that the obtained values matches the theoretical values for a constant electric field E:

d e,h = E e,h th Eq . (14)
The extracted values are shown in figure 3, one can observe the good agreement in both cases with the theoretical value for a constant electric field. The difference in behavior of the impact ionization models is underlined by the probability density function of the distance traveled by a carrier before it suffered its first impact ionization, represented in figure 4. The local model enables impact ionization to occur immediately after the carrier was injected, while dead-space models give results where a certain distance is needed for the carrier to gain a threshold energy, necessary to trigger an impact ionization. The hard dead-space model results in a zone where the probability of triggering an impact ionization is strictly zero, while the soft model slightly relax this condition, with a low, but positive ionization probability.

C. Device simulation

Device simulation procedure begins with the design of the desired device followed by simulation of the manufacturing process with a commercial process simulator to obtain the final doping profile of the device. The electrostatics of the device is then calculated with standard commercial TCAD simulation software [START_REF] Selberherr | Analysis and Simulation of Semiconductor Devices[END_REF]. This preliminary workflow results in doping and electric field profiles in two or three dimensions. Our methodology then uses these profiles whenever it is necessary for the simulation, by interpolating the doping and the electric field within the mesh of the simulation. In what follows, the frozen field approximation is used [START_REF] Higman | Coupled Monte Carlodrift diffusion analysis of hot-electron effects in MOSFETs[END_REF], which means that the electrostatics of the device is not updated at each iteration of the Monte Carlo simulation, with respect to the new charge distribution. Our simulation could be extended to a so-called self-consistent device simulation, at the cost of computing, with a solver of the Poisson equation, the charge densities and its corresponding electric field after every given time step [START_REF] Hess | Monte Carlo Device Simulation[END_REF]. In Algorithm 1, the overall design of our methodology is represented, non-local approximations are not considered for clarity purposes. The algorithms with hard and soft deadspace models are shown in the Algorithms section of the supplementary material file. Figure 5 represents the trajectory of a single electron injected in the volume of the SPAD, which drifts and diffuses following the electric field, and finally

time ← 0 X 0 ← x 0 p ii ← 0 Generate r → U(0, 1)
while time ≤ T max and no EndCondition do for each particle do Generate W n → N (0, dt)

X n+1 = X n + v(X n , n • dt)dt + σ(X n , n • dt)W n p ii = p ii + v(X n , n • dt) • dt • α e (X n ) if p ii ≥ r then
Impact Ionization Occurs : Inject a new electron and hole at position X n+1 reset p ii = 0 reset r → U(0, 1) if particle hits an Ohmic contact then Remove particle from the simulation

III. RESULTS AND DISCUSSION

A. Device architecture and simulation parameters 1) Device architecture: Our method was tested against characterization results from several SPAD devices. The considered SPADs are made of a main N on top of P regions that form the main junction, also hereafter called the avalanche region.

A guard ring formed of lower p-type and n-type annular areas is also used to prevent lateral breakdown [START_REF] Wang | High-Performance Structure of Guard Ring in Avalanche Diode for Single Photon Detection[END_REF]. These SPADs are backside illuminated. The optical stack ends with an antireflective coating layer, the interface between the optical stack and the silicon is structured to enhance the light absorption [START_REF] Liu | Advances in near-infrared avalanche diode single-photon detectors[END_REF]. This interface is passivated by the means of an additional P+ layer, to reduce the dark count rate due to the interfacial defects [START_REF] Van Sieleghem | A Near-Infrared Enhanced Silicon Single-Photon Avalanche Diode With a Spherically Uniform Electric Field Peak[END_REF]. A schematic representation of a two-dimensional sectional view of the SPAD architecture is given in figure 6.

The dimensions of the different regions as well as the doping level in each of them can be modulated in order to optimize the device performances.

2) Simulation parameters: The simulations were performed with a time step of 10 fs and a maximum simulation time of 5 ns. Simulations are typically made by injecting electrons in the volume of the SPAD, with starting points distributed on a grid of 50 × 50 points × 150 points, the values are taken to approximately match the SPAD's dimension ratios. On each of these points, 250 simulations are performed, leading to a total of around 10 8 simulations. During the simulation, if more than 200 particles are active (i.e. not absorbed by a contact), we consider that the simulation ends in an avalanche breakdown event. All the simulations are independent from each other, which allows us to parallelize the simulation process. One may note that, unlike in standard "free-flight"-based Monte Carlo, the time step in our method is constant and equal for every particle. This implies that all the particles are always synchronized in time, which is a strong advantage for parallelization over the particles, without performing costly synchronization. Simulations on three-dimensional SPADs, using the above set of parameters, performed with parallelization over the simulations, using 40 cores, take around 15 hours. Although comparison with standard BTE Monte Carlo highly depends on implementation, we found that our method is at least 200 times faster than typical Full-Band Monte Carlo simulation. The temperature is set to 300 K, to match the experimental conditions, for comparison purposes.

The only considered wavelength is 940 nm, which is the wavelength of the light emitted by the laser diode used in the experiments. At such large wavelengths, the absorption coefficient of silicon is very low, and the optical absorption can thus be considered uniform in the silicon volume thanks to the optical structuration that scatters the incoming light in all directions and the metal and oxide reflectors that reflect the light back into the silicon volume. For shorter wavelengths, the optical absorption might become non-uniform and one would have to couple the optical simulation with the ADMC simulation. This can easily be done by setting the number of simulations starting from a given position to a number proportional to the optical absorption at this position.

B. Avalanche breakdown probability

The detection efficiency of a SPAD is determined by, on one hand, its optical absorption (OA), and on the other hand its avalanche breakdown probability (BrP), that is, the probability that a carrier injected within the device will eventually trigger an avalanche that will be detected by the readout circuit. The probability for an incoming photon to trigger a count is called the Photon Detection Efficiency (PDE). One has also to take into account the fact that only a fraction of the sensor surface is photosensitive, this fraction is called the optical fill factor (FF). The optical absorption as much as the optical fill factor are ruled by the optical stack and the back-end circuit. Hence, only the avalanche breakdown probability can be optimized at a Silicon level. Within our model, the breakdown probability is defined as the ratio of the number of simulations that turned into an avalanche over the total number of simulations performed. For a given simulation, if a certain threshold number of active particles is reached, then we consider that the avalanche mode is attained. This threshold number of particles can either be defined as the number for which a high enough current will be measured, and a count triggered, either an arbitrary threshold can be used, set to a number after which it is highly unlikely that the process will stop by itself. Figure 7 represents the number of active particles versus time for 20 000 simulations in a single SPAD.

With this figure, it appears that if a simulation reaches 1000 active particles, it is then very unlikely that it will not lead to an avalanche (one event over 20 000).

In modern SPAD architecture the electrostatics is neither homogeneous nor isotropic [START_REF] Morimoto | 3.2 Megapixel 3D-Stacked Charge Focusing SPAD for Low-Light Imaging and Depth Sensing[END_REF], thus the avalanche breakdown probability has the same properties. In order to determine the local BrP, we performed several simulations with the same starting point r 0 , so we can determine the BrP of this location. By repeating this operation on many points within the device simulation domain, we create a map of avalanche breakdown probability. Such a map is shown in figure 8 this location to eventually trigger an avalanche. These kinds of maps were already produced by means of deterministic methods, such as solving the McIntyre equations over electric field lines [START_REF] Helleboid | Comprehensive modeling and characterization of photon detection efficiency and jitter tail in advanced SPAD devices[END_REF]. Those methods could not take into account the statistical intrinsic behavior of the particles. With our method, a particle can, for example, overcome a potential barrier thanks to its thermal energy. One can appreciate that, even though the breakdown probability maps are overall similar (see figure 8) with both methodologies, the one made with the ADMC method shows that some zones that seems to have a zero-BrP with the deterministic methods, have actually a low, but positive probability to be active.

On figure 9, many simulations from two single starting points are shown. For an electron injected near the top right anode (white star), either it enters directly in the peripheral junction, where the electric field is low, and cannot trigger any avalanche, or it diffuses downward enough to get around the ring region, to eventually be collected in the main junction and trigger an avalanche. This second pathway cannot be appreciated with straight line or field lines method. Considering the random nature of transport is mandatory to fully capture the detailed breakdown mechanisms and understand the device operation.

Being able to accurately predict low breakdown probability areas is as well mandatory to simulate dark count rate precisely. Indeed, small areas can have a massive contribution to dark count due to localized clusters of defects with high Shockley-Read-Hall (SRH) generation rate [START_REF] Jay | Clusters of Defects as a Possible Origin of Random Telegraph Signal in Imager Devices: a DFT based Study[END_REF]. This can be particularly the case in the near contact region (top left of figure 6) that are processed by means of heavily doped implant, that can lead to clusters of extending defects likely to enhance SHR generation rate. Dark count rate can be computed by integrating the result of the multiplication of the local breakdown probability by the local SRH generation rate [START_REF] Sicre | Dark Count Rate in Single-Photon Avalanche Diodes: Characterization and Modeling study[END_REF]. Photon detection efficiency can also be computed by integrating the result of the multiplication of the local breakdown probability by the local photon absorption rate, that can be obtained from optical simulation.

C. Jitter

The main application of SPAD devices is ToF for threedimensional imaging. In order to accurately determine the distance to an object, the time for a photon between its emission and its reception after being reflected by the object must be measured precisely. This measure is perturbed by different mechanisms: the limited time resolution of the light source, the transport time from the absorption location within the SPAD volume and the avalanche region, the time for the avalanche to grow up to measurable level of induced current and finally, the finite precision of the readout circuit.

The statistical distribution of the detection times for a given distance is called jitter. Among the different sources of jitter, we focus on those directly induced by the silicon device : the transport time and the avalanche build-up time.

The statistical aspect of the jitter is key, as a high discrepancy in the timing will be interpreted as a wrong distance, for example, an additional 1 ns leads to an error of 0.15 m in the distance hand. Figure 10 illustrates the jitter due to the avalanche build-up. The distribution of the time interval between the first impact ionization and several numbers of active particles thresholds is shown. The spread of the distribution indicates that the time required to have a detectable current is a significant source of jitter. Building those statistical distributions requires an important number of results, which are hardly attainable with standard Monte Carlo device simulation. On figure 10, each curve gathers 10 6 simulations. While recent works explored analytically the statistics of the avalanche build-up under constant one-dimensional electric fields [START_REF] Windischhofer | The statistics of electron-hole avalanches[END_REF] [50], our method enables further investigations on three-dimensional arbitrary fields. Figure 11 exhibits some properties of the build-up jitter. The average plot highlights the exponential growth of the avalanche process, with a very stable rate after 1000 particles. The plot of the standard deviation emphasizes the fact that an important proportion of the build-up time fluctuation occurs during the beginning of the avalanche. Indeed after around 1000 particles, the standard deviation reaches a plateau. As a matter of fact, when the multiplication process is made of many particles, its timing properties are averaged, and the random history of a single particle has a lesser importance in the overall avalanche process.

In order to verify the accuracy of our model, we compared the jitter distribution obtained from our simulation with jitter characterization performed on several SPADs, manufactured within a CMOS technology process. Figure 12 gathers three comparisons between simulation and measurement of the jitter statistical distribution. We picked three diodes with significantly different electrostatics, to underline the robustness of our methodology against different architectures. Experimental and simulation parameters were kept identical for the three The exponential growth of the avalanche process is clearly visible. After around 1000 particles, the standard deviation reaches a plateau, this indicates that the build-up time standard deviation is mostly determined by the timing of the first impact ionization events. devices. All three SPADs have the same epi-layer thickness, SPAD A and B have a lighter P-SPAD implant, which helps increase the depletion region. SPAD C has a heavier P implant, leading to a narrower depletion region, which causes more diffusion of the electrons in the SPAD volume. This is reflected in the jitter distribution, as the SPAD C has a higher jitter tail than the other two SPADs. The results show that for all devices, the simulation well reproduces the experimental results, despite very different distribution shapes. The beginning of the peak cannot be fully recovered because of the two extrinsic causes of jitter: the finite timing resolution of the laser pulse and of the readout circuit. To assess the performance of a given device, it is often useful to reduce the full jitter distribution into a single number figure. We use a very common figure that is the time for which a given proportion of the population has a smaller jitter. For example, the 50% population time is the median of the timing jitter data. We compared the 50% population time and the 90% population time for 16 different layouts. The correlation plots between values extracted from measurements and simulations are represented in figure 13. For both metrics, we found a very significant correlation between the simulation and the characterization. This corroborates the fact that the part of the jitter our model cannot fully assess is the same for all diodes, which is consistent with its origin (laser pulse and readout circuit), which is also the same among all silicon layouts.

IV. CONCLUSION

Since SPAD operation from carrier transport to multiplication process is intrinsically stochastic, an accurate modeling of such a device must include this randomness behavior. We introduced a new kind of Monte Carlo simulation that couples transport and multiplication mechanisms to achieve full multiparticle simulation. Our methodology, by its flexibility, makes it possible to carefully tune the transport as well as the impact ionization models. The simulations are utterly faster than any standard BTE Monte Carlo solver, which allows us to simulate the statistical behavior of the device with large amounts of data within reasonable computational time. The new methodology was tested against previous methods and characterization results to ensure its accuracy.

Fig. 1 .

 1 Fig.1. Comparison of the carriers' drift velocity between input model[START_REF] Arora | Electron and hole mobilities in silicon as a function of concentration and temperature[END_REF],[START_REF] Canali | Electron and hole drift velocity measurements in silicon and their empirical relation to electric field and temperature[END_REF] and the Monte-Carlo results. Error bars represent the 99% confidence interval. The velocity is extracted by computing the distance over time ratio for 1000 electrons at different electric fields. The results show that our algorithm reproduces perfectly the input mobility model at any electric field.

1 /

 1 Electric field (cm • V -1 )

Fig. 3 .

 3 Fig.3. Extracted dead space at various electric fields for hard and soft threshold models. Error bars represent the 99% confidence intervals. Plain blue and red lines represent theoretical values from formula (14) for electron and hole, respectively. Both models give values close to the theoretical calculation with an expected higher variability for the soft dead space model.

Fig. 4 .

 4 Fig.[START_REF] Bronzi | SPAD Figures of Merit for Photon-Counting, Photon-Timing, and Imaging Applications: A Review[END_REF]. Probability density function of the position of the first impact ionization for a cold electron injected in a constant electric field of 600kV/cm with three different impact ionization models. Each curve is made out of 1 000 000 single-particle simulation results.

Fig. 5 .

 5 Fig.5. Illustration of the device multi-particle simulation. The background map represents the strength of the electric field, the red region corresponds to the SPAD main junction. An electron is injected in the volume of the SPAD (yellow cross), it drifts and diffuses towards the junction and triggers several impact ionization events, leading to a self-sustained avalanche. The color of the trajectories' dot represents the intrinsic time of the particles. A zoomed view of the avalanche is shown in the inset. The third dimension was dismissed for visualization purposes.

Fig. 6 .

 6 Fig. 6. Schematic sectional view of the considered SPAD architecture.

Fig. 7 .

 7 Fig. 7. Number of particles versus time for 20 000 simulations of a SPAD device. The simulations are synchronized via the time of the first impact ionization event. Particles are created through the impact ionization process and removed when collected at Ohmic contacts. The log scale allows us to appreciate the exponential growth of the number of particles when the avalanche mode is reached.

Fig. 8 .

 8 Fig. 8. Two-dimensional map of breakdown probability obtained by the method of Helleboid and al. using the McIntyre model over electric field lines [28] (left). Two-dimensional map of breakdown probability obtained from the ADMC method (right).

Fig. 9 .

 9 Fig. 9. Many electron trajectories before their first impact ionization. Only two starting points are considered, represented by the star and the cross. The important diffusion leads to very different paths, thus the particles enter the avalanche region at very diverse locations.

Fig. 11 .

 11 Fig.11. Mean (A) and standard deviation (B) time to reach different numbers of particles in the avalanche process. The exponential growth of the avalanche process is clearly visible. After around 1000 particles, the standard deviation reaches a plateau, this indicates that the build-up time standard deviation is mostly determined by the timing of the first impact ionization events.

Fig. 12 .Fig. 13 .

 1213 Fig.[START_REF] Gao | 3D Electro-optical Simulations for Improving the Photon Detection Probability of SPAD Implemented in FD-SOI CMOS Technology[END_REF]. Probability density function of the jitter for three different SPAD layouts. The densities extracted from simulations are compared to characterization results. The comparison shows an excellent agreement between simulation and measurement. All simulations were done in three-dimensional meshes, at the same excess bias as the characterization and each curve gathers around one million simulations, from starting points evenly distributed within the device volume.

  Fig.10. Distribution of the time between the first impact ionization event and different number of particles thresholds. This figure illustrates the timing jitter due to the avalanche process itself. There is a not negligible spread of time for the avalanche to be established at measurable levels of currents.

					Number of particles reached
	Probability density	function (arb. unit)					10 100 1000 10000 50000
		0.00	0.01	0.02	0.03	0.04	0.05	0.06
					time (ns)

This work was supported by the French "Association Nationale Recherche Technologie" ANRT and the French GeSPAD Project of the Agence Nationale de la Recherche (ANR-20-CE24-0004).


Françoise Lecaignard
email: francoise.lecaignard@inserm.fr

Jérémie Mattout
email: jeremie.mattout@inserm.fr

Forward models for EEG/MEG

published or not. The documents may come

Synopsys

Like in most inverse problems, where causes have to be inferred from consequences, estimating brain activity from scalp electrophysiological data requires a model. Many EEG/MEG biophysical models have been developed over the last few decades, mostly since the inception of whole head MEG and high density EEG systems. These so-called forward models map neuronal activity onto sensor signals. In this chapter, we briefly cover the main variants and components of these models. They mostly differ by the way they represent neuronal populations, and the way they account for individual anatomy and physical properties. They range from simple (analytical) models to more realistic (numerical) ones. While there is still room for improvement, the question of whether more complex forward models are useful can now be addressed empirically using Bayesian model comparison.

This chapter is made of three sections. The first section describes our current knowledge of what most contributes to EEG/MEG signals and how laws of physics enable us to quantify this contribution. The second section covers the main steps and assumptions that yield the computation of EEG/MEG forward models. It also stresses the main factors of uncertainty in this computation. Finally, the third section emphasizes the difference between the complexity of forward model computation and (probabilistic) model complexity. This notion becomes particularly important when one wants to evaluate different forward models of EEG/MEG data.

I.

The origin of EEG and MEG signals a. From microscopic (single neuron) activity to macroscopic (population level) currents

Neurons have the property to be electrically excitable, to produce and to propagate nerve impulses. Neuronal activity expresses as transmembrane electric currents, which contribute to the macroscopic electric potential differences and magnetic fields measurable on scalp and its vicinity. Pyramidal and stellate cells are the two main types of neurons. They mostly differ in shape and orientation relative to the cortical surface. Whatever their functional specialization, their typical structure consists of a cell body (soma) surrounded by dendritic branches and extending in a nerve fiber, the axon, which projects onto other neurons through synaptic junctions (Fig. 1). The neuron's membrane contains pores that are selectively permeable to ions. In the absence of excitation, continuous ionic exchanges across the pores maintain the membrane potential V m at a resting value:

V m = V r = V int ra -V extra ≈ -70mV .
Typical dendrites receive inputs from pre-synaptic neurons by way of a chemical exchange called neurotransmission. Fixation of released neurotransmitters on post-synaptic receptors triggers the opening of specific pores, allowing for ionic movements across the membrane. These currents cause the post-synaptic (membrane) potentials (PSP). Synaptic integration of PSP over time and space results in an action potential (AP) propagating along the axon towards the apical dendrites, upon condition that a depolarization threshold has been reached (

mV V V t m 50 - ≈ >
). Excitatory synapses tend to evoke a local membrane depolarization ( V m > V t ), whereas inhibitory synapses cause local membrane hyperpolarization ( V m < V t ).

From the extracellular medium, a site along the membrane where ions move into the cell corresponds to a sink current (current disappearance), whereas a site where ions go out of the cell defines a source current (current appearance). Transmembrane and intracellular currents are grouped together at the micro-(neuronal) scale and more generally at the meso-(cell assembly) scale to form the primary currents p J , often described as currents related to postsynaptic activity. They produce an electric field in the medium, which in turn gives rise to conduction currents E J e σ = , also referred to as volume or secondary currents. In every point within the head, the electric current can be described as

E J J p σ + = .
Source and sink current distributions on the membrane can be modeled by multipolar (Taylor) developments [START_REF] Mosher | MEG Source Imaging Using Multipolar Expansions[END_REF]. Simple dipole models capture well the contribution of PSPs [START_REF] Koles | Trends in EEG source localization[END_REF], whereas APs are well described by quadrupole models whose contribution decreases rapidly with distance and can be neglected in the context of EEG and MEG [START_REF] Crouzeix | Méthodes de localisation des générateurs de l'activité électrique cérébrale à partir de signaux électro-et magnétoencéphalographiques[END_REF]) (Fig. 2). Furthermore, approximately 50,000 neurons with simultaneous activity and similar orientation are required to be picked up by MEG and EEG sensors [START_REF] Pernier | Electro et magnéto encéphalographie Biophysique, techniques et méthodes[END_REF]. Consequently, it is usually assumed that only the primary currents produced by a synchronous population of neurons with a preferred orientation such as the pyramidal cells contribute to the electric potentials and magnetic fields measurable on scalp [START_REF] Crouzeix | Méthodes de localisation des générateurs de l'activité électrique cérébrale à partir de signaux électro-et magnétoencéphalographiques[END_REF][START_REF] Nunez | Electric fields of the brain: The neurophysics of EEG[END_REF]. b. From macroscopic currents to observed electric and magnetic fields Predicting the electric ( E ) and magnetic ( B ) fields produced by neuronal activity on EEG and MEG sensors requires solving Maxwell's equations in head tissues. In their general form, Maxwell's equations write:

∇.E  = ρ ε t B E ∂ ∂ - = ∧ ∇ ! 0 . = ∇ B ∇ ∧ B  = µ J   +ε ∂  E ∂t $ % & ' ( ) (1.1)
where J and ρ indicate the volume density of current and volume density of charge, respectively. ε and µ are the electric permittivity and magnetic permeability of the medium, respectively.

The head is composed of several tissues with various conductivities and therefore can be described as a finite inhomogeneous conducting volume. It is commonly accepted that ε and µ are equal to free space value (ε 0 and µ 0 respectively) (Malmivuo and Plonsey, 1995).

Given the frequency spectra of EEG and MEG signals is much below 1 kHz, capacitive effects can be neglected [START_REF] Schwan | Capacitive properties of body tissues[END_REF]. The duration of the electromagnetic wave propagation (from neuronal sources to sensors) is negligible compared to physiological time constants, hence scalp measures appear as instantaneous and synchronous, and the quasistatic approximation holds [START_REF] Plonsey | Considerations of quasi-stationarity in electrophysiological systems[END_REF]. Under this regime, Maxwell's equations simplify as follows:

0 . ε ρ = ∇ E 0 = ∧ ∇ E 0 . = ∇ B J B 0 µ = ∧ ∇ (1.2)
The electric and magnetic fields are now decoupled and importantly, the forward computation has become independent of time. The later means that only the location, orientation and amplitude of the neuronal sources need to be known to compute the sensor signals.

Since the electric potential simply relates to the electric field by

 E r, t ( ) = -grad    (V (r, t)) (1.3)
and accounting for Ohm's law, the law of conservation of charge and Biot-Savart's law which relates magnetic fields to the underlying currents, equations (1.2) yield the following formulations for the electric potential V(r,t) and the magnetic field B(r,t) distributions, given a time-varying source distribution Jp(r,t) and tissue conductivity σ(r):

) , ( . ) , ( ( . ) ( t r J t r V grad r p ∇ = ∇ σ (1.4) [ ] ∫ - - ∧ - = ' ' ) ' ( )) , ' ( ( ) ' ( ) , ' ( 4 ) , ( 3 0 dr r r r r t r V grad r t r Jp t r B σ π µ (1.5)
at spatial location r and time t. Difference between EEG and MEG data generated by the same underlying sources are illustrated in Fig. 3.

Fig.3. a). Current dipole with electric field lines (red) and magnetic field lines (green). Because of different conductivities in different head tissues, current lines are distorted when they cross tissue boundaries (figure adapted from [2]).b) Scalp topography of real EEG data (auditory evoked response N100). c) Corresponding scalp topography of MEG data (from simultaneously recorded EEG and MEG responses to auditory tones).

II.

The main forward model assumptions

A forward model of EEG and MEG data calls for a description of the neuronal current distribution (referred to as the source distribution), a specification of the conducting properties of the head tissues and some information about the sensors (locations and orientations in the case of MEG). A large number of models have been proposed in the literature. They differ by the realism with which they account for the physical and geometrical properties of head tissues. In other words, the practicality of deriving EEG and MEG forward models often boils down to trading off between the complexity of computing a fine individual structural and physical model, and the accuracy of the ensuing forward prediction. In addition, contrary to simple models, complex ones will not have analytical solutions and numerical approximations will be needed.

a. The source model

At the macro-scale, the coherent activity of a neuronal assembly is most commonly modeled by an equivalent current dipole (ECD) (Fig. 3a). Higher order models have also been explored [START_REF] Nolte | Current multipole expansion to estimate lateral extent of neuronal activity: a theoretical analysis[END_REF], mostly as an attempt to better capture the extent of active cortical areas. However, such models introduce a higher number of unknown parameters that are also difficult to interpret. Nevertheless, quadrupoles may be of interest, particularly in MEG, to model the complex fields generated by extended and deep cortical sources [START_REF] Jerbi | Localization of realistic cortical activity in MEG using current multipoles[END_REF].

Two main approaches can be distinguished when considering the ECD as a unitary source model:

•

The dipolar approach relies on the assumption that a fairly small number of ECDs (fixed a priori, typically less than 8) contribute simultaneously and significantly to the scalp data. Each ECD is fully described by 6 parameters, 3 for its location and 3 for its orientation and magnitude. This yields a well-conditioned system where the number of unknown parameters to be fitted is smaller than the number of independent data points (which at each point in time roughly corresponds to the number of sensors).

•

The distributed or imaging approach relaxes the strong constraint on the number of active regions. A few thousands ECDs are typically used to model the entire source space with fixed positions, either distributed over a 3D regular grid covering the whole brain volume or limited to the cortical sheet with a possible normal orientation constraint. For the latter, only each ECD magnitude is left to be estimated, yielding a linear but ill-conditioned system to be solved. A unique solution to the ill-posed EEG/MEG inverse problem will be obtained here by incorporating additional constraints about the source configuration (see Chapter 331). This source model enables one to produce images of cortical or brain activity, which makes statistical inference at the subject or group level more tractable and more sensitive [START_REF] Mattout | Canonical source reconstruction for MEG[END_REF][START_REF] Litvak | Electromagnetic source reconstruction for group studies[END_REF].

b. The head model

This is another critical part of the forward model. It embeds our knowledge and assumptions on head tissue geometry and conducting properties.

• Head geometry models Spherical models are the simplest ones, consisting of concentric spheres with homogeneous and isotropic conductivity in each compartment or layer [START_REF] Rush | Current distribution in the brain from surface electrodes[END_REF][START_REF] De Munck | Mathematical dipoles are adequate to describe realistic generators of human brain activity[END_REF]. Although the head is not spherical, these models are attractive because of the ensuing exact analytic expressions for the electric potential and magnetic field on the head surface. They have been extensively evaluated empirically.

The "3 shell" model distinguishes between the scalp, skull and brain layers with radius ratio and isotropic conductivities as proposed in [START_REF] Rush | EEG electrode sensitivity--an application of reciprocity[END_REF]. It is largely used for EEG and available in most routine softwares. Since the magnetic permeability is homogenous over tissues and volume currents barely contribute to the external magnetic field, spherical models appear more suitable for MEG than EEG. MEG spherical models are blind to volume currents contributions but present the advantage that neither tissue conductivity knowledge nor the radius of the sphere is needed.

The "overlapping spheres" model, which refines the single sphere model by fitting a sphere to each sensor location, provides a better MEG forward solution [START_REF] Huang | A sensor-weighted overlappingsphere head model and exhaustive head model comparison for MEG[END_REF].

Realistic models are numerical models that have been mostly developed for EEG, to better account for the shape and electrical properties of the tissues, namely their electrical conductivity as with the Boundery Element Method (BEM) [START_REF] Hämäläinen | Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data[END_REF] and their anisotropy as with the Finite Element Method (FEM) [START_REF] Marino | A finite difference model of electric field propagation in the human head : implementation and validation[END_REF] or the Finite Difference Method (FDM) [START_REF] Lemieux | Calculation of electrical potentials on the surface of a realistic head model by finite differences[END_REF]. The BEM relies on surface meshes derived from MRI segmentation and assigns each layer with homogenous and isotropic conductivity. Studies comparing spherical and BEM models in EEG and MEG obtained better source estimates using BEM for dipoles below the supratemporal plane [START_REF] Yvert | A systematic evaluation of the spherical model accuracy in EEG dipole localization[END_REF][START_REF] Crouzeix | An evaluation of dipole reconstruction accuracy with spherical and realistic head models in MEG[END_REF].FEM and FDM rely on 3D meshes where each finite element can be ascribed with a different, anisotropic conductivity tensor matrix. Models then differ in their number of compartments, conductivity values, anisotropic ratios and tensor orientations. The latter can be derived from Diffusion Tensor Imaging (DTI) [START_REF] Tuch | Conductivity tensor mapping of the human brain using diffusion tensor MRI[END_REF].

We refer the reader to [START_REF] Rush | EEG electrode sensitivity--an application of reciprocity[END_REF][START_REF] Sarvas | Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem[END_REF][START_REF] Meijs | On the numerical accuracy of the boundary element method[END_REF][START_REF] Hämäläinen | Magnetoencephalography -theory, instrumentation, and applications to noninvasive studies of the working human brain[END_REF][START_REF] Mosher | EEG and MEG: forward solutions for inverse methods[END_REF] for a detailed description of the forward model computation, under the spherical assumption (analytic form) and the more realistic assumption (numerical form).

• Head tissue conductivities

Individual conductivity values are of high importance, particularly for realistic models. Estimations from dead tissues have proved very different from in vivo values and first in vivo observations have come from anaesthetized animals [START_REF] Robillard | Specific-impedance measurements of brain tissues[END_REF]. More recently, in vivo measures have been made possible in humans thanks to the advent of Electric Impedance Tomography (EIT) [START_REF] Ferree | Regional head tissue conductivity estimation for improved EEG analysis[END_REF][START_REF] Gonçalves | The application of electrical impedance tomography to reduce systematic errors in the EEG inverse problem--a simulation study[END_REF] as well as Diffusion Tensor Imaging (DTI) [START_REF] Tuch | Conductivity mapping of biological tissue using diffusion MRI[END_REF] but these techniques can still hardly be used routinely. Hence empirical values reported in experimental studies are largely used as a first approximation [START_REF] Rush | Current distribution in the brain from surface electrodes[END_REF]. Several studies based on simulated data have compared the relative sensitivity of EEG and MEG forward solutions to conductivities. Critically, EEG is highly sensitive to the brain/skull conductivity ratio [START_REF] Vallaghé | A global sensitivity analysis of three-and four-layer EEG conductivity models[END_REF] as well as to white matter anisotropy [START_REF] Güllmar | Influence of anisotropic electrical conductivity in white matter tissue on the EEG/MEG forward and inverse solution. A high-resolution whole head simulation study[END_REF]. MEG is particularly sensitive to brain tissue conductivity [START_REF] Gencer | Sensitivity of EEG and MEG measurements to tissue conductivity[END_REF][START_REF] Van Uitert | Influence of Brain Conductivity on Magnetoencephalographic Simulations in Realistic Head Models[END_REF] and white matter anisotropy [START_REF] Güllmar | Influence of anisotropic electrical conductivity in white matter tissue on the EEG/MEG forward and inverse solution. A high-resolution whole head simulation study[END_REF].

c. Sensor registration

Sensor description relative to the head model is achieved by means of a spatial transformation based on head landmarks (least-square fitting) or head surface (surfacematching methods) or both, identified in both the MRI and the electrophysiological coordinate systems.

Various sources of errors are associated with sensor coregistration, particularly landmark identification on MR images, electrode and landmark digitization and head movements during MEG acquisition. Typically, coregistration errors range between 5 and 10 mm [START_REF] Whalen | Validation of a method for coregistering scalp recording locations with 3D structural MR images[END_REF][START_REF] Hillebrand | Practical constraints on estimation of source extent with MEG beamformers[END_REF], with moderate consequences on EEG inverse solutions [START_REF] Wang | The influence of electrode location errors on EEG dipole source localization with a realistic head model[END_REF][START_REF] Acar | Effects of Forward Model Errors on EEG Source Localization[END_REF] but potentially dramatic effects on MEG ones [START_REF] Hillebrand | The use of anatomical constraints with MEG beamformers[END_REF]. Interestingly, uncertainty about the forward model, due to coregistration, could be accounted for in the source reconstruction process thanks to probabilistic or Bayesian methods [START_REF] López | A general Bayesian treatment for MEG source reconstruction incorporating lead field uncertainty[END_REF].

III.

Empirical evaluation of forward model assumptions a. From forward to generative models

The forward relationship between source parameters 𝜃 and observed EEG or MEG data 𝑌 is of the general form

𝑌 = 𝐿(𝜃) (3.1)
where 𝐿 indicates the lead-field operator and embodies all the pre-cited anatomical and biophysical assumptions one needs to account for in the forward model. Data 𝑌 is a 𝑁 * 𝑇 matrix, where 𝑁 is the number of sensors and 𝑇 the number of time samples. 𝜃 is a 𝑃-long vector made of all source location, orientation and amplitude parameters.

Forward models have in themselves barely no interest; they are only useful and even mandatory when one aims at reconstructing brain activity from scalp recordings, that is inverting equation (3.1) to estimate 𝜃. Contrary to the forward computation of 𝐿, this inverse problem is ill-posed and requires additional (prior) information or constraints to ensure a unique solution (see chapter 331). Two types of additional assumptions can be specified: assumptions about measurement noise and prior knowledge about parameters 𝜃 . In particular, when source locations are fixed as in distributed approaches, equation (3.1) becomes linear:

𝑌 = 𝐿. 𝜃 + 𝜀 (3.2)
where 𝐿 is the 𝑁 * 𝑃 lead-field (or gain) matrix operator and 𝜀 models an additive measurement noise which is usually assumed to follow a Gaussian distribution with zero mean and a fully known or parameterized variance structure [START_REF] Mattout | MEG source localization under multiple constraints: an extended Bayesian framework[END_REF].

This highlights the fact that solving the inverse problem requires the specification of not only the lead-field operator 𝐿 but also of the prior distributions over noise 𝜀 and parameters 𝜃.

Altogether, those assumptions make a full generative model, which could be used to simulate realistic EEG or MEG data.

The probabilistic or Bayesian framework is very much appropriate to define and invert generative models. Indeed, probabilistic distributions can flexibly describe our knowledge or uncertainty about a phenomenon. Moreover, advanced inference techniques have been developed to invert complex probabilistic models (see chapter 325).

b. Bayesian (forward) model comparison

Importantly, the Bayesian framework enables formal model comparison given empirical observations. Since the forward model is part of the generative model, Bayesian model comparison offers a principled way to compare forward model assumptions, as long as all other assumptions (namely priors over noise and source parameters) are kept the same for each compared generative model, and provided that all models are fitted to the same dataset. This is a recent and important extension to previous evaluation approaches of EEG and MEG forward models which mostly rested on numerical simulations [START_REF] Crouzeix | An evaluation of dipole reconstruction accuracy with spherical and realistic head models in MEG[END_REF][START_REF] Vatta | Realistic and Spherical Head Modeling for EEG Forward Problem Solution: A Comparative Cortex-Based Analysis[END_REF][START_REF] Acar | Effects of Forward Model Errors on EEG Source Localization[END_REF] and on a few empirical measures using biophysical phantoms [START_REF] Leahy | A study of dipole localization accuracy for MEG and EEG using a human skull phantom[END_REF][START_REF] Baillet | Evaluation of inverse methods and head models for EEG source localization using a human skull phantom[END_REF].

Bayesian model comparison rests on computing the model evidence 𝑝 𝑌 𝑀 (see chapter 328). The higher the model evidence, the better the model. A useful approximation to the logevidence is the Free energy ( 𝐹 ) [START_REF] Penny | Comparing Dynamic Causal Models using AIC, BIC and Free Energy[END_REF]. It can be obtained using variational techniques (see Chapter 327) and has already been used to compare forward models of EEG and MEG data. Namely, it could show that canonical cortical meshes may carry sufficient structural information to solve the MEG inverse problem [START_REF] Mattout | Canonical source reconstruction for MEG[END_REF][START_REF] Henson | Selecting forward models for MEG source-reconstruction using model-evidence[END_REF].

c. A note on model complexity

Given data 𝑌 and model 𝑀 with parameters 𝜃, the free energy writes

𝐹 = ln 𝑝 𝑌 𝜃, 𝑀 ! -𝐾𝐿 𝑞 𝜃 𝑝 𝜃 𝑀 (3.3)
where • 𝑝 𝑌 𝜃, 𝑀 and 𝑝 𝜃 𝑀 are the likelihood and prior distributions, respectively. They fully define the generative model 𝑀. • 𝑞 𝜃 is the approximate posterior distribution over model parameters (the outcome of the inverse inference process). • KL is the Kullback-Leibler divergence which can be interpreted as a statistical distance between two distributions. Here, it quantifies the distance between the posterior and the prior distribution over 𝜃.

Importantly, the first term in equation (3.3) corresponds to model accuracy, while the second term quantifies model complexity. In the general case of Gaussian distributions, this term writes

𝐾𝐿 𝑞 𝜃 𝑝 𝜃 𝑀 = 1 2 ln 𝐶 𝜃 - 1 2 ln 𝐶 𝜃 𝑌 + 1 2 𝜇 𝜃 𝑌 -𝜇 𝜃 𝑇 𝐶 𝜃 -1 𝜇 𝜃 𝑌 -𝜇 𝜃 + 𝑇𝑟𝑎𝑐𝑒 𝐶 𝜃 -1 𝐶 𝜃 𝑌 + 𝑐𝑠𝑡 (3.4)
where 𝜇 ! , 𝐶 ! and 𝜇 ! ! , 𝐶 ! ! are the mean and variance of the prior and posterior distributions, respectively.

Given those equations, changing the lead-field operator by moving from a simple spherical head model to a more realistic one might increase the free energy in two ways:

by improving the fit of the data (increasing model accuracy); -by reducing model complexity through a posterior distribution that would decrease the above Kullback-Leibler divergence. Namely, this could be the case if the realistic model would yield a smaller posterior correlation between parameters [START_REF] Penny | Comparing Dynamic Causal Models using AIC, BIC and Free Energy[END_REF].

Importantly, this means that a more realistic model, although more complex in a computing sense (because it requires the fine extraction of individual anatomical and biophysical features) might yield a significantly higher free energy.

However, this will be the case only if such a model offers a more realistic and higher spatial resolution that the data can accommodate. In other words, whether it is worth deriving a fine and realistic head model for source reconstruction depends on the spatial precision that the data can offer.

In [START_REF] Henson | Selecting forward models for MEG source-reconstruction using model-evidence[END_REF], using data from a face perception MEG experiment; it was shown that a BEM model should be preferred to a spherical one, provided individually-defined inner skull and scalp meshes were used.

Finally, beside head models, Bayesian model comparison can also be used to evaluate the ability of EEG and MEG data to inform advanced source models based on neural masses. Recent developments of Dynamic Causal Models to study brain effective connectivity have led to more biologically plausible models of neuronal populations. As an example, a recent study suggests that EEG data can distinguish between the dynamics of local neuronal excitatory and inhibitory subpopulations [START_REF] Moran | Free energy, precision and learning: the role of cholinergic neuromodulation[END_REF].

Related chapters

Fig. 1 .

 1 Fig.1. Schematic view of the neuron's structure.

Fig. 2 .

 2 Fig.2. Schematic view of neuronal current distributions. An excitatory synapse creates a local sink current, which induces primary currents across the membrane and within the cell. The conducting currents close the loop in the extracellular space, hence creating a distribution of source currents along the membrane. In the same way, an inhibitory synapse creates a source current.


S Taileb

E Rougon

V Robin

V Rodriguez

P Vidal

S.-M Lau-Chapdelaine

J Melguizo-Gavilanes

A Chinnayya

A three-step, three-gamma model for the numerical modeling of the critical height of the propagation of semi-confined detonation waves

Introduction

Gaseous detonation waves are supersonic combustion waves, propagation of which depends of its surroundings [START_REF] Lee | The detonation phenomenon[END_REF]. When the gaseous mixture is confined by an inert layer, lateral expansion of the products and streamline divergence are responsible of the velocity deficit, as compared to the ideal Chapman-Jouguet velocity, and ultimately can quench the detonation wave. The determination of this critical height is of importance, for propulsion applications, in rotating detonation engines [START_REF] Shepherd | Analytical models for the thrust of a rotating detonation engine[END_REF] and for safety [START_REF] Rudy | Critical conditions of hydrogen-air detonation in partially confined geometry[END_REF] in industrial and nuclear plants. This canonical configuration has also been heavily studied in the context of explosives (see review of [START_REF] Short | High explosive detonation-confiner interactions[END_REF]). Reynaud et al. [START_REF] Reynaud | A computational study of the interaction of gaseous detonations with a compressible layer[END_REF][START_REF] Reynaud | Computation of the mean hydrodynamic structure of gaseous detonations with losses[END_REF], Mi et al. [START_REF] Mi | Effect of spatial inhomogeneities on detonation propagation with yielding confinement[END_REF], Houim and Fievisohn [START_REF] Houim | The influence of acoustic impedance on gaseous layered detonations bounded by an inert gas[END_REF] performed numerical simulations using a one-step, one-γ model. The numerical critical heights from Reynaud were not consistent with the experimental results. Taileb et al. [START_REF] Taileb | Influence of the chemical modeling on the quenching limits of gaseous detonation waves confined by an inert layer[END_REF] and Shigeoka et al. [START_REF] Shigeoka | The Numerical Investigation of Hydrogen Detonation Propagating in Semi-confined Layers, ICDERS 2019[END_REF] found that the use of detailed chemistry enabled to recover the experimental trends, at least for hydrogen-oxygen at ambient conditions. However, the computational cost required to the use of detailed chemistry remains quite prohibitive. Calibration processes have been proposed by Kaplan and al. [START_REF] Lu | A chemical-diffusive model for simulating detonative combustion with constrained detonation cell sizes[END_REF], where fewer calibration coefficients (isentropic coefficient and heat of reaction) have been used in a one-step global reaction. Houim and Taylor [START_REF] Houim | Detonation initiation from shock and material interface interactions in hydrogen-air mixtures[END_REF] have shown that they could reproduce the shock to detonation transition. Veiga et al. [START_REF] Lopez | A methodology to develop simplified kinetic schemes for detonation simulation[END_REF] have tried to match the D -κ curve. Nevetheless, chemical conversion of reactants does not always take place through a global step, but rather through a sequence of intermediate steps. The three-step model [START_REF] Dold | Comparison between shock initiations of detonation using thermally-sensitive and chain-branching chemical models[END_REF][START_REF] Sharpe | Homogeneous explosion and shock initiation for a three-step chain-branching reaction model[END_REF][START_REF] Melguizo-Gavilanes | A comparison between constant volume induction times and results from spatially resolved simulation of ignition behind reflected shocks: implications for shock tube experiments[END_REF] enjoys some of these properties, through a sequence consisting of initiation, branching and termination steps. Some of the initial constraints of this model, which were probably initially present for analytic traceability can be loosened, namely different molecular weights and heat capacities. This numerical study reports some of these calibration efforts used to recover the critical height from a three-step, three-γ model, the reference being the one obtained from the detailed chemistry.

Mathematical model and numerical methods

The 2-D reactive Euler equations are used for this study. Full details of the high-fidelity numerical schemes and the computational implementation can be found in Reynaud et al. [START_REF] Reynaud | A computational study of the interaction of gaseous detonations with a compressible layer[END_REF][START_REF] Reynaud | Computation of the mean hydrodynamic structure of gaseous detonations with losses[END_REF], Taileb et al. [START_REF] Taileb | Influence of the chemical modeling on the quenching limits of gaseous detonation waves confined by an inert layer[END_REF]. A three-step chain branching model is used to model the kinetics of 2 H 2 + O 2 at atmospheric pressure and ambient temperature. The chemical decomposition is represented by:

Initiation F -→ Y (1)
Branching

F + Y -→ 2Y (2)
Termination Y -→ νP (3)
where ν = W Y /W P is the stoichiometric coefficient, which is the ratio of the molecular weights of the intermediate and the products species. The enthalpy reads h = Y F h F + Y Y h Y + Y P h P such as each specific enthalpy is assumed to be a linear function of the temperature

h k = h 0 k + c p,k T .
The mixture capacity heat is given by c p,m = c p,k and the molecular weight is

1/W = (Y k /W k ).
Only the third step releases heat. From the ZND model [START_REF] Higgins | Steady one-dimensional detonations[END_REF], the Master equation for the velocity can be obtained

(1 -M 2 ) du dx = W W F - c p,F c p,m dY F dt + W W Y - c p,Y c p,m dY Y dt + W W P - c p,P c p,m dY P dt (4)
+ q c p,m T dY P dt ,
the right-hand side being the thermicity. The terms related to variation of the heat capacities and molar masses are represented in magenta, which are not present in the original three-step model. The consumption rates for the different species are

r I = k I k I (ρY F ) exp(-E I /R u T ) r B = k B k B (ρY F )(ρY Y )/W F exp(-E B /R u T ) (5) r T = k T k C (ρY Y )
The corresponding pre-exponential factors for the initiation, branching and termination are

k I = k C exp(E I /R u T I ), k B = k C exp(E B /R u T B ) and k T = k C ,
respectively. E I /R u and E B /R u are the initiation and branching activation temperatures, and T I and T B are the crossover temperatures at which the chain-initiation and branching rates become equal to the chain-termination rate. k I , k T , and k B are calibration coefficients, in green in Eq. 5.

The parameters are: [START_REF] Taileb | Influence of the chemical modeling on the quenching limits of gaseous detonation waves confined by an inert layer[END_REF] after a calibration against constant volume combustion induction times). From Thompson [START_REF] Thompson | Compressible fluid dynamics[END_REF], the CJ detonation velocity can be found to be

k C = 2 × 10 7 s -1 , E I /R u = 25 000 K, E B /R u = 8500 K, T I = 8500 K, T B = 1350 K (from
M CJ = H + (γ P -1)(γ F + γ P ) 2γ F (γ F -1) 1/2 + H + (γ P -γ F )(γ P + 1) 2γ F (γ F -1) 1/2 , H = (γ P -1)(γ P + 1)q 2γ F R F T 0 (6)
3 Numerical results

Results from the ZND model

The ZND reference profiles have been obtained from the SD-Toolbox [START_REF] Browne | [END_REF]. The heat capacity and the molecular weight are related through the Meyer relation R u /W k = γ k -1 γ k c p,k . The molecular weights were determined at first, W F = 12 g/mol, W P = 14.4 g/mol at the von Neumann (vN) and CJ states, respectively. The latter is different from that of the water vapor, due to dissociation at the high temperatures near stochiometry. This enabled to obtain temperatures that were close to the initial and end states of the ZND profile. Then, the isentropic coefficients were determined. Low values are indeed known to enhance cell irregularity through cell bifurcation [START_REF] Sow | The effect of the polytropic index γ on the structure of gaseous detonations[END_REF][START_REF] Taileb | The influence of the equation of state on the cellular structure of gaseous detonations[END_REF] and potential hot spots. The isentropic coefficient for F was determined in order to match as much as possible the temperature at the vN state. Then, γ P was taken from the results of the SD-Toolbox. Then the heat of reaction is deduced from Eq. 6. All the thermodynamic parameters for the intermediate species Y were taken equal to that of the species. Another approach is to calculate the heat of reaction q from the equilibrium detonation velocity D CJ after supplementing Eq. 6 with the exact expression Eq. ( 72) for the equilibrium γ P derived in [START_REF] Vidal | A Velocity-Entropy Invariance theorem for the Chapman-Jouguet detonation[END_REF].

In From Fig. 1, the end state for the temperatures is almost the same for both configurations. However, the main difference was that thermicity was much better captured in configuration 2 with a rather slight modification of the pre-exponential factors. This also impacted the rate at which the other variables reached the end states. This also means that the χ parameter which measures the relative time scales of induction and heat deposition, is also better recovered. We recall χ = T vN τ ind • ∂τ ind ∂T • ∆ ind ∆reac , ∆ reac = u CJ / σmax . Again, this parameter has been used as a measure of the detonation cell regularity [START_REF] Radulescu | The propagation and failure mechanism of gaseous detonations: experiments in porous-walled tubes[END_REF][START_REF] Ng | Detonation instability[END_REF].

Results from 2-D simulations

From the original parameters [START_REF] Taileb | Influence of the chemical modeling on the quenching limits of gaseous detonation waves confined by an inert layer[END_REF] (all molecular weights equal to 12 g/mol, and constant isentropic coefficient equal to 1.33), first simulations have shown that the critical height for the transmission and propagation of the detonation was about 22 mm, that is much greater from the experimental value of 6 mm. First computations, not shown here, have highlighted that the critical height for propagation was different from that of transmission and propagation. Progressive decrease of the reactive height has shown that the critical height for the detonation propagation was less than 15 mm. This is reminiscent from the work of Murray and Lee [START_REF] Murray | The influence of physical boundaries on gaseous detonation waves[END_REF], that showed that the rate of expansion imposed on the detonation wave was important for the determination of the critical height. Figure 2 shows numerical soot foil for the transmission and extinction of the detonation wave in a 15 mm height reactive mixture, for the initial parameters. The red dotted line delineates the reactive mixture from the inert one. Hence, the reactive height is taken to be 15 mm in the following study. The resolution was 10 pts per ∆ ind ≈ 34 µm. As for the initial configuration, L x × y = 10 mm × 25 mm and for the two other configurations, L x × L y = 6 mm × 26 mm. We use 24 processors in the x direction and 30 in the y direction, with a typical mesh size of 13.5 million computational cells.

Discussions

The initial three-step model has been at first modified in order to take into account thermodynamic changes of the products (molecular weight and γ). Then, a second modification has been made to the values of the pre-exponential factors in order to recover the thermicity profile of the ideal ZND model from detailed chemistry. To some degree, matching the thermicity profile, the calibration has also relied on targeting the γ and the χ parameters, which have been used to qualify the degree of irregularity of the detonation cellular structure. Preliminary 2-D results have shown that these two successive modifications enabled to come closer to the reference critical height, determined from detailed chemistry. The results above are promising, but further computations should assess to which extent this modeling approach can reproduce the critical height determined from detailed chemistry. A more general goal should be to widen its thermodynamic and spatial representativeness.

Figure 1 :

 1 Figure 1: Comparaison between the thermicities (top left), temperatures (top right) and polytropic coefficient (bottom). Reference in black, configuration 1 in red, and configuration 2 in green.

Figure 2 :

 2 Figure 2: Numerical soot foil for the transmission and extinction of detonation wave in a 15 mm height reactive mixture. The total length is 65.24 mm.

Fig. 3

 3 Fig. 3 shows the three different configurations. The left figure shows the extinction of the detonation with the original parameters. The middle figure shows a re-initiation process with the calibration of the thermodynamic parameters. The right figure shows that a steady propagation detonation with the calibration of the thermodynamic and kinetics calibrations.

Figure 3 :

 3 Figure 3: Comparaison between the 2-D simulations with three different chemical model parameters. Left: original [9], 56.87 mm-65.24 mm. Middle: thermodynamic parameters calibrated (configuration 1), 132.583 mm-137.625 mm. Right: thermodymamics and kinetics calibrated (configuration 2), 90.324 mm-95.356 mm. The reactive height is 15 mm.

th ICDERS -June 19 -24, 2022 -Napoli, Italy

Acknowledgements

The computations were performed using HPC resources from GENCI-CINES (Grant 2020-A0092B07735). This work was supported by the CPER FEDER Project of Région Nouvelle Aquitaine. Part of this research was supported by the French National Research Agency (ANR) under the project ANR-21-CE05-0002-01.