domain: spi phys

text:
Nicolas Dellinger
email: nicolas.dellinger@onera.fr

Nicolas Bertier

Francis Dupoirieux

Guillaume Legros

Hybrid Eulerian-Lagrangian method for soot modelling applied to ethylene-air premixed flames

Keywords: CFD, Combustion, PAH, Sectional model, Soot, Lagrangian description

Soot particles dynamics is described through a detailed Lagrangian tracking method.

• Soot particles growth and oxidation processes are taken into account.

• Soot particles interactions are limited with success to coalescence as a first guess.

• Particle size distributions are extracted efficiently through population reduction.

• Targeted flames for validation are laminar premixed ethylene-air flames.

Introduction

Soot particles created in gas turbine combustors contribute to local air pollution, particularly in airport areas, and enhance the radiative transfer in the combustion chamber, increasing the gas temperature and the thermal load on the chamber walls and turbine blades. The production of other pollutants, such as nitrogen oxides which are particularly sensitive to temperature, might be influenced too [START_REF] Dorey | Modélisation des phénomènes couplés combustion-formation des suies-transferts radiatifs dans les chambres de combustion de turbine à gaz[END_REF]. As a consequence of both concerns, efficient and accurate modelling of soot formation and evolution is an important issue in the design of aircraft engines, but remains an open field due to the complexity of the phenomena involved. The chemical kinetics of gaseous soot precursors formation from the fuel, as well as of soot particles formation from the precursors, has in particular to be described. The interactions of the soot particles with the surrounding gas, leading to their growth or oxidation, as well as the interactions between the particles must also be considered. Soot models are generally classified in three categories. Empirical models use correlations from experimental measurements to compute the soot volume fraction and number density, for instance as a function of the fuel concentration [START_REF] Moss | Modelling Soot Formation for Turbulent Flame Prediction[END_REF], the concentration of an intermediate species [START_REF] Magnussen | On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion[END_REF], or a location in the flame [START_REF] Legros | Phenomenological model of soot production inside a non-buoyant laminar diffusion flame[END_REF]. Suitable with global reaction schemes, and thus allowing the calculation of complex configurations, these models are however weakly predictive over wide ranges of parameters.

Semi-empirical models, as the Leung one [START_REF] Leung | A simplified reaction mechanism for soot formation in nonpremixed flames[END_REF], are based on a simplified description of soot nucleation, coagulation, growth, and oxidation kinetics. Source terms for the soot volume fraction and number density are calculated from kinetic constants, determined from experimental measurements, or derived from a population balance equation for coagulation. More predictive, these models generally require the concentration of acetylene, as light soot precursor, computed by means of expensive reaction mechanisms or cheaper tabulated chemistry [START_REF] Dorey | Soot and radiation modeling in laminar ethylene flames with tabulated detailed chemistry[END_REF]. Application of semi-empirical models to the simulation of industrial burners is therefore usually based on the tabulation approach.

Detailed models, of particular interest for the present study, include a complex description of both soot chemical kinetics and particles dynamics. Reactions for surface growth by addition of acetylene and condensation of polycyclic aromatic hydrocarbons (PAHs), as well as oxidation by O 2 and OH are generally included in the reaction mechanism [START_REF] Kennedy | Models of soot formation and oxidation[END_REF]. Following the requirement to correctly predict these species and accurately describe soot nucleation, a large number of intermediate species, from benzene up to large PAHs, has to be added to the kinetic scheme, resulting in mechanisms of tens of species and hundreds of reactions [START_REF] Wang | A detailed kinetic modeling study of aromatics formation in laminar premixed acetylene and ethylene flames[END_REF]. Particles dynamics, which includes soot particles evolution via the aforementionned chemical processes associated with coalescence and aggregation, can be modelled using sectional methods [START_REF] D'anna | A model of particulate and species formation applied to laminar, nonpremixed flames for three aliphatic-hydrocarbon fuels[END_REF][START_REF] Saggese | Kinetic modeling of particle size distribution of soot in a premixed burner-stabilized stagnation ethylene flame[END_REF][START_REF] Rodrigues | Coupling an LES approach and a soot sectional model for the study of sooting turbulent non-premixed flames[END_REF] or methods of moments [START_REF] Frenklach | Detailed Mechanism and Modeling of Soot Particle Formation[END_REF][START_REF] Mueller | Hybrid Method of Moments for modeling soot formation and growth[END_REF]. Unlike empirical and semi-empirical models, the two approaches account for both the polydispersion in size and the non-sphericity of the particles, which can be described as fractal aggregates [START_REF] Saggese | Kinetic modeling of particle size distribution of soot in a premixed burner-stabilized stagnation ethylene flame[END_REF][START_REF] Mueller | A joint volume-surface model of soot aggregation with the method of moments[END_REF][START_REF] Wen | Study of soot growth in a plug flow reactor using a moving sectional model[END_REF][START_REF] Park | A novel fixed-sectional model for the formation and growth of aerosol agglomerates[END_REF][START_REF] Eberle | Large Eddy Simulations of a Sooting Lifted Turbulent Jet-Flame[END_REF]. Methods of moments, for which a restricted set of moments can be sufficient to characterize the soot particles population, are yet preferable to limit the CPU cost in industrial applications, sectionals methods requiring a high number of sections and, eventually, variables per section.

Because of the large chemical mechanism required to predict the evolution of the soot precursors and the need to properly describe the particles dynamics, detailed soot models can be prohibitely expensive for the computation of aircraft engines combustion chambers. Tabulated chemistry [START_REF] Gicquel | Laminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ILDM with differential diffusion[END_REF][START_REF] Van Oijen | Modeling of complex premixed burner systems by using flamelet-generated manifolds[END_REF][START_REF] Fiorina | A filtered tabulated chemistry model for LES of premixed combustion[END_REF]] might be a solution to limit the number of transported variables relative to the fuel combustion. However, the modelling of heat losses as well as the consideration of combustion regimes and multiphase flow effects increase the complexity and thus the computational cost of such approaches. Also, the extension of tabulated chemistry to species having long evolution time scales, such as NO or soot, is an issue [START_REF] Boucher | A method to extend flamelet manifolds for prediction of NOx and long time scale species with tabulated chemistry[END_REF]. Alternatively, sectional methods can be used to characterize PAH growth [START_REF] Saggese | Kinetic modeling of particle size distribution of soot in a premixed burner-stabilized stagnation ethylene flame[END_REF][START_REF] Zamuner | Numerical Simulation of Soot Formation in a Turbulent Flame with a Monte-Carlo PDF Approach and Detailed Chemistry[END_REF][START_REF] Richter | Detailed modeling of PAH and soot formation in a laminar premixed benzene/oxygen/argon low-pressure flame[END_REF] and combined with a reduced chemical mechanism for lighter species [START_REF] Di Domenico | Development and validation of a new soot formation model for gas turbine combustor simulations[END_REF][START_REF] Blacha | Soot predictions in premixed and non-premixed laminar flames using a sectional approach for PAHs and soot[END_REF]. The present study follows the latter strategy and focuses on the Lagrangian description of soot particles as a dispersed solid phase. The Lagrangian tracking method is a substitute to sectional methods and methods of moments for computing particle dynamics. It gives the possibility to follow a numerical particle continuously along its trajectory, and to monitor its individual interaction with both the gas phase and the other particles. No classes or sections of particles have here to be prescribed, and the calculation of the particle size distribution, of first interest for soot applications, is straightforward.

Nevertheless, the computational cost could be a limiting factor, the statistical convergence of a Lagrangian method being inversely proportional to the square root of the number of particles.

The Lagrangian tracking approach has been recently proposed by Gallen et al. [START_REF]Lagrangian tracking of soot particles in LES of gas turbines[END_REF] and Ong et al. [START_REF] Ong | Evaluation of a Lagrangian Soot Tracking Method for the prediction of primary soot particle size under engine-like conditions[END_REF]. While acetylene is used in these two studies, the soot particles are nucleated from PAH molecules and radicals in the present one. More detailed soot evolution models, which allow the two-way coupling of the gaseous and solid phases, are additionally used. An innovative algorithm to control the number of tracked particles is finally introduced in order to keep the computational cost moderate.

The model was implemented in the aerothermochemistry CEDRE code, developed at ONERA [START_REF] Refloch | Cedre software[END_REF]. The approach being unusual compared to more common Eulerian methods, the subsequent simulations must be validated by comparison with academic flame configurations. While numerical simulations of steady twodimensional laminar sooting flames have already been conducted with some success [START_REF] Consalvi | Numerical study of soot formation in laminar coflow diffusion flames of methane doped with primary reference fuels[END_REF], the discrimination between failures of the soot model and of the flow description that may be responsible for the discrepancies between experimental measurements and numerical fields especially remains an issue. For that reason, the present study is specifically focused on the validation based on steady one-dimensional laminar ethylene-air flames at three fuel equivalence ratios.

Gas phase models and methods

Through the coupling of dedicated solvers, the CEDRE code allows the simulation of multiphase compressible reacting flows on unstructured meshes. In particular, the conservation equations for the reacting gas phase in the present study are resolved through the finite-volume Eulerian solver CHARME, able to deal with future targeted applications to industrial configurations.

Conservation equations

For compressible laminar reacting flows, the mass, momentum, and energy conservation equations write:

∂ρY k ∂t + ∇ • (ρY k V) = -∇ • J k + ωk + S d k , (1) Mass
∂ρV ∂t + ∇ • (ρV ⊗ V) = ∇ • (τ -pI) + ρg + S d q , (2)
Momentum

∂ρe t ∂t + ∇ • (ρe t V) = -∇ • J et + ∇ • [(τ -pI) V] + ρg • V + S d et + S exp et (3) Energy
where t is the time, ρ the density of the ideal gas mixture, composed of N s species, Y k the mass fraction of species k, p the pressure, V the velocity vector, g the gravity vector, and I the identity matrix.

Mass conservation equation (1)

J k is the molecular diffusion flux of species k, calculated using Hirschfelder and Curtiss' approximation [START_REF] Hirschfelder | Molecular theory of gases and liquids[END_REF] as:

J k = ρV k Y k = -ρ D k M k M ∇X k + V c Y k , (4)
where V k , M k and X k are respectively the diffusion velocity, the molecular weight and the mole fraction of species k, and M the molecular weight of the mixture. D k is the diffusion coefficient of species k in the mixture, computed by equation ( 5) as a function of the binary diffusion coefficients D jk of each species j in species k.

D k = (1 -Y k )   j =k X j D jk   -1 (5)
V c is a correction velocity, given by equation ( 6), ensuring global mass conservation.

V c = - Ns k=1 D k M k M ∇X k (6)
ωk is the chemical source term of species k, calculated from the reaction mechanism [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. S d k , the feedback source term of the dispersed phase -i.e. the soot particles -on species k, will be discussed in Subsection 4.3.

Momentum conservation equation (2)

τ is the viscous tensor, calculated for Newtonian fluids considered in the present study as:

τ = µ ∇V + (∇V) T + κ - 2 3 µ (∇ • V) I, ( 7
)
where µ is the dynamic viscosity of the mixture, defined as the arithmetic mean of the species viscosities µ k weighted by the species mass fractions, and κ the bulk viscosity, commonly assumed to be negligible.

S d q , the feedback source term of the dispersed phase, will be discussed in Section 3.

Energy conservation equation (3)

The energy conservation equation is written using the total energy, defined by equation ( 8), where T is the temperature, T 0 the reference temperature, h 0 f k the specific enthalpy of formation of species k at T 0 , c v the specific heat capacity at constant volume of the mixture and R the gas constant:

e t = Ns k=1 Y k h 0 f k + T T0 c v dT - RT 0 M + 1 2 V • V. (8)
J et is the heat diffusion flux, computed under the assumption of negligible Dufour effect as:

J et = -λ∇T + Ns k=1 h 0 f k J k , ( 9
)
where λ is the thermal conductivity of the mixture. λ/µ is defined as the arithmetic mean of the ratios λ k /µ k weighted by the species mass fractions, with λ k the thermal conductivity of species k.

S d et , the feedback source term of the dispersed phase, will be discussed in Section 3. Radiative heat transfers are not taken into account in equation (3), because, in the validation cases hereinafter, experimental temperature profiles are imposed through the corrective source term S exp et [START_REF] Dorey | Soot and radiation modeling in laminar ethylene flames with tabulated detailed chemistry[END_REF][START_REF] Xu | Soot formation in laminar premixed ethylene/air flames at atmospheric pressure[END_REF]. It allows to make the calculation of the targeted burner-stabilised flames easier.

Numerical methods

A second order MUSCL (Monotonic Upwind Scheme for Conservation Laws) type scheme is used for space discretisation and combined with a HLLC (Harten-Lax-van Leer-Contact) scheme for hyperbolic fluxes. A second order Gear implicit scheme associated with a GMRES method for the resolution of the linear system is applied for time integration.

Lagrangian dispersed phase solver

Soot particles dynamics is described in the present study through the newly developed Lagrangian module SPARTACCUS of the CEDRE code.

Preliminary hypothesis

Soot particles, whose diameter is in the range 1-80 nm hereinafter, are supposed to behave as tracers embedded in the gas flow. Their velocity is therefore taken equal to the gas velocity and their temperature, assumed uniform, to the gas temperature. Soot particles are furthermore modelled as solid spheres, with a constant density ρ s .

Two collision regimes should be considered, coalescence and aggregation, depending on the age of the particles. Young soot particles, considered as liquid-like spheres [START_REF] Dobbins | Soot inception temperature and the carbonization rate of precursor particles[END_REF], are indeed deformable and will preferably coalesce with their collision partner to form a larger spherical particle. Older carbonised particles, on the contrary, are the elementary solid bricks for the building of soot aggregates. The description of the latter is however not of primary interest for the comparison of numerical simulations and experimental measurements in the present study. Also computationally1 challenging, it will have to be addressed separately in the future. Aggregation is here replaced by "non-coalescence" to bound coalescence within its physical domain.

To this end, the ageing of the soot particles is assumed to depend on a single parameter, their diameter, and a critical diameter d crit is introduced. Coalescence only occurs if one of the collision partners has a diameter smaller than d crit . If not, "non-coalescence" takes place and the kinematic parameters of the collisions partners are equalised by analogy to aggregation.

The present formulation allows to directly compare the particle diameter to the primary particle diameter measured by Xu et al. [START_REF] Xu | Soot formation in laminar premixed ethylene/air flames at atmospheric pressure[END_REF] as part of the validation process in Section 5. Comparison of the soot volume fraction -the quantity usually available for complex configurations -is also straightforward, since not modified through aggregation at constant volume, at least as a first-order approximation. For the same reason, the soot radiative heat transfer might additionally be modelled as a function of the volume fraction [START_REF] Dorey | Soot and radiation modeling in laminar ethylene flames with tabulated detailed chemistry[END_REF]. Note that, in a flame with a mean temperature of 1500 K as in the following, the peak wavelength λ involved in the radiative heat transfer would be of the order of 2 µm. This means that the size parameter πd/ λ of a particle, with d its diameter, would be much less than one and implies that light absorption would exceed light scattering. The sensitivity of the radiative heat transfer to the shape of the soot aggregates -which are not considered -would therefore not be of first order.

Numerical methods

The particle number density function, governed by the Williams kinetic equation [START_REF] Williams | Spray Combustion and Atomization[END_REF], is discretised in the Lagrangian dispersed phase solver by a finite number of tracked Dirac distributions or numerical particles, each having a numerical weight w p . A numerical particle p thereby represents w p physical particles with identical state vectors composed of the position x p , the velocity V p , the temperature T p and the radius r p .

For soot applications, the state vector of the numerical particles is updated at each iteration of duration ∆t through a collision and a transport step, detailed below.

Collision step

Collisions are described through a four-step algorithm, presented as part of the Supplementary materials, designed to reproduce the real physics of collisions from the tracking of numerical particles, the number of which being necessarily limited compared to the real number of physical particles. The underlying theory of the algorithm is out of the scope of this study and the interested reader is referred to the original work of Villedieu and Hylkema [START_REF] Villedieu | Une méthode particulaire aléatoire reposant sur une équation cinétique pour la simulation numérique des sprays denses de gouttelettes liquides[END_REF][START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes : application aux propulseurs à poudre[END_REF] for details.

According to the collision regime, the state vectors of particles 1 and 2 are modified as follows:

• For coalescence, the total volume of the particles is conserved. The numerical weight of particle 1 is therefore set to w 1w 2 and the new state vector of particle 2 is computed through the system of equations [START_REF] Saggese | Kinetic modeling of particle size distribution of soot in a premixed burner-stabilized stagnation ethylene flame[END_REF]:

                 x coa 2 = a 1 x 1 + a 2 x 2 , V coa 2 = a 1 V 1 + a 2 V 2 , T coa 2 = a 1 T 1 + a 2 T 2 , r coa 2 = 3V coa 4π 1/3 , (10)
where a 1 and a 2 are the ratios of the volumes V 1 and V 2 of particles 1 and 2, respectively, to the one of the coalesced particle

V coa = V 1 + V 2 .
It means that the result of the coalescence process is a spherical particle containing all the mass of the collision partners.

Note that, if both particles have originally the same numerical weight, their state vectors are taken equal to the one of the coalesced particle, and their numerical weights are halved.

• For "non-coalescence", the numerical weight of particle 1 is set to w 2 , and its position and velocity are taken equal to the one of particle 2. It means that the result of the non-coalescence process is a set of spherical particles located at the same point with the same velocity, as if they were stuck to each other by a single contact point by analogy to aggregation.

Note that, if both particles have initially different numerical weights, a new particle, with a numerical weight equal to w 1w 2 and a state vector identical to the one of particle 1, has to be created first.

Transport step

The transport step describes the evolution, due to interactions with the gas phase, of the particles velocity, position and temperature on the one side, which is straightforward for tracers, and of the particles radius through chemical processes on the other side. At each time step, the new state vector -with superscript n + 1 -of a numerical particle p located in cell c is calculated as a function of the former one -with superscript n -through an explicit scheme as follows:

             V n+1 p = V n c , x n+1 p = x n p + V n c ∆t, T n+1 p = T n c , m n+1 p = m n p + ωn p V c ∆t, (11)
where V c and T c are the gas velocity and temperature in cell c, and m p the mass of the particle. The mass equation replaces here the radius equation for the sake of clarity, the new radius being inferred from equation 12 for spherical particles. The mass source term ω will be discussed in Subsection 4.3.

r n+1 p = 3 3m n+1 p 4πρ s (12)

Optimisation

To reduce the computational cost of the Lagrangian solver, a quadtree-like procedure [START_REF] Berg | Computational Geometry: Algorithms and Applications[END_REF] grouping particles by radius r and residence time t r -defined as the time required by a particle to leave the cell where it is located -is applied every N ts time steps. Its goal is to identify and suppress the numerical particles carrying redundant informations in each cell c where the number of numerical particles N np,c is higher than a constant N np . A dispersion criterion is used to that end in the two-dimensional (r, t r ) -space, of primary interest here since the evolution of a particle in a cell is driven -among other parameters (cf. 4.3) -by its size and its interaction time with both the gas phase and the other particles.

In the reduction procedure described as part of the Supplementary materials, the sum of the particles volumes and the sum of the particles surfaces are conserved at the expense of their number, the reactive surface being of particular interest for the coupling with the gas phase. It is however worth noting that the conservation of all quantities is achieved when the maximum number of particles per cell N np is sufficiently high (namely when σ q tends towards zero in the Supplementary materials). The impact of N np on the particle size distribution will therefore be discussed in Subsection 5.3.1. The particles in a cell are furthermore assumed to have close enough temperatures -even though particles transported from a neighbouring cell might introduce some discrepancies -so that their specific heat capacity is considered identical over the range of particle temperatures. The latter hypothesis allows to compute the temperature of each mean particle without the use of an iterative variable change.

Momentum and energy feedback source terms

For each time step ∆t and grid cell c, S d q,c and S d et,c are calculated as:

S d q,c = - 1 V c ∆t p∈Pc(t) w p t+∆t t d (m p V p ) dt Variation of the particle momentum dt, (13)
S d et,c = - 1 V c ∆t p∈Pc(t)
where P c (t) is the set of particles located in cell c at time t and c p the specific heat capacity of particle p.

Note that, even if the velocity and temperature of a particle are set to the ones of the gas phase at the transport step, their variations are to be known for the computation of the feedback source terms. Hence, both variables are modified during the collision and population reduction steps.

Modelling of soot nucleation and evolution

Gas-phase chemistry

Gas-phase chemistry is described in the CHARME solver using a reduced mechanism for ethylene combustion, successfully applied in previous studies [START_REF] Blacha | Soot predictions in premixed and non-premixed laminar flames using a sectional approach for PAHs and soot[END_REF][START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF], which includes 43 species and 304 reactions.

PAH 0 PAH 1 PAH 2 PAH * 0 PAH * 1 PAH * 2 SOOT 0 0.1 0.2 0.4 0.8 1.6
M (kg/mol)

PAH modelling

Polycyclic aromatic hydrocarbons are modelled in the CHARME solver using the sectional method proposed by Di Domenico et al. [START_REF] Di Domenico | Development and validation of a new soot formation model for gas turbine combustor simulations[END_REF], Blacha et al. [START_REF] Blacha | Soot predictions in premixed and non-premixed laminar flames using a sectional approach for PAHs and soot[END_REF], and recently updated by Eberle et al. [START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF]. Three sections of PAH and PAH radicals (respectively noted PAH i and PAH * i , 0 ≤ i ≤ 2), whose molecular weights range from 0.1 kg/mol to 0.8 kg/mol, are considered (see Fig. 1). A soot section SOOT 0 is added for both the derivation of the reaction mechanism and the modelling of the nucleated soot particles.

The atomic composition, the thermodynamic properties, as well as the viscosity for each section are calculated from reference polycyclic aromatic hydrocarbons [START_REF] Blacha | Soot predictions in premixed and non-premixed laminar flames using a sectional approach for PAHs and soot[END_REF][START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF]. For the sake of simplicity, the binary diffusion coefficients of PAH i and PAH * i are here assumed to be equal to the ones of the smallest PAH (benzene) and of its radical (phenyl), respectively. Their diffusion coefficients are computed through equation ( 5).

PAH i -PAH *

i reaction mechanism The reaction mechanism is composed of five sub-mechanisms [START_REF] Blacha | Soot predictions in premixed and non-premixed laminar flames using a sectional approach for PAHs and soot[END_REF][START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF]: PAH-gas phase interaction, PAH collisions, PAH growth by addition of C 2 H 2 , PAH dehydrogenation and PAH oxidation by O, OH and O 2 . PAH i and PAH * i being minor species, their third-body efficiencies are additionally set to zero in the gas phase reaction mechanism for the present study.

Soot nucleation rate

Soot nucleation is modelled through the collision of PAH 2 and PAH * 2 with the other sections [START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF]:

PAH * i≤2 + PAH * 2 → ν 1 PAH 2 + ν 2 SOOT 0 + ν 3 H 2 , PAH * i≤2 + PAH 2 → ν 1 PAH 2 + ν 2 SOOT 0 + ν 3 H 2 + H, PAH i≤1 + PAH * 2 → ν 1 PAH 2 + ν 2 SOOT 0 + ν 3 H 2 + H, (15)
where the stoichiometric coefficients ν 1 , ν 2 and ν 3 , calculated using the work of Pope et al. [START_REF] Pope | Simultaneous Particle and Molecule Modeling (SPAMM): An Approach for Combining Sectional Aerosol Equations and Elementary Gas-Phase Reactions[END_REF], ensure the conservation of both mass and elements. The reaction rates are computed similarly to the study by Blacha et al. [START_REF] Blacha | Soot predictions in premixed and non-premixed laminar flames using a sectional approach for PAHs and soot[END_REF] and the reaction efficiencies from the work of Eberle et al. [START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF] are used without modification.

The mass production rate of SOOT 0 -the soot nucleation rate ω0 -is computed as for a standard gas species [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF] assuming a molecular weight M 0 equal to 1.2 kg/mol and transmitted to the dispersed phase solver SPARTACCUS at each time step, so that no transport equation is required for the latter species.

Also note that the soot volume fraction, lower than 1 × 10 -6 in the simulation cases hereinafter (see Fig. 5), is far below usual validity limits of two-phase dispersed flow hypotheses [START_REF] Murrone | Numerical modeling of dispersed two-phase flows[END_REF] and therefore neglected in the CHARME solver.

Soot particles modelling

Soot particles are modelled in the SPARTACCUS solver by solid spheres, as mentioned in Section 3.

Their density ρ s is assumed constant -typically between 1700 and 2000 kg/m 3 in the litterature -and chosen equal to 1800 kg/m 3 [START_REF] Zamuner | Numerical Simulation of Soot Formation in a Turbulent Flame with a Monte-Carlo PDF Approach and Detailed Chemistry[END_REF] hereinafter. From the preliminary work of Blacha et al. [START_REF] Blacha | Soot predictions in premixed and non-premixed laminar flames using a sectional approach for PAHs and soot[END_REF], their specific heat capacity is taken equal to the one of acetylene.

Particles nucleation

For each time step ∆t and grid cell c, the mass of soot m c and the number of physical particles N pp,c to be injected are calculated as follows:

m c = ω0,c V c ∆t, (16)
N pp,c = m c N a M 0 , (17)
where N a is the Avogadro number.

If m c has a non-zero value, a single numerical particle is generally injected with a numerical weight equal to N pp,c . Its radius r 0 is given as a function of ρ s and M 0 by [START_REF] Gicquel | Laminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ILDM with differential diffusion[END_REF]:

r 0 = 3M 0 4πN a ρ s 1/3 = 0.64 nm (18)
Its velocity and temperature are initialised to the values in cell c, and its position is determined randomly from the equivalent cell diameter. In the present study, the numerical weight w p of each particle p is nevertheless bounded at injection to ensure the representativeness of the particles population and limit the number of non-significant subjects. Considering an initial mass m c , a numerical particle of maximum numerical weight is therefore injected until m c reaches zero or w p stands below the lower bound. Eventually, the remaining non-injectable mass is stored for a use in the next time step.

r 0 is additionally used as a removal criterion for the non-significant particles at the end of the transport step. Hence, particles are suppressed if their radius becomes lower than 2 -1/3 r 0 , namely if their volume is at least twice as small as the one of the nucleated soot particles.

Chemically-driven evolution

General formulation

The generic reaction of a soot particle SOOT with a species R is first considered:

SOOT + ν R R → SOOT # + νH 2 + ν P P, (19)
where SOOT # is the state of SOOT after reaction with R, and P the associated product. The term νH 2 ensures mass conservation.

The molecular weight of SOOT # and the stoichiometric coefficient ν are determined iteratively at each time step as follows:

M # = n # C M C + H/C M # M H , (20)
and

ν = M -M # + ν R M R -ν P M P 2M H , (21)
where n # C is the number of carbon atoms in SOOT # , calculated from the one in SOOT, R and P, and M the molecular weight of SOOT. H/C is the ratio of hydrogen to carbon atoms in a soot particle, computed in the present study through the correlation introduced by Blacha et al. [START_REF] Blacha | Soot predictions in premixed and non-premixed laminar flames using a sectional approach for PAHs and soot[END_REF] from the work of Richter et al. [START_REF] Richter | Detailed modeling of PAH and soot formation in a laminar premixed benzene/oxygen/argon low-pressure flame[END_REF]:

H/C = 0.4405M -0.10524 , (22)
with M in kg/mol.

For a particle located in cell c at time t, the progress rate Q of reaction ( 19) is given by equation ( 23),

Q = K [R] ν R c 1 N a V c Apparent concentration of SOOT , (23)
with K the reaction rate constant and [R] c the concentration of species R in c.

The evaluation of the reaction rate constants for each evolution model is detailed hereinafter, with [Y] c the set of gas species concentrations in c. The oxidation models are presented for the sake of completeness, however it is worth noting that their contribution to soot particles evolution in the following test cases is not significant.

Surface growth by addition of C 2 H 2

For surface growth by addition of C 2 H 2 , generic reaction [START_REF] Van Oijen | Modeling of complex premixed burner systems by using flamelet-generated manifolds[END_REF] writes:

SOOT + C 2 H 2 → SOOT # + νH 2 . ( 24
)
The reaction rate constant is calculated using the model of Frenklach and Wang [START_REF] Frenklach | Detailed modeling of soot particle nucleation and growth[END_REF][START_REF] Kazakov | Detailed modeling of soot formation in laminar premixed ethylene flames at a pressure of 10 bar[END_REF] as:

K C2H2 = K C2H2,s (T c ) αχ s (T c , [Y] c ) S, (25)
where K C2H2,s is a per-site reaction rate constant, α the fraction of reactive surface sites, χ s the number density of surface sites and S the particle surface.

Oxidation by O 2

For oxidation by O 2 , generic reaction (19) reads:

SOOT + O 2 → SOOT # + νH 2 + 2CO. ( 26
)
The reaction rate constant, also computed using the model of Frenklach and Wang, is given by equation [START_REF] Ong | Evaluation of a Lagrangian Soot Tracking Method for the prediction of primary soot particle size under engine-like conditions[END_REF].

K O2 = K O2,s (T c ) αχ s (T c , [Y] c ) S (27)
Oxidation by OH For oxidation by OH, generic reaction [START_REF] Van Oijen | Modeling of complex premixed burner systems by using flamelet-generated manifolds[END_REF] writes:

SOOT + OH → SOOT # + νH 2 + CO. ( 28
)
The reaction rate constant is calculated using the model of Neoh et al. [START_REF] Neoh | Soot Oxidation in Flames[END_REF] as:

K OH = N a γ OH RT c 2πM OH 0.5 S, (29)
where the collision efficiency γ OH is equal to 0.13.

Condensation of PAH i and PAH *

i By analogy to the nucleation reactions [START_REF] Wen | Study of soot growth in a plug flow reactor using a moving sectional model[END_REF], the reaction product P of a soot particle with a PAH radical in ( 19) is chosen to be H. Therefore, the generic reaction reads for the condensation of PAH * i :

SOOT + PAH * i≤2 → SOOT # + ν * H 2 + H, (30)
and can be written for the condensation of PAH i as:

SOOT + PAH i≤2 → SOOT # + νH 2 . ( 31
)
The reaction rate constant is computed for both as follows:

K PAH ( * ) = k vdw N a γ PAH ( * ) β * i,s (d c,i , S, T c , l c , µ c ) , (32)
where k vdw is the van der Waals enhancement factor, corresponding to the increase of the collision rate due to the attractive van der Waals forces between the collision partners, equal to 2.2 [START_REF] Harris | The Coagulation of Soot Particles with van der Waals Forces[END_REF]. β * i,s is the collision frequency and γ PAH ( * ) the collision efficiency, initially fixed to 0.3 [START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF] and discussed in Subsection 5.4.

The collision frequency between PAH i or PAH *

i and a soot particle is determined using the general Fuchs formulation [START_REF] Fuchs | The mechanics of aerosols[END_REF], including transition from the free molecular regime to the continuum regime with increasing particle diameter. Apart from S, β * i,s is a function of the PAH or PAH * collision diameter d c,i , as well as of the gas temperature T c , mean free path l c and viscosity µ c . d c,i , in Ångström, is calculated for both PAH i and PAH * i through the correlation of Frenklach and Wang [START_REF] Frenklach | Detailed Mechanism and Modeling of Soot Particle Formation[END_REF]:

d c,i = 1.395 √ 2 √ n C i Å , (33)
where n C i is the number of carbon atoms of the species PAH ( * ) i .

Particle mass source term

The mass source term of a particle is calculated as a function of the mass reaction rates associated to each evolution model M as:

ω = M M # M -M Q M . (34)
The mass of a particle being positive, the term ωV c ∆t in equation 11 has to be higher than -m, with m the original mass of the particle. The reaction rate of all evolution models is therefore corrected by a factor defined by equation ( 35),

= m/ | ωV c ∆t| if ωV c ∆t < -m, 1 otherwise, ( 35
)
so that the apparent reaction time associated to the particle for the feedback to the gas phase becomes ∆t.

Feedback source terms

For a given evolution model M and a given particle p, the mass source term of species R, H 2 and P is expressed as follows:

  ωR ωH2 ωP   p,M = w p p Q p,M   -M R ν R M H2 ν p M P ν P   M , (36)
and the feedback source term S d k,c of the dispersed phase on the gas phase in cell c as:

S d k,c = M p∈Pc(t) ωk,p,M , k = R, H 2 , P. (37)
S d k,c , along with the aforementioned source terms, allows a complete two-way coupling of the gas-phase and dispersed phase solvers.

Collision parameters

While soot particles are assumed to behave as tracers at the macroscopic scale, their trajectory at the mesoscopic one is driven by Brownian motion, namely by their interaction with the surrounding molecules of the gas phase. Hence, the collision frequency β 12 (cf. Supplementary materials) of two soot particles is calculated as:

β 12 = k vdw β * 12 , (38)
where β * 12 is the collision frequency function of Brownian motion introduced in the work of Fuchs [START_REF] Fuchs | The mechanics of aerosols[END_REF].

d crit is also fixed to 13 nm as a first guess, one order of magnitude higher than the diameter of nucleated particles.

Validation and discussion

Experimental test cases

The ethylene-air flames simulated in the present study are burner-stabilised laminar premixed flames at atmospheric pressure investigated by Xu et al. [START_REF] Xu | Soot formation in laminar premixed ethylene/air flames at atmospheric pressure[END_REF] on a McKenna burner. Experimental measurements are advantageously available at three fuel equivalence ratios φ, especially delivered as profiles of the major gas species mole fraction, the primary particle diameter, and the soot volume fraction. The operating conditions for the three test cases are reported in Tab. 1. Soot volume fraction profiles measured by intrusive sampling are given in Fig. 5 for indicative purposes only.

Also note that the diameter of the transported spherical particles will be compared to the measured primary particle diameter in the following, which is consistent with the collision model used in the present study.

Numerical setup

The computational domain is 2D, 40 mm long and 1.5 mm wide. The mesh grid contains a total of 210 cells and is refined within the flame front, with a single cell in the transverse direction. The inlet boundary conditions are deduced from Tab. 1 for each test case. Atmospheric pressure combined with non-reflective acoustic conditions are imposed at the outlet, symmetry conditions are used elsewhere. Temperature profiles suggested by the International Sooting Flame Workshop [START_REF]International Sooting Flame (ISF) Workshop[END_REF] are additionally prescribed at each iteration of the CHARME solver. The time step is fixed to 1 µs for both the gas phase and the dispersed phase solver, and the simulated physical time is chosen large enough to reach a steady-state solution.

A preliminary computation with Cantera [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF] was carried out at each equivalence ratio. The solution obtained by only considering the gas phase reaction mechanism allows here to define a proper initial state for the CEDRE calculation.

Evaluation of the optimisation algorithm

Reduction intensity of the optimisation algorithm

In order to evaluate the statistical convergence of the quadtree-like algorithm, six cases are simulated at the intermediate equivalence ratio φ = 2.64. All are identical, except for the maximum number N np of numerical particles per cell, successively limited to 64, 128, 256, 512, 1024 and 2048 every 100 time steps. In each case, the properties of the particles are measured 20 mm and 30 mm above the burner for comparison purposes, these two locations being denoted as measurement points in the following. The gap between d 10 and its reference value remains lower than 1 %. In contrast, the gap reaches 14 % for σ d at the point located 30 mm above the burner for the lowest value of N np , i.e. 64. Thus, while the algorithm allows to preserve the mean of the particle size distribution, even for low values of the parameter N np , the gap for the standard deviation increases with the decrease of N np . This behaviour, intrinsically related to the reduction of the particles population, is not surprising. N np should just be set as a compromise between the computational cost and the accuracy in predicting the standard deviation, depending on the application. In the following, it will be taken equal to 256, which is assumed sufficient for the comparison with the experimental measurements, only available upstream of the second point.

Mean and standard deviation

Particle size distribution

Figure 3 shows the particle size distribution at both measurement points and for all cases. As expected, the deviation from the reference case varies inversely with the parameter N np . Nonetheless, even for strong limitations of the maximum number of particles per cell, the shape of the distribution -representative of the individual history of the particles -is preserved. Two peaks are visible at each point. The first one matches the diameter of nucleated soot particles d 0 . The second one is lower bounded by the diameter d 1 = 2 1/3 d crit , which corresponds to the diameter of a particle formed through the coalescence of two particles with a diameter equal to d crit . It indicates that coalescence stops very early downstream of the flame. The evolution of the soot particles, evidenced by the weakening of the first peak and the shift of the second peak towards higher diameters with increasing height above the burner (HAB), is associated to surface growth and condensation, as well as to coalescence of the remaining or newly nucleated small particles with the bigger ones.

Reduction frequency of the optimisation algorithm

The quadtree-like algorithm is now successively applied every N ts time steps, with N ts = 1, 10 and 1000, N np being fixed to 256. The three cases are simulated at the intermediate equivalence ratio φ = 2.64 to assess the impact of this reduction frequency on the first and second order statistical moments of the particle size distribution. The gap between d 10 or σ d and their reference values always remains smaller than 2 %, so that the influence of the reduction frequency is considered to be negligible within the explored range of values. The parameter N ts will therefore be taken equal to 100 in Subsection 5.4 to keep the computational cost moderate as discussed below.

Remark on the computational cost

As a guide for comparison to sectional methods or methods of moments, the computational cost of the Lagrangian solver represents 38 % of the total cost for N np = 256 and N ts = 100, 53 transport equations being solved for the gas-phase. For N ts = 100, the mean computational cost of the Lagrangian solver decreases from 4.5 s/time step for N np = 2048 to 0.64 s/time step for N np = 64, due to the reduced size of the particles population. For N np = 256, the mean computational cost of the Lagrangian solver increases from 0.89 s/time step for N ts = 100 to 2.3 s/time step for N ts = 1000, due to the growth of the particles population between two reductions. From N ts = 100 to N ts = 1, the mean computational cost also increases following the high number of reduction events, for which monitoring files are additionally written.

Note that the performance of the Lagrangian solver could be improved by dynamic load balancing in the case of massively parallel simulations. From a global point of view, the model proves its ability to reproduce the growth of both quantities with increasing values of φ. Fairly good agreement with the measurements is especially noticed for the soot volume fraction at φ = 2.64 and φ = 2.94. However, while the mean particle diameter is correctly predicted in the first case, it is underestimated in the second one. The sensitivity of the model to the limitation of coalescence processes will therefore be discussed below. At φ = 2.34, the soot volume fraction and the mean particle diameter are both underrated. The temperature profile being fixed, the influence of two parameters has here to be investigated, namely the collision efficiency γ PAH ( * ) , introduced in Subsection 4.3, and the kinetic constants for C 2 H 2 addition. Note that similar discrepancies are also observed by Eberle et al. [START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF] at low fuel equivalence ratios for a comparable configuration and an identical nucleation mechanism, showing that an upgrade of the nucleation model could be considered. Such reflections are however beyond the scope of the present study. (reference simulation) and 1.0. Both quantities increase with γ PAH ( * ) due to improved condensation. The consumption of the PAH molecules and radicals also modifies the distribution of these species above the burner, which reduces the soot nucleation rate, displayed in Fig. 7, and the particle number density. The growth of the particles however overcompensates this phenomenon, as evidenced by the elevation of the soot volume fraction. For the sake of completeness, Fig. 7 shows the aforementioned profiles for two additional combinations of collision efficiencies, i.e. 1.0 for PAH molecules with 0.3 for radical species, and vice versa.

Limitation of coalescence processes

The introduction of such an asymmetry demonstrates that PAH condensation onto soot particles overrides that attributed to PAH * . In the present study, α is calculated through the correlation of Kazakov et al. [START_REF] Kazakov | Detailed modeling of soot formation in laminar premixed ethylene flames at a pressure of 10 bar[END_REF] from the maximum flame temperature assimilated to the adiabatic one T ad . This is reported in Tab. 2 for the three fuel equivalence ratios investigated. At φ = 2.34, no significant sensitivity of the volume fraction and mean particle diameter to variations of the temperature value used could however be highlighted. Hence, the present study focuses on their sensitivity to the pre-exponential factor A sg of the per-site reaction rate constant.

Figure 8 shows the profiles of soot volume fraction and mean particle diameter at φ = 2.34 for A sg /A 0 sg equal to 1 (reference simulation) and 10. As expected, both f v and d 10 increase with the pre-exponential factor. Nevertheless, the higher consumption of acetylene due to surface growth has no discernible impact on its distribution above the burner. The reduction of the soot nucleation rate, although limited, is therefore directly linked to the stronger consumption of PAH molecules and radicals with increasing particle size in the flame zone.

Discussion

According to the sensitivity analysis, a critical diameter equal to 18 nm offers the best compromise for the three fuel equivalence ratios, at least as a first approximation. As a result, this value is retained, although its validity for different fuels, combustion and turbulence regimes should be investigated in the future.

Concerning the collision efficiencies for PAH/PAH * condensation and the kinetic constants for C 2 H 2 addition, no change in the standard values can be justified without further validation. The reported results however highlight the major role of condensation on both particle nucleation and growth.

Temperature profiles

In order to emphasize the sensitivity of the soot formation and evolution processes to the temperature level, the profile imposed at φ = 2.34 is rescaled through an homothetic transformation, so that profiles with maximum values varying of ±50 K compared to the initial one ( T ) can be tested. Strictly speaking, for a burner-stabilised flame, the temperature profile should result from the coupling between injection, combustion and heat transfer to the gas phase, the dispersed phase and the burner itself. Such a simulation is however beyond the scope of the present study.

Figures 5 and6 compare the profiles of soot volume fraction and mean particle diameter for increasing temperature levels, d crit being fixed to its default value, i.e. 13 nm. As discussed in the work of Mueller et al. [START_REF] Mueller | Hybrid Method of Moments for modeling soot formation and growth[END_REF], the soot volume fraction is higher with lower flame temperatures, which is here directly linked to enhanced nucleation. The mean particle diameter also increases following the higher number of particles and the induced multiplication of coalescence events. Coalescence finally stops earlier, hence a ceiling of the particle diameter, already visible at φ = 2.64 and φ = 2.94, is observed.

Major gas species

Conclusion

A Lagrangian description of soot particles evolution and dynamics, coupled with a reduced gas phase reaction mechanism and a sectional model for polycyclic aromatic hydrocarbons, is proposed. It was applied with some success to the simulation of laminar premixed ethylene-air flames at three fuel equivalence ratios. The statistical convergence of the population reduction algorithm -as a function of the reduction intensity and frequency -was first highlighted. The evolution model was then contrasted with experimental measurements, showing its ability to reproduce the growth of the soot volume fraction and the mean particle diameter with increasing equivalence ratio and decreasing temperature. In particular, the suggested description of coalescence processes allows fairly good agreement of the mean particle diameter with the measurements, in terms of both curve shape and level. Finally, the strong sensitivity of soot nucleation and growth to PAH condensation was evidenced.

The present study, restricted to laminar flames, constitutes a prerequisite towards the use of the Lagrangian tracking method in a realistic turbulent combustion chamber. In addition, the modelling of soot aggregates, relevant to radiative heat transfer, has to be addressed, as well as the sensitivity of the critical coalescence diameter over which coalescence is prohibited to the fuel composition and the combustion regime.

  particle kinetic energy + m p c p dT p dt ↑ Variation of the particle thermal energy dt,

Figure 1 :

 1 Figure 1: Definition of PAH i and PAH * i sections [39].

Figure 2

 2 Figure 2 compares the mean particle diameter d 10 and its standard deviation σ d at both measurement points and for all cases. For the sake of simplicity, the variables are normalised by the values reached at N np = 2048, this simulation being considered as a reference given the very good convergence of d 10 and σ d with increasing N np .

Figure 2 :

 2 Figure 2: Normalised mean particle diameter (I) and standard deviation(G) as a function of N np at the measurement points located 20 mm and 30 mm above the burner. For N np = 2048, d 10 = 25.4 nm and σ d = 2.64 nm (a), d 10 = 26.9 nm and σ d = 2.60 nm (b).

Figure 3 :

 3 Figure 3: Particle size distribution as a function of N np at the measurement points located 20 mm and 30 mm above the burner (G = 64, = 128, G = 256, = 512, G = 1024, = 2048). The bottom figures are a zoom of the top ones in the region delimited by the black dashed lines.

Figure 4

 4 Figure 4 compares the mean particle diameter d 10 and its standard deviation σ d at both measurement points and for all cases. The variables are normalised by the values previously obtained for N ts = 100 for the sake of clarity.

Figure 4 :

 4 Figure 4: Normalised mean particle diameter (I) and standard deviation (G) as a function of N ts at the measurement points located 20 mm and 30 mm above the burner. For N ts = 100, d 10 = 25.4 nm and σ d = 2.57 nm (a), d 10 = 26.9 nm and σ d = 2.50 nm (b).

5. 4 .

 4 Figures 5 and 6 present the profiles of soot volume fraction f v and mean particle diameter d 10 at the three fuel equivalence ratios. From a global point of view, the model proves its ability to reproduce the

Figure 6 Figure 7

 67 Figure 6 show the profile of d 10 at the three fuel equivalence ratios for two additional values of the critical diameter d crit , 9.1 and 18 expected, the increase of d crit promotes which augments the size of the particles. At φ = 2.34, coalescence stops far downstream of the flame zone, so that the impact of the critical diameter is only discernible close to the outlet. Note that the profiles of f v are not presented since not deeply affected by the critical diameter. Collision efficiency γ PAH ( * ) Figure 7 compare the profiles of f v and d 10 at φ = 2.34 for two values of the collision efficiency, i.e. 0.3

Figure 5 :

 5 Figure 5: Soot volume fraction along height above the burner: Xu -laser extinction (M), Xu -sampling (G), CE-DRE at d crit = 13 nm ( ), CEDRE at T -50 K ( ) and T +50 K ().

Figure 6 :

 6 Figure 6: Mean particle diameter along height above the burner: Xu (G), CEDRE at d crit = 9.1 nm ( ), 13 nm ( ), 18 nm ( ), CEDRE at T -50 K ( ) and T +50 K (). The critical diameters are reminded by the black dashed lines.

Figures 9 and 10

 10 Figures 9 and 10 present the mole fraction profiles of the major gas species at the three fuel equivalence ratios for the retained value of d crit , i.e. 18 nm. Fairly good agreement with the measurements is noticed for the main combustion products (CO, CO 2 and H 2 O) and for the intermediate species C 2 H 2 . The ethylene mole fraction is overestimated downstream the flame front in all cases. However, at the measured levels, experimental uncertainties are known to exceed 15 %[START_REF] Xu | Soot formation in laminar premixed ethylene/air flames at atmospheric pressure[END_REF].

Figure 7 :

 7 Figure 7: Influence of collision efficiency for condensation: γ PAH /γ PAH * = 0.3/1.0 (M), 1.0/0.3 (G), 0.3/0.3 ( ), 1.0/1.0 ( ). Variables of interest along height above the burner.

  • ωsg [kg/(s. cm 3 )] (d) Acetylene mass reaction rate (surface growth)

Figure 8 :

 8 Figure 8: Influence of the chemical kinetics for C 2 H 2 addition: A sg = 8.4 × 10 5 m 3 /(mol. s) (reference, ), 8.4 × 10 6 m 3 /(mol. s) (). Variables of interest along height above the burner.

Figure 9 :

 9 Figure 9: Mole fractions of the main combustion products: CO ( / M), CO 2 ( / G) and H 2 O ( / I).Lines are for CEDRE, symbols for experimental measurements.

Figure 10 :

 10 Figure 10: C 2 H 2 C H 4 mole fractions: 5 × C 2 H 2 ( /M) and 50 × C 2 H 4 ( / G). Lines are for CE-DRE, symbols for experimental measurements.

Table 1 :

 1 Test cases inlet conditions.

	Equivalence ratio	[-]	2.34	2.64 2.94
	Cold gas velocity	[cm/s] 6.73	6.73 6.73
	Cold gas temperature	[K]	298	298	298
	X C2H4	[%]	14.08 15.60 17.0
	X O2	[%]	18.05 17.73 17.4
	X N2	[%]	67.84 66.67 65.6

Table 2 :

 2 Adiabatic flame temperatures (estimated from a preliminary equilibrium calculation carried out with Cantera and only considering the gas phase reaction mechanism).Kinetic of C 2 H 2 additionAs indicated in Subsection 4.3, the reaction rate constant for C 2 H 2 addition is proportional to the fraction of active surface sites α.

	Equivalence ratio	[-] 2.34 2.64 2.94
	Estimated adiabatic [K] 1728 1583 1449 flame temperature

In terms of code architecture

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:


Jingwei Liang
email: jingwei.liang@ensicaen.fr

Jalal M Fadili
email: jalal.fadili@ensicaen.fr

Gabriel Peyré
email: gabriel.peyre@ceremade.dauphine.fr

Local Linear Convergence of Inertial Forward-Backward Splitting for Low Complexity Regularization

In this abstract, we consider the inertial Forward-Backward (iFB) splitting method and its special cases (Forward-Backward/ISTA and FISTA). Under the assumption that the non-smooth part of the objective is partly smooth relative to an active smooth manifold, we show that iFB-type methods (i) identify the active manifold in finite time, then (ii) enter a local linear convergence regime that we characterize precisely. This gives a grounded and unified explanation to the typical behaviour that has been observed numerically for many low-complexity regularizers, including 1 , 1,2 -norms, total variation (TV) and nuclear norm to name a few. The obtained results are illustrated by concrete examples.

I. INTRODUCTION

Consider the following structured optimization problem

min x∈R n Φ(x) def = F (x) + J(x) , (P)
where J ∈ Γ0(R n ), the set of proper, lower semi-continuous and convex functions, F is convex, C 1,1 (R n ) with ∇F being β-Lipschitz continuous. We assume that Argmin Φ = ∅.

In this paper, we consider a generic form of inertial Forward-Backward for solving (P) which reads,

y k a = x k + a k (x k -x k-1 ), y k b = x k + b k (x k -x k-1 ), x k+1 = Proxγ k J y k a -γ k ∇F (y k b ) , (I.1)
where

a k ∈ [0, ā] and b k ∈ [0, b], (ā, b) ∈ [0, 1] 2
, and the step-size 0 < γ ≤ γ k ≤ γ < 2/β. For γ > 0, the proximity operator is defined as ProxγJ (x) = argmin z∈R n 1 2γ ||z -x|| 2 + J(z). iFB (I.1) covers various special cases in the literature, including the (unrelaxed) Forward-Backward (FB) [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF] and FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. In the original FISTA, only convergence of the objective function is guaranteed. Recently in [START_REF] Chambolle | How to make sure the iterates of FISTA converge[END_REF], the iterates are proved to be convergent under a k = b k = (t k-1 -1)/t k where t k = (k + p -1)/p, p ≥ 2.

II. PARTLY SMOOTH FUNCTIONS AND FINITE IDENTIFICATION

The class of partly smooth functions [START_REF] Lewis | Active sets, nonsmoothness, and sensitivity[END_REF], is specialized here to functions in Γ0(R n ). Denote par(C) the linear subspace parallel to the non-empty convex set C ⊂ R n , and ri(C) its relative interior.

Definition II.1. Let J ∈ Γ0(R n ) and x ∈ R n such that ∂J(x) = ∅. J is partly smooth at x relative to a set M containing x if (Smoothness) M is a C 2 -manifold, J|M is C 2 around x; (Sharpness) The tangent space TM(x) = Tx def = par(∂J(x)) ⊥ ; (Continuity) The ∂J is continuous at x relative to M.
Examples of such functions are given in Section IV, see also [START_REF] Vaiter | Partly smooth regularization of inverse problems[END_REF].

Theorem II.2 (Finite activity identification). Suppose x k converges to a minimizer x of (P) such that J is partly smooth at x relative to Mx , and

-∇F (x ) ∈ ri ∂J(x ) , (II.1)
then there exists a K > 0 such that for all k ≥ K, x k ∈ Mx . If moreover Mx is affine/linear, then y k a , y k b ∈ Mx for k > K.

Condition (II.1) can be viewed as a geometric generalization of the strict complementarity of non-linear programming, and is almost necessary for the finite identification of Mx [START_REF] Lewis | Active sets, nonsmoothness, and sensitivity[END_REF].

III. LOCAL LINEAR CONVERGENCE

We now turn to the local linear convergence of the iFB-type methods with partly smooth functions. For space limitations, we mainly focus on the case where a k = b k , and denote

d k+1 = x k+1 -x x k -x .
Theorem III.1. We assume the conditions of Theorem II.2 hold. If moreover F is C 2 near x and there exists α ≥ 0 such that

PT x ∇ 2 F (x )PT x αId.
Then for all k large enough, we have

1) Q-linear rate: if 0 < γ ≤ γ k ≤ γ < min(2αβ -2 , 2β -1 ), then
given any ρ such that 1 > ρ ≥ ρ k , the iterates satisfy

||x k+1 -x || 2 ≤ ||d k+1 || 2 ≤ ρ||d k || 2 ,
where η = max q(γ), q(γ)

∈ [0, 1[, q(γ) = 1-2αγ +β 2 γ 2 , ρ k =    (1+a k )η+ √ (1+a k ) 2 η 2 -4a k η 2 , η ∈ [ 4a k (1+a k ) 2 , 1[, √ a k η, η ∈ [0, 4a k (1+a k ) 2 ]. 2) R-linear rate: if Mx is affine/linear, then ||x k+1 -x || 2 ≤ ||d k+1 || 2 ≤ ρ k ||d k || 2 , where ρ k ∈ [0, 1[ ρ k =    |(1+a k )η k |+ (1+a k ) 2 η 2 k -4a k η k 2 , η k ∈] -1, 0] ∪ [ 4a k (1+a k ) 2 , 1[, √ a k η k , η k ∈ [0, 4a k (1+a k ) 2 ],
and Example IV.4 (Nuclear norm). The nuclear norm is partly smooth relative to the manifold of fixed rank matrices, M = z ∈ R n 1 ×n 2 : rank(z) = r .

η k ∈] -1, 1[ is an eigenvalue of Id -γ k P T 1 0 ∇ 2 F (x + t(y k a -x ))dtP T . IV. NUMERICAL EXPERIMENTS Example IV.1 ( 1-norm). The 1-norm is partly smooth relative to M = Tx = {u ∈ R n : supp(u) ⊆ supp(x)}. Example IV.2 ( 1,2-norm). 1,2-norm is partly smooth relative to M = Tx = {u ∈ R n : supp B (u) ⊆ supp B (x)},
We now consider the problem min x∈R n 1 2 ||y -Ax|| 2 + λJ(x), where y ∈ R m is the observation, A : R n → R m is drawn from the standard Gaussian ensemble, and λ > 0 is the regularization parameter. The convergence profiles are depicted in Figure 1. ). All algorithms were tested with γ k ≡ 1/||A|| 2 . The solid lines are the practical observed profiles the dashed ones the theoretical predictions. The beginning of the dashed lines are the points when x k identifies the manifold Mx . As one can observe, FISTA has the fastest manifold identification, however, locally it is the slowest for all tested examples. Indeed, when the manifold is affine, it can be shown from Theorem III.1 that ρ k ∈]η k , √ η k ] for a k > η k , i.e. FISTA is locally slower than FB.

  where supp B (x) = {b : x b = 0}, and b∈B b = {1, . . . , n}. Example IV.3 (TV semi-norm). The TV semi-norm ||x|| TV = ||∇x|| 1 is partly smooth relative the subspace M = Tx = u ∈ R n : supp(∇u) ⊆ I , I = supp(∇x).

Fig. 1 :

 1 Fig. 1: Local linear convergence of iFB-type methods in terms of ||x k -x ||. The forward model of the problem of interests reads y= Ax0 + ε, ε ∼ N (0, δ 2 ). (a) 1-norm, (m, n) = (48, 128), x0 is 8-sparse; (b) 1,2-norm, (m, n) = (60, 128), x0 has 3 non-zero blocks with block-size 4; (c) 1D TV semi-norm, (m, n) = (48, 128), ∇x0 is 8-sparse; (d) Nuclear norm, (m, n) = (1425, 2500), x0 ∈ R 50×50and rank(x0) = 5. The red, black and blue lines are respectively the results of FB, FISTA[START_REF] Chambolle | How to make sure the iterates of FISTA converge[END_REF] and iFB (with a k = b k ≡ √ 5 -2.01). All algorithms were tested with γ k ≡ 1/||A|| 2 . The solid lines are the practical observed profiles the dashed ones the theoretical predictions. The beginning of the dashed lines are the points when x k identifies the manifold Mx . As one can observe, FISTA has the fastest manifold identification, however, locally it is the slowest for all tested examples. Indeed, when the manifold is affine, it can be shown from Theorem III.1 that ρ k ∈]η k , √ η k ] for a k > η k , i.e. FISTA is locally slower than FB.

ACKNOWLEDGMENT

This work has been partly supported by the European Research Council (ERC project SIGMA-Vision). JF was partly supported by Institut Universitaire de France.


Jingwei Liang
email: jingwei.liang@ensicaen.fr

Jalal M Fadili
email: jalal.fadili@ensicaen.fr

Gabriel Peyré
email: gabriel.peyre@ceremade.dauphine.fr

Russell Luke
email: r.luke@math.uni-goettingen.de

Local Linear Convergence of Douglas-Rachford/ADMM for Low Complexity Regularization

The Douglas-Rachford (DR) and ADMM algorithms have become popular to solve sparse recovery problems and beyond. The goal of this work is to understand the local convergence behaviour of DR/ADMM which have been observed in practice to exhibit local linear convergence. We show that when the involved functions (resp. their Legendre-Fenchel conjugates) are partly smooth, the DR (resp. ADMM) method identifies their associated active manifolds in finite time. Moreover, when these functions are partly polyhedral, we prove that DR (resp. ADMM) is locally linearly convergent with a rate in terms of the cosine of the Friedrichs angle between the tangent spaces of the two active manifolds. This is illustrated by several concrete examples and supported by numerical experiments.

I. INTRODUCTION

In this work, we consider the problem of solving

min x∈R n J(x) + G(x), (1)
where both J and G are in Γ0(R n ), the class of proper, lower semi-continuous and convex functions. We assume that ri(domJ) ∩ ri(domG) = ∅, where ri(C) is the relative interior of the nonempty convex set C, and domF is the domain of the function F . We also assume that the set of minimizers of (1) is non-empty, and the two functions are simple (prox γJ , prox γG , γ > 0, are easy to compute), where the proximity operator is defined, for γ > 0, as prox γJ (z) = argmin x∈R n 1 2 ||x -z|| 2 + γJ(x). An efficient and provably convergent algorithm for solving (1) is the DR method [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF], which reads, in its relaxed form, (

z k+1 = (1 -λ k )z k + λ k BDRz k , x k+1 = prox γJ z k+1 , (2)
where

BDR def. = 1 2 `(2prox γJ -Id) • (2prox γG -Id) + Id ´, for γ > 0, λ k ∈]0, 2] with P k∈N λ k (2 -λ k ) = +∞.
Since the set of minimizers of (1) is non-empty, so is the set of fixed points fix(BDR) since the former is nothing but prox γJ (fix(BDR)). See [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF] for a more detailed account on DR. Though we focus in the following on DR, all our results readily apply to ADMM since it is the DR applied to the Fenchel dual problem of (1).

II. PARTLY SMOOTH FUNCTIONS AND FINITE IDENTIFICATION

Beside J, G ∈ Γ0(R n ), our central assumption is that J, G are partly smooth functions. Partial smoothness was originally defined in [START_REF] Lewis | Active sets, nonsmoothness, and sensitivity[END_REF]. Here we specialize it to the case of functions in Γ0(R n ). Denote par(C) the subspace parallel to the non-empty convex set C ⊂ R n .

Definition II.1. Let J ∈ Γ0(R n ) and x ∈ R n such that ∂J(x) = ∅. J is partly smooth at x relative to a set M containing x if (Smoothness) M is a C 2 -manifold, J|M is C 2 around x;
(Sharpness) The tangent space Tx(M) = Tx def.

= par(∂J(x)); (Continuity) The ∂J is continuous at x relative to M. The class of partly smooth functions at x relative to M is denoted as PSFx(M).

The class of PSF is very large. Popular examples in signal processing and machine learning include 1, 1,2, ∞ norms, TV semi-norm and nuclear norm, see also [START_REF] Vaiter | Partly smooth regularization of inverse problems[END_REF]. Now define the variable

v k = prox γG (2prox γJ -Id)z k .
Theorem II.2 (Finite activity identification). Let the DR scheme (2) be used to create a sequence

(z k , x k , v k ). Then (z k , x k , v k ) → (z , x , x )
, where z ∈ fix(BDR) and x is a global minimizer of (1). Assume that

J ∈ PSFx (M J x ), G ∈ PSFx (M G x ), and
z ∈ x + γ `ri(∂J(x )) ∩ ri(-∂G(x )) ´. (3)
Then, the DR scheme has the finite activity identification property, i.e. for all k sufficiently large,

(x k , v k ) ∈ M J x × M G x . Condition (3) implies that 0 ∈ ri(∂J(x ) + ∂G(x ))
, which can be viewed as a geometric generalization of the strict complementarity of non-linear programming. In a compressed sensing scenario, it can be guaranteed for a sufficiently large number of measurements.

III. LOCAL LINEAR CONVERGENCE OF DR

We now turn to local linear convergence properties of DR for the case of locally polyhedral functions. This is a subclass of partly smooth functions, whose epigraphs look locally like a polyhedron. In the following, we will refer to the Friedrichs angle between two subspaces V and W , denoted θF (V, W ). In fact, θF (V, W ) is the (d + 1)-th principal angle between V and W , where d = dim(V ∩ W ), see also [START_REF] Bauschke | Optimal rates of convergence of matrices with applications[END_REF].

Theorem III.1. Assume that J and G are locally polyhedral, and the conditions of Theorem II.2 hold with λ k ≡ λ. Then there exists K > 0 such that for all k K,

||z k -z || ρ k ||z0 -z ||, (4)
where

ρ = p (1 -λ) 2 + λ(2 -λ) cos 2 θF (T J x , T G x ) ∈ [0, 1
[. This rate is optimal. For the special case of basis pursuit, we recover the result of [START_REF] Demanet | Eventual linear convergence of the Douglas-Rachford iteration for basis pursuit[END_REF], but with less stringent assumptions.

IV. NUMERICAL EXPERIMENTS

As examples, we consider the 1, ∞ norms and the anisotropic TV semi-norm which are all polyhedral, hence partly smooth relative the following subspaces

1 : Tx = {u ∈ R n : supp(u) ⊆ supp(x)}, ∞ : Tx = ˘u : uI = rsI , r ∈ R ¯, s = sign(x), I = {i : |xi| = ||x||∞}, TV : Tx = ˘u ∈ R n : supp(∇u) ⊆ I ¯, I = supp(∇x),
where ∇ is the gradient operator.

Figure1 displays the observed and predicted convergence profiles of DR when solving several problem instances, including compressed sensing, denoising and inpainting.

Fig. 1 .

 1 Fig. 1. Observed (solid) and predicted (dashed) convergence profiles of DR (2) in terms of ||z k -z ||. For the first 4 subfigures, we solve a problem of the form min x∈R n J(x) s.t. Ax = y, where A is either drawn randomly from the standard Gaussian ensemble (CS), or random binary (inpainting). (a) CS with J = || • || 1 , A ∈ R 48×128 . (b) CS with J = || • ||∞, A ∈ R 120×128 . (c) CS with J = || • || TV , A ∈ R 48×128 . (d) TV image inpainting, A ∈ R 512×1024 . (e) Uniform noise removal by solving min x∈R 128 ||x|| TV s.t. ||x -y||∞ δ. (f) Outliers removal by solving min x∈R 128 λ||x|| TV + ||x -y|| 1 . The starting points of the dashed lines are the iteration at which the active subspaces are identified.

ACKNOWLEDGMENT

This work has been partly supported by the European Research Council (ERC project SIGMA-Vision). JF was partly supported by Institut Universitaire de France.


Jean-Yves Fortin
email: jean-yves.fortin@univ-lorraine.fr

Charge oscillations in a simple model of interacting magnetic orbits

Keywords: Magnetic field, Landau levels, Tunneling

Exact eigenstates for a set of two or more interacting electronic orbits in a magnetic field are studied for a class of factorized Hamiltonians with coupled Fermi surfaces.

We study the condition for the existence of annihilation-creation operators that allows for the construction of eigenstates. For the case of two interacting cyclotronic orbits, we consider the oscillations of the overlap function and the transfer of charge density between the orbits as function of the inverse field. The expressions of the Fourier frequencies are given in the semiclassical regime and they depend on the geometrical structure of the electronic bands. A generalization of this construction is provided for a chain of several interacting orbits with exact eigenfunctions.

I. INTRODUCTION

The two-state problem in quantum mechanics, such as the Rosen-Zener two-level model, gives precise information on how the wavefunction propagates through a junction and and has important applications in time-driven quantum systems [1][START_REF] Hänggi | Quantum Transport and Dissipation[END_REF][START_REF] Torosov | [END_REF][4], where resonance and phase shifts are studied in details for exact solvable cases of potentials. This general problem can be applied in the case of magnetic breakdown as well, where a quasi-particle orbits a Fermi surface which can be composed of multiple sheets connected by junctions through which the particle can tunnel. A realization of such Fermi surfaces, resembling a linear chain of coupled orbits, is presented in figure 1 for the organic conductor (BEDO-TTF) 5 [CsHg(SCN) 4 ] 2 [5] (BEDO-TTF is the abbreviation for bis-ethylenedioxi-tetrathiafulvalene molecule): an incoming wave-packet (a) on a giant orbit β is transmitted to the small cyclotronic orbit α (b) and reflected onto the same orbit β (c) within the chain. Fermi surfaces with a finite number of interacting orbits also exist and the Fourier spectrum of their magnetic oscillations have been studied in details. For example, compensated Fermi structures with only three bands, made of one hole and two electron pockets, can be found in compound α-'pseudo-κ'-(ET) 4 H 3 O[Fe(C 2 O 4 ) 3 ]•(C 6 H 4 Br 2 ), where the α-type and 'pseudo-κ'-type are conducting and insulating layers respectively [6]. At the magnetic breakdown junction, the Hamiltonian can be linearized and the two sheets hybridized with some energy coupling g. This simplest form of two-level Hamiltonian was solved by Rosen and Zener in a different context [7] using this approximation around the tunneling region. The probability of tunnelling is exponentially small in the ratio between a breakdown field and the magnetic field [8], and a field dependence of the Onsager phase [9,10], in addition to the π/2 phase each time the particle is reflected at the junction.

In this paper we would like to construct simple and exact eigenstates for multi-band conductors that incorporate magnetic tunneling on such Fermi surfaces made for example of a chain of two or more coupled orbits [11].

In general there exist several possible techniques for solving Hamiltonians with complex potentials, in addition to the semi-classical techniques provided by WKB method [12,13] dealing with scattering matrices for interacting orbits. For example, the factorization approach, or Darboux method, [14][15][16][17][START_REF] Cooper | Supersymmetry and quantum mechanics[END_REF][START_REF] Fellows | [END_REF] consists in considering an Hamiltonian which can be factorized as

H 1 = D † D = -∂ 2 x + V 1 (x), where D = ∂ x + Q(x) and D † = -∂ x + Q(x)
H 2 = DD † = -∂ 2 x + V 2 (x). V 1 (
x) and V 2 (x) are supersymmetric partner potentials expressed with Q(x) and its derivative. It follows that the eigenfunctions and eigenvalues of H 2 are closely related to those of H 1 . If H 1 is exactly solvable, then H 2 is also exactly solvable [15,20,21]. This allows to construct exact solutions of the Schrödinger equation with non-trivial potentials V 2 (x), for example for classes of potentials with two wells [22].

Another possibility that we explore in this paper is to consider the same form of factorized Hamiltonian H = D † D, where Q(x) is now a generalized complex matrix but whose structure is constrained by the condition that the commutator between the operators D and D † satisfies [D, D † ] = 1. In this case D and D † are called ladder or annihilation/creation operators and it is straightforward to construct the eigenstates from the knowledge of the ground state, and the energy spectrum is simply discrete. In the case of a single band, the constraint leads simply to Q(x) being a linear function of x (harmonic oscillator), but for two-band or multi-band systems Q(x) can have a more complicated structure. In addition, we will see that the matrix Q(x) can be chosen such that the two sheets of the Fermi surface of the Hamiltonian present a magnetic gap, which is the purpose of this paper. We will also show that a simple extension to multi-band Hamiltonians leads to a modelization of the linear chain of coupled orbits discussed above and in figure 1.

II. MODEL

The problem of treating Bloch electrons in a magnetic field has been generally approached using the Peierls substitution [24,25]. The Bloch Hamiltonian for the energy dispersion H(k) is replaced by H(-i∇ + eA/ ), where -e is the charge of the electron. This approximation is usually valid near the Fermi surface, or for tight-binding Hamiltonians [26]. The quantification of the quasi-momentum -i∇ + eA/ in the Landau gauge A = (0, Bx 1 , 0) leads to an analogy with a one-dimensional particle motion in a potential [START_REF] Landau | Statistical Physics: Part 2 : Theory of Condensed State, Volume 9 of Course of theoretical physics[END_REF] if we identify the dimensionless operators x = a 2 (-i∂ 2 + x 1 /l 2 B ) and p = -ia 1 ∂ 1 with effective position and momentum, where (a 1 , a 2 ) are the typical sizes of the Brillouin zone unit cell and l -2 B = eB/ the inverse squared of the magnetic length. These operators satisfy the commutation rela- It can also be identified as the ratio between the magnetic flux though the unit cell and the quantum flux h/e. Usually it is considered as a small control parameter in the semiclassical theories if the cell units in the Brillouin zone are of the order of the Angström, or if the magnetic field is of the order of the Tesla. We should notice that in the case of positive charges, e is instead negative and to keep the effective Planck constant h positive we can invert the role of x and p to keep the same commutation relation. Starting with the Hamiltonian expression for an electron in a uniform magnetic field B, where H(k) is quadratic in momenta k and describing a simple harmonic oscillator near the Fermi surface which can be therefore factorized in terms of ladder or creation and annihilation operators [START_REF] Yoshioka | The Quantum Hall Effect[END_REF],

tion [x, p] = ia 1 a 2 eB/ ≡ ih, where h = a 1 a 2 /l 2 B is a dimensionless effective Planck constant.
we introduce the representation p = -ih∂ x and consider a class of differential operators for multiband systems defined by

D † = 1 √ 2h -hσ 0 ∂ x + Q † (x) , D = 1 √ 2h [hσ 0 ∂ x + Q(x)] (1)
from which is constructed the factorized Hermitian Hamiltonian Ĥ

Ĥ = 2hD † D = -σ 0 h 2 ∂ 2 x + [Q † (x) -Q(x)]h∂ x -hQ ′ (x) + Q † (x)Q(x) (2)
Q(x) is a general matrix depending on x, and σ i=0,••• ,3 are the Pauli matrices. From a supersymmetric point of view [15], it is known that the product Q † (x)Q(x) describes interaction between bosons, whereas the term hQ ′ (x) deals with fermion-boson interaction, both being dependent of the same function Q(x). The fermion-boson interaction is absent if Q ′ (x) is constant, as for the bosonic harmonic oscillator. The classical version H c of the Hamiltonian in the limit h → 0 is obtained by replacing the operator p by momentum variable p

H c (x, p) = (Q † (x) -iσ 0 p)(Q(x) + iσ 0 p) = σ 0 p 2 + Q † (x)Q(x) + ip[Q † (x) -Q(x)] (3)
The Fermi surface H x (x, p) = Eσ 0 in the plane (x, p) is given by the eigenvalues of H c at constant energy E, and is generally composed of two closed surfaces. In the following we consider examples of matrices Q(x) that lead to magnetic tunneling between the two sheets of the Fermi surface. The commutation relation between operators D and D † is given by matrix

P (x) [D, D † ] = 1 2 ∂ x Q(x) + ∂ x Q † (x) + 1 2h [Q(x), Q † (x)] = P (x) (4)
If P (x) is the unity matrix, this relation is the usual bosonic commutator for creation and annihilation operators D † and D respectively. And the eigenstates Ψ n = (φ n , φn ) ⊤ are constructed from the ground state Ψ 0 such that DΨ 0 = 0, by successively acting the ladder operator D † on Ψ 0 . This is, for example, the case when the system is composed of two independent harmonic oscillators, centered at positions ±x c

Q(x) = Q † (x) =   x + x c 0 0 x -x c   , P =   1 0 0 1   (5)
We would like the study the condition for which tunneling can occur when P (x) = 1 for general matrices Q(x). We first decompose Q(x) in the base of the Pauli matrices Q(x) = 3 i=0 α i (x)σ i , with α i complex functions, and for which the commutator P (x) is equal to

P (x) = 1 2 ∂ x (α 0 (x) + ᾱ0 (x))σ 0 + 1 2 ∂ x (α + ᾱ).σ + i h (α × ᾱ).σ (6)
Functions α i are then chosen such as P (x) = σ 0 . This leads to a set of differential equations

∂ x Re(α 0 ) = 1, ∂ x Re(α) + 2 h Re(α) × Im(α) = 0, ∂ x Re(α) 2 = 0 (7)
This is similar to a rotational equation of motion of a particle located at Re(α) around the angular vector Ω = (2/h) Im(α) on the sphere of constant radius Re(α) , where the dynamical time variable is played by x. The solution for α 0 is given by α

0 (x) = x + iw(x)
up to a constant, and where w(x) is any real function of x. We will restrict in the following this solution to α 0 (x) = x, which corresponds to the displacement of the oscillator in the harmonic case. We will see further below that the presence of w(x) is equivalent to add a gauge term e ih -1 x w(y)dy to the wavefunctions. For the other components, let us consider the solution α = (u(x), v(x), iω(x) + x c ) ⊤ where u, v, and ω are functions with real values. The constant x c gives the location of the center of the two oscillators at x = ±x c . The system of equations satisfying ( 7) is given by

u ′ (x) + 2 h ω(x)v(x) = 0, v ′ (x) - 2 h ω(x)u(x) = 0 (8)
whose solutions are chosen such that u(x) = g cos θ(x) and v(x) = g sin θ(x), with θ(x) = 2h -1 x ω(y)dy and g a constant (coupling) parameter. Therefore α = (g cos θ(x), g sin θ(x), iω(x)

+ x c ) ⊤ (9)
The matrix Q(x) is then expressed as

Q(x) =   x + x c + iω(x) ge -iθ(x) ge iθ(x) x -x c -iω(x)   (10)
and the classical Hamiltonian is given by

H c (x, p) =   (p + ω) 2 + (x + x c ) 2 + g 2 2ge -iθ (x -iω) 2ge iθ (x + iω) (p -ω) 2 + (x -x c ) 2 + g 2   (11)
The Fermi surface is given by two sheets, solutions of the equation det(H c -σ 0 E) = 0

E = p 2 + ω(x) 2 + x 2 + x 2 c + g 2 ± 2 (pω(x) + xx c ) 2 + g 2 (x 2 + ω(x) 2 ) ( 12
)
As example, we have plotted in figure 2 12) for parameters E = 3, ω(x) = 1, x c = 1, and g = 0.05. The red surface is given in equation ( 12) by the minus sign, and the blue small pocket by the plus sign.

FIG. 3: Fermi surface of equation (12) for parameters E = 5, ω(x) = x, x c = 1, and g = 0.5. There is no gap in this case as the two sheets intersect at x = 0. a gap controlled by g. In figure 3 is another example of section with ω(x) = x, x c = 1, and g = 0.5, at energy E = 5. In the latter case, the two surfaces intersect at x = 0 because ω(0) = 0 and the square root in equation ( 12) is zero. In general if ω(0) = 0 the two surfaces do not intersect and a gap appears.

We should notice that H c depends on h through the non-diagonal elements e ±iθ(x) which become highly oscillating in the classical limit when h → 0, but the equation for the Fermi surface ( 12) in this limit does not depend on h since the two exponentials in equation (11) cancel each other. The form of the matrix Q in equation (10) ensures that the commutator P (x) = σ 0 and that the eigenfunctions Ψ n = (φ n , φn ) ⊤ can be constructed exactly by applying D † successively on the ground state Ψ 0 :

Ψ n = D † n Ψ 0 / √ n!, with Landau level energies E n = 2hn.
The ground state energy is E 0 = 0, contrary to the usual harmonic oscillator for which the spectrum is E n = (2n + 1)h and the zero point energy E 0 = h.

The zero point energy depends on the term hQ ′ (x) in equation ( 2) which generates the interaction between bosons and fermions in the supersymmetry theories [15].

To solve for the eigenfunctions, we introduce the transformation

Q(x) = U † (x)Q 0 (x)U (x) + hU † (x)∂ x U (x), Q 0 (x) =   x + δ 0 0 x -δ   , U (x) = e i 2 ϕσ 2 e i 2 θ(x)σ 3 (13)
where U (x) is a unitary matrix, and (ϕ, δ ≥ 0) two parameters satisfying the following relations:

g = δ sin(ϕ), x c = δ cos(ϕ) (14)
The Hamiltonian can then be expressed as

H = U † H 0 U , with H 0 = 2hD † 0 D 0 , D 0 = (hσ 0 ∂ x + Q 0 (x))/ √ 2h.
The eigenfunctions are therefore related to the elementary pair of harmonic oscillators with Hermite polynomials

Ψ n =   φ n φn   = 1 √ n!2 n (πh) 1/4 U †   AH n x+δ √ h e -1 2 (x+δ) 2 /h BH n x-δ √ h e -1 2 (x-δ) 2 /h   (15)
where A and B are complex constants. U † takes explicitly the form

U † =   cos ϕ 2 e -i 2 θ(x) -sin ϕ 2 e -i 2 θ(x) sin ϕ 2 e i 2 θ(x) cos ϕ 2 e i 2 θ(x)   (16)
and the values of the components of Ψ n are simply given by

φ n (x) = e -i 2 θ(x) √ n!2 n (πh) 1/4 A cos ϕ 2 H n x+δ √ h e -1 2h (x+δ) 2 -B sin ϕ 2 H n x-δ √ h e -1 2h (x-δ) 2 , φn (x) = e i 2 θ(x) √ n!2 n (πh) 1/4 A sin ϕ 2 H n x+δ √ h e -1 2h (x+δ) 2 + B cos ϕ 2 H n x-δ √ h e -1 2h (x-δ) 2 (17)
We should notice that if we take into account the general function w(x) from the solutions of equation ( 7), the matrix Q 0 (x) would be instead

Q 0 (x) =   x + iw(x) + δ 0 0 x + iw(x) -δ   ( 18
)
and it is easy to show that the solutions (φ n , φn ) transform into (φ n , φn )e ih -1 x w(y)dy . We also define the currents

J n (x) = 1 2 h Im(φ * n (x)∂ x φ n (x)) and Jn (x) = 1 2 h Im( φ * n (x)∂ x φn (x)
). The sum of the two current satisfies

J n (x) + Jn (x) = 1 2 ω(x) | φn (x)| 2 -|φ n (x)| 2
The function ω(x) is physically equivalent to an internal current, proportional to the imbalance between the two orbitals occupancy. A and B are chosen such that for each orbital there are ρ 0 and ρ0 electrons in the ground state: dx|φ 0 | 2 = ρ 0 and dx| φ0 | 2 = ρ0 . We could also eventually choose ρ 0 and ρ0 as function of A and B. A and B satisfy the system of equations

|A| 2 cos 2 ϕ 2 + |B| 2 sin 2 ϕ 2 -cos ϕ 2 sin ϕ 2 (A * B + AB * )e -δ 2 h = ρ 0 , |A| 2 sin 2 ϕ 2 + |B| 2 cos 2 ϕ 2 + sin ϕ 2 cos ϕ 2 (A * B + AB * )e -δ 2 h = ρ0 (19)
We can parametrize coefficients A and B such that

A = ρ 0 + ρ0 cos(φ)e iθ A , B = ρ 0 + ρ0 sin(φ)e iθ B ( 20
)
where φ satisfies the following equation

cos(2φ) - g x c cos(θ A -θ B ) sin(2φ)e -δ 2 /h = ρ 0 -ρ0 ρ 0 + ρ0 δ x c (21)
We can choose θ A and θ B such that cos(θ A -θ B ) = 0 to remove the dependence of φ with h. The resulting equation for cos(2φ) has solutions provided that |ρ 0 -ρ0 |/(ρ 0 + ρ0 ) < x c /δ.

This excludes configurations such that ρ 0 = 1 and ρ0 = 0 for example, since x c /δ < 1 for g > 0, in which case the two components of the wavefunction are "intricated": each quasiparticle in one band has a non-zero component in the other band when g > 0, due to the tunnel effect induced by the unitary rotation in equation ( 13) when ϕ = 0, even in absence of magnetic field. The critical ratio between the two densities when cos(2φ) = ±1 is ρ 0 /ρ 0 = (δ ± x c )/(δ ∓ x c )). This means that for a density ρ 0 in the first orbital, there is a transfer of a small component proportional to g 2 ρ 0 to the second orbital when g ≪ 1. For example, when B = 0, the ratio between the two densities is simply

ρ0 ρ 0 = tan 2 ϕ 2 = g 2 x c + δ (22)
Unless g = 0, in which case ρ 0 = |A| 2 and ρ0 = 0, the transmitted part ρ0 is proportional to g 2 ρ 0 . When ρ 0 = ρ0 , we will consider in the following the real solutions of the system (20) with θ A = θ B = 0 and for which φ is explicitly given by tan(2φ) = x c g e

δ 2 h (23)
If ω(x) is constant, the total integrated current vanishes in this case. In the limit of small coupling g, we obtain |A| 2 ≃ ρ 0 , |B| 2 ≃ ρ0 , and φ ≃ π/4: the two functions φ 0 and φ0 decouple as each electron is localized in his own orbital.

A. Overlap function

As a first application, we consider the case in figure 2 for which ω(x) = 1 and define the overlap or transmission coefficient between the two wavefunctions (17) to evidence and study the oscillatory nature of the system as function of the inverse magnetic field. For this, we define the inner product or transfer function T mn = ∞ -∞ dxφ m (x) φn (x)dx between the Landau states m and n. After computing the product using integral formula for Hermite polynomials [START_REF] Prudnikov | Special Functions[END_REF], we obtain

T m≤n = m!2 n n!2 m (i √ h) m-n ge -1/h δ cos(φ) 2 e 2iδ/h -sin(φ) 2 e -2iδ/h L n-m m ( 2 h ) -e -δ 2 /h sin(φ) cos(φ) g x x + δ (1 + iδ) n-m - x x + δ g (1 -iδ) n-m L n-m m ( 2 h (1 + δ 2 )) , T m≥n = n!2 m m!2 n (i √ h) n-m ge -1/h δ cos(φ) 2 e 2iδ/h -sin(φ) 2 e -2iδ/h L m-n n ( 2 h ) -e -δ 2 /h sin(φ) cos(φ) g x x + δ (1 -iδ) m-n - x x + δ g (1 + iδ) m-n L m-n n ( 2 h (1 + δ 2 )) , T nn = ge -1/h δ cos(φ) 2 e 2iδ/h -sin(φ) 2 e -2iδ/h L n ( 2 h ) + x c e -δ 2 /h g sin(2φ)L n ( 2 h (1 + δ 2 )) (24)
where L n-m m are generalized Laguerre polynomials. We define the transmission function

T (h) at a given energy E T (h) = n≥0 δ n - E 2h T nn (h) (25)
In the numerical situation of figure 5, we will use a broadening of the Landau levels with a Gaussian function of width Γ ≪ 1 instead of the delta function δ

(x) → e -x 2 /2Γ 2 (2πΓ 2 ) -1/2 .
Oscillations as function of the inverse field B -1 are characterized by their frequencies F in the Fourier spectrum. For example, in the de Haas-van Alphen semi-classical theory, the Fourier spectrum of the magnetization oscillations are given by a set of discrete fundamental frequencies F and their harmonics that are proportional to the surface area A of all the closed cyclotronic orbits in the Fermi surface (x, p) at a given energy: they are equal to F = A /2πe. In term of inverse effective field variable h -1 , they are expressed instead as F = A /2π. In absence of coupling g = 0 and with ω constant, we can estimate the geometrical expressions of the surface areas A α enclosed by the blue curve and A β > A α in red on figure 2(b) which are given respectively by

A α = 2E arcsin E -ω 2 -x 2 c E -2 (E -ω 2 -x 2 c )(ω 2 + x 2 c ), A β = 2πE -A α = 2A 0 -A α (26)
where A 0 = πE is the area of the two circular cyclotronic orbits in the non interacting case, with radius √ E. In figure 5 is plotted the real part of the overlap T (h) as function of 1/h

with parameters E = 3, x c = 1, and g = 0.05. We observe that the amplitude of T (h) decreases with h in the semiclassical limit. The presence of two frequencies, one small and one large, is observed with beating effect in field regime considered. With these parameters and g = 0, we have the de Haas-van Alphen frequency F 0 = E/2 = 1.5 [START_REF]spectrum E n = 2hn gives the dominant frequency F 0 = A 0 /2π[END_REF], in addition to

F α ≃ 0.138, F β ≃ 2.
862. We also assume that the areas are still approximately given by these values when g is small. We then compute the Fourier amplitude of T (h) with respect

to h -1 A(F ) = ∞ 0 dh -1 T (h)e -2iπF/h = 1 F 0 n≥1 T nn (F 0 /n)e -2iπnF/F 0 (27)
In inset of figure 5, we have plotted |A(F )| which presents double peaks at regular locations we identify to be equal nF 0 + (F 0 ± F α )/2, with (F 0 -F α )/2 ≃ 0.68 and (F 0 + F α )/2 ≃ 0.82.

Each pair of peaks is separated by F 0 , and the beating frequency is F α . F α can be considered physically as the overlap area between the two orbits and therefore characterizes the function T nn (h) at small coupling g, whereas F 0 represents the periodicity of the energy density.

In the semiclassical regime h ≪ 1, T (h) can be found analytically with the saddle point method. Using the representation integral of the delta function, one may consider the following function relevant to compute T (h):

G(h, w) = e -w/h ∞ -∞ dx 2π n≥0 L n ( 2w h ) e ix(n-E/2h) (28)
with w > 0. Using the generating function for the Laguerre polynomials, we can perform the series and obtain formally [31]

G(h, w) = e -w/h ∞ -∞ dx 2π e -iϕ(x)/h 1 -e ix , ϕ(x) = 1 2 Ex - 2iw e -ix -1 (29)
We then apply the saddle point method to search for the solutions of ϕ ′ (x) = 0. We find that the solutions are given by the equation sin(x/2) 2 = w/E or x = ±2 arcsin( w/E) + 2πk, with k integer. Using the second derivative ϕ ′′ (x) = 1 2 w cos(x/2)/ sin(x/2) 3 , we obtain, after some algebra, the following behavior

G(h, w) ≃ h π 1 (w(E -w)) 1/4 sin 2π F (w) h + π 4 + γ w Ω(h), Ω(h) = n δ n - E 2h , tan γ w = w E -w (30)
The frequency F (w) is given by

2πF (w) = E arcsin w E + w(E -w) (31)
From this expression, T (h) is given in this limit by

T (h) ≃ g δ cos(φ) 2 e 2iδ/h -sin(φ) 2 e -2iδ/h G(h, 1) + x c δ sin(2φ)G(h, 1 + δ 2 ) Ω(h) (32)
The Fourier spectrum should therefore present a series of peaks at harmonics of F = F 0 due to the periodicity of Ω(h), and two beating frequencies F (1) and F (1 + δ 2 ) ≃ (F 0 -F α )/2, the latter being dominant in amplitude when g is small. The red curve in figure 5 represents the expression (32) which is in good agreement with the solutions given by equation (25) for the field range.

B. Charge transfer oscillations

The second example of application concerns the computation of the electron density or charge density for each orbital ρ n = dx|φ n | 2 and ρn = dx| φn | 2 , as function of the Landau level n. We can, as before, integrate over Hermite polynomials and obtain

ρ n = ρ 0 + g (ρ 0 + ρ0 ) 2δ sin(2φ) 1 -L n 2δ 2 h e -δ 2 /h ρn = ρ0 -g (ρ 0 + ρ0 ) 2δ sin(2φ) 1 -L n 2δ 2 h e -δ 2 /h (33)
The sum of the densities is conserved but there is a transfer of charge as function of index n. We can consider the transfer charge function (we set ρ 0 = ρ0 = 1 in the following)

∆ρ(h) = n≥0 δ n - E 2h (ρ n -ρn ) (34)
and study as before the semi-classical limit h ≪ 1, using function G(h, δ 2 ), equation [START_REF] Yoshioka | The Quantum Hall Effect[END_REF].

We obtain the expression for ∆ρ(h) in this limit

∆ρ(h) ≃ 2g δ sin(2φ) e -δ 2 /h - h π 1 √ δ(E -δ 2 ) 1/4 sin 2π F (δ 2 ) h + π 4 + γ δ 2 Ω(h) (35) ∆ρ (

III. GENERAL CASE

We generalize in this section the model of two oscillators by considering the following extension of the matrix Q(x) in equation (10) representing a system of 2(N + 1) interacting orbits, where N ≥ 0,

Q(x) =        x + (2N + 1)x c + iω(x) ge -iθ(x) 0 • • • ge iθ(x) x + (2N -1)x c -iω(x) ge iθ(x) 0 • • • • • • • • • 0 ge iθ(x) x -(2N + 1)x c -iω(x)        (36)
This extended matrix satisfies the commutation relation P (x) = 1 provided that θ(x) = 2h -1 x ω(y)dy as before. There are alternating internal current generated by ±ω(x) in each orbital. For N = 1, we have plotted in figure 7 the classical Fermi surface composed of a linear chain containing 4 orbits. We can compute the ground state by solving the differential equation DΨ 0 = 0. Noting the (2N + 2)-component vector Ψ 0 = (φ

(N ) 0 , • • • , φ (-N -1) 0
) ⊤ , the system of equations h∂ x Ψ 0 + Q(x)Ψ 0 = 0 can be reduced by the transformation Ψ 0 → Ψ0 where

Ψ0 = ( φ(N) 0 , • • • , φ(-N-1) 0 ) ⊤ , Ψ 0 = e -x 2 2h diag(e -i 2 θ(x) , e i 2 θ(x) , • • • , e i 2 θ(x) ) Ψ0
This simplifies the system of differential equations -h∂ x Ψ0 = Q Ψ0 , where Q is the constant matrix:

Q =        (2N + 1)x c g 0 • • • g (2N -1)x c g 0 • • • • • • • • • 0 g -(2N + 1)x c        (37)
We can diagonalize this matrix by searching for the eigenvectors

v λ k , with k = -N -1, • • • , N v λ k-1 + v λ k+1 + (2k + 1 -λ)x c g v λ k = 0 (38)
whose solutions are given by Bessel functions

v λ k = aJ µ-k (g/x c ) + bY µ-k (g/x c ) (39)
where we have defined λ = 2µ + 1. We impose as boundary conditions v λ N +1 = v λ -N -2 = 0. This leads to the equation for the discrete eigenvalues µ

J µ+N +2 (g/x c )Y µ-N -1 (g/x c ) = Y µ+N +2 (g/x c )J µ-N -1 (g/x c ) (40)
There are 2(N + 1) solutions λ k = 2µ k + 1 to this equation which is a polynomial in µ: In this paper we have shown that we can construct exact solutions of interacting bands in a magnetic field, from a factorized Hamiltonian and using a unitary transformation acting on independent band orbitals, see equation (13). The simple extension of the model leads to the matrix structure, equation (36), which is a Fermi surface composed of several coupled orbits and allows us to determine the eigenvectors. In particular for the simple case of two coupled orbits, we have evaluated the overlap function, equation (25), whose oscillating behavior contains the geometrical characteristics of the cyclotronic trajectories on the Fermi surface. We also considered the charge transfer between the two orbitals, and showed that it oscillates as function of the inverse field with a beating frequency F (δ 2 ), in addition to the dominant de Haas-van Alpen frequency of the Landau energy spectrum.

FIG. 1 :

 1 FIG. 1: Fermi surface of the organic conductor (BEDO-TTF) 5 [CsHg(SCN) 4 ] 2 , from [5, 23]. It is composed of a chain of small orbits α and magnetic breakdown orbits β. Electron trajectories are represented by incoming waves (a) that are scattered in (c) and transmitted in (b).

  FIG.2: Fermi surface of equation(12) for parameters E = 3, ω(x) = 1, x c = 1, and

FIG. 4 :

 4 FIG. 4: Overlap function modulus |T nn | as function of the inverse field h -1 for parameters g = 0.5, x c = 1, and for n = 5, 10, 15.

FIG. 5 :

 5 FIG. 5: Real part of the overlap T (h) as function of 1/h with parameters E = 3, g = 0.2, ω(x) = 1, and x c = 1. Broadening parameter Γ = 0.1. Red curve is the semiclassical expression (32). Inset: Fourier transform (FT) of Re(T (h). The first peak is estimated to be located at frequency F = 0.68 and the second at F = 0.82.

FIG. 6 :

 6 FIG. 6: Density fluctuations ∆ρ(h), equation (34), as function of inverse fields h -1 (black line) for the parameters ρ 0 = ρ0 = 1, E = 3, g = 0.2, and x c = 1. Broadening parameter Γ = 0.1. Red line is the semiclassical result (35). Inset: Fourier amplitude of the black curve. Fundamental frequency is located at F = 0.527, which corresponds to F (δ 2 ) (see text), and harmonics at nF 0 ± F (δ 2 ) with n integer.

FIG. 7 :

 7 FIG. 7: Fermi surface of equation (36) with N = 1, representing 4 coupled orbits with parameters E = 2, g = 0.1, ω(x) = 1/2, and x c = 1.

µ- 1 C l v λ l k e -(x+λ l xc) 2 2h- 1 C 2 2h

 1212 N > • • • > µ -N -1 (for small g, we have the approximation µ k ≃ k). x = -λ k represents the position of the k th oscillator on the chain. Then we can choose the following eigenvectors:v λ k = Y µ-N -1 (g/x c )J µ-k (g/x c ) -J µ-N -1 (g/x c )Y µ-k (g/x c ) Y µ-N -1 (g/x c )J µ-N (g/x c ) -J µ-N -1 (g/x c )Y µ-N (g/x c )(41)For example, for N = 1, equation (40) leads to the quartic equation for µ16µ 4 + 32µ 3 -4{3(g/x c ) 2 + 4}µ 2 -4{3(g/x c ) 2 + 8}µ + (g/x c ) 4 = 0 (42)It is then easy to show that the components of Ψ 0 can be expressed as φ The sign of the wavefunction phase θ(x) is alternating depending on the oscillator position, with a minus sign for the first oscillator located at x = -λ -N -1 .The excited states φ (k) n≥0 (x) with energy E n = 2hn are then constructed by identification to equation (17) with Hermite polynomials φ (k) n (x) = e ± i 2 θ(x) l v λ l k H n x+λ l xc √ h e -(x+λ l xc)

ACKNOWLEDGMENTS

This work was supported by the Brain Pool Program through the National Research

Foundation of Korea (NRF-2018H1D3A2065321).