Wouter J T Bos

Jean-Pierre Bertoglio

Louis Gostiaux

J.-P Bertoglio

Modeling the role of clusters and diffusion in the evolution of COVID-19 infections during lock-down

The dynamics of the spread of epidemics, such as the recent outbreak of the SARS-CoV-2 virus, is highly nonlinear and therefore difficult to predict. As time evolves in the present pandemic, it appears more and more clearly that a clustered dynamics is a key element of the description. This means that the disease rapidly evolves within spatially localized networks, that diffuse and eventually create new clusters. We improve upon the simplest possible compartmental model, the SIR model, by adding an additional compartment associated with the clustered individuals. This sophistication is compatible with more advanced compartmental models and allows, at the lowest level of complexity, to leverage the well-mixedness assumption. The so-obtained SBIR model takes into account the effect of inhomogeneity on epidemic spreading, and compares satisfactorily with results on the pandemic propagation in a number of European countries, during and immediately after lock-down. Especially, the decay exponent of the number of new cases after the first peak of the epidemic is captured without the need to vary the coefficients of the model with time. We show that this decay exponent is directly determined by the diffusion of the ensemble of clustered individuals and can be related to a global reproduction number, that overrides the classical, local reproduction number. 1 Lock-down and modeling of the COVID-19 epidemic Since the beginning of the COVID-19 epidemic in December 2019, many countries have chosen severe restrictions on mobility and transport. It has become a controversial ques-

tion in public debate whether these decisions were necessary and inevitable, and whether the lock-down was maintained for an unnecessary long period or, on the contrary, was relaxed too soon in an unsafe way. In order to evaluate retrospectively such decisions, and to anticipate similar measures in the future, an in-depth understanding of epidemic propagation is needed. In the present investigation we focus on the seemingly universal way in which the SARS-CoV-2 disease is controlled by lock-down.

Indeed, in a number of countries in which lock-down was applied, the number of new infections decays exponentially after the peak of the epidemic with a very similar timescale (see Fig. 4). In the present investigation we shed light on the role of the finite size of spatially localized clusters on this phenomenon, by improving upon a simple paradigm model, the classical SIR model, and by assessing the resulting SBIR model against realistic data.

The SIR model in its simplest form contains three compartments, corresponding to susceptible, infected and recovered or removed individuals, and the dynamics are governed by two model parameters γ (recovery rate) and β (infection rate). A wide variety of SIR modifications have seen the light since its original formulation [START_REF] Kermack | A contribution to the mathematical theory of epidemics. Proceedings of the royal society of london Series A[END_REF]. Many of them propose additional compartments (see for instance Massonis, Banga, and Villaverde [START_REF] Massonis | Structural identifiability and observability of compartmental models of the covid-19 pandemic[END_REF] for an analysis of a great number of compartmental models), nonlinear incidence rates [START_REF] Capasso | A generalization of the kermack-mckendrick deterministic epidemic model[END_REF], and non-constant model coefficients [START_REF] Bacaër | Un modèle mathématique des débuts de l'épidémie de coronavirus en france[END_REF][START_REF] Derrode | Piecewise estimation of R0 by a simple SEIR model. Application to COVID-19[END_REF]. An overview of these models is beyond the scope of the present work and we refer to the review articles on the subject [START_REF] Hethcote | The mathematics of infectious diseases[END_REF][START_REF] Brauer | Mathematical models in population biology and epidemiology[END_REF].

In a large number of these epidemiologic studies using compartmental models, the dynamics are extremely sensitive on the temporal variations of the model parameters, and saturation of the epidemic at a realistic level is only obtained if the main parameters, such as the infection rate (or reproduction number R 0 = β /γ) are changed during the analysis. Sophisticated approaches have been developed to esti-mate the temporal evolution of β by assuming for instance a prescribed shape of the time-evolution of the reproduction number [START_REF] Linka | The reproduction number of covid-19 and its correlation with public health interventions[END_REF] or by adjusting the value of β instantaneaously to the observations [START_REF] Derrode | Piecewise estimation of R0 by a simple SEIR model. Application to COVID-19[END_REF]. We show that even using constant coefficients, saturation of the epidemic is obtained when the local nature of interactions within clusters is taken into account.

Indeed, perhaps the largest known flaw of the SIR model is the assumption of well mixedness. By assimilating all susceptible individuals into one single compartment, it is tacitly assumed that all individuals are in contact with all others, and the spatial structure of the epidemic is thereby ignored. An approach to overcome this weakness is to use either networks, self-propelled interacting agents, nested or coupled local SIR-type models [START_REF] Mideo | Linking within-and between-host dynamics in the evolutionary epidemiology of infectious diseases[END_REF] or diffusion equations [START_REF] Murray | On the spatial spread of rabies among foxes[END_REF][START_REF] Kuperman | Small world effect in an epidemiological model[END_REF][START_REF] Martcheva | An introduction to mathematical epidemiology[END_REF]. Such modifications change the level of complexity of the model considerably requiring either the solution of a two-dimensional reaction-diffusion problem, tracking a large number of propelled particles, or the solution of integrodifferential equations. We refer to Rahmandad and Sterman [START_REF] Rahmandad | Heterogeneity and network structure in the dynamics of diffusion: Comparing agent-based and differential equation models[END_REF] for a discussion of the differences between agent-based and compartmental models, and to Zohdi [START_REF] Zohdi | An agent-based computational framework for simulation of global pandemic and social response on planet x[END_REF] for an agentbased model of epidemic spreading applied to a planetary domain.

The model we introduce in the present work is much simpler and uses a two-scale approach: either individuals are far enough away from the individuals not to be in contact, or they pertain to a cluster of individuals in contact with the infected. This approach allows to take into account spatial structuring at the lowest level of complexity, i.e., by adding one single compartment. We sub-divide the group of susceptibles into B, those in a so-called "blob", who are or have been in contact with the infected, and S, those who are (spatially) distant from the infected. We will show, by comparison with data for a variety of European countries, that this model reproduces the exponential decay of the number of newly infected after its peak-value, with an exponent determined by the spatial diffusion coefficient of the blob.

An important simplification of reality which we will assume throughout this investigation is the constancy of the model-coefficients. All countries introduced different measures which change the effective transmission rate by physical effects such as quarantine-measures [START_REF] Berger | An seir infectious disease model with testing and conditional quarantine[END_REF], or by changing the awareness of the individuals [START_REF] Wu | The impact of awareness on epidemic spreading in networks[END_REF]. The present investigation aims primarily at the assessment of a novel approach to take into account inhomogeneity. The lock-down period for a number of similar occidental countries seems the correct starting point, since despite differences in detailed measures, the parameters governing the dynamics of the epidemic are expected to become constant during lock-down. The universal features captured by the model can therefore be disentangled from other effects, and the results can be compared to the data within the limits of the quality of the measurements. In SIR and CSIR three compartments interact, with the only difference in CSIR that those three compartments do not contain the entire population, since part of the population is supposed to be protected by cross-immunity. In the SBIR model, the susceptible part of the population is subdivided into individuals B inside the blob, and S, outside the blob, and the four compartments interact.

In the following we will first present the SBIR model and its difference with the classical SIR approach (section 2). We will analytically determine the exponential decay time-scale for the SIR model with and without cross-immunity, as well as for the SBIR model. We show, by comparison with data from the Johns Hopkins University data-base [START_REF] Dong | An interactive webbased dashboard to track covid-19 in real time[END_REF], that the SBIR model allows to simultaneously fit the data in the increasing and decreasing phase of the epidemic without the need of adjusting the model-parameters during the epidemic (section 3). We further estimate the reproduction number and show that it reaches a stable minimal value during lockdown. The appendices contain the mathematical procedure to analytically estimate the decay exponents of the different models and the initial and final reproduction numbers.

A cluster-based compartmental model

Spatially resolved epidemiological models are elegant and physically better justified than the SIR model but have the disadvantage that they are harder to analyze analytically, which makes insights on the effect of the different parameters more complicated. Indeed the complexity increases in general the descriptive capacities of a model, but decreases the insights which can be obtained.

We propose here a new model, that does only introduce one new compartment, representing the individuals who are or have been spatially close to the infected people. Thereby our model remains based on global parameters and does I 0 e λ t I 1 e -λ t not introduce the need to solve a two-dimensional diffusion problem. This simplicity allows the analytical derivation of the long-time behavior and thereby elucidates the role of local cluster-saturation on global decay-rates. We will in this section first recall the SIR model, give an analytical estimate of the decay-exponent in the final phase of the epidemic, and show the influence of cross-immunity in a simple CSIR model. After that we introduce the SBIR model and discuss its physical and some of its mathematical properties.

The decay exponent in SIR and cross-immunity

The simplest model for macroscopic modeling of an epidemy is the so-called SIR model [START_REF] Kermack | A contribution to the mathematical theory of epidemics. Proceedings of the royal society of london Series A[END_REF][START_REF] Diekmann | The legacy of kermack and mckendrick[END_REF]. The model describes a population divided in three categories, Susceptible S, Infected I and Removed or Recovered R (Fig. 1, top). The total population N = S + I + R is supposed constant on the timescale of the epidemic. The evolution of, and the interaction between the three compartments is given by a set of three coupled ordinary differential equations,

dS dt = -β IS N dI dt = β IS N -γI dR dt = γI. (1)
In the following it is convenient to normalize these equations by N, so that R, S, I all take values between one and zero and S + R + I = N = 1. Initially R(0) = 0, I(0) = I 0 1 and S(0) = 1 -I(0) ≈ 1. The dynamics of the SIR model are characterized by an exponential growth with exponent λ + = βγ, upto a time τ where saturation of the infected is obtained, by collective, or herd immunity, obtained when S < γ/β . Note that in the present epidemic of SARS-CoV-2, we remain very far below such collective immunity at the time of the first peak. Long enough after this t = τ, exponential decay is observed until all infected have disappeared. The behavior is sketched in Fig. 2. We show in appendix A a simple and straightforward way to determine the approximate exponential growth and decay rate of the epidemic. This analysis, assuming continuity of two exponentials (as in Fig. 2), yields the decay-exponent

λ -= γλ + β ≡ γλ + λ + + γ , (2)
which gives a direct relation between the exponential during the initial, growing phase of the epidemic, and the decaying phase. We note that an exact analytical solution of the system exists [START_REF] Harko | Exact analytical solutions of the susceptible-infected-recovered (sir) epidemic model and of the sir model with equal death and birth rates[END_REF], but the extension of the method to more complicated systems might not be straightforward. The present approach yields less accurate but readily interpretable results, allowing direct insights in the effect of the parameters on the evolution of the epidemic. Furthermore, the peak-value of the number of infections is given by

I(τ) = λ + β 2 , (3)
and at this time τ, the rate of infection P = β SI is close to its maximum and is given by (see appendix A)

P(τ) = γ λ + β 2 . ( 4
)
An important question is whether cross-immunity could explain the saturation of an epidemic such as the ongoing one at the observed levels. We consider thereto a population where a part of the individuals has attained cross-immunity -the precise source of it, past contact with similar viruses or any other hidden mechanism being beyond the scope of this study. This part of the individuals is called C (Fig. 1, center). To obtain the CSIR equations, the SIR-equations are now only changed on the level of the global community size, where we replace

N = C + S + I + R = 1 → N * = S + I + R = 1 -C. (5)
The growth-exponent λ + is not changed by this modification for given γ, β , but the peak-values of I and P are multiplied by a factor (1 -C) as explained in appendix B. For instance

P CSIR (τ) = γ(1 -C) λ + β 2 = (1 -C)P SIR (τ). (6)
In the long-time however, the decay-exponent λ -is the same for the SIR and the CSIR model.

SBIR: SIR in a blob

We consider a cluster-based compartmental model, where only a sub-ensemble of individuals N * < N will interact with infected individuals. Via social interaction, a set of clusters diffusing around infected individuals will emerge, this diffusion process being distinct from virus transmission, and having its own dynamics. We define the number of persons belonging to one of these diffusing clusters as N * = B+I +R ≡ 1 -S, where I and R have the same interpretation as before, i.e., infected and recovered/removed, respectively, and B is referred to as the blob, namely the people who joined a diffusing cluster through social interaction and who are not yet infected or recovered. These individuals can be seen as local susceptibles, and in the following the B, I and R compartments are considered globally, independently of the number of clusters they correspond to. The persons in S are susceptible individuals, that need first to join the blob before being possibly infected. The 4 compartments, as illustrated in Fig. 1, are now related by the following system of ordinary differential equations,

dS dt = -κ SI N + ζ SB N (7)
dB dt = κ SI N -ζ SB N -β BI N * (8)
dI dt = β BI N * -γI (9)
dR dt = γI. ( 10
)
If the growth exponent λ + = βγ has the same expression as in the SIR model, the final decay exponent, if β > κ, is given by (see appendix C),

λ -= γ -κ (11)
which is different from the expressions for SIR and CSIR, and now depends on the spatial diffusion-coefficient of the blob related to κ.

In the limit of very fast diffusion, B will approach S, and the SIR-dynamics is observed subsequently. Furthermore, mathematical robustness of the formulation is ensured, since none of the quantities S, B, I, R can become negative.

The physics behind SBIR

Before integrating the SBIR model, it is insightful to show what the physical picture behind the model is. This will help to interpret the results and to guide the choice of the parameters. The parameters of the model are γ, β , κ and ζ . Whereas γ is an intrinsic parameter characterizing the disease, β , κ and ζ depend on physical, biological and psychological influences. The β parameter, representing the infection rate is already present in the classical SIR model. It will be obtained here from an estimation of λ + in the early stage of the epidemics, which is known to be barely affected by the spatial structure of the population [START_REF] Trapman | Inferring R 0 in emerging epidemics-the effect of common population structure is small[END_REF].

In order to understand the dynamics of the blob, let us consider the coupling between the S and B compartments, and first neglect the term proportional to ζ . We have then

dN * dt = κ SI N , (12)
If we assume S/N ≈ 1 (as expected in early stages) and that I varies slowly, we have

N * ≈ κIt. ( 13
)
When a constant population density σ P is introduced, the surface of the blob N * /σ P corresponds to the increase of a surface explored by a random walk. Indeed, the separation of a random walker from its initial position is

r(t) ≡ (x(t) -x(0)) 2 ∼ √ Dt, ( 14
)
where D is the diffusion coefficient, and the surface explored is proportional to the square of this length. The interaction between the two compartments is therefore governed by a Brownian-like diffusion process, with a diffusion coefficient D = κI/σ P . In most epidemics, I rapidly grows so that the time-varying effective diffusion coefficient increases as well, and super-diffusive growth of the blob is observed. We can ask how this diffusion is related to the movement of individuals. In the kinetic theory of gases, the macroscopic diffusion coefficient is related to microscopic motion by

D m ∼ ul, (15)
where in this expression u is the velocity magnitude of the gas-molecules and l the mean-free-path. In our description, where the blob diffuses with a typical timescale κ -1 , the effective mean free path of individuals in contact with infected is proportional to the number of infected I, assuming that their velocity, which is related to the displacement of individuals, is constant (and reduced) during lock-down. Even though this is obviously a very rough approximation, we could use this idea to estimate the change of diffusion after lock-down using smartphone mobility data (as for instance [1, 2]). In the present investigation we will use constant model-coefficients and show that this approximation, during confinement, gives good agreement with data. The influence of non-constant coefficients is left for future investigation.

We still have the parameter ζ to discuss. The evolution equation of S can be written as

dS dt = -κ SI N 1 -Γ -1 B I , (16)
where The ratio B/I indicates how many individuals in the cluster are typically interacting with one infected. Indeed, if we set ζ to zero, there is no bound on the blob-size, and a blob of one-million individuals in contact with a single infected case would not make sense physically. This is precisely the flaw in the SIR model (the well-mixedness hypothesis) that we want to improve upon in the SBIR model. The combination of the κ and ζ term will allow the system to relax to an equilibrium value Γ of the typical infection density per cluster. In the present investigation we will use Γ not only to set the ratio of κ and ζ , but also to set the initial condition for the ratio B/I.

Γ = κ/ζ . (17)

Dependence of epidemic spreading on the diffusion rate

We have thus given physical interpretations of the model parameters κ and ζ . Before showing how the SBIR model compares to realistic data for constant coefficients, we will determine the long-time evolution of the system for well defined initial conditions. The values we used are

β = 3γ, ζ = κ/10, γ = 0.1, B(0) = R(0) = 0, I(0) = 1 • 10 -5
. The temporal evolution of the different quantities is documented in Appendix D. We have observed that at very long times, the only stable solution of the system is I = 0, B = 0, S = 1 -R. The final state of the population is shown in Fig. 3, where the values of S(∞) and R(∞) are plotted as a function of κ/γ. Values of R(∞) can range from R(∞) = I(0) for κ = 0 to R(∞) = O(1) for κ 1, where the SIR limit is attained. It is clear from this plot that the final state of the SBIRpopulation is very importantly determined by the diffusion rate. This allows to model the evolution of the COVID-19 epidemic during its first months, without the need to change the model coefficients during a simulation.

Comparison with data on confirmed infections

In this section we compare the results of our model to recent data on the epidemic. The number of daily confirmed new infected cases is obtained from the Johns Hopkins University data-base [START_REF] Dong | An interactive webbased dashboard to track covid-19 in real time[END_REF]. This quantity is related to the source term in equation ( 9),

P SBIR ∆t = β BI N * ∆t. (18)
With N * = B + I + R and ∆t = 1 day. We have chosen to compare to several European countries where the lock-down is applied and Sweden, a country who has only applied "mild" restrictions on its population. Our procedure is the following. We use for the parameter γ = 0.1 and determine in the data the time t = 0 where the prevalence of daily confirmed new infected cases exceeds P(0) = 1 • 10 -6 . The value γ -1 determines the typical time-interval during which an individual is infectious and can thereby transmit the disease to other individuals. The exact length of this time interval can vary between 3 and 18 days [START_REF] Linka | The reproduction number of covid-19 and its correlation with public health interventions[END_REF], and our choice of 10 days is in agreement with recent investigations [START_REF] Bar-On | A quantitative compendium of covid-19 epidemiology[END_REF][START_REF] Derrode | Piecewise estimation of R0 by a simple SEIR model. Application to COVID-19[END_REF][START_REF] Tian | An investigation of transmission control measures during the first 50 days of the covid-19 epidemic in china[END_REF]. The value γ = 0.1 corresponds to the case where from a group of isolated infected individuals, 50% will recover after approximately 7 days. Choosing a different value for this parameter would not dramatically change the results, but would directly change the choice of β used in the remainder. This latter parameter is determined as follows. During the early growth phase of the epidemic, the data is fitted by an exponential, thereby determining λ + = βγ, and thus also β . We use I(0) = P(0)/β and the value of B(0) is adjusted to match the height of the peak P(τ). The value of the cluster diffusion κ is varied to get a best fit on the decay phase. We note that in the phase of the epidemic we focus on, the influence of the parameter ζ on the results is weak compared to the influence of the choice of κ. It is chosen so that the initial ratio B(0)/I(0) is equal to κ/ζ . Finally, we impose R(0) = 0.

The results of the integration of the SBIR model are compared in Fig. 4 to the data, and the parameters are indicated on each plot. The general agreement is rather encouraging. It is observed that for all countries the value of β is of order 0.3. The value of the diffusivity varies within the subset of considered countries in the range κ ∈ [4.5 • 10 -2 , 7.5 • 10 -2 ] and the ζ parameter is typically 2 to 3 orders of magnitude smaller. For all countries a revival of the epidemic starts some time after the lock-down period, which is not reproduced by the model in its present form, and the reason for this is discussed in the conclusion. The lock-down period is indicated in Fig. 4 by a shaded area. We have used an objective definition for the dates of this lock-down period using mobility data from Google [1], given in appendix E.

Also shown in the figure for comparison are the CSIR results. We have again fitted the initial growth exponent, and we have adjusted the value of the cross-immunity C to reproduce the peak of infection. To obtain a correct fit we need to use values for the part C of individuals having cross-immunity in the interval C ∈ [0.997, 0.999]. The cross-immunity needed to attain a saturation of the number of newly infected is therefore unrealistic. Furthermore, the decay exponent λ -in the CSIR model λ - SIR ≈ 0.07 is too large, which largely underestimates the duration of the epidemic. The SBIR model, with the introduction of the diffusivity parameter κ, gives fairly good agreement in the decay regime with λ - SBIR ≈ 0.04. This ability to correctly represent and explain the fast increase and the slow decrease of the epidemic is the major advantage of the SBIR model.

We also compare fhe SBIR model to a country where lock-down restrictions were very light, Sweden. It is observed that the model also describes the behaviour of the data in this country, but with a value of κ more than twice larger than in countries with severe lock-down restrictions. The spatial diffusion of the blob in which the ensemble of cluster evolves is in Sweden thus, according to our model, twice as large as in the other considered countries.

Estimating the reproduction number

The evolution of the number of new infected individuals can also be quantified by the reproduction number R t . We recall that the initial value of the reproduction number within the SIR or SBIR models is defined by R 0 = β /γ.

We can rewrite the I-equation of the SBIR model [Eq. ( 9)] as 1

I dI dt = γ (R t -1) (19)
where the instantaneous reproduction number is given by the expression

R t = β γ B(t) N * (t) . (20)
Since both B and N * evolve in time, this quantity will also vary. Initially, since R, I B, we have R t ≈ R 0 . In Fig. 5, we show the evolution of R t obtained from the SBIR-integration. To compare with data we use the local logarithmic derivative of the newly infected cases P to estimate R t

R t = 1 + 1 γP dP dt . (21)
The use of P instead of I in Eq. ( 19) is justified as long as B/N * is a slowly varying quantity. In order to reduce the noise in the value of the derivative, data is preprocessed using a two-week moving average. We see that the values of R t from SBIR do agree relatively well with the data. Around t = 0 the scatter is very large, which is in part caused by the quality of the data which depends among others on the number of tests. During lockdown, the value of R t drops from its initial value which varies for the considered countries around 3 -4, to a value around 0.5. The rapid decay of R t is present in both the data and the model estimates. The shape of the temporal evolution is close to a hyperbolic tan shape, as used in a recent modeling investigation of the COVID-19 epidemic [START_REF] Linka | The reproduction number of covid-19 and its correlation with public health interventions[END_REF]. However, there seems to be a time-delay between the two values. This might stem from the modeling approach, as no delay is introduced here between exposition and infection. The SBIR model, which is a direct improvement of the SIR model could possibly reproduce this when an Ecompartment is added as in classical SEIR models.

Finally, the relation between the value of the decay exponent λ -observed during lockdown in Fig. 4 and R can be derived when a pure exponential is assumed,

λ -= γ(1 -R t ). (22)
Therefore, using Eq. ( 11), we have for SBIR when the blobsize is determined by the diffusion of B,

R t = κ γ ≡ R κ . (23)
We could call R κ the macroscopic or global reproduction number, whereas the classical quantity R 0 = β /γ is the microscopic or local reproduction number. When R 0 > R κ , the number of infected individuals in the blob will increase faster than new susceptibles can enter the blob, so that diffusion will sooner or later become the limiting factor in the spreading of the epidemic. During lock-down it is this factor which will determine the decay of the number of new infections.

Discussion

We have shown that the introduction of an additional compartment in the SIR model allows to take into account the presence of clusters in the evolution of an epidemic. We have shown that the SBIR model reproduces the decay-exponent of new infections during lockdown, without the need to vary in time the model parameters over a duration exceeding 100 days.

The results in the present paper have to be understood as an illustration of the potential of the approach. A more comprehensive study should take into account the increase in the number of tests and its consequences on the number of new infected cases [START_REF] Sauter | Testing informed sir based epidemiological model for covid-19 in luxembourg[END_REF]. Also, we have deliberately chosen in the present work to consider the simplest model, SIR, and to add the effect of clustering upon this model. This allows to assess most easily the refinement. However, the nature of the present subdivision of compartments is so simple that most existing compartmental models (such as the SEIR model), can be modified in a similar manner. Effects such as vaccination, a finite incubation period, structuring of the compartments by age, and incomplete immunity could be taken into account in this manner.

The SBIR model shows that the close-to-universal decay exponent of the number of newly infected during lockdown, as observed in the data, can be linked to the cluster-diffusion. We also show that, during lockdown, a highly infectuous disease can be characterized by a global reproduction number R κ = κ/γ, which overrides the classical local reproduction number R 0 = β /γ.

We also show that the SBIR model in its present form reproduces observations in countries during lock-down, but no longer reproduces observations long time after restrictions are weakened. We trace back the origin of this different behavior to two reasons. The main reason is that we have assumed a constant diffusion coefficient κ during lock-down and we have not increased this coefficient at the end of the lock-down period. Using realistic data on mobility from mobile network providers (as used here to determine the period of lockdown) to modulate κ, would allow the blob to increase more rapidly after lockdown, and to allow the epidemic to rapidly generate new clusters. Another reason is the fact that in the SBIR model, as formulated here, all the recovered individuals R are counted in the blob-population N * . This is plausible in countries during lock-down, where clusters are fixed in space, hence will saturate and suffocate when R increases. In countries where the clusters not only expand, but also relocate, the epicenter of the epidemic will move spatially and leave behind regions where a large part of the population has been touched by the epidemic, but where the disease is less active. The way in which we can model these two distinct phenomena will be discussed and assessed in a future investigation.

A Analytical determination of the decay exponent in SIR

The decay exponent in the SIR model can be approximated analytically. We define t = τ the time where the number of infected I peaks (dI/dt = 0). We observe exponential increase and decay before and after this time, respectively (see Fig. 2),

I(t τ) = I 0 exp(λ + t), I(t τ) = I 1 exp(-λ -t). (24)
For short times S ≈ 1 and for very long times, I R, S. Therefore

λ + = β -γ, ( 25
) -λ -= β S ∞ -γ (26) = λ + -β R ∞ . (27)
In order to determine λ -, we need therefore the total number of recovered, R ∞ , given by

R ∞ = γ ∞ 0 I(t)dt. (28)
Since we assume exponential growth for short time, and exponential decay for long times, respectively, continuity at t = τ allows to write,

R ∞ = γ τ 0 I(t)dt + γ ∞ τ I(t)dt (29)
= γ 1

λ + + 1 λ -I(τ), (30)
which expresses R ∞ as a function of a new unknown, I(τ), the peakvalue of I. Considering the I-equation,

d t I = (Sβ -γ) I, (31)
the term in brackets vanishes for S = γ/β . We know that S = 1 -I -R, we also have that in the initial exponential phase

R = γ λ + I (32)
and therefore

γ/β = 1 -(1 + γ/λ + )I(τ) (33) = 1 -(β /λ + )I(τ). (34)
so that

I(τ) = λ + β 2 (35)
Which yields an expression for I(τ). We have the relations ( 27), (30), (34) for the unknown quantities, I τ , R ∞ and λ -, and solving this system yields

λ -= γλ + β = γ(β -γ) β . (36)
Note that for typical values, for instance γ = 0.2; β = 0.4, this yields λ -= 0.1. Numerical integration yields λ -= 0.119. This error of 20% is due to the assumptions of the two exponentials near the peak, where the rounding of the integrals significantly contributes to the error. The quantity that we will consider in the comparison with data is the daily confirmed new infected P∆t, which is associated with the source term in the I-equation, In Fig. 6 (c) we display the temporal behavior of the number of infected individuals for different values of κ. It is observed that for κ < γ, the number of infected decreases almost immediately: the epidemic has no time to spread, since the infected individuals remain confined locally around the source of the epidemic outbreak. On the contrary, when diffusion is strong, the epidemic will spread and rapidly all susceptible individuals will be part of the blob. The system behaves in this limit of strong diffusion as a well-mixed system, and the classic SIRbehavior is observed for the highest values of κ/γ considered here.

P = β SI, (37)
Finally, Fig. 6 (d) shows all recovered and removed individuals as a function of time. At the end of the evaluated time-interval the value of R is a continuously increasing function, ranging from R ≈ I(0) for κ ↓ 0 to a value of order (but somewhat below) unity, corresponding to the asymptotic SIR value for κ → ∞.

E Use of Google mobility-data to determine lock-down dates

We use the mobility data provided by Google [1] to determine the dates of the lock-down in the different countries. Fig. 7 shows the variations of mobility compared to a reference-value measured mid-january for the six countries considered in this work. The data used here corresponds to the retail and recreation data. Several other quantities are reported by Google, but most reproduce the same trend. Since we only use this data to illustrate the duration of the lockdown and do not use this to change the model-coefficients, the precise choice does not seem very important.

We have smoothened the data to get rid of high-frequency variations. We call the smoothened variable (t). The general trend shows that during the epidemic the mobility drops to a value depending on the restrictions in the different countries. An approximate plateau value is observed. We take for each country the minimum value of this plateau as the reference mobility-level during lock-down, indicated by  0 . The beginning of lock-down is then determined for the time when for the first time (t) < 1.15  0 and the end, the first later date where (t) > 1.5  0 , as can be seen on Fig. 7. These begin and end dates are used to indicate the lock-down periods as a shaded region in Fig. 4.

We note that this approach could allow to modulate the model coefficients. Since this will hinder the analysis of the present approach and since the assessment of the SBIR model is the main goal of the present work, this modulation is not implemented here. However, if a second wave of infections is to be reproduced, changing the model coefficients, or the dynamics of the R-individuals will become com-pulsory. Realistic mobility data can allow to get a handel on the way in which to modify the coefficients.

Fig. 1

 1 Fig.1Diagrammatic representation of three compartmental models, SIR, CSIR and SBIR. In SIR and CSIR three compartments interact, with the only difference in CSIR that those three compartments do not contain the entire population, since part of the population is supposed to be protected by cross-immunity. In the SBIR model, the susceptible part of the population is subdivided into individuals B inside the blob, and S, outside the blob, and the four compartments interact.

Fig. 2

 2 Fig. 2 Schematic representation of the evolution of the number of infected individuals in SIR and SBIR under the assumption of constant model-parameters. Time τ represents the peak-value of the number of infected.

Fig. 3

 3 Fig.3Influence of the diffusion rate κ on the final number of susceptible individuals S(∞) = 1 -R(∞) and recovered or removed individuals R(∞) for given constant coefficients.

Fig. 4

 4 Fig.4Comparison of SBIR and CSIR results, with data from JHU[START_REF] Dong | An interactive webbased dashboard to track covid-19 in real time[END_REF] for 8 European countries. The data represent the reported daily new infections. The shaded region corresponds to the lock-down period in the different countries.

Fig. 5

 5 Fig.5Comparison of the SBIR results for R t for 8 European countries. The shaded region is the superposition in transparency of the lock-down periods in the considered countries.

Fig. 6

 6 Fig. 6 Temporal evolution of (a) susceptible individuals S, (b) individuals in the blob B, (c) the number of infected individuals I, (d) recovered/removed individuals R. The different colors indicate different values of the diffusion rate in the range κ ∈ [0, 10γ].

Fig. 7

 7 Fig. 7 Estimation of the lock-down dates from Google mobility reports. Thick lines represent the lock-down period determined for each country.

Acknowledgments

The authors thank Samuel Alizon (MIGEVEC -CNRS) and Loïc Méès (LMFA -CNRS) for stimulating discussions, as well as Benjamin Pillot (UMR Espace-Dev 228) for developing the Pycovid-19 package [20] used to access data of the JHU-database and for useful advices. We also thank Pascal Roy (LBBE UMR 5558) and Vitaly Volpert (ICJ UMR 5208) for their helpful advices.

Data Availability Statements

The data that support the findings of this study are openly available from the Johns Hopkins University data-base [START_REF] Dong | An interactive webbased dashboard to track covid-19 in real time[END_REF].

Mobility data used in appendix E are openly available from [1].

where in our system, where the time-unity is days, ∆t = 1.

Interesting relations, when substituting this value for λ -in the equations are

B Exponent in SIR with Cross immunity (CSIR)

Let us now consider the case of a partially immune population,

where I(0)/S(0) and R(0)/S(0) are small (as in the previous case), so that S(0) ≈ 1 -C. The I-equation reads

For short times, we have thus, since S ≈ 1 -C,

as for SIR. The growth of I halts at time τ when dI/dt = 0, so that we find for S(τ)

We find then for the different quantities, following the same reasoning as in the previous section,

and for the exponential decay-rate,

as for SIR.

What we learn from this is that, cross-immunity or not, if we measure λ + , initially, λ -is not influenced for a given γ. Only the height of the peak of the number of infected I, and of the number of daily confirmed cases P will be multiplied by 1 -C.

C Exponent in SBIR

Let us now consider the SBIR model. The I equation of SBIR reads,

Let us consider the situation where the number of infected individuals grows faster that the number of individuals in the blob. The growth of the number of infections is then bounded by the source-term in the blob-evolution equation, and β IB/N * ≈ κIS/N. In this saturated phase, the I-equation becomes then,

If, as during the first wave of the COVID-19 epidemic, the global infection rate is still low and S/N ≈ 1, it is immediately found that

so that exponential decay is observed with a decay exponent λ -= γκ for κ < γ, or a slow increase if β < κ < γ. The reproduction number is in this phase therefore determined by the spatial spreading, so that the value of R t evolves from its local, microscopic value R 0 = β /γ to its diffusion determined value R κ = κ/γ,

Let us evaluate when this is the case. Let us rewrite the I-equation as

We will again assume exponentials as in Fig. 2. At short times, since B I + R, we have the same increasing exponential as in SIR and CSIR, with λ + = βγ. This increase will continue as long as I + R B. Since R is in the exponential phase proportional to I (as for SIR), we have

Let us therefore focus on the quantity I/B to see whether this quantity will be negligible for long times,

Considering the initial phase, where I/B 1, and estimating the order of magnitude of the different terms, it is observed that the leading order contributions of this equation are the last two terms in brackets, yielding,

Within a finite time, the blob will thus contain a considerable number of infected and recovered individuals, so that there is no more room for additional infected.

D Numerical integration of the model

We present here an analysis of the behavior of the model by numerical integration of the ODEs ( 7)-( 10) constituting the SBIR system. Since the principal difference of SBIR with respect to SIR is the diffusion from the S to the B compartment, we will assess the sensitivity of the model on a change of the diffusion coefficient. For this, we fix the parameters as follows β = 3γ, ζ = κ/10, γ = 0.1, B(0) = R(0) = 0, I(0) = 1 • 10 -5 and S(0) = 1 -I(0) ≈ 1. Using these parameters we vary the value of the diffusion rate in the interval κ ∈ [0, 10γ]. In Fig. 6 (a) we show the susceptible individuals as a function of time for different values of κ/γ. As expected, for small values of κ, S remains close to unity, whereas for large values the S-compartment is emptied almost completely. For long times the S value is slowly increasing, an effect associated with the flux of individuals from the B to the S compartment when the number of infected becomes very small. We have observed that this latter effect can be very slow for certain sets of parameters.

Fig. 6 (b) shows the temporal behavior of the individuals contained in the blob. For small values of the diffusion, B reaches rapidly a local maximum, and the number of individuals in the blob does not exceed the value 10 -5 . The epidemic is therefore rapidly limited in size. The value of the local maximum increases when κ increases and the blob will contain almost all susceptible individuals after around 100 days for the largest values of the diffusivity. For these values, the SIR behavior is recovered.


A Pires Da Cruz
email: antonio.pires-da-cruz@ifp.fr

C Pera

J Anderlohr

R Bounaceur

F Battin-Leclerc

A complex chemical kinetic mechanism for the oxidation of gasoline surrogate fuels: n-heptane, iso-octane and toluene -Mechanism development and validation

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

Introduction

Fundamental research in kinetic modeling of mixtures representative of commercial engine fuels (gasoline, diesel, jet...) has been strongly increasing in the last decade. Examples of such efforts can be found in references [1]- [START_REF] Naik | Society of Automotive Engineers[END_REF]. Motivation for this research has come from an increasing demand on understanding and modeling engine combustion towards a reduction in fuel consumption, CO 2 and pollutant emissions. However, commercial engine fuels may contain hundreds of different molecules and their composition is not completely well known. For modeling purposes, it is therefore necessary to define model fuels composed of molecules able to represent the most important classes of organic compounds present in real fuels as well as the most important features of the combustion process.

In this paper, the development and test of a kinetic mechanism for a surrogate fuel representative of gasoline used in spark ignition as well as gasoline Homogeneous Charge Compression Ignition (HCCI) internal combustion engines are presented.

Most past efforts in this area were concentrated in modeling gasoline PRF fuel mixtures composed of n-heptane and iso-octane in variable proportions. For this purpose, detailed mechanisms [2], [3] as well as reduced mechanisms have been developed [START_REF] Tanaka | [END_REF], [11]. PRF mixtures are used for octane number determination and are representative of the spontaneous auto-ignition resistance in the conditions of the RON (Research Octane Number -procedure D-2699 [START_REF]Annual Book of ASTM Standards. Test methods for rating motor, diesel and aviation fuels[END_REF]) and MON (Motor Octane Number -procedure D-2699 [START_REF]Annual Book of ASTM Standards. Test methods for rating motor, diesel and aviation fuels[END_REF]) experiments. However, RON and MON are not representative of gasoline auto-ignition resistance in all possible engine conditions [START_REF] Heywood | Internal combustion Engine Fundamentals[END_REF]. Also, laminar flame speeds of PRF mixtures are different from gasoline laminar flame speeds [START_REF] Metgalchi | [END_REF], [START_REF] Farrel | Society of Automotive Engineers[END_REF] and may thus induce errors during flame propagation simulations.

Since kinetic modeling of large molecules includes hundreds of species involved in thousands of reactions, surrogate fuels cannot contain a large number of species in order to keep the chemical mechanisms manageable and computations possible. Gasoline is mostly composed of saturated hydrocarbons (normal and branched alkanes), olefins and aromatics. A three components gasoline surrogate including PRF mixtures and toluene has then been chosen here as a gasoline surrogate fuel. Toluene should allow better predictions of auto-ignition delays in various conditions since aromatics have a large octane number sensitivity (differences between RON and MON above 10) and different laminar flame speeds compared to alkanes. Olefins have not been considered in this work. The main reasons for this were: The surrogate mechanism described here is based on already well validated available individual component mechanisms and such a mechanism for olefins is still under development [START_REF] Touchard | [END_REF]; chemical behavior of unsaturated hydrocarbons is strongly dependant on the degree of unstauration making it difficult to choose one particular olefin; very few experimental results were found in the literature to allow for a thorough validation of a more complex surrogate mechanism. Olefin addition to the gasoline surrogate mechanism is however part of work in progress.

The full mechanism is described in the first part of the paper. Special attention is given to the coupling between the toluene and PRF sub-mechanisms. The second part of the paper deals with mechanism test and validation against a wide number of literature experimental configurations and thermodynamic conditions and surrogate fuel composition.

Specific Objectives

The target of the surrogate mechanism is to simulate the auto-ignition of three components mixtures with different percentages of each fuel in a variety of thermodynamic conditions representative of spark ignition and gasoline HCCI internal combustion engines. The proposed mechanism is validated against available experimental results from the literature. These include HCCI engine experiments, rapid compression machines and a shock tube. Jet stirred reactor simulations have also been compared to available experimental results. All simulations were performed with the Chemkin 4.1 combustion modeling package [START_REF] Kee | CHEMKIN Release 4.1, Reaction Design[END_REF].

Surrogate mechanism development

The modeling approach followed in this paper is similar to that of Andrae et al. [4], who have also started from two different validated literature mechanisms: The PRF model from Lawrence Livermore National Laboratory [2] and the toluene mechanism from Dagaut et al. [START_REF] Dagaut | [END_REF]. Validation of their mechanism was performed against Homogeneous Charge Compression Ignition (HCCI) engine experiments.

The surrogate mechanism developed is the present study is based on the available PRF mixtures (n-heptane and iso-octane) mechanism of Buda et al. [19] and on the toluene mechanism of Bounaceur et al. [20], with the addition of the needed co-oxidation reactions.

PRF and toluene mechanisms

The PRF mechanism of Buda et al. [19] and the toluene model of Bounaceur et al. [20] have the advantage of being issued from the same research group, meaning that they include a common C 0 -C 2 reaction base. Thermochemical data of molecular and radical species have been calculated with the THERGAS software [21].

The PRF mechanism has been automatically generated by the EXGAS-ALKANE software [19] for the n-heptane/iso-octane mixture and includes all the usual low temperature reactions of alkanes (additions to oxygen, isomerization of peroxy radicals, formation of conjugated alkenes, cyclic ethers, aldehydes, ketones and hydroperoxides species).

The toluene model includes a sub-mechanism for the oxidation of benzene and of unsaturated C 3 -C 4 species, as well as reactions of toluene and benzyl, tolyl (methylphenyl), peroxybenzyl, alcoxybenzyl and cresoxy free radicals in the primary mechanism and those of benzaldehyde, benzyl hydroperoxyde, cresol, benzylalcohol, ethylbenzene, styrene and bibenzyl in the secondary mechanism [20]. The simulations of experiments at pressures above 40 bar have required multiplying the rate constant of the decomposition of H 2 O 2 by a factor 4 showing a not yet identified problem in the mechanism at high pressures.

These two mechanisms have already been validated in different experimental configurations under a wide range of thermodynamic conditions. They behave well in terms of auto-ignition delays measured in shock tubes and a rapid compression machine, although in some conditions, the agreement is not satisfactory for toluene. Mechanism results also agree well in terms of species evolution as a function of temperature measured in a jet stirred reactor.

Co-oxidation reactions

According to Andrae et al. [4], who have already proposed a coupled mechanism for n-heptane, iso-octane and toluene, co-oxidation reactions play an important role in both the auto-ignition of PRF mixtures as well as of n-heptane/toluene mixtures. These authors have proposed a set of coupling reactions which includes metathesis of primary fuels (RH) with R• radicals (heptyl, octyl and benzyl) and reactions between peroxy radicals (ROO•), toluene and PRF fuels.

Our previous work on PRF fuels has shown that co-oxidation reactions and species deriving from n-heptane and iso-octane were negligible. The coupling effect was only due to the pool of small radicals [22]. Concerning toluene, due to the strong electronic delocalization in benzyl radicals, co-oxidation involving these radicals has to be taken into account.

In summary, the following co-oxidation reactions have been written:

• Metathesis of benzyl radicals with n-heptane and iso-octane leading to toluene and alkyl radicals (respectively heptyl and octyl) with rate parameters given in Table 1; • Terminations between benzyl and alkyl radicals, with a rate constant of 1.0x10 13 cm 3 .mol -1 .s -1 according to collision theory; • Metathesis of secondary allylic radicals (iso-butenyl, iso-octenyl, heptenyl) with toluene with a rate constant of 1.6x10 12 exp(-7600/T) cm 3 .mol -1 .s -1 as for allyl radicals [20]; • Reactions between toluene and alkyl peroxy radicals. A sensitivity analysis has shown that among all possible alkyl peroxy reactions, only methyl peroxy (CH 3 O 2 ), due to its high concentration, needed to be considered, with a rate constant of 4.0x10 13 exp(-6040/T) cm 3 .mol -1 .s -1 [4].

Type of H-atom

Results and Discussion

The gasoline surrogate mechanism has been tested over a wide range of thermodynamic conditions for different surrogate fuel compositions. Table 2 summarizes the different gasoline three component fuel surrogates ran by each author and tested in this work. Experiments cover HCCI engine applications, rapid compression machines, a shock tube and a jet stirred reactor.

HCCI engine, Andrae et al. [4]

Andrae et al. [4] have performed experiments in a HCCI engine with PRF and n-heptane/toluene mixtures in the proportions given in Table 2 and operating conditions given in Table 3.

Simulations were ran from 99° Before Top Dead Center (BTDC) with initial pressure and temperature fit from experimental data. Since the auto-ignition delay calculation is the main objective of this work and since in an internal combustion engine, most heat losses occur after Top Dead Center, these have been neglected. The air/fuel mixtures were considered homogeneous. Figure 1 and Figure 2 show that for both PRF fuel and n-heptane/toluene mixtures, the mechanism reproduces the sensitivity of the ignition delay to the fuel composition and the initial engine conditions. Over-prediction of the maximum pressure is due to both adiabatic and homogeneous simulations whereas in the engine, the load is never really homogeneous and heat losses play an important role during piston expansion.

Authors

Rapid Compression Machine, Tanaka et al. [24]

Tanaka et al. [24] have studied the auto-ignition behavior of different fuels in a Rapid Compression Machine (RCM) over a variety of thermodynamic conditions. We have chosen here a fixed temperature and pressure (318 K and 0.1 MPa respectively) and three n-heptane, iso-octane and toluene mixtures (cf. Table 2) in order to assert the sensitivity of the mechanism to small variations of the iso-octane to n-heptane ratio with a fixed amount of toluene. Results presented in Figure 3 show this sensitivity is well captured by the mechanism although the simulated auto-ignition delay differences between fuels are smaller than the experimental ones.

Rapid Compression Machine, Vanhove et al. [25]

Vanhove et al. [25] have studied the auto-ignition as a function of temperature of a 50/50 (volume) n-heptane and toluene mixture and a 65/35 iso-octane and toluene mixture in another RCM. They have varied the initial temperatures and pressures and measured both cool flame (t cf ) and main ignition delays (t ig ). As a reference, auto-ignition delays of pure n-heptane and pure iso-octane have been measured in similar thermodynamic conditions and simulated with the same kinetic mechanism.

Figure 4 shows that both ignition delays are well captured by the mechanism and that the retarding influence of toluene over n-heptane is retrieved. Nevertheless, for the latter case, the Negative Temperature Coefficient (NTC) zone is not well reproduced. The mechanism does not show enough reactivity between lower and intermediate temperatures. Figure 5 represents the same ignition data in the case of the iso-octane/toluene and pure iso-octane mixtures. Again, the main observed experimental trends are captured by the mechanism. However, the strong inhibiting effect of toluene on iso-octane observed by Vanhove et al. is underestimated by the mechanism. The authors have attributed this effect to a deactivation of the radical pool caused by the presence of toluene. Toluene can react either with OH or HO 2 to form unreactive species such as benzyl, H 2 O or H 2 O 2 respectively. This has been confirmed by studying the sensitivity of the mechanism to the kinetic constant of reaction toluene+OH leading to benzyl+H 2 O which was shown to be very high. The interactions between iso-octane and toluene at the radical pool level are thus an open subject and are the subject of work in progress.

Shock tube, Gauthier et al. [26] Gauthier et al. [26] have studied the ignition of n-heptane, iso-octane and toluene mixtures at high temperature and pressure in a shock tube. Their experiences also covered a broad range of equivalence ratios and dilution gas proportions. We have chosen to present their results at stoichiometry and zero dilution although all reported experimental conditions have been tested.

Figure 6 shows that the mechanism is again able to capture the influence of coupled temperature, pressure and fuel differences in terms of auto-ignition delays in conditions very close to those inside an internal combustion engine. Results of simulations at different equivalence ratios and dilution gas concentrations have also shown to agree well with experimental results. Jet stirred reactor, Dubreuil et al. [27] In spite of the fact that our first objective is to conveniently simulate auto-ignition delays of representative engine fuels, it is also interesting to test the species evolution predicted by the mechanism against available Jet Stirred Reactor (JSR) results. In this context, we have simulated the Dubreuil et al. [27] n-heptane/toluene mixture experiment where the authors have measured reactant as well as main intermediate species molar concentrations as a function of temperature.

Figure 7 shows that both n-heptane and toluene evolutions as a function of temperature are well captured by the mechanism although around 650 K to 700 K, the mechanism is not reactive enough with a slower consumption of n-heptane compared to experiments. This agrees with the discrepancies obtained in the same temperature range by Vanhove et al. [25], as it can be seen in Figure 5. It is however interesting to note that one main combustion product (CO 2 ) as well as two main intermediates (CO and formaldehyde) evolutions as a function of temperature are well reproduced by the mechanism (see Figure 8).

Conclusions

A three components chemical kinetic mechanism developed to model gasoline surrogate fuel ignition and combustion was developed. The full mechanism is based on existing n-heptane/iso-octane and toluene sub-mechanisms. The main objective of this paper was to present the coupling strategy between the sub-mechanisms and to test the final version against experimental results available from the literature.

Simulations with the three components gasoline surrogate mechanism show that its behavior is in agreement with experimental results for most of the conditions tested. We emphasize the model predictivity in terms of simulating auto-ignition delays of different fuel mixtures for different temperatures, pressures and equivalence ratios covering most internal combustion engines operating conditions as well as different experimental setups.

The auto-ignition delays measured by Andrae et al. [4] in a HCCI engine with n-heptane/iso-octanne and n-heptane/toluene mixtures are reproduced by the mechanism. The simulations of the Tanaka et al. [24] experiments in a rapid compression machine show that the mechanism is also capable of retrieving the sensitivity of auto-ignition delay variations to very small iso-octane addition to toluene and n-heptane mixtures. The major trends of the results obtained by Vanhove et al. [25] in a rapid compression machine have been captured as well as the Gauthier et al. [26] shock tube auto-ignition measurements with two different three components surrogate fuels. Finally, the simulation of the Dubreuil et al. [27] experiments in a jet stirred reactor shows that major species evolutions as a function of temperature of a binary toluene/n-heptane mixture are also conveniently captured by the three components mechanism.

Progress topics have been identified and are the subject of work in progress. These include a better coupling between n-heptane and toluene and iso-octane and toluene respectively in the NTC region and on the radical pool level and also the development of a gasoline surrogate including olefins, a major component of this fuel.

A

  complex chemical kinetic mechanism for the oxidation of gasoline surrogate fuels: n heptane, iso octane and toluene -Mechanism development and validation A. Pires da Cruz, C. Pera, Jörg Anderlohr, Roda Bounaceur, Frédérique Battin-Leclerc To cite this version: A. Pires da Cruz, C. Pera, Jörg Anderlohr, Roda Bounaceur, Frédérique Battin-Leclerc. A complex chemical kinetic mechanism for the oxidation of gasoline surrogate fuels: n heptane, iso octane and toluene -Mechanism development and validation. Proceedings of the European Combustion Meeting 2007, Apr 2007, Chania, Greece. pp.CD-ROM du congrès. hal-00370855

Figure 1 :

 1 Figure 1: Experimental and simulated pressure curves for fuels A (Run 1) and B (Run 1 and Run 2). Initial temperature and pressure at start of simulations (99° BTDC) -Fuel A/Run 1: 460 K, 1.2 bar, Fuel B/Run 1: 470 K, 1.7 bar, Fuel B/Run 2: 420 K, 3.2 bar. HCCI engine of Andrae et al. [4].

Figure 2 :

 2 Figure 2: Experimental and simulated pressure curves for fuels C (Run 1) and D (Run 1 and Run 2). Initial temperature and pressure at start of simulations (99° BTDC) -Fuel C/Run 1: 470 K, 1.6 bar, Fuel D/Run 1: 470 K, 1.6 bar, Fuel D/Run 2: 420 K, 3.1 bar. HCCI engine of Andrae et al. [4].

Figure 3 :

 3 Figure 3: Experimental and simulated pressure change ( p) curves for fuels A, B and C. Equivalence ratio: 0.4, Initial temperature: 318 K, Initial pressure: 0.1 MPa, Compression ratio: 16. Rapid Compression Machine of Tanaka et al. [24].

Figure 4 :

 4 Figure 4: Experimental and simulated main (tig) and cool flame (tcf) auto-ignition delays for fuel A and n-heptane. Pressures between 3.7 bar and 4.9 bar and Equivalence ratio: 1. Rapid Compression Machine of Vanhove et al. [25].

Figure 5 :

 5 Figure 5: Experimental and simulated main (tig) and cool flame (tcf) auto-ignition delays for fuel B and iso-octane for pressures between 12 bar and 16.1 bar and Equivalence ratio: 1. Rapid Compression Machine of Vanhove et al. [25].

Figure 6 :

 6 Figure 6: Experimental and simulated auto-ignition delays for surrogate fuels A and B for pressures between 12 atm and 60 atm, Equivalence ratio: 1 and initial temperatures between 860 K and 1022 K. Shock tube of Gauthier et al. [26].

Figure 7 :

 7 Figure 7: Experimental and simulated mole fractions of n-heptane and toluene as of function of temperature for fuel A. Pressure: 10 bar, Equivalence ratio: 0.2, Residence time=0.5 s. Jet Stirred Reactor of Dubreil et al. [27].

Figure 8 :

 8 Figure 8: Experimental and simulated mole fractions of CO2, CO and CH2O as of function of temperature for fuel A. Pressure: 10 bar, Equivalence ratio: 0.2, Residence time=0.5 s. Jet Stirred Reactor of Dubreil et al. [27].

Table 1 :

 1 Rate expressions for metathesis with benzyl radicals (in s -1 , kcal, mol units), taken from the values proposed by Tsang for allyl radicals[23].

		A	n	E a
	Primary	1.6	3.3	19.84
	Secondary	1.6	3.3	18.17
	Tertiary	1.6	3.3	17.17

Table 2 :

 2 Three components fuel mixtures tested in the different authors experiments and used in this work.

	Run 1	900	120	1	0.29
	Run 2	900	40	2	0.25

Op. Cond. RPM T init (°C) p init (bar)

Table 3 :

 3

HCCI engine operating conditions of Andrae et al.

[4]

. Compression ratio: 16.7.

Acknowledgements

This research was partially funded by the GSM (Groupement Scientifique Moteur: PSA Peugeot Citroen, Renault and IFP).


Brigitte Bidégaray-Fesquet

Positiveness and Pauli exception principle in raw Bloch equations for quantum boxes

HAL is

Introduction

Bloch equations are a very common model to describe the time evolution of a system of electrons in different contexts such as gas of electrons, glasses, crystals. Mathematical results have been obtained in these contexts which ensure in particular that the variable of Bloch equations, namely the density matrix, keeps through the time evolution some properties, such as hermicity and positiveness (of diagonal elements and of the matrix as an operator) [START_REF] Bidégaray | Contributions à l'électromagnétisme dans le domaine temporel. Modélisation classique et quantique en optique non linéaire[END_REF][START_REF] Bidégaray | Introducing physical relaxation terms in Bloch equations[END_REF]. In the above contexts, the indices in the density matrix are integers which distinguish between the different electron levels. Diagonal entries model the population of the levels, while off-diagonal entries model coherences between two levels.

Bloch equations have also been derived for quantum wells [START_REF] Hess | Maxwell-Bloch equations for spatially inhomogeneous semiconductor lasers. I. Theoretical formulation[END_REF][START_REF] Kuhn | Monte Carlo simulation of ultrafast processes in photoexcited semiconductors: Coherent and incoherent dynamics[END_REF]. They look very similar but the variables are different. In the quantum well context, we distinguish between two species of electrons (valence and conduction electrons) and variables are indexed by (two-dimensional) wave vectors.

Electrons can interact directly only if they have the same wave vector, which induces intra-band coherences to be zero and valence and conduction electrons to be coupled two-by-two. This is not a mathematical problem since at first order, i.e. without considering relaxation effects, such a system is a juxtaposition of two-level models.

The model for quantum wells has been extended to quantum boxes [START_REF] Gehrig | Mesoscopic spatiotemporal theory for quantumdot lasers[END_REF]. Due to the three-dimensional confinement of electrons, variables are now again indexed by integers. Valence and conduction electron levels are not coupled two-by-two any more and any valence electron can interact with any conduction electron. But in the propounded model intra-band coherences are still considered as zero. This destroys the Liouville structure of the system which in other respects resembles much the gas model. Another puzzling point is the fact that levels are supposed to have any degree of degeneracy. This is not much coherent with the Pauli exception principle that electrons (which are fermions) are supposed to follow. This induces also mathematical and modelling problems which we also address here.

In this paper we are interested in the derivation of a raw Bloch model for a quantum box system. This model will be derived carefully from computations involving creation and annihilation operators for both species of electrons. Our goal is to obtain a model with appropriate mathematical properties in order to be able to prove the conservation of certain physical properties through time. Our study will be restricted to the leading order terms in the Hamiltonian, leaving extra contributions such as Coulomb effect or electron-phonon interaction to future works. More precisely this is not a restricted study but the first step towards a broader study, since the properties proved on the first order terms are necessary to study the higher terms.

The outline of this paper is as follows. In Section 2 we derive carefully a one-species Bloch model such as the one for gazes. The structure of this model is analysed. Positiveness results are recalled and boundedness of the populations is discussed. This motivates Section 3 where the problem of degeneracy is thoroughly addressed. Section 4 is devoted to our central problem, i.e. the derivation of a two-species Bloch model. This model has the desired mathematical properties. The details of the computations and an analogy with a boson model are postponed to various appendices.

2 The case of one species of electrons For an operator A, we will define the associated observable A A = Tr(S 0 A) by averaging with respect to the initial state density S 0 of the system. If the system is described by an Hamiltonian H, the time-evolution for this observable is given by the Heisenberg equation

ih∂ t A = [A, H] .
Bloch equations are the Heisenberg equation of motion when the observable is the density matrix.

Operators and commutation rules

We suppose there is only one species of electrons, like in the case of gas atoms, and that levels are indexed by integers i ∈ I. For the ith level, we define the creation and annihilation operators c † i and c i . Electrons are fermions and should respect the Pauli exception principle. The corresponding skew-commutation rules are

{c i , c j } = {c † i , c † j } = 0, {c i , c † j } = δ i,j .

This implies in particular that

c i c † i c i = c i , c † i c i c † i = c † i , c i c i = c † i c † i = 0
, this last equality meaning clearly that it is impossible to create twice or annihilate twice the same electron. This is the Pauli exception rule. Let A and B be the products of n A and n B such operators which all skew-commute two-by-two, then

[A, B] = 2AB if n A or n B is odd, 0 else.
Since we only use two-operator observables, we will mostly always be in the second case.

Hamiltonian and observables

We define the Hamiltonian as the sum of a part due to the electrons and only, and a part due to their interaction with an external electric field E(t):

H = H e + H L with H e = k∈I ǫ k c † k c k
where ǫ k is the energy of an electron in the kth level and

H L = 1 2 (k,l)∈I 2 (E(t) • M kl c † k c l + E * (t) • M * kl c † l c k ),
where M kl is an element of the dipolar moment matrix M. The 1/2 coefficient is due to the fact that interactions are counted twice in this sum. The dipolar moment M kl can be expressed as ψ * l |r|ψ k in terms of the wave functions associated to each level and the local position vector r. This implies in particular that the case k = l does not contribute to the sum since M kk = 0. We will also use the fact that M is hermitic : M * kl = M lk . The observable we are interested in is the density matrix, which elements are the ρ ij = c † j c i . It is clear that this matrix is Hermitian. In the computations, we distinguish between diagonal terms of this matrix also called populations ρ ii = c † i c i and off-diagonal terms, ρ ij , i = j also called coherences.

Computation of the raw Bloch equations

We have to compute the commutators of the density matrix with both Hamiltonians. The details of the computations are given in Appendix A. We obtain

[c † j c i , H e ] = (ǫ i -ǫ j )ρ ij , and 2 [c † j c i , H L ] = Re E(t) • k∈I (M ik ρ kj -M kj ρ ik ).
Gathering these results and the Heisenberg equation, raw Bloch equations

read ih∂ t ρ ij = (ǫ i -ǫ j )ρ ij + Re E(t) • k∈I (M ik ρ kj -M kj ρ ik ).
If we denote E = diag({ǫ i } i∈I ), this can be cast as a Liouville equation

ih∂ t ρ = [V (t), ρ],
where

V (t) = E + Re E(t) • M. (1)
In Appendix B we show that considering electrons as bosons leads to exactly the same equation, although the intermediate computations are different. This is not true any more if you consider extra contributions to the Hamiltonian such as Coulomb effect or electron-phonon interaction, but this is an other story.

Positiveness and trace results

Equation ( 1) clearly preserves the Hermitian structure of ρ. We can give an exact solution to equation ( 1), namely

ρ(t) = exp - i h t 0 V (τ )dτ ρ(0) exp i h t 0 V (τ )dτ .
This evolution equation preserves positiveness. Indeed, if ρ(0) is a non negative operator, that is for all X ∈ l 2 (I), X * ρ(0)X ≥ 0, then

X * ρ(t)X = X * exp - i h t 0 V (τ )dτ ρ(0) exp i h t 0 V (τ )dτ X = Y * ρ(0)Y, where Y = exp i h t 0 V (τ )dτ X.
Therefore X * ρ(t)X ≥ 0. A corollary of this is that for all time t, and all level j, ρ jj (t) ≥ 0. The total population is the sum of populations, which can be described by the trace of the density matrix. Since commutation is a trace preserving operation, we have ih∂ t Tr ρ = Tr[V (t), ρ] = 0.

Hence Tr ρ(t) = Tr ρ(0) ≡ N tot . The raw Bloch equations are trace preserving. In other words the total number of electrons is preserved through the time evolution. This together with non negativeness ensures that for all time t 0 ≤ ρ jj (t) ≤ N tot .

We want also to know whether each level remains bounded by 1, which would be the expression of the Pauli exception principle. This is indeed the case. If we set ρ = Id -ρ, where Id is the identity matrix with the same dimensions as ρ, then

[V, Id -ρ] = [V, Id] -[V, ρ] = -[V, ρ], ∂ t (Id -ρ) = -∂ t ρ.
Therefore ρ is also solution to the raw Bloch equations, and in particular ρii (t) ≥ 0, which means ρ ii (t) ≤ 1.

Degenerate levels

In [START_REF] Gehrig | Mesoscopic spatiotemporal theory for quantumdot lasers[END_REF], levels are supposed to be degenerate, with a degeneracy order that depends on the levels. In a fermion description levels cannot be degenerate since each level population is bounded by 1. This would a priori be possible in a boson description. Let us first investigate this situation.

The same proof as above shows that if ρ is solution to the raw Bloch equations and the level populations are initially bounded by the same constant (say d), then ρ ii (t) ≤ d. To derive a d-degenerate raw Bloch model, it is therefore possible to consider electrons as bosons and impose a condition on the initial data. But this solution is not very satisfactory from two points of view. First this cannot be extended to the case when the degeneracy depends on the level, i.e. the ith level is d i -degenerate. Second if we want to look at other contributions in the Hamiltonian, bosonic and fermionic commutation rules do not lead to the same model.

The question is therefore: can we derive a degenerate model from fermion calculations with level dependent degeneracies?

Deriving a degenerate model

For i ∈ I and n = 1, . . . , d i , we denote by c n † i and c n i the creation and annihilation operators associated to the nth degenerate sub-level of level i. To ensure a coherent degenerate model we have to assume the following commutation rules

{c n i , c m j } = {c n † i , c m † j } = 0, {c n i , c m † j } = δ i,j δ n,m .
The energy ǫ i does not depend on the sub-level index n. Since the levels are exactly the same (and not just have the same energy), the associated wave functions are the same and therefore the dipolar moment entries only depend on the levels and not on the sub-level. We still denote them by M ij . The Bloch equations governing

ρ nm ij = c m † j c n i is the Full Degenerate Bloch (FDB) model ih∂ t ρ nm ij = (ǫ i -ǫ j )ρ nm ij + Re E(t) • k∈I d k p=1 (M ik ρ pm kj -M kj ρ np ik ). ( 2
)
We notice that natural quantities that arise in the right handside are

ρ +m kj = d k p=1 ρ pm kj , ρ n+ ik = d k p=1 ρ np ik ,
and we write the equation for ρ

n+ ij ih∂ t ρ n+ ij = ih∂ t d j m=1 ρ nm ij = (ǫ i -ǫ j ) d j m=1 ρ nm ij + Re E(t) • k∈I (M ik d j m=1 ρ +m kj -M kj d j m=1 ρ n+ ik ) = (ǫ i -ǫ j )ρ n+ ij + Re E(t) • k∈I (M ik ρ ++ kj -M kj d j ρ n+ ik ),
where we have introduced

ρ ++ ij = d i n=1 d j m=1 ρ nm ij
and we obtain a closed set of equations for the density matrix ρ ++ , which is the Condensed Degenerate Bloch (CDB) model

ih∂ t ρ ++ ij = ih∂ t d i n=1 ρ n+ ij = (ǫ i -ǫ j ) d i n=1 ρ n+ ij + Re E(t) • k∈I (M ik d i n=1 ρ ++ kj -M kj d j d i n=1 ρ n+ ik ) = (ǫ i -ǫ j )ρ ++ ij + Re E(t) • k∈I (M ik d i ρ ++ kj -M kj d j ρ ++ ik ).
It is not exactly the same equation as for non degenerate levels since the degeneracies occur in the equation coefficients. If we set

σ ij = ρ ++ ij d i d j , N ij = M ij d i d j ,
we recover the usual Bloch equation

ih∂ t σ ij = (ǫ i -ǫ j )σ ij + Re E(t) • k∈I (N ik σ kj -N kj σ ik ). (3)
This ensures in particular that σ (and therefore ρ ++ ) defines a non negative operator if this is valid at the initial time.

Boundedness of degenerate levels

The diagonal elements are σ ii = ρ ++ ii /d i and we would expect this model to ensure σ ii (t) ≤ 1 (i.e. ρ ++ ii ≤ d i ) for all time. By the same arguments as above applied on equation (3), this is true if this hold at the initial time. The problem is that this condition is not natural when dealing with the variables of the FDB model [START_REF] Bidégaray | Introducing physical relaxation terms in Bloch equations[END_REF]. Indeed the diagonal entries of the CDB model are not the sum of the populations of a given level but also include intra-level coherences. An other consequence of the non negativeness of ρ is that |ρ ij | ≤ √ ρ ii ρ jj (for the variables of equation [START_REF] Bidégaray | Contributions à l'électromagnétisme dans le domaine temporel. Modélisation classique et quantique en optique non linéaire[END_REF], see [START_REF] Bidégaray | Introducing physical relaxation terms in Bloch equations[END_REF]). If plugged in the definition of σ ii this yields

σ ii ≤ 1 d i d i n=1 d i m=1 ρ nm ii ≤ 1 d i d i n=1 d i m=1 ρ nn ii ρ mm ii ≤ 1 d i d i n=1 ρ nn ii d i m=1 ρ mm ii ≤ d i .
Since there can be configurations where σ ii = d i , the only way to ensure that σ ii (0) ≤ 1 is to impose that intra-level coherences are initially zero. This certainly is valid in most experimental situations. Under this vanishing condition, the CDB model preserves the natural property ρ ++ ii ≤ d i . In the opposite situation, it is not possible to preserve the property ρ ++ ii ≤ d 2 i which reverts to preserve σ ii ≤ d i and is the same problem as deriving a degenerate model based on a boson derivation of Bloch equations, the initial problem that we have eluded. [START_REF] Gehrig | Mesoscopic spatiotemporal theory for quantumdot lasers[END_REF] The case of two species of electrons

Conduction-valence notations

To be closer to the above computations, we will first present the model of two species of fermions considering electron operators (and not electron and hole operators, as in [START_REF] Gehrig | Mesoscopic spatiotemporal theory for quantumdot lasers[END_REF]). For i ∈ I c (a set of integers, indexing conduction electrons), we keep the notations c † i and c i for the creation and annihilation operators associated to conduction electrons. For i ∈ I v (indexing valence electrons), we will use the non conventional notation v † i and v i for the creation and annihilation operators associated to valence electrons. Along with the already defined commutation rules for conduction electrons

{c i , c j } = {c † i , c † j } = 0, {c i , c † j } = δ i,j ,
we have the same rules for valence electrons

{v i , v j } = {v † i , v † j } = 0, {v i , v † j } = δ i,j
and commutation rules between the two species

[c i , v j ] = [c † i , v † j ] = [c i , v † j ] = [c † i , v j ] = 0.
We now consider the Hamiltonians

H c = k∈I c ǫ c k c † k c k , H v = k∈I v ǫ v k v † k v k , H Lc = 1 2 (k,l)∈(I c ) 2 (E(t) • M c kl c † k c l + E * (t) • M c * kl c † l c k ), H Lv = 1 2 (k,l)∈(I v ) 2 (E(t) • M v kl v † k v l + E * (t) • M v * kl v † l v k ), H Lcv = (k,l)∈I c ×I v (E(t) • M cv kl c † k v l + E * (t) • M cv * kl v † l c k ).
The density matrices we are interested in are intra-band densities

ρ c ij = c † j c i , ρ v ij = v † j v i ,
and inter-band densities

ρ cv ij = v † j c i , ρ vc ij = c † j v i .
We will not write the equations for ρ vc but an appropriate definition of M vc ensures that

ρ tot = ρ c ρ cv ρ cv * ρ v
is Hermitian. We now have to compute the commutators of these matrices with the three above Hamiltonians. This is done in Appendix C.

Towards a two-species Bloch equation

Gathering all the computations of Appendix C we obtain

ih∂ t ρ c ij = (ǫ c i -ǫ c j )ρ c ij + Re E(t) • k∈I c (M c ik ρ c kj -M c kj ρ c ik ) + E(t) • k∈I v M cv ik ρ vc kj -E * (t) • k∈I v M cv * jk ρ cv ik , ih∂ t ρ v ij = (ǫ v i -ǫ v j )ρ v ij + Re E(t) • k∈I v (M v ik ρ v kj -M v kj ρ v ik ) -E(t) • k∈I c M cv kj ρ vc ik -E * (t) • k∈I c M cv * ki ρ cv kj , ih∂ t ρ cv ij = (ǫ c i -ǫ v j )ρ cv ij + Re E(t) • ( k∈I c M c ik ρ cv kj - k∈I v M v kj ρ cv ik ) + E(t) • ( k∈I v M cv ik ρ v kj - k∈I c M cv kj ρ c ik ). If we denote E c = diag({ǫ c i } i∈I c ), E v = diag({ǫ v i } i∈I v ) and V (t) = E c + Re E(t) • M c E(t) • M cv E * (t) • M cv * E v + Re E(t) • M v ,
we have the Bloch equation in Liouville form

ih∂ t ρ tot = [V (t), ρ tot ].
We can in a straightforward way apply our discussion about the degeneracy of levels to this two-species context.

Electron-hole formulation

Since all valence electrons do not play a role in the interaction with conduction electrons, only those who interact are described and not by their presence but their absence: a hole in the valence band. We formulate anew the former notations and results, denoting for i ∈ I h = I v by d † i = v i the creation operator of a hole (annihilation of the corresponding electron) and

d i = v †
i its annihilation operator. The hole energy is the opposite of energy of the corresponding electron

ǫ h i = -ǫ v i . Moreover d † i d i = 1 -d i d † i = 1 -v † i v i , hence the
Hamiltonian for holes is derived from the Hamiltonian for valence electrons via

H v = k∈I v ǫ v k v † k v k = k∈I h -ǫ h k (1 -d † k d k ) = k∈I h -ǫ h k + k∈I h ǫ h k d † k d k .
The first term is a constant that will not play any rôle in the commutators, hence we define

H h = k∈I h ǫ h k d † k d k .
The intra-band and inter-band density matrices involving holes are defined by

ρ h ij = d † j d i = δ ij -ρ v ji and ρ ch ij = d j c i = ρ cv ij .
The system for these variables read

ih∂ t ρ c ij = (ǫ c i -ǫ c j )ρ c ij + Re E(t) • k∈I c (M c ik ρ c kj -M c kj ρ c ik ) + E(t) • k∈I h M ch ik ρ hc kj -E * (t) • k∈I h M ch * jk ρ ch ik , ih∂ t ρ h ij = (ǫ h i -ǫ h j )ρ h ij + Re E(t) • k∈I h (M v jk ρ h ik -M v ki ρ h kj ) + E(t) • k∈I c M ch ki ρ hc jk -E * (t) • k∈I c M ch * kj ρ ch ki , ih∂ t ρ ch ij = (ǫ c i + ǫ h j )ρ ch ij + Re E(t) • ( k∈I c M c ik ρ ch kj - k∈I h M h kj ρ ch ik ) + E(t) • ( k∈I h M ch ik (δ jk -ρ h jk ) - k∈I c M ch kj ρ c ik ). (4)
which has lost the Liouville structure. Here we have chosen to denote

M h ij = M v ij and M ch ij = M cv ij but
no other choice would help to recover this structure. Mathematical results have to be obtained from the electronelectron formulation. The electron-hole formulation should be only kept for simulation or intuition in introducing new terms such as electron-hole recombination.

Vanishing intra-band coherences

In comparison with the model proposed in [START_REF] Gehrig | Mesoscopic spatiotemporal theory for quantumdot lasers[END_REF], a natural question is now whether we recover their model by imposing vanishing intra-band coherences, i.e. ρ c ij = ρ v ij = 0 for i = j, . In this reference the variables are n e i = ρ c ii , n h j = ρ h jj and p ji = ρ ch ij (with a subtlety about coupling with forward and backward propagating optical fields, which we do not separate here). With this set of variables and vanishing intra-band coherences, the system (4) now reads

ih∂ t n e i = E(t) • k∈I h M ch ik p * ki -E * (t) • k∈I h M ch * jk p ki , ih∂ t n h j = E(t) • k∈I c M ch kj p * jk -E * (t) • k∈I c M ch * kj p jk , ih∂ t p ji = (ǫ c i + ǫ h j )p ji + Re E(t) • ( k∈I c M c ik p jk - k∈I h M h kj p ki ) + E(t) • M ch ij (1 -n h j -n e i ).
The difference with the equations announced in [START_REF] Gehrig | Mesoscopic spatiotemporal theory for quantumdot lasers[END_REF] (apart from extra terms describing phenomena out of the scope of the present paper) are the terms involving M c and M h . Instead of assuming that intra-band coherences are zero, we could assume that M c and M h are zero (or, since this is certainly not true, that the Hamiltonians H Lc and H Lh should not be taken into account in the derivation), we also then do not recover the equations in [START_REF] Gehrig | Mesoscopic spatiotemporal theory for quantumdot lasers[END_REF] because of extra terms in the equation for p ji . Under this only assumption intra-band coherences are not zero, even if taken to be zero at the initial time.

The equations in [START_REF] Gehrig | Mesoscopic spatiotemporal theory for quantumdot lasers[END_REF] have therefore lost the Liouville structure and are not suitable for the mathematical analysis. The model we propose is far better in this respect. If assuming zero intra-band coherences is not coherent with the Liouville structure, we can easily have them very small by simply assuming them to be initially zero and adding in their equations strong damping terms. As shown in [START_REF] Bidégaray | Introducing physical relaxation terms in Bloch equations[END_REF], positiveness results can still be obtained in this context. Besides the frequencies of E are chosen to be close to the gap frequency, in order to match valence-conduction transition frequencies. The intra-band coherences are hence not directly excited by such a wave.

Conclusion

We have derived raw Bloch equations for quantum boxes, i.e. two species of electrons. These equations have the Liouville form which allows to prove positiveness results. We have also discussed the problem of level degeneracy and derived a condensed degenerate Bloch model in which diagonal entries remain bounded by the corresponding level degeneracy through time evolution, even when the degeneracy order is level-dependent.

Positiveness and boundedness results are necessary to prove similar properties when extra terms are added in these raw equations. They are also required to prove sharp existence results when coupled with a model for wave propagation (e.g. Maxwell equations) as already done for one species of electrons [START_REF] Dumas | Global existence for Maxwell-Bloch systems[END_REF].

Future work will consist in introducing with the same care extra con-tributions to the Hamiltonian to model Coulomb effect, electron-phonon interactions, electron-hole recombination, . . . Our goal is to obtain a model close to that of [START_REF] Gehrig | Mesoscopic spatiotemporal theory for quantumdot lasers[END_REF] but with anew possible other contributions, and full equations for intra-band coherences.

A Detail of one species fermion computations

Let us first compute the commutator with the free electron Hamiltonian. We compute separately the cases i = j and i = j but notice afterwards that results reads the same in both cases namely

[c † i c i , H e ] = M kj c † j c i c † k c j -c † k (1 -c † j c j )c i +E * (t) • k =j M * ki c † j (1 -c † i c i )c k -c † i c k c † j c i +E(t) • M ij c † j (1 -c † i c i )c j -c † i (1 -c † j c j )c i +E * (t) • M * ji c † j (1 -c † i c i )c j -c † i (1 -c † j c j )c i = E(t) • l =j M il c † j c l -E * (t) • l =i M * jl c † l c i -E(t) • k =i M kj c † k c i + E * (t) • k =j M * ki c † j c k +E(t) • M ij c † j c j -c † i c i + E * (t) • M * ji c † j c j -c † i c i [c † j c i , H L ] = Re E(t) • k M ik ρ kj -Re E(t) • k M kj ρ ik .

B Boson computations

If electrons are considered as bosons, or in other words without Pauli exception rule, the definitions for the operators, Hamiltonian and density matrix are the same as for fermions, but the computations differ because of com-mutation rules between creation and annihilation operators which are now

[c i , c j ] = [c † i , c † j ] = 0, [c i , c † j ] = δ i,j .
We still have [c † i c i , H e ] = 0. We first compute the commutators of the density matrix elements with the electron Hamiltonian

[c † j c i , H e ] = k ǫ k c † j c i c † k c k -c † k c k c † j c i = ǫ i c † j c i c † i c i -c † i c i c † j c i + ǫ j c † j c i c † j c j -c † j c j c † j c i = ǫ i c † j (1 + c † i c i )c i -c † i c i c † j c i + ǫ j c † j c i c † j c j -c † j (1 + c † j c j )c i = ǫ i c † j c i -ǫ j c † j c i = (ǫ i -ǫ j )ρ ij .
Then we compute the commutators of the density matrix elements with the laser Hamiltonian

2 [c † i c i , H L ] = E(t) • l =i M il c † i c i c † i c l -c † i c l c † i c i +E * (t) • l =i M * il c † i c i c † l c i -c † l c i c † i c i +E(t) • k =i M ki c † i c i c † k c i -c † k c i c † i c i +E * (t) • k =i M * ki c † i c i c † i c k -c † i c k c † i c i = E(t) • l =i M il c † i (1 + c † i c i )c l -c † i c l c † i c i +E * (t) • l =i M * il c † i c i c † l c i -c † l (1 + c † i c i )c i +E(t) • k =i M ki c † i c i c † k c i -c † k (1 + c † i c i )c i +E * (t) • k =i M * ki c † i (1 + c † i c i )c k -c † i c k c † i c i = E(t) • l =i M il c † i c l -E * (t) • = E(t) • k =j M ik ρ kj -E * (t) • k =i M kj ρ ik -E(t) • k =i M kj ρ ik + E * (t) • k =j M ik ρ kj +E(t) • M ij (ρ jj -ρ ii ) + E * (t) • M ij (ρ jj -ρ ii ), [c † j c i , H L ] = Re E(t) • k M ik ρ kj -Re E(t) • k M kj ρ ik .
Although calculations are different, we find the same results as in the boson case. If we consider the Pauli exception rule or not we obtain the same raw Bloch equations.

C Detail of two-species fermion computations

Because conduction and valence operators commute, we clearly have

[c † j c i , H v ] = [v † j v i , H c ] = [c † j c i , H Lv ] = [v † j v i ,
H Lc ] = 0. Besides we have already computed in the one-species case

[c † j c i , H c ] = (ǫ c i -ǫ c j )ρ c ij , [v † j v i , H v ] = (ǫ v i -ǫ v j )ρ v ij , [c † j c i , H Lc ] = Re E(t) • k∈I c (M c ik ρ c kj -M c kj ρ c ik ), [v † j v i , H Lv ] = Re E(t) • k∈I v (M v ik ρ v kj -M v kj ρ v ik ).
We now compute

[v † j c i , H c ] = k∈I c ǫ c k v † j c i c † k c k -c † k c k v † j c i = ǫ c i v † j c i c † i c i -c † i c i v † j c i = ǫ c i v † j c i = ǫ c i ρ cv ij , [v † j c i , H v ] = k∈I v ǫ v k v † j c i v † k v k -v † k v k v † j c i = ǫ v j v † j c i v † j v j -v † j v j v † j c i = -ǫ v j v † j c i = -ǫ v j ρ cv ij .
These computations are valid also for i = j, which is a case which has to be treated apart for the following computations.

[c † i c i , H Lcv ] = E(t) • (k,l)∈I c ×I v M cv kl c † i c i c † k v l -c † k v l c † i c i +E * (t) • (k,l)∈I c ×I v M cv * kl c † i c i v † l c k -v † l c k c † i c i

2. 1

 1 Commutators and Heisenberg equationLet A and B be two operators, we define their commutator by [A, B] = AB -BA and their skew-commutator by {A, B} = AB + BA.