M Perin

C Chandre

E Tassi

Hamiltonian fluid reductions of drift-kinetic equations and the correspondence with water-bag distribution functions

Hamiltonian models for the first three moments of the drift-kinetic distribution function, namely the density, the fluid velocity and the parallel pressure, are derived from the Hamiltonian structure of the drift-kinetic equations. The link with the water-bag closure is established, showing that, unlike the one-dimensional Vlasov equations, these solutions are the only Hamiltonian fluid reductions for the drift-kinetic equation. These models are discussed through their equations of motion and their Casimir invariants.

I. INTRODUCTION

An effective model for describing the dynamics of strongly magnetized plasmas is provided by the drift-kinetic equations (see, e.g. Ref. [1]). These describe the evolution of distribution functions of the guiding centers of the plasma particles and apply to phenomena where the characteristic frequencies are much lower than the cyclotron frequency and where the characteristic wave lengths are much longer than the Larmor radius of the particle species under consideration. The advantage of the drift-kinetic description with respect to the more complete Vlasov description, is that in the former, the distribution functions, at the leading order, are defined on a reduced phase space, when compared to the latter. Indeed for drift-kinetic distribution functions the dependence on the coordinate corresponding to the gyration angle of the particles around the magnetic field, is ignored, because its evolution takes place on time scales too short for the phenomena of interest.

In this article we consider a simple drift-kinetic system in slab geometry. We assume that the evolution of a particle species (ions or electrons) in the plasma is described by the drift-kinetic equation

∂ f ∂t = -v ∂ f ∂z + q m ∂φ ∂z ∂ f ∂v - c B [φ, f ] + µc ql [x, f ]. (1)
Such equation can be obtained, for instance, from the gyrokinetic equation of Ref. [2], in the limit of small Larmor radii.

In Eq. ( 1), we indicate with f (r, v, µ, t) the guiding center distribution function of the particle species under consideration, where r = (x, y, z) denotes the position in configuration space, v corresponds to the velocity of the particles along the direction of the magnetic field, µ is the particle magnetic moment and t denotes time. The magnetic field is assumed to be straight and weakly inhomogeneous in space. Its expression is given by B(x, y) = B(1-x/l)ẑ, where B is a constant and l is a characteristic length of variation of the magnetic field, such that x/l ≪ 1. As the magnetic field is along z, we will refer to z as the longitudinal direction and x and y as the transverse directions. The constants c, q and m indicate the speed of light, for all functions g and h. We also assume that the distribution function f be of the form

f (r, v, µ, t) = mf (r, v, t)δ(µ -µ 0 )/(2πB),
that is we assume a beam-like distribution with respect to the magnetic moment. Making use of this assumption, from Eq. ( 1) one easily obtains the following evolution equation for the function f (r, v, t):

∂f ∂t = -v ∂f ∂z + q ∂φ ∂z ∂f ∂v -[φ, f ] + q ∂f ∂y , (2)
where the value of the physical constants c, m, B, µ 0 and l is set equal to unity and where, to keep track of the particle charge, we introduced the symbol q, whose value corresponds to 1 or -1 depending on whether the drift-kinetic species consists of ions or electrons, respectively.

Equation ( 2) is complemented by an equation relating the electrostatic potential φ to f .

The results presented in the paper hold for a generic equation of the type φ = L dvf , where L is some self-adjoint linear operator. This allows the treatment of different specific examples, just by choosing different expressions for L. For instance, the case of a Poisson's equation for drift-kinetic electrons with a neutralizing ion background corresponds to

L f dv = ∆ -1 f dv -1 .
On the other hand, the case of a quasi-neutrality relation with drift-kinetic ions and adiabatic electrons, as given in Ref. [3], corresponds to

L f dv = f dv -1.
With regard to the physical interpretation of the terms appearing in Eq. ( 2), we remark that the first three terms correspond to the one-dimensional Vlasov-Poisson equation which describes the dynamics of the distribution function along z and v. In particular, the second and third term account for the advection and for the action of the electric force, respectively, along the longitudinal direction. The last two terms, on the other hand, account for the E × B and the ∇B drifts respectively, which affect the evolution of the distribution function in the direction transverse to the strong magnetic field.

A way of solving the drift-kinetic equation is to use fluid reductions, that describe the time evolution of P i (r, t), namely the i-th fluid moment of the distribution function with respect to the velocity v, such that

P i (r, t) = v i f (r, v, t) dv. ( 3
)
The use of fluid reductions of kinetic models is ubiquitous in plasma physics. Indeed, whereas the drift-kinetic equation involves the distribution function of particles f defined in the fourdimensional phase space (r, v), fluid models rely on quantities such as the density ρ = P 0 , the velocity u = P 1 /P 0 and the pressure p = P 2 -P 2 1 /P 0 which are defined solely in the three-dimensional configuration space r. In addition of being more tangible than their kinetic counterpart, fluid models require less numerical resources to be solved. On the other hand, fluid models derived from kinetic equations lead to an infinite set of moment equations. Specific closures have to be performed in order to truncate this infinite hierarchy of moments resulting from the reduction procedure. A way of closing the system is to compute the closure in order to recover some kinetic effects such as the growth rate of kinetic instabilities (see e.g.

Refs. [4] and [5]). Our purpose in this article is to find closures that preserve the geometrical features of the drift-kinetic equation. Indeed, Eq. ( 2) can be cast into a Hamiltonian form.

Not preserving this Hamiltonian structure during the fluid reduction may lead to models that exhibit some non-physical dissipation. A Hamiltonian fluid model for the first two moments of the distribution function derived from the drift-kinetic equation has been found in Ref. [6]. Indeed, if we consider, e.g., functionals of the first two fluid moments, namely functionals of the kind F [P 0 , P 1 ], we have to perform a closure on P 2 by expressing it as a function of P 0 and P 1 such that P 2 = P 2 (P 0 , P 1 ). Moreover, it has been shown that in order to preserve the Hamiltonian structure of the parent model, namely the drift-kinetic equation, one has to impose P 2 = P 2 1 /P 0 + AP 3 0 , where A is a constant [7]. This closure corresponds to a single water-bag distribution function of height A, namely a distribution function of the type

f (r, v, t) = A Θ v - P 1 (r, t) P 0 (r, t) + P 0 (r, t) 12A -Θ v - P 1 (r, t) P 0 (r, t) - P 0 (r, t) 12A ,
where Θ is the Heaviside distribution (see e.g. Ref. [START_REF] Besse | [END_REF]). A closure preserving the Hamiltonian structure for an arbitrary number of moments in the "δf " approximation of the drift-kinetic equation has been presented in Ref. [9]. So far it has not been generalized to the fully non-linear case.

In this article, we derive a fluid model for the first three moments of the fully non-linear drift-kinetic distribution function. We consider the density ρ, the velocity u and the second reduced moment S 2 as dynamical variables (see Ref. [10]). Compared to the two field model, the three field model has the advantage of including S 2 (or equivalently P 2 ) as a dynamical variable. This is significant since the Hamiltonian depends explicitly on P 2 through the kinetic energy. We construct a closure on the third and fourth moments respectively such that the resulting fluid model preserves the Hamiltonian structure of the parent kinetic one. Furthermore, we investigate the role of the dimensionality of the system in the closure procedure. Indeed, the drift-kinetic model, which accounts for the evolution of a distribution function defined in a four-dimensional phase space, represents a useful step between the onedimensional Vlasov-Ampère equations (see Ref. [11]) and the full Vlasov-Maxwell system.

The method developed hereby makes use of some Mathematica c code which is provided in the appendix in order for the reader to be able to reproduce the calculations.

II. HAMILTONIAN FLUID REDUCTION OF THE DRIFT-KINETIC EQUA-TION

The drift-kinetic model described by Eq. ( 2) can be cast in a Hamiltonian form by introducing the bracket given in Ref. [6] {F,

G} = f ∂F f ∂z ∂G f ∂v - ∂F f ∂v ∂G f ∂z -q[F f , G f ] drdv, (4)
where F f = δF/δf denotes the functional derivative of F with respect to f such that

F [f + δf ] -F [f ] = δF δf δf drdv + O δf 2 .
The Hamiltonian writes

H[f ] = 1 2 f v 2 + 2[1 -x] + qL f dv ′ drdv, (5)
in which the first term accounts for the kinetic energy of the particles in the longitudinal direction while the other terms account for the kinetic energy of the particles in the direction perpendicular to the magnetic field. Equation ( 2) is then obtained by using ∂ t f = {f, H}.

Using Eq. ( 3) to perform the change of variables from the distribution function f to the fluid moments P i , Bracket (4) and Hamiltonian (5) write respectively

{F, G} = n,m∈N nP n+m-1 G n ∂F m ∂z -F n ∂G m ∂z -qP n+m [F m , G n ] dr, (6)
and

H[P 0 , P 1 , P 2 ] = 1 2 (P 2 + P 0 [2(1 -x) + qLP 0 ]) dr, (7)
where F m = δF/δP m denotes the functional derivative of F with respect to P m . Bracket (6) can be rewritten by introducing the operator J such that

{F, G} = n,m∈N F n J nm G m dr, (8)
where

J nm = -m∂ z P n+m-1 -(n + m)P n+m-1 ∂ z + q[P n+m , •].
The fluid reduction process described above allows us to replace the distribution function of particles f defined in phase space with an infinite number of fluid moments P i defined in configuration space. In order to be able to compute the dynamics of these variables, it is mandatory to truncate the infinite hierarchy of moments. Due to the terms P n+m-1

and P n+m in Bracket ( 8), functionals of the kind F [P 0 , . . . , P N ] for any N ≥ 1 do not constitute a Poisson subalgebra of the bracket. In other words, the bracket between two functionals of the type F [P 0 , . . . , P N ] not only depends on (P 0 , . . . , P N ) but also on N higher order moments (P N +1 , . . . , P 2N ). As a consequence, one has to express these higher order moments P N +1≤n≤2N with respect to the variables P n≤N in order to be able to close the infinite hierarchy of moments.

When considering functionals of the type F [P 0 , P 1 , P 2 ], Bracket [START_REF] Besse | [END_REF] becomes

{F, G} = 2 n,m=0 F n J nm G m dr. ( 9
)
Denoting F χ = (F 0 , F 1 , F 2 ) and defining the inner product

a, b = 2 n=0
a n b n dr, Bracket (9) reads {F, G} = F χ , J G χ , and can be split into a longitudinal and a transverse (with respect to the magnetic field) part such that {F, G} = {F, G} + {F, G} ⊥ , where

{F, G} = F χ , J G χ with J = -      0 ∂ z P 0 2∂ z P 1 0 ∂ z P 1 2∂ z P 2 0 ∂ z P 2 2∂ z P 3      -      0 P 0 2P 1 P 0 2P 1 3P 2 2P 1 3P 2 4P 3      ∂ z ,
corresponds to the one-dimensional Vlasov bracket found in Refs. [12,13] that accounts for the dynamics along the direction of the magnetic field, namely z. The bracket corresponding to the transverse part reads {F, G} ⊥ = F χ , J ⊥ G χ , where

J ⊥ = q      [P 0 , •] [P 1 , •] [P 2 , •] [P 1 , •] [P 2 , •] [P 3 , •] [P 2 , •] [P 3 , •] [P 4 , •]     
, and describes the dynamics in the transverse direction to the magnetic field, namely the xy-plane. Like in the one-dimensional Vlasov case, J introduces a contribution of the third order moment P 3 which does not belong to the subset of functional of the form F [P 0 , P 1 , P 2 ].

In addition, J ⊥ introduces a contribution of the fourth order moment P 4 which constitutes an important difference with the one-dimensional Vlasov equation. In this case, the closure has to be performed on both P 3 and P 4 as functions of P 0 , P 1 and P 2 such that P 3 = P 3 (P 0 , P 1 , P 2 )

and P 4 = P 4 (P 0 , P 1 , P 2 ). We want to compute these closures such that the resulting fluid model is Hamiltonian, i.e., in particular that Bracket (9) satisfies the Jacobi identity, namely

J = {F, {G, H}} + {H, {F, G}} + {G, {H, F }} = 0,
for any functionals F , G and H. J = 0 provides the necessary and sufficient condition on

J ⊥ = {F, {G, H} ⊥ } ⊥ + = 0.
Indeed, J and J ⊥ are both bilinear differential operators of order one. However, while J differentiates along the z direction, J ⊥ only involves derivatives in the transverse, respectively x and y, directions. As these directions are independent from one another, it is required that both brackets satisfy the Jacobi identity. Furthermore, it has been shown in Ref. [11] that in order for the bracket {•, •} to satisfy the Jacobi identity (or equivalently to make J vanish), we must have

P 3 = 3 P 1 P 2 P 0 -2 P 3 1 P 2 0 + P 4 0 S 3 P 2 P 3 0 - P 2 1 P 4 0 , (10)
where S 3 is an arbitrary function of its argument. This result suggests the introduction of the reduced moments, defined in Ref. [10] as

ρ = f dv, u = 1 ρ vf dv, S 2 = 1 ρ 3 (v -u) 2 f dv.
The change from the variable χ = (P 0 , P 1 , P 2 ) to ψ = (ρ, u, S 2 ) is given by

ρ = P 0 , u = P 1 P 0 , S 2 = P 2 P 3 0 - P 2 1 P 4 0 .
The inverse of this transformation is given by

P 0 = ρ, P 1 = ρu, P 2 = ρu 2 + ρ 3 S 2 .
By using the chain rule, noting that

δF δP 0 = δF δρ - u ρ δF δu + u 2 ρ 3 -3 S 2 ρ δF δS 2 δF δP 1 = 1 ρ δF δu -2 u ρ 3 δF δS 2 δF δP 2 = 1 ρ 3 δF δS 2 , Bracket (9) becomes {F, G} = F ψ , KG ψ , (11)
where K can be decomposed such that K = K + K ⊥ and where

K =      0 0 0 0 0 ρ -1 ∂ z S 2 0 -ρ -1 ∂ z S 2 -2∂ z (ρ -2 S 3 )      -      0 1 0 1 0 0 0 0 4ρ -2 S 3 (S 2 )      ∂ z ,
is associated with the bracket in the longitudinal direction and

K ⊥ = -q      0 0 0 0 [S 2 , ρ] 3ρ -1 [S 3 , ρ] + ρ -1 [u, S 2 ] 0 ρ -1 [S 3 , ρ] -ρ -1 [u, S 2 ] 6[S 2 2 , ρ -1 ] -3 [S 4 , ρ -1 ] -2 [u, ρ -2 S 3 ]      + q      [ρ, •] [u, •] [S 2 , •] [u, •] ρ -2 [ρ 3 S 2 , •] ρ -4 [ρ 4 S 3 , •] [S 2 , •] ρ -4 [ρ 4 S 3 , •] ρ -6 [ρ 5 S 4 , •] + 4ρ -2 S 3 [u, •] -3ρ -4 [ρ 3 S 2 2 , •]     
, leads to the definition of the bracket in the transverse direction. The Hamiltonian of the system writes

H[ρ, u, S 2 ] = 1 2 ρu 2 + ρ 3 S 2 + ρ [2(1 -x) + qLρ] dr. ( 12
)
Expressed in terms of the variables, namely (ρ, u, S 2 ), the closure on P 3 given by Eq. ( 10) eventually writes S 3 (S 2 ). The generic closure on P 4 is now replaced by the corresponding reduced moment S 4 which reads

S 4 = P 4 (ρ, u, S 2 ) ρ 5 - u 4 ρ 4 -6 u 2 ρ 2 S 2 -4 u ρ S 3 (S 2 ) = S 4 (ρ, u, S 2 ).
As stated previously, in order for Bracket (11) to satisfy the Jacobi identity we must have J = J ⊥ = {F, {G, H} ⊥ } + {F, {G, H} } ⊥ + = 0 where {F, G} ⊥ = F ψ , K ⊥ G ψ and {F, G} = F ψ , K G ψ . Due to the constraint S 3 = S 3 (S 2 ), we already have J = 0. We now look for constraints in order to make J ⊥ = {F, {G, H} ⊥ } ⊥ + vanish. This is done by computing the Jacobi identity for particular choices of the functionals F , G and H. The calculations are lengthy and, consequently, are done with Mathematica [14]. In Appendix, we provide the details that allow us to compute the necessary and sufficient constraints on S 3 and S 4 in order for Bracket (11) to satisfy the Jacobi identity. These constraints write respectively

3S 2 -(S ′ 3 ) 2 + 4S 3 S ′′ 3 = 0, (13)
and

S 4 = 1 5 4S 3 S ′ 3 + 9S 2 2 + C ρ 5 , ( 14
)
where C is a constant. Equations ( 13) and ( 14) are necessary and sufficient conditions in order to cancel J ⊥ . In addition, it is shown in Appendix that they are sufficient conditions in order for Bracket (11) to satisfy the Jacobi identity. We notice that the fluid closures for the drift-kinetic equation are more constrained than in the one-dimensional Vlasov equation where any function S 3 = S 3 (S 2 ) leads to a Hamiltonian closure.

III. CORRESPONDENCE WITH WATER-BAGS

As shown in the previous section, Hamiltonian models with three moments derived from the drift-kinetic equation are those whose third and fourth order reduced moments, namely S 3 and S 4 , satisfy Eqs. ( 13) and ( 14) respectively. However, in order for S 4 , which corresponds to the kurtosis of the distribution function, to stay bounded as the particle density ρ tends to 0, one has to impose C = 0 in Eq. ( 14).

Assuming that there exists an open set in which the solution of Eq. ( 13) is of constant sign, it can be transformed into an Emden-Fowler equation (see e.g. Ref. [START_REF] Polyanin | Handbook of Exact Solutions for Ordinary Differential Equations[END_REF]) whose solution is given in parametric form

S 2 (λ) = 4 3 b 2 + cλ + b 2 λ 2 - 1 64 λ 4 , (15)
S 3 (λ) = ± c + bλ - 1 16 λ 3 2 , ( 16
)
where λ parametrizes the curve (S 2 (λ), S 3 (λ)) and where b and c are constant. This implies that the solution S 3 (S 2 ) is given by parts such that S 3 is of constant sign. S 3 (λ) has either one or three real roots. The solutions S 3 (S 2 ) to Eq. ( 13) can be sorted out accordingly depending on the sign of ∆ = 64b 3 -27c 2 . Moreover, the roots of S 3 (λ) are such that S 2 ≥ 0, which is a necessary condition as S 2 corresponds to the variance of the distribution function. Indeed, let λ 0 be such that S 3 (λ 0 ) = 0. We thus have

S 2 (λ 0 ) = 3 64 λ 2 0 - 16 3 b 2 ≥ 0.
The case S 2 (λ 0 ) = 0 corresponds to ∆ = 0, otherwise S 2 (λ 0 ) > 0. For ∆ < 0, the solution S 3 (S 2 ) to Eq. ( 13) has two branches with a single root. For ∆ > 0, it has four branches and three roots. This is summarized in Fig. 1 which represents the nature of the solutions depending on the position in the parameter space (b, c). It is worth noting that ∆ is even in c such that Fig. 1 is symmetric for c < 0 through the change of parameter λ → -λ.

Moreover, for all b < 0, we have ∆ < 0 such that the corresponding solutions have only one real root (thus belong to the dark pink region in Fig. 1).

f (r, v, t) = a 1 [Θ (v -v 1 (r, t)) -Θ (v -v 2 (r, t))] + (a 1 + a 2 ) [Θ (v -v 2 (r, t)) -Θ (v -v 3 (r, t))] ,
where a 1 and (a 1 + a 2 ) are the heights of the two water-bags respectively and v 1 (r, t), v 2 (r, t) and v 3 (r, t) are the contours of the bags (see e.g. Refs. [3,[START_REF] Besse | [END_REF][START_REF] Gravier | [END_REF]). The associated closure

S W B 3 = S W B 3 (S W B 2
) is given in Ref. [17] in parametric form, namely

S W B 2 (n 1 ) = a 2 1 + 6a 2 a 1 n 2 1 + 4a 2 (a 2 -a 1 )n 3 1 -3a 2 2 n 4 1 12a 2 1 (a 1 + a 2 ) 2 , S W B 3 (n 1 ) = - a 2 (n 1 -1) 2 n 2 1 (a 1 + a 2 n 1 ) 2 4a 3 1 (a 1 + a 2 ) 3
= a 2 1 + a 1 a 2 + a 2 2 12a 1 |a 2 |(a 1 + a 2 ) , c = (a 2 -a 1 )(2a 1 + a 2 )(2a 2 + a 1 ) 54[a 1 |a 2 |(a 1 + a 2 )] 3/2 , λ = 2(a 1 -a 2 ) + 6a 2 n 1 3 a 1 |a 2 |(a 1 + a 2 ) .
It is shown, e.g., in Ref. [18] that the water-bag distribution function is preserved by the dynamics given by the drift-kinetic equation. Hence, any fluid model derived from it is Hamiltonian for an arbitrary number of fields. Furthermore, as it has been shown, the water-bag closure corresponds to the only physically relevant closure with three moments derived from the drift-kinetic equation given by Eq. ( 2) as it is the case with two moments.

Whether this result extends to higher numbers of moments is an open question.

For ∆ = 0, there is a one-parameter family of solutions which do not exactly correspond to the water-bags. More precisely, these solutions, which are given by Eqs. ( 15) and ( 16)

with c = 8(b/3) 3/2 , constitute the limiting case a 2 /a 1 → +∞ and the corresponding closures are represented on Fig. 2. These solutions have two branches, one physically irrelevant depicted in dashed line and one physically relevant for which S 3 = 0 when S 2 = 0.

IV. THE HAMILTONIAN THREE MOMENT DRIFT-KINETIC MODEL

The bracket and the Hamiltonian of the three moment model are given respectively by Eq. ( 11) and Eq. ( 12). The corresponding equations of motion are then obtained by using 16)) as a function of S 2 (given by Eq. ( 15)) for b = 0.4 and c = 8(b/3) 3/2 (which corresponds to ∆ = 0). All the quantities are expressed in arbitrary units. Ḟ = {F, H} which leads to

∂ρ ∂t = - ∂(ρu) ∂z -[Lρ, ρ] + q ρ y , ∂u ∂t = -u ∂u ∂z - 1 ρ ∂(ρ 3 S 2 ) ∂z -q ∂(Lρ) ∂z -[Lρ, u] + q u y , ∂S 2 ∂t = -u ∂S 2 ∂z - 1 ρ 3 ∂(ρ 4 S 3 ) ∂z -[Lρ, S 2 ] + q S 2 y .
It is more common to consider the pressure p = ρ 3 S 2 as dynamical variable instead of S 2 .

The corresponding equation of motion writes

∂p ∂t = -u ∂p ∂z -3p ∂u ∂z - ∂(ρ 4 S 3 ) ∂z -[Lρ, p] + q p y .
The equations of motion are very much like the ones for the one-dimensional Vlasov equation.

The only difference comes from the advection term associated with the transverse drifts (two rightmost terms in the above equations). As a consequence, if one just adds these terms in the equations of motion while keeping the closure used with the one-dimensional fluid model resulting from the Vlasov equation (see Ref. [11]), it is very much likely that the Hamiltonian structure is lost as the drift-kinetic closure is much more constrained. This would result in a dissipative system. The consequence of this can be seen by considering the Casimir invariants as we show below.

Casimir invariants C of a Poisson bracket {•, •} are quantities that commute with all functionals F , namely {F, C} = 0 for all F . In the three field model, one can show that the system possesses three Casimir invariants given by One can show that if S 3 = 0 (the case S 3 = 0 is trivial as it is not Hamiltonian), the previous quantities satisfy {F, C} = 0 for all F if and only if Eqs. ( 13) and ( 14) are satisfied or equivalently if and only if the associated bracket satisfies the Jacobi identity.

C 1 = ρ dr, C 2 = ρκ 0 dr, C 3 = u - 1
Therefore, there is an equivalence between the model being Hamiltonian and the fact that it possesses as many Casimir invariants as dynamical field variables. This feature is strongly linked to the fact that the possible Hamiltonian closures are the water-bag ones. In other words, if advection terms associated with transverse drifts are added to the equations of motion without ensuring that the closure used is Hamiltonian (or equivalently water-bag), two conserved quantities are lost and non-physical diffusion across surfaces defined by these quantities occur. The brackets introduce the closures S 3 (S 2 ) and S 4 (ρ, u, S 2 ) which we define by S3=Subscript [S,3][psi[ [3]]];

S4=Subscript [S,4][psi[ [1]],psi[ [2]],psi[ [3]]];

It is worth noting that the definition of the function bracket, while matching the definition given by Eq. ( 11), does not include the integration. This is done on purpose in order to simplify the code. As shown later, the results will not be affected. J is computed by using of the integral in the definition of the bracket as stated previously. Furthermore, the contributions of second order derivatives of the functionals can be omitted due to the lemma given in Ref. [START_REF] Morrison | AIP Conf. Proc[END_REF].

The Jacobi identity has to be satisfied for any triplet of functionals (F, G, H). We choose in particular (F, G, H) = (u, S 2 dr, u dr). This is done by defining dG={0,0,1}; dH={0,1,0};

computing J and then computing

VariationalD[J,dF[[2]],r]

the latter being the translation of F ψ = (0, δ, 0). However, as the dynamical variables are independent from one another, the terms proportional to (∂ y ρ) 2 ∂ 2 x u, which are found by using

  the particle charge and mass, respectively. The canonical Poisson bracket [•, •] is defined by[g, h]

P 3 and P 4

 4 for the closure to be Hamiltonian. J can be decomposed in the following way asJ = {F, {G, H} } + {F, {G, H} ⊥ } + {F, {G, H} } ⊥ + {F, {G, H} ⊥ } ⊥ + ,where means that the summation is performed over all the circular permutations of the functionals F , G and H. The two brackets {•, •} and {•, •} ⊥ have to satisfy the Jacobi identity independently from each other, i.e., we must have J = {F, {G, H} } + = 0 and

FIG. 1 .

 1 FIG.1. Diagram representing the nature of the parametric solutions given by Eqs. (15) and (16) as a function of b and c. In each regions, the insets represent the solutions S 3 (S 2 ). All the quantities are expressed in arbitrary units.

, where n 1 ∈(S W B 2 (n 1 ), S W B 3 (n 1 ) 2 (n 1 3 (n 1 )

 121312131 [0; 1] is the density of particles in the first bag which parametrizes the curve ). One can verify that for any value of the parameters b and c within the light blue region in Fig.1 (i.e. the region where ∆ > 0) there exists a change of parameters that leads to the water-bag closure. Indeed, this is done by performing the invertible change λ → αn 1 + β and solving the equations S 2 (αn 1 + β) = S W B ) and S 3 (αn 1 + β) = S W B which leads to b

3 FIG. 2 .

 32 FIG.2. S 3 (given by Eq. (16)) as a function of S 2 (given by Eq. (15)) for b = 0.4 and c = 8(b/3) 3/2

3 .

 3 C 1 and C 2 correspond respectively to the total mass and the total entropy of the system. C 3 is referred to as the generalized velocity. If the Hamiltonian structure is lost, C 2 and C 3 are no longer conserved in general which is at the origin of the fake dissipation introduced in the model. For water-bags, these Casimir invariants can be expressed as functions of the water-bags dynamical variables v i . Indeed, they are linear combinations of v 1 dr, v 2 dr and v 3 dr which are the Casimir invariants of the two water-bag model.

  In summary, starting from the drift-kinetic equation, we derive a Hamiltonian model for the first three moments of the distribution function, namely the density, the fluid velocity and the pressure. The computation of the Jacobi identity provides constraints on the fluid closures in order for the resulting model to be Hamiltonian. The solutions of these equations are given in parametric form, which allow us to identify them with the water-bag closure which, in this case, corresponds to the only Hamiltonian closures. Even though the resulting equations of motion differ only by the terms associated with the transverse drifts from those derived from the one-dimensional Vlasov equation, we show that the closure is more constrained in the drift-kinetic case. The canonical bracket bracketXY is defined such that bracketXY[a_,b_]:=D[a,x]D[b,y]-D[a,y]D[b,x];

  the command line J=bracket[dF,VariationalD[bracket[dG,dH],psi,r]]+ bracket[dH,VariationalD[bracket[dF,dG],psi,r]]+ bracket[dG,VariationalD[bracket[dH,dF],psi,r]]; The function VariationalD[a,b,c] computes the functional derivative of a(b) dc, performing operations such as integration by parts by assuming vanishing boundary terms. As the function VariationalD acts on the argument of the integral, this justifies the omission

ACKNOWLEDGMENTS

This work was supported by the Agence Nationale de la Recherche (ANR GYPSI) and the French Research Federation for fusion study.

Appendix: Computation of J

The computation of the constraints given by Eqs. ( 13) and (14), which are necessary and sufficient conditions in order for Bracket (11) to satisfy the Jacobi identity, is done by using Mathematica [14]. In this appendix, we provide the detailed code along with the procedure required in order to obtain these constraints.

We need to include the "VariationalMethods" package in order to be able to compute functional derivatives. This is done with the command-line

Needs["VariationalMethods'"];

We then define the coordinate we are considering, which in this article correspond to x, y and z. This is done with the command-line r={x,y,z};

We define then the field variables of the system, ψ = (ρ, u, S 2 ):

Then, we proceed to the definition of a set of three functional derivatives, respectively F ψ , G ψ and H ψ which we, for sake of simplicity, denote dF, dG and dH respectively.

have to vanish. This leads to the necessary constraint ∂S 4 /∂u = 0. Choosing (F, G, H) = (S 2 , S 2 dr, u dr) and requiring that the Jacobi identity is satisfied leads to an additional constraint. However, in order to make the term proportional to (∂ x ρ) 2 ∂ 2 y u vanish, we have to impose

whose solution writes

where T 4 (S 2 ) is an arbitrary function. We choose (F, G, H) = (S 2 , xS 2 dr, S 2 dr). In order to make the term proportional to ∂ x ρ∂ y ρ∂ y u cancel, we must have

which, as ρ and S 2 are independent variables, leads to the necessary constraints T ′ 4 = 0 and

Computing J, one can verify that these constraints are sufficient to ensure that the terms J ⊥ and {F, {G, H} ⊥ } + {F, {G, H} } ⊥ + vanish such that Bracket (11) satisfies the Jacobi identity and hence is of Poisson type.


Thanh Nhan Vo

Karine Amis

Thierry Chonavel

Pierre Siohan

Bit loading in MIMO-PLC systems with the presence of interference

Keywords: Bit-loading, Power allocation, MIMO, Interference, Power Line Communication, Greedy based approach

Power line has become attractive for data transmission thanks to its simple implementation that does not require any additional wires. Windowed orthogonal frequency division multiplex (OFDM) with adaptive bit loading has been selected to satisfy spectrum mask constraints of current regulations. In broadband indoor power line communication (PLC) systems, multiple input multiple output (MIMO) techniques have been introduced to address the increasing demand for high data rates under the constraint of limited allocated bandwidth. Whereas the self inter-antenna interference can be dealt with on each subcarrier, both inter-carrier and inter-symbol interference can occur yielding to sub-optimal bit loading if not considered. In this paper, we extend to the MIMO case the low-complexity bit/power allocation algorithm based on the Greedy principle and taking interference into account, that was first applied in the SISO case in our previous work, called Reduced Complexity Algorithm (RCA). Moreover, we also consider this problem for two cases: optimum eigen beamforming and spatial multiplexing. Simulation results show the efficiency of the RCA in terms of throughput and computation cost for both cases.

I. INTRODUCTION

MIMO technique has been considered as a major key to increase the data rate in the next generation of broadband indoor power line communication systems. The HomePlug AV2 technology exploits the MIMO scheme to achieve higher data rates as well as a larger coverage [START_REF] Yong | An overview of the HomePlug AV 2 technology[END_REF]. MIMO-PLC is feasible since a protective earth (PE) wire is available in addition to phase (P) and neutral (N) wires. Hence, a scheme of MIMO with two transmitting antennas can be applied in the context of PLC. In interference-free OFDM systems, the bit loading problem is solved by the Water-filling algorithm [START_REF] Baccarelli | Novel efficient bit-loading algorithms for peak-energy-limited ADSL-Type multicarrier systems[END_REF], [START_REF] Papandreou | Bit and power allocation in constrained multicarrier systems: The single-user case[END_REF] for the continuous case and by Greedy based algorithms such as bit-adding and bit-removal [START_REF] Al-Hanafy | A new low-cost discrete bit loading using greedy power allocation[END_REF], [START_REF] Papandreou | A new computationally efficient discrete bit-loading algorithm for DMT applications[END_REF] for the discrete case. In SISO-OFDM systems with the presence of interference, the bit loading problem has been considered in [START_REF] Wen | Resource management techniques for OFDM systems with the presence of inter-carrier interference[END_REF], [START_REF] Chan | Multiuser spectrum optimization for discrete multitone systems with asynchronous crosstalk[END_REF]. Recently, an efficient Greedy based approach has been proposed in [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF] to solve the optimal bit/power allocation in SISO-PLC systems with the presence of interference resulting from an insufficient guard interval. Its achievable throughput is almost the same as the one obtained by the conventional Greedy approach. However, the complexity is significantly reduced.

In this paper, we consider the bit/power allocation problem for MIMO-PLC systems with the presence of inter-symbol and inter-carrier interference. The self inter-antenna interference is dealt with on each subcarrier. The key point is the transformation into an equivalent bit loading problem in SISO-PLC systems and then the application of the reduced complexity approach in [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF] with some small modifications involving power constraints to solve this problem.

The paper is organized as follows. In Section II, the bit loading problem in a 2x2 MIMO-PLC system with the presence of interference is introduced. In Section III, the transformation of this problem into an equivalent problem in SISO-PLC systems is shown. We then introduce the method of [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF] to solve the problem for both cases: optimum eigen beamforming and spatial mutiplexing. The simulation results and the complexity study are reported in Section IV. Finally, Section V is dedicated to conclusions and perspectives.

II. MIMO-PLC AND INTERFERNCE MODELS

In this paper, a 2x2 MIMO-PLC system is considered for the sake of simplicity. Let us consider a 2x2 MIMO-PLC system with L used subcarriers. On subcarrier m 0 , we have

y 1 (m 0 ) y 2 (m 0 ) = h 11 (m 0 ) h 12 (m 0 ) h 21 (m 0 ) h 22 (m 0 ) x 1 (m 0 ) x 2 (m 0 ) + I 11 (m 0 ) + I 12 (m 0 ) I 21 (m 0 ) + I 22 (m 0 ) + b 1 (m 0 ) b 2 (m 0 ) (1)
where y j (m 0 ), x j (m 0 ) and b j (m 0 ) are the received signal, the transmitted signal and the Gaussian noise with variance σ 2 j (m 0 ) at antenna j on subcarrier m 0 ; h ji (m 0 ) and I ji (m 0 ) are the channel frequency response and interference term from antenna i to antenna j on subcarrier m 0 . Rewritting (1) under matrix form, we obtain

Y(m 0 ) = H(m 0 )X(m 0 ) + I(m 0 ) + B(m 0 ) (2)
In SISO-PLC systems, the interference power on a given subcarrier depends on the signal power on all subcarriers [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF], [START_REF] Achaichia | OFDM/OQAM: A solution to efficiently increase the capacity of future PLC networks[END_REF]. Similarily, the interference power in MIMO-PLC at a given receiver antenna and on a given subcarrier depends on the signal power at all antennas and on all subcarriers [START_REF] Vo | Influence of Interference in MIMO Power Line Communication systems[END_REF].

In [START_REF] Vo | Influence of Interference in MIMO Power Line Communication systems[END_REF], the covariance matrix of the noise plus interference C I+B (m 0 ) has been written in the following form:

C I+B (m 0 ) = [W 1 P](m 0 ) + σ 2 1 (m 0 ) [W 2 P](m 0 ) [W * 2 P](m 0 ) [W 3 P](m 0 ) + σ 2 2 (m 0 ) (3
) where P = [P 1 (1), ..., P 1 (L), P 2 (1), ..., P 2 (L)] T denotes the power allocation vector; W 1 , W 2 and W 3 are the matrices of interference contributions, with size Lx2L (see [START_REF] Vo | Influence of Interference in MIMO Power Line Communication systems[END_REF] fore more details about these matrices) and [W i P] is a vector with size Lx1 resulting from the product of W i and P.

In the following, to alleviate the notation, the subcarrier index m 0 is omitted. If the MIMO equalizer is based on zeroforcing (ZF), the detection matrix is given by

A = H -1 (4)
Then, the signal of interest is

Z = AY = X + A(I + B) (5)
Let us denote A = a 11 a 12 a 21 a 22 . The SINR level after ZF equalizer is

s 1 = P 1 |a 11 | 2 [W 1 P] + 2 (a * 12 a 11 [W 2 P]) + |a 12 | 2 [W 3 P] + N 1 s 2 = P 2 |a 21 | 2 [W 1 P] + 2 (a * 22 a 21 [W 2 P]) + |a 22 | 2 [W 3 P] + N 2 where N i = |a i1 | 2 σ 2 1 + |a i2 | 2 σ 2 2 .
The bit loading is done based on the SINR level:

b i = log 2 1 + s i Γ , b i ∈ A, (i = 1, 2) (6)
where A denotes the set of allowable number of bits specified in IEEE P1901 standard [START_REF] Society | IEEE Standard for Broadband over Power Line Networks: Medium Access Control and Physical Layer Specication[END_REF], Γ is the SNR gap that models the practical modulation and coding scheme for a targeted symbol error rate:

Γ = 1 3 Q -1 ( SER 4 )
2 , where SER is the target Symbol Error Rate and Q -1 (x) is the inverse tail probability of the standard normal distribution [START_REF] Cioffi | A multicarrier prime[END_REF].

A. Bit-power relation

The bit-power relations are in the form

P1 = s1 |a11| 2 [W1P] + 2 (a * 12 a11[W2P]) + |a12| 2 [W3P] + N1 (7)
P2 = s2 |a21| 2 [W1P] + 2 (a * 22 a21[W2P]) + |a22| 2 [W3P] + N2 (8)
where, using [START_REF] Wen | Resource management techniques for OFDM systems with the presence of inter-carrier interference[END_REF], s i = (2 bi -1)Γ. We can see that the number of bits allocated to a given subcarrier depends not only on its power level but also on the power allocated to other subcarriers. Then, the bit-power relation problem is as follows:

Given a vector of number of bits b allocated to the subcarriers and on the antennas, we have to find the corresponding vector of allocated power P so that Eq.( 6) is satisfied for all subcarriers.

B. Bit loading problem

The bit loading problem consists in finding the vector of allocated power that optimizes the achievable throughput. For 2x2 MIMO-PLC systems, it can be written as follows:

                     max P L m=1 2 i=1 b i (m) L m=1 2 i=1 P i (m) ≤ P tot 2 i=1 P i (m) ≤ P max (m) (9)
This problem is quite difficult to solve because the number of bits allocated on any subcarrier depends on the power allocated on all subcarriers. Thus, we have to take into account the power allocated on all subcarriers to allocate the number of bits on any subcarrier. Fortunately, we addressed a similar problem in [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF] for SISO-PLC systems. Hence, in this paper, we aim to transform the bit loading problem (9) into a problem equivalent to the case of SISO-PLC systems. Then, the reduced complexity approach in [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF] can be reused with some small modifications to solve this problem.

III. BIT LOADING PROBLEM IN MIMO-PLC WITH THE PRESENCE OF INTERFERENCE

In this section, a solution of problem ( 9) is given. To this end, we transform it into an equivalent SISO problem. Let us consider the following notations

Λ = diag s 1 (1), ..., s 1 (L), s 2 (1), ..., s 2 (L) (10)
A 1 = diag |a 11 (1)| 2 , ..., |a 11 (L)| 2 , |a 21 (1)| 2 , ..., |a 21 (L)| 2 (11)
A 2 = diag a * 12 (
A 3 = diag |a 12 (1)| 2 , ..., |a 12 (L)| 2 , |a 22 (1)| 2 , ..., |a 22 (L)| 2 (12)
W = A 1 W 1 W 1 + 2 A 2 W 2 W 2 + A 3 W 3 W 3 (13)
N = N 1 (1), ..., N 1 (L), N 2 (1), ..., N 2 (L) T (14)
Then, rewriting Eq. ( 7), ( 8) for all L subcarriers in a vector form, we obtain

P = Λ WP + N ⇒ P = I 2Lx2L -ΛW -1 ΛN (16)
Eq.( 16) is the same as the one in [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF] (c.f. Eq.( 15)) except that the matrices are twice as big. Thus, we can extend the method obtained in the case of SISO-PLC in [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF] with some slight modifications upon power constraints to solve the bit-loading problem for MIMO-PLC in the presence of interference. We recall that in the conventional Greedy algorithm, which is a popular tool for discrete optimization problems [START_REF] Cormen | Introduction to Algorithms[END_REF], at every iteration and for each subcarrier, we have to calculate the required power to increase throughput by one bit. Then, we choose the subcarrier that requires minimum power to increase after verifying the constraints of the total power allocated on every subcarrier and on all subcarriers remain satisfied. Eq. [START_REF] Hashmat | A channel model for multi input multi output in-home power line networks[END_REF] shows that the computational complexity to calculate the power after increasing one bit on a subcarrier is that of a 2Lx2L matrix inversion. Thus, the total complexity for one iteration in the conventional Greedy algorithm is the product of the number of subcarriers that can be loaded up and the complexity of a matrix inversion. Due to its important complexity, the conventional Greedy algorithm cannot be used in practical systems.

To reduce complexity, in [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF] we have proposed the Reduced Complexity Algorithm (RCA). Its principle consists in the simplification of the cost function (which is used to determine the subcarriers to increase number of bits) in the Greedy procedure and the use of a judicious initial bit vector instead of a null one as in the conventional Greedy approach and the efficient selection of K subcarriers to load up at every iteration. The main advantage of the RCA is that it almost converges to the same achievable throughput as the one obtained with the conventional Greedy algorithm. In addition, the RCA significantly outperforms the conventional Greedy algorithm in terms of complexity. For the sake of further complexity reduction, for the bit vector initialization, we use the Equal Power Allocation (EPA), i.e. the total power P tot is identically allocated between the subcarriers and the antennas, instead of the Constant Power Water-Filling (CPWF) as in the original RCA in [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF]. More details about the RCA can be found in [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF].

Complexity

) GR + Init ≈ O(4β 1 U 3 ) N 1 << N s N 1 * O(4β 1 ) RCA ≈ O(U 3 ) N K << N s N K * O(1)
Table 1. Complexity per iteration, number of iterations and equivalent total number of matrix inversions.

The complexity comparison between the conventional Greedy, the Greedy with EPA initialized bit vector and the reduced complexity algorithms is illustrated in Table 1 where 1 << β 1 < β, U is the length of the power vector P and N s , N 1 , N K are the total number of iterations used in the standard Greedy algorithm, the Greedy algorithm with the EPA initialized bit vector and RCA respectively. In SISO systems, U is equal to L, the number of active subcarriers. In 2x2 MIMO systems, U is equal to 2L. For more details about how to calculate the complexity of the algorithms, interested readers can refer to [START_REF] Vo | Achievable Throughput Optimization in OFDM Systems in the Presence of Interference and its Application to Power Line Networks[END_REF]. We can see that the RCA reduces not only the number of iterations but also the computation cost per iteration, due to the decrease of the average number of subcarriers for which we have to calculate the value of the cost function in the Greedy procedure.

In the following, we consider the bit/power allocation problem in MIMO-PLC systems with the presence of interference for the cases of optimum eigen beamforming and of spatial multiplexing. Indeed, we can use the proposed solution in section III to solve the problem in both cases. However, we have to change the detection matrix A. We give here the formula to calculate A in both cases.

A. Optimum eigen beamforming

In this case, the eigen beamforming is used as it offers a good performance by adapting the transmission in an optimum way to the underlying eigen modes of the MIMO-PLC channel. The underlying MIMO streams are obtained by a Singlar Value Decomposition (SVD) of the channel matrix [START_REF] Paulraj | Introduction to Space-Time Wireless Communications[END_REF] where U, V are unitary matrices, i.e. U -1 = U H and V -1 = V H , H denotes the Hermitian operator and D is a diagonal matrix containing the single values of H.

H = UDV H (
In the case of eigen beamforming with precoding matrix V, the channel matrix H is replaced by HV in (2). Then, the detection matrix A can be expressed as [START_REF] Yong | An overview of the HomePlug AV 2 technology[END_REF]

A = HV -1 = D -1 U H . (18)

B. Spatial multiplexing

In this case, no precoding is done and then the detection matrix is A = H -1 .

IV. SIMULATION RESULTS

The simulation parameters are the following:

• Exploitable subcarriers defined in the IEEE P1901 standard [START_REF] Society | IEEE Standard for Broadband over Power Line Networks: Medium Access Control and Physical Layer Specication[END_REF] with sampling rate f s = 100 MHz, frequency shift ∆f = 24.414 kHz. • Channel model: SISO-PLC Class 2 channel of Tonello's model [START_REF] Tonello | A fitting algorithm for random modeling the PLC channel[END_REF] and 2x2 MIMO-PLC (same circuit using P-N, N-PE) of Hashmat model [START_REF] Hashmat | A channel model for multi input multi output in-home power line networks[END_REF]. • Noise model: 2x2 MIMO-PLC independent noise of Esmailian model [START_REF] Hashmat | Analysis and modeling of background noise for inhome MIMO-PLC channels[END_REF]. • Spectral mask constraint: P 1 (m) + P 2 (m) ≤ P 0 , ∀m.

The total power allocated on any subcarrier is less than P 0 , where P 0 = -50 dBm/Hz [START_REF] Society | IEEE Standard for Broadband over Power Line Networks: Medium Access Control and Physical Layer Specication[END_REF].

• A = {1, 2, 3, 4, 6, 8, 10, 12} (bits)
.

• Number of channel realizations: 1000. For the sake of simplicity, we only take into account the first 100 useful subcarriers out of 917 defined in IEEE 1901 standard in our simulation, i.e. L = 100. The calculation of covariance matrix C I+B is completely based on the calculation in [START_REF] Vo | Influence of Interference in MIMO Power Line Communication systems[END_REF] with the value of guard interval of 5.56 µs, the minimum allowable value defined in the IEEE P1901 standard. In this case, the effect of interference on the system performance is significant. We illustrate in Figures 1 and3 the achieved throughput and the power used with the EPA (in both cases: continuous bit-loading and discrete bitloading), the conventional Greedy algorithm (GR), the Greedy with judicious initial bit vector (GR + Init) and the RCA, respectively in the case of spatial multiplexing. Figures 2 and4 supply the same informations in the case of optimum eigen beamforming. In the RCA, the number of subcarriers chosen to simultanously increase the number of bits is set to 5. achievable throughput of the proposed RCA is little degraded when compared to that obtained with the GR algorithms and significantly improved as compared to the EPA algorithm.

Algorithm

The maximal shift of achievable throughput between RCA and the conventional GR is about 0.2% and in average, the RCA increases the achievable throughput by 17% w.r.t the EPA algorithm in both cases. In addition, we can see that the achievable throughput obtained by the RCA with optimum eigen beamforming is increased by about 4% (in average) as compared to the one of spatial multiplexing. Hence, as compared to the spatial multiplexing, the achievable throughput gain obtained by the optimum eigen beamforming is not really significant in 2x2 MIMO-PLC systems with significant interference resulting from an insufficient guard interval. The achievable throughput comparison when the total normalized power P tot = 50 and with spatial multiplexing is shown in Table 2.

Figure 5 shows the complexity of different algorithms that is represented by the total number of matrix inversions used to find the bit/power allocation in the case of spatial multiplexing. Note that the vertical axis is in logarithmic scale. The complexity of different algorithms with optimum eigen beamforming is almost the same when compared to the one illustrated in Figure 5. For both cases, the RCA has reduced the complexity by 97% as compared to the conventional GR algorithm. However, as mentioned above, the achievable throughput for both algorithms is almost the same. The runtime comparison when P tot = 50 is also shown in Table 2. As compared to the conventional Greedy algorithm, the RCA only degrades the achievable throughput by 0.1%. However, the run-time is reduced by about 35 times. It is shown that the RCA is efficient to solve the bit/power allocation in MIMO-PLC systems with the presence of interference. in presence of interference. To this end, we have transformed the MIMO-PLC problem into an equivalent SISO-PLC problem. Simulation results have clearly shown that the reduced complexity algorithm is also efficient for the problem of achievable throughput maximization in MIMO-PLC systems in presence of interference, i.e. the achievable throughput loss is negligible and the complexity is significantly reduced as compared to the Greedy solutions. It is also shown that the RCA algorithm outperforms the EPA algorithm. Therefore, it is a good candidate to solve resource management problems for single-user Windowed-OFDM systems in the presence of significant interference for both SISO and MIMO systems. Our future work aims to solve the bit/power allocation problem for the multi-user with/without the presence of interference in SISO/MIMO-OFDM systems.

Fig. 3 :

 3 Fig. 3: Total consumed power of different algorithms with spatial multiplexing.

Fig. 4 :

 4 Fig. 4: Total consumed power of different algorithms with optimum eigen beamforming.

Fig. 5 :

 5 Fig. 5: Complexity of different algorithms with spatial multiplexing.

Table 2 .

 2 Throughput and run-time comparison between different algorithms with spatial multiplexing.For both cases, the performance of GR, GR + Init algorithms in terms of achievable throughput is almost the same. The

		Throughput Number of Run-time
		(Mbits)	inversions	(s)
	EPA (continuous)	29.42	-	-
	Standard GR	29.12	1070	17.2
	GR + Init	29.13	207	2.7
	RCA	29.09	38	0.5
	EPA (discrete)	24.90	-	-


Roger Collet
email: roger.collet@ifsttar.fr

Matthieu De Lapparent
email: matthieu.de.lapparent@ifsttar.fr

Laurent Hivert
email: laurent.hivert@ifsttar.fr

Are French Households Car-Use Addicts? A microeconomic perspective

Keywords: Household car use, Rational Addiction, GMM, Elasticities

This article presents a microeconometric analysis of household car use in France. To feature car use dependence, the myopic and rational addiction models are estimated using panel data drawn from the French "Car Fleet" survey. Significantly, the assumption of rational addiction to car use cannot be rejected, and is even supported by a plausible estimate of the intertemporal rate of substitution (17%). Furthermore, the rational model yields realistic estimates of the fuel cost-and income-elasticities of household car use, respectively -0.23 and +0.10 for the short run, and -0.37 and +0.16 for the long run.

Introduction

Currently, the contribution of human activity to climate change and the related ecological issues are recognized, and a large international consensus has emerged about the urgency for making economic growth cleaner. In France, policymakers have set environmental targets to converge towards a more sustainable development, as shown by the energy law voted in 2005 by the French Parliament, which aims at dividing the 1990 national level of greenhouse gas emissions by 4 by the year 2050. To reach this objective, a particular effort in the transportation sector is deemed essential. Not only would attention be paid to the volume of emissions, but the perspectives of global peak oil and oil scarcity plead for the necessity to find new and/or more efficient ways to travel, as they should entail an increase in fuel prices in the future. Thus, the question of car use intensity emerges, given that current technology is mainly based on fossil resources. In a context of highly volatile and increasing fuel prices, it is of great importance to understand and to quantify to what extent households adapt their behaviors in terms of car use. Moreover, drivers can be viewed as constituting a "club" whose advantages have grown with their number, leading to a worsening of the social marginalization of non-driving persons. Indeed, a car provides users with a larger control of space and time, and allows a better accessibility to jobs, leisure places, health care, public amenities, etc. However, car users might develop some kind of automobile dependence in order to achieve their mobility, insomuch that several authors go so far as to compare car use to the consumption of tobacco or drugs. 1In this article, we propose a microeconometric analysis of the annual mileage traveled by French households with their personal cars. In the literature, this mileage is sometimes referred to as household "automobility", and we will use both expressions interchangeably throughout this paper. The analysis is achieved using the sample of households participating in each of the annual waves of the French "Car Fleet" survey from 1999 to 2001. This period is characterized by peak fuel prices in 2000. It has led households to adjust their car use behavior substantially, thus making fuel cost-elasticities of household automobility easier to measure accurately. Moreover, the panel layout of the data enables us to use dynamic specifications and to derive short and long-run effects. To feature the dependence of households on car use, the myopic and rational addiction models [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF]) are investigated. By applying these models, generally used to describe the consumption of cigarettes, alcohol or drugs, the addiction hypothesis regarding car use can be tested from the microeconomic point of view. To our knowledge, this article provides the first application of the rational addiction model to describe car use behavior, while "automobile dependence" is a major topic of research in transportation economics. The rest of the article is fashioned as follows. Section 2 presents the theoretical framework: the rational addiction model and its properties are examined. Section 3 discusses the econometric specifications and the estimation method. Section 4 provides a description of the 1999-2001 French "Car Fleet" panel dataset. Section 5 reports and debates the estimation results. The last section concludes.

Microeconomic Framework

The rational addiction model

In the literature dealing with habit formation, pioneering contributions considered individual preferences to be endogenously determined, assuming that the current utility level depends on past consumption levels [START_REF] Gorman | Tastes, habits and choices[END_REF][START_REF] Pollak | Habit formation and dynamic demand functions[END_REF][START_REF] Ryder | Optimum growth with intertemporally dependent preferences[END_REF]. In this framework, individuals are said to be myopic, since they totally ignore the effect of their current consumption on their future utility. In other words, individuals are making their choices period after period, disregarding the habit-forming property of consumption. [START_REF] Becker | A theory of rational addiction[END_REF] argue against the myopia hypothesis of individual behaviors, and have come up with the rational addiction theory, in which individuals are forward-looking agents and aware that consumption is habit forming. In line with [START_REF] Iannaccone | Addiction and satiation[END_REF], this theory describes an individual as "addicted" to a good if, all else equal, an increase in their past consumption yields a rise in their current consumption. The rational addiction model derived in [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF] examines the behavior of an individual, whose current utility U depends on the consumption amounts of two goods: a quantity X t of a composite good X, and a quantity C t of an addictive good C. The current utility also depends on a set of individual characteristics related to life-cycle and summarized by e t . The addictive and composite goods are different in that current utility also depends on the previously consumed quantities of C. According to the authors, these quantities are accumulated into an addictive capital stock S t , which is supplied at each current period by the past level of consumption C t-1 , so that S t = C t-1 + (1τ)S t-1 . The most commonly used expression to define this stock of habits assumes a complete decay from a period to the next one (τ = 1), thus writing S t = C t-1 . At period t = 1, let the rational individual look to maximize their discounted utility over an infinite lifetime horizon. Their program is given by: max

∞ ∑ t=1 B t-1 •U(C t ,C t-1 , X t , e t ), (1)
where B = (1 + ρ) -1 is the discount factor, and is the intertemporal rate of substitution. In addition, it is assumed that the interest rate of the economy is equal to ρ, and that the composite X is the equivalent of the money with a price conventionally set to 1. Maximization of the individual's utility is subject to an initial condition regarding the addictive good, and the individual is bound by the intertemporal budget equilibrium:

C 0 = C 0 ; A 1 = ∞ ∑ t=1 B t-1 (X t + P t C t ), (2)
where A 1 is the net present value of the individual's wealth, and where P t is the price of the addictive good at date t. Let the utility function be concave and quadratic, and such as:

U(C t ,C t-1 , X t , e t ) = α C C t + α S C t-1 + α X X t + α e e t + α CC 2 C 2 t + α SS 2 C 2 t-1 + α XX 2 X 2 t + α ee 2 e 2 t (3) +α CS C t C t-1 + α CX C t X t + α Ce C t e t + α Xe X t e t .
Under the previous hypotheses, the optimal demand function for the addictive good is derived from solving the maximization program of the individual. It is found to be a function of the nearest past and future consumptions, the current nominal price P t and the characteristics e t :

C t = δC t-1 + δ BC t+1 + δ 1 P t + δ 2 e t , (4)
where

θ = -((α CC α XX -α 2 CX ) + B(α SS α XX )) -1 (α XX α CS ). (5)
For α CS strictly positive, the past and current consumptions of C are said to be complementary. In that case, the current marginal utility of the addictive good, say U C , is an increasing function of C t-1 :

U C = dU dC t = α C + α CC C t + α CS C t-1 + α CX X t + α Ce e t . (6)
The higher the level of C t-1 and the value of α CS , the higher the marginal utility of the addictive good. By analogy with the learning-by-doing, the individual draws all the more utility from the consumption of the addictive good as they "practiced" it in the near past, and as the "learning speed" α CS is high.

The intertemporal complementarity of the consumptions of C is the origin of addiction. It implies θ > 0 in (5). The larger the value of θ , the greater the level of addiction. It can be noticed that the static and autoregressive consumption models are particular cases of the dynamic demand function (4). For θ = 0, the demand function neither depends on past nor on future levels of consumption, and the static case emerges. When θ > 0 and B = (1 + ρ) -1 drops to 0 (that is, for an infinite preference for the present), the demand function (4) takes the form of a first-order autoregressive process. In that case, the individual is not a forward-looking agent. [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF] defines such a behavior as myopic addiction. In any other situation where θ and B are strictly positive, the addiction behavior is said to be rational. As both these parameters can be estimated, it is possible to determine which formulation is the most relevant from the empirical perspective.

Elasticities in the rational addiction model

Without loss of generality, let model (4) be rewritten as:

C it = θC it-1 + θ BC it+1 + X it β + ε it , (7)
where the subscript i identifies the individual, and X is now a set of exogenous covariates which includes the price for the addictive good. The error term ε summarizes all the unobserved factors from the modeler's point of view. The impacts on the current consumption resulting from a variation in the past and future consumptions can be deduced from the characteristic roots of the homogeneous equation of (7). These are given by:

ϕ 1 = 1 - √ 1 -4θ 2 B 2θ ; ϕ 2 = 1 + √ 1 -4θ 2 B 2θ . (8)
In ( 8), ϕ 1 and ϕ -1 2 measure the effect on C t induced by a shift in C t+1 and in C t-1 respectively. The elasticities in the rational addiction model can be expressed as functions of these characteristic roots. Let X k denote the k rd continuous covariate of X, and β k the related parameter in model ( 7). The short-and long-run elasticities of the demand for C with respect to X k evaluated at the sample averages X k and C, and respectively denoted e SR C/X k and e LR C/X k , are given by (Becker 1996, p.113):

e SR C/X k = β k (θ (1 -ϕ 1 )ϕ 2 ) × X k C , (9)
e LR C/X k = β k (θ (1 -ϕ 1 )(ϕ 2 -1)) × X k C . ( 10
)
It can be noted that the elasticities resulting from the myopic addiction model logically emerge as particular cases of ( 9) and ( 10), for ϕ 1 = 0 and ϕ -1 2 = θ . The elasticity that stems from the static demand equation, which is the same in both the short and the long run, arises for θ = 0.

Testing for rational addiction

Testing the rational addiction theory has been performed on various data that pertain to different topics of research. Many applications dealt with drug consumption to tackle public health issues. For instance, [START_REF] Baltagi | The econometrics of rational addiction: The case of cigarettes[END_REF], [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF], [START_REF] Gardes | Evidence on addiction effects from households expenditure surveys: the case of the polish panel[END_REF][START_REF] Gardes | Evidence on addiction effects from households expenditure surveys: the case of the polish panel[END_REF][START_REF] Tiezzi | An empirical analysis of tobacco addiction in italy[END_REF] used addiction models to explain the consumption of tobacco. Grossman et al. (1998), [START_REF] Bentzen | Rational addiction and alcohol consumption: evidence from the nordic countries[END_REF] and [START_REF] Baltagi | Rational addiction to alcohol: panel data analysis of liquor consumption[END_REF] applied them to model the demand for alcohol. Van Ours (1995), and Grossman and Chaloupka (1998) tested the relevance of addiction models to explain the consumption of opium and cocaine. However, as mentioned by [START_REF] Becker | A theory of rational addiction[END_REF], there is no need to express a biological dependence to be considered as addicted to a good. There are other areas of application than tobacco or drug consumption. For instance, [START_REF] Mobilia | Gambling as a rational addiction[END_REF] tested for addiction to gambling, [START_REF] Cawley | Rational addiction, the consumption of calories, and body weight[END_REF] focused on the consumption of calories, [START_REF] Villani | Rational addiction in the arts[END_REF] dealt with addiction to art, [START_REF] Cameron | Rational addiction and the demand for cinema[END_REF] and [START_REF] Sisto | Rationally addicted to cinema and tv? an empirical investigation of italian consumers[END_REF] applied addiction models to the demand for cinema. To our knowledge, the rational addiction model has not been applied so far to describe car use, while "automobile dependence" is a major topic of transportation research. However, it is advisable to be cautious when interpreting the results of the rational addiction model. Indeed, while the effects of past and future consumptions are often found to be significantly positive in the literature, the resulting intertemporal rates of substitution may take unlikely values. All in all, a rule of thumb would be to accept the assumption of rational addiction also based on a "reasonable" rate.2

Modeling and estimation

Notations

Unless explicitly stated, the following notations will be used all throughout the rest of our study. In addition to the subscripts i and t which identify respectively the household and the period, we also consider a subscript v which is related to a specific car owned by i at date t. These three subscripts are used all together to index the following variables:

• KM: the annual mileage converted into kilometers;

• FE: the average fuel efficiency in kilometers per liter;

• DP: the price for one liter of diesel oil in constant (2000) Euros;

• PP: the price for one liter of premium-petrol in constant (2000) Euros;

• FOC: the fuel operating cost per 100 kilometers. It is defined as

FOC ivt = 100 × (FP ivt /FE ivt )
, where FP stands here for the price of the appropriate type of fuel (DP or PP);

• AD: the age of the car if it is diesel-powered;

• AP: the age of the car if it is petrol-powered.

Additional variables that describe the characteristics of the household will also be used. They are detailed in the next subsection.

The econometric model

The annual mileage of a car v that is owned by household i at date t is modeled as:

KM ivt = ∑ s β 0s 1(R it = s) + β 1 FOC ivt + β 2 AD ivt + β 3 AP ivt + ε ivt . (11)
The error term ε ivt is assumed to be drawn from a zero-mean normal distribution. The model intercept β 0 is here differentiated according to the household residential location: either Paris city (R it = 1), or the inner suburbs of Paris (R it = 2), or the outer suburbs of Paris (R it = 3), or the Provinces (R it = 4). The automobility of household i at date t is obtained by totalizing the annual mileages of the cars it owns on this date. Thus aggregated, the following model emerges:

KM it = ∑ v KM ivt = ∑ s β 0s 1(R it = s)NC it + β 1 ∑ v FOC ivt ( 12
)
+β 2 ∑ v AD ivt + β 3 ∑ v AP ivt + ∑ v ε ivt ,
where NC refers to the number of cars. The previous specification is enlarged by introducing additional variables to allow for a better control of household heterogeneity: the household annual income in constant (2000) Euros, three dummy variables describing the age class of the household head ([18-39], [40-65], >65), the number of adults (except the head), the number of adults in employment, the number of adult women, the number of driving-license holders, and the number of children in the household. Finally, according to the type of addiction model -myopic or rational -to be estimated, the specification also includes the past and future household automobility (KM t±1 ) as covariates.

Selectivity, heteroscedasticity and endogeneity

The addiction models require the dependent variable not be censored. In fact, the survey from which the data are drawn contains households that do not own cars. The annual mileage for these non-motorized households is zero, corresponding in microeconomics to a corner solution. The automobility models have been estimated using the subsample of households that declared to own at least one car in 2000. However, excluding the non-motorized households may result in a sample selection problem. This issue has been controlled by applying a two-step estimation procedure [START_REF] Heckman | Sample selection bias as a specification error[END_REF]. In a first step, a dichotomous Probit model has been estimated to explain household car ownership. Then, the results from this so-called "selection model" have been used to estimate the current individual inverse Mills ratios λ it , say λit . The latter has in turn been introduced as covariate into the household mileage equation ( 12), which has been estimated using the subsample of car-owner households. This corrects for the potential selection bias by capturing the correlation in the error terms between the ownership-based selection and the automobility models. Testing for selection bias is easily carried out by checking whether the estimated coefficient that weighs λit , say βλ , is statistically different from zero.

This correction method of sample selection is not without some difficulties regarding inference. Actually, the introduction of λit into the automobility model generates heteroscedasticity, as it makes the variance of errors in (12) depend on the individual covariates used in the ownership Probit model [START_REF] Heckman | Sample selection bias as a specification error[END_REF]. Correcting the sample selection is not the only reason why heteroscedasticity has to be taken into account. The structure of the proposed model is also a source. Indeed, the summation of car mileages to compute household automobility is likely to produce heteroscedasticity, related to the car ownership level. Thus, a robust estimation method allowing for a general form of heteroscedasticity should be applied.

By involving simultaneously both the lagged and the forwarded dependent variables as covariates, the dynamic specification of the rational addiction model makes them necessarily endogenous, even assuming the serial independence of individual errors. In addition, the errors are likely to be serially correlated due to an unobservable time-invariant household-specific factor η i assuming that ∑ v ε ivt = ε it = η i +u it . Therefore, KM it±1 and ε it are correlated variables, and estimating the automobility models by means of ordinary least squares would yield biased estimates.

Estimation strategy and tests

A solution resides in turning to estimators that use instrumental variables. Among those that exist, the 2SLS estimator has almost but not all the desired properties. Although convergent, it is not consistent in the presence of heteroscedastic error terms.

It may however be used to test for existence of heteroscedasticity [START_REF] Breusch | A simple test for heteroscedasticity and random coefficient variation[END_REF]. If the test rejects the homoscedasticity assumption, the Generalized Method of Moments (GMM hereafter) should then be applied [START_REF] Hansen | Large sample properties of generalized method of moments estimators[END_REF]. However, a condition for applying the GMM is to use a set of "good" instruments. They are required to be orthogonal to the estimation residuals and sufficiently correlated with the endogenous variables to be instrumented. Both properties have therefore to be examined. To that extent, the tests proposed by [START_REF] Hansen | Large sample properties of generalized method of moments estimators[END_REF] and Bound et al. (1995, BJB hereafter) are implemented. The readers are referred to [START_REF] Baum | Instrumental variables and gmm: Estimation and testing[END_REF] for a detailed presentation of these methods.

Data and descriptive statistics

Data source

Data are drawn from the French "Car Fleet"3 panel survey, which is achieved annually since the mid 1980s by the private pooling institute TNS-Sofres4 . The survey aspires to a better knowledge of several dimensions of automobile demand in France, especially car ownership and use. It depicts with a great level of detail many of the attributes of the cars owned by the households, as well as many of the household characteristics. A nationally representative sample of 10,000 households is surveyed each year. The panel is rotating: about one third of the sample is renewed each year. Such a methodology allows for a longitudinal follow-up of some households for at least 3 years. In the present paper, we focus on households that were surveyed over the period 1999-2001, and we estimate the addiction models for t = 2000.

Descriptive statistics

We have identified 3010 households continuously present in the 1999, 2000 and 2001 waves of the survey. On annual averages, slightly less than 20 percent of them have no car, about 50 percent have one car, slightly more than 25 percent have two cars and about 5 percent have three cars or more.

The average automobility of households is monotonically decreasing over the covered period, from 15,610 kilometers in 1999 down to 14,826 kilometers in 2001. Excluding non-motorized households to cancel out the decision related to car ownership, we observe the same decreasing trend: the average household mileage also declines monotonically from 19,279 kilometers to 18,189 kilometers over the period. In the Appendix, Table A1 reports the annual descriptive statistics related to the characteristics of the sampled households.

In 1999, the households in the panel report a total of 3552 cars. There is a total of 3576 cars in 2000 and a total of 3605 cars in 2001. Based on this sample, Table A2 in the Appendix provides the descriptive statistics of the car attributes and mileages. In 2000, the average car is almost 6.8 years old. According to the type of engine, petrolpowered cars are about two years older on average than diesel-powered cars. This difference results from the very active dieselization trend in the French car fleet [START_REF] Hivert | Dieselisation and the «new dieselists» behaviour: recent developments in the french car fleet[END_REF]. In 1980, diesel cars represented less than 5 percent of the total fleet in France. This share has continuously increased to reach about 15 percent in 1990, 30 percent in 1995, 35 percent in 1999 and 40 percent in 2001 5 . In accordance with these figures, the proportion of diesel cars is also increasing in our data, from 35 percent in 1999 to 38 percent in 2001. Dieselization partially explains the improvement over time in the average energy efficiency of vehicles, since diesel cars consume less fuel on average than petrol cars of about 0.9 liters per 100 kilometers (Table A2). It also derives from the improvement in the fuel efficiency for both types of car over time: globally in the data, the average vehicle consumed 7.33 liters of fuel per 100 kilometers in 2001 to 7.44 liters initially in 1999.

After a period of low fuel prices during the 1990s in France, the year 2000 marked an episode of significant rise: in 1999, due to the decision of various oil-producing countries (including OPEC) to limit the production, the price of the crude oil barrel has soared. On annual averages, the price per liter of premium-petrol in France raised from A C0.98 in 1999 to A C1.14 in 2000 (from A C0.69 to A C0.85 for diesel-oil). But this increase has been short-lived, because during 2000, the production of oil increased again, resulting in a decline in fuel prices for the following year: in 2001, the price for one liter of premium-petrol dropped to A C1.09 (A C0.80 for diesel-oil). However, the fuel price trend over 1999-2001 still represents an increase for both premium-petrol (+11%) and diesel-oil (+16%) (Table A3, in the Appendix).

These differences in fuel price and energy efficiency partially explain the more intensive use of diesel cars. These have covered 17,085 km on average in 2000 to "only" 10,412 km for petrol cars. Between 1999 and 2001, average mileage decreased by 1600 km and 840 km for diesel and petrol cars respectively (Table A2). Thus, the decrease has been higher for diesel-powered cars than for petrol-powered cars. One reason is that the diesel-oil price per liter has increased faster than the premium petrol price.

Results

In the Appendix, Table A4 shows the estimates of the selection model, which is the first stage of our modeling approach. It simply models the household probability to own at least one car in 2000. As mentioned earlier, these estimates are used to compute the vector of individual correction factors λ , which is introduced into the automobility model to control for selection. The estimates for addiction models are reported in Table 1. In this section, the myopic and the rational addiction models are compared so as to draw conclusions about the intertemporal dynamics of household automobility demand. Then, we mainly focus on the fuel cost and income effects, and derive the related elasticities of household car use. But first of all, we examine how the GMM estimations has performed.

Fitting properties and selectivity

Both myopic and rational addiction models present good fitting properties on our disaggregate data. The R 2 statistics are respectively 0.58 and 0.67, and the Fisher statistics show that the set of explanatory variables is relevant to explain the dynamics of household automobility. The Breusch-Pagan statistic after estimating the models by 2SLS rejects the homoscedasticity assumption, and justifies resorting to the GMM estimator as IV technique.

The instruments that were used to implement the GMM are the set of current exogenous covariates (the "included" instruments) and a set of past and future household characteristics (the "excluded" instruments). In both models, the Hansen test cannot reject the null hypothesis that residuals and instruments are orthogonal. Moreover, the BJB test concludes to acceptance of the alternative hypothesis that excluded instruments are jointly significant in explaining the endogenous covariates. This is supported by the Anderson-Rubin test,6 which suggests that the endogenous covariates are not too weakly explained by the "excluded" instruments.

Both addiction models agree not to conclude to significance of the correction factor λ . Indeed, the null hypothesis H0 : β λ = 0 cannot be rejected at the 5 percent level, meaning that the selection of car owner households in 2000 to estimate the automobility models has not been a significant source of bias. Finding no evidence for selection bias should not be surprising here. It only means that the ownership and automobility models are specified in such a way that their errors are not significantly correlated. In other words, given that the explanatory variables used in the household ownership model are also included in the automobility models, the unobserved factors impacting ownership are found not to influence significantly household car use.7

Addiction and intertemporal rate of substitution

While the results substantiate the addiction hypothesis regarding household automobility behavior, the rational addiction model emerges as the most relevant to describe the type of car use dependence captured by the data. Indeed, the myopic model confirms the significant effect of households' past annual mileage on their current automobility. But on the other hand, the rational model rejects the behavioral myopia assumption, as the parameter that weighs the forwarded annual mileage is also statistically significant. Therefore, households are forward-looking agents in setting their automobility demand. Furthermore, the rational addiction model yields an intertemporal rate of substitution of about 17 percent, which is a rather plausible value. 2

Cost and income elasticities of automobility

The estimates related to the fuel operating cost have the expected negative signs in both models, but their levels of significance slightly differ. It is significant at the 1 percent level in the myopic model while it is just significant at the 5 percent level in the rational model.

Elasticity measures of household automobility with respect to fuel operating cost are reported in Table 2. The myopic and rational addiction models agree upon measuring the short-run cost-elasticity at -0.22. Because of a larger intertemporal perspective in the rational addiction model, the long-run cost-elasticity is higher (-0.37) than in the myopic model (-0.31). A raise by (constant 2000) A C1 in the cost to achieve 100 km entails an estimated decrease in the annual mileage per car by 380 km in the short run, and 623 km in the long run (Table 3).

Estimating cost-sensitivities of car use has been the topic of many earlier works. To cite a few, such figures were estimated and discussed in [START_REF] Hensher | The demand for vehicle use in the urban household sector: theory and empirical evidence[END_REF], [START_REF] Oum | Concepts of price elasticities of transport demand and recent empirical estimates: an interpretative survey[END_REF], [START_REF] Eltony | Transport gasoline demand in canada[END_REF], [START_REF] Rouwendal | An economic analysis of fuel use per kilometre by private cars[END_REF], [START_REF] Johansson | Measuring the long run fuel demand of cars: separate estimations of vehicle stock, mean fuel intensity, and mean annual driving distance[END_REF], and opposite, household income is shown to determine significantly both car ownership and use. If it has not been included in the explanatory variables, the selection correction factor would have been significant. 1). 95% confidence intervals in brackets. [START_REF] Berri | Étalement urbain et croissance de la circulation automobile des ménages en france. projection aux horizons 2010 et[END_REF]. [START_REF] Graham | The demand for automobile fuel: a survey of elasticities[END_REF] collected many existing results in the literature dealing with car use sensitivities to fuel price, which can be expected to be close to our fuel cost-elasticities. In their survey, the work of [START_REF] Goodwin | A review of new demand elasticities with special reference to short and long run effects of price changes[END_REF] is mentioned: based on four elasticities drawn from empirical studies in the 1980s, the author reported an average fuel price sensitivity of automobile traffic of -0.16 for the short run, and -0.33 for the long run. [START_REF] Goodwin | Elasticities of road traffic and fuel consumption with respect to price and income: a review[END_REF] updated these metaanalytic results using empirical works published in the 1990s and in the early 2000s. They reported an average sensitivity of car use with respect to fuel price of -0.10 for the short-term and -0.30 for the long-term.

In our application, the most relevant model to be compared with the literature is probably the myopic addiction model, since short-and long-run elasticities are usually derived from first-order autoregressive specifications. In this model, the long-run fuel cost-elasticity of household automobility, estimated at -0.31, thus emerges as a very plausible value. Relatively, the short-run elasticity obtained from this model (-0.22) could seem high, meaning a fast convergence to the long-term equilibrium. In other words, French households would adapt quickly their automobility to a change in the fuel operating cost. This conclusion is supported by Graham and Glaister (2002, p.22, fig.1). Indeed, the authors reported price-elasticities of the demand for gasoline for a set of western countries: comparatively, France presents one of the highest sensitivities in the short run and one of the lowest in the long run. The short-long ratio is about 3/4, while it is clearly below 0.5 for the other comparable countries (including Germany, United Kingdom, Austria and Canada).

Table 4 reports the sensitivity of car annual mileage to the fuel price, by type of engine. The elasticity of annual mileage for petrol-powered cars with respect to premiumpetrol price is estimated at -0.32 in the short run and at -0.52 in the long run. Regarding diesel-powered cars, these elasticities with respect to the diesel-oil price are -0.13 and -0.21 respectively.8 Thus, the use of petrol cars is about 2.5 times more sensitive to fuel price variations than diesel cars, both in the short and the long runs. These Notes: Elasticities evaluated at the average petrol and diesel cars in 2000, using the estimates of the rational addiction model and the fuel prices in 2000 (see Footnote 8 for details). 95% confidence intervals in brackets.

short-run fuel price elasticities of car use are very close to those recently published by Calvet and Marical (2012, Tab. 1). Based on the 2006 French family budget survey, the figures reported by the authors are -0.11 for diesel-powered cars, and -0.35 for petrol-powered cars. However, these estimates are found not to be significant, whereas they are in our study.

The estimated coefficient that pertains to the household annual income is significant and has the expected positive sign in both addiction models. Table 2 also reports the short-and long-run elasticities of households' automobility with respect to their income. Whichever addiction model is considered, these elasticities take similar values: about +0.10 in the short run and about +0.16 in the long run. Using the results from the rational addiction model, an increase by (constant 2000) A C1000 in the household annual income yields an increase in its automobility by about 79 km in the short run, and 131 km in the long run (Table 3). In [START_REF] Hensher | The demand for vehicle use in the urban household sector: theory and empirical evidence[END_REF], the income-elasticities of car use for households living in the Sydney urban area ranged from +0.05 to +0.14 according to the car ownership level. Mentioned for comparison in their article, [START_REF] Greene | The influence of the price of gasoline on vehicle use in multivehicle households[END_REF] and [START_REF] Mannering | A dynamic empirical analysis of household vehicle ownership and utilization[END_REF] respectively estimated this elasticity at +0.13 and +0.11 for the United States. Thus, our elasticity figures do not differ much from these references. Nevertheless, since income is a key factor of household car ownership, income-elasticities of automobility should vary widely depending on whether the level of car equipment is held constant, as in our study, or not.9

Conclusions

This article focuses on modeling the annual mileage covered by households with their personal cars, that is, their "automobility". It sheds new light on the car dependence issue. The results that are presented put the emphasis on the rational addiction model proposed by [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF], which had not been applied to automobile data so far. Uncommon within our context, this model lives up to expectations when applied to describe the empirical automobility behavior of French households in 2000. Indeed, the assumption of addiction to car use is not statistically rejected. Therefore, assertions that relate car use to an addictive consumption in the transportation literature are reinforced by the microeconomic point of view. Then, the rational version of the addiction model proves to be more relevant than the myopic version. Indeed, both past and future household car use are significant in explaining current automobility in the rational addiction model. Moreover, the intertemporal rate of substitution is estimated at 17 percent, which is a likely value. Such results show that the household behavior is consistent with a theoretical intertemporal optimization scheme. This conclusion stands our work apart from earlier dynamic studies, which may have missed an important point in explaining car use demand. Indeed, models based on a first-order autoregressive specification are useful for deriving short and long-run elasticities, but they also require the individuals to be myopic as it regards the future. Nonetheless, the myopic model is also reported to make comparisons with other studies.

In France, the fuel operating cost-and income-elasticities of household automobility derived from the rational addiction model coincide with expectations. The respective estimations are -0.23 and +0.10 for the short run, -0.37 and +0.16 for the long run. These figures do not diverge from existing results that were reported in the transportation literature. According to the type, petrol cars are more sensitive than diesel cars to a change in fuel price. For diesel cars, the elasticity of annual mileage with respect to the diesel-oil price is estimated at -0.13 in the short run, while that for petrol cars with respect to the premium-petrol price is measured at -0.32. The long-run elasticities are found to be about 1.6 times higher. Our elasticity results are sensible whichever addiction model is applied, myopic or rational. As it regards the latter, it strengthens our recommendation to apply it when data are available, in order to draw further conclusions related to car dependence and intertemporal substitution. 610 13,752 15,193 13,240 14,826 13,371 Source: 1999-2001 French Car Fleet panel (3010 households).

Table 1 :

 1 Estimates of the myopic and rational addiction models applied to household automobility

	Addiction Model	Myopic	Rational
	Covariates	Coefficients	Coefficients
	Dynamics
	Past automobility (KM it-1 )	0.307 * * *	0.346 * * *
	Future automobility (KM it+1 )	-	0.295 * * *
	Economic factors
	Fuel operating cost per 100 km (∑ v FOC ivt )	-362.68 * * *	-223.29 * *
	Annual income (thousands of cst. 2000A C)	87.37 * * *	46.37 * *
	Number of cars (NC it ) for households living in:
	Paris-city (R it = 1)	8813.36 * * *	5514.74 * * *
	the inner suburbs of Paris (R it = 2)	8540.45 * * *	4888.76 * * *
	the outer suburbs of Paris (R it = 3)	10810.79 * * *	6310.44 * * *
	the Provinces (R it = 4)	11038.71 * * *	6189.06 * * *
	Age of cars:
	Diesel cars (∑ v AD ivt ) Petrol cars (∑ v AP ivt )	-24.75 -191.68 * * *	-12.39 -104.05 * *
	Household characteristics
	Number of driving license holders	2138.45 * * *	1324.65 * *
	Number of adults (except the household head)	510.27	409.90
	Number of employed persons	815.36 * *	286.32
	Number of women	-2041.94 * * *	-1158.26 * *
	Number of children	-174.64	-297.28
	Age of the household head (ref: [18-40[)
	[40-65[	-1617.14 * * *	-1376.77 * *
	≥65	-2648.16 * * *	-1608.01 * *
	Selection correction factor λ	2430.34	1555.84
	Intercept	-1210.78	-954.01
	Significance levels:

* * * 1%; * * 5%; * 10%. Table continued on next page.

Table 1 :

 1 Estimates of the myopic and rational addiction models applied to household automobility (continued from previous page) Notes: p-values in parentheses. Estimation using the subsample of 2631 households which described at least one car in 2000. Dependent variable: KM in 2000 (Section 3). GMM estimation (2 covariates related to age-and type-ambiguous cars not reported).

	Addiction Model	Myopic	Rational
	Intertemporal rate of substitution (ρ)	+∞	17.04%
	R 2	0.58	0.67
	Fisher	F(19, 2611) : 82.6 (0.00)	F(20, 2610) : 115.5 (0.00)
	Breusch-Pagan, after 2SLS	χ 2 (40) : 777.0 (0.00)	χ 2 (40) : 878.8 (0.00)
	Hansen	χ 2 (21) : 30.46 (0.08)	χ 2 (20) : 23.22 (0.28)
	BJB	F(22, 2590) : KM t-1 : 3.80 (0.00)	F(22, 2590) : KM t-1 : 3.80 (0.00) KM t+1 : 4.51 (0.00)
	Anderson-Rubin	χ 2 (22) : 40.47 (0.01)	χ 2 (22) : 40.47 (0.01)

Table 2 :

 2 Short and long-run elasticities of household automobility with respect to the fuel operating cost and income

	Addiction model	Myopic	Rational
	Fuel operating cost
	Short run	-0.22 [-0.34; -0.10]	-0.23 [-0.41; -0.05]
	Long run	-0.31 [-0.47; -0.16]	-0.37 [-0.72; -0.08]
	Income
	Short run	+0.11 [+0.06; +0.16]	+0.10 [+0.03; +0.16]
	Long run	+0.16 [+0.09; +0.23]	+0.16 [+0.06; +0.29]
	Notes: 95% confidence intervals in brackets.

Table 3 :

 3 Short and long-run marginal effects on household automobility

	Horizon	Short run	Long run
	Effect on household automobility induced by:
	• a (constant 2000) A C1000 raise in the	+79 km	+131 km
	household annual income	[+20; +125]	[+39; +223]
	• a (constant 2000) A C1 raise in the fuel	-380 km	-623 km
	cost per 100 km, per car	[-676; -79]	[-1194; -138]
	Notes: Evaluation from the rational addiction model (Table

Table 4 :

 4 Elasticities of annual mileage of cars with respect to fuel prices

	Horizon	Car type (fuel used)	Elasticity
		Diesel car	-0.13
	Short-run	(diesel-oil) Petrol car	[-0.19; -0.10] -0.32
		(premium-petrol)	[-0.46; -0.25]
		Diesel car	-0.21
	Short-run	(diesel-oil) Petrol car	[-0.40; -0.14] -0.52
		(premium-petrol)	[-1.03; -0.36]

Table A2 :

 A2 Descriptive statistics of car characteristics -2001 French Car Fleet panel. Notes: All the personal cars described by the households of the panel. Statistics for fuel-ambiguous cars (2% of the car sample for 1999, 1% for 2000 and 2001) not reported.

	Year	1999		2000		2001
		Mean	St. dev.	Mean	St. dev.	Mean	St. dev.
	Repartition
	Diesel cars	0.34	0.23	0.36	0.23	0.38	0.24
	Petrol cars	0.64	0.23	0.63	0.23	0.61	0.24
	Mileage (in km)
	Diesel cars	17786	9014	17085	9099	16187	8662
	Petrol cars	11001	6564	10412	6138	10157	6172
	Age (in years)
	Diesel cars	5.19	4.01	5.52	4.28	5.48	4.38
	Petrol cars	7.22	5.70	7.50	5.86	7.56	5.99
	Energy consumption (L/100 km)
	Diesel cars	6.82	1.29	6.80	1.49	6.74	1.35
	Petrol cars	7.74	1.63	7.72	1.58	7.67	1.56
	Harmonic mean of the vehicle fuel efficiency (km/L)
	Diesel cars	14.66	14.71	14.84
	Petrol cars	12.91	12.95	13.03
	Number of observa-tions	3552 cars	3576 cars	3605 cars
	Source: 1999

Table A3 :

 A3 Fuel prices at filling stations in France Calculations from the yearbooks of the French professional comity of oils.

	Year	1999	2000	2001
	Diesel-oil (A C/Liter)	0.6890	0.8461	0.7958
	Premium-petrol (A C/Liter)	0.9825	1.1380	1.0877
	Source:

Table A4 :

 A4 Selection model estimates -Car ownership Probit model Notes: Probit estimation for wave 2000 of the panel, 3010 households. Dependent variable: Y it = 1 if the household owned at least one car (2631 cases), 0 otherwise (379 cases).

	Covariates	Coefficients	t-values
	Household annual income (thousands of constant 2000A C)	1.42	4.91
	Age of the household head (reference: [40 ; 65[)
	[18; 40[	-0.10	-0.80
	> 65	0.22	2.19
	Household location (reference: Paris-city)
	Inner suburbs of Paris	0.68	3.81
	Outer suburbs of Paris	1.64	8.10
	The Provinces	1.65	11.13
	Number of driving license owners	1.50	20.22
	Number of adults (except the head)	0.07	0.91
	Number of employed persons	0.11	1.43
	Number of women	-0.36	-4.05
	Number of children	0.15	2.44
	Intercept	-2.63	-12.99

[START_REF] Dupuy | From the "magic circle" to "automobile dependence": measurements and political implications[END_REF] and[START_REF] Wickham | Scenarios for a sustainable society: Car transport and the sociology of embedded technologies, final report of the project[END_REF], for example.

Remember the assumption in the rational addiction model that the intertemporal rate of substitution is equal to the interest rate of the economy. Arbitrarily, one can consider that a plausible rate should range between 0% and 25%. For instance,Becker (1996, p.103) reports a rate of 15% from applying the model to the consumption of cigarettes, arguing that it is a "quite reasonable" value.

"Car Fleet" is the literal translation from French of the original name of the survey, "Parc Automobile".

"Sofres" is the acronym for "Société Française d'Enquête par Sondages".

This phenomenon has even continued at the same pace during the 2000s: this proportion was about 50% in 2005 and 55% in 2008.

[START_REF] Anderson | Estimation of the parameters of a single equation in a complete system of stochastic equations[END_REF].

For example, variables as car parking facilities are likely to impact household car ownership. But while they are captured by the error term in the ownership model, there is no reason that such variables also determine household car use, resulting in the non significance of the selection correction factor. At the

Example: in 2000, the fuel efficiency of the average diesel car was 14.71 km per liter of diesel-oil, the average mileage for this type of car was 17,085 km (TableA2), and the price for one liter of diesel-oil was A C0.8461 (TableA3). Moreover, the long-term marginal effect of the fuel operating cost per 100 km on automobility is estimated at -623 km by car (Table3). For diesel cars, the long-run elasticity of mileage to the diesel-oil price is evaluated at: (-623/0.1471) × (0.8461/17, 085) =-0.21.

For example,[START_REF] Collet | Household car use in france: a demographic and economic analysis[END_REF] reports an income-elasticity of household automobility in France of about +0.5, but applying a static model that does not control for the household car ownership level.

Appendix