M Idouakass

M Faganello

H L Berk

X Garbet

S Benkadda

Linear study of the precessional fishbone instability

Keywords: numbers: 52.35.Py, 52.35.-g, 52.25.Dg, 52.35.Qz

The precessional fishbone instability is an m = n = 1 internal kink mode destabilized by a population of trapped energetic particles. The linear phase of this instability is studied here, analytically and numerically, with a simplified model. This model uses the reduced magneto-hydrodynamics (MHD) equations for the bulk plasma and the Vlasov equation for a population of energetic particles with a radially decreasing density. A threshold condition for the instability is found, as well as a linear growth rate and frequency. It is shown that the mode frequency is given by the precession frequency of the deeply trapped energetic particles at the position of strongest radial gradient. The growth rate is shown to scale with the energetic particle density and particles energy while it is decreased by continuum damping.

I. INTRODUCTION

The resonant interaction between a distribution of particles and waves evolving in a plasma is one of the peculiarities of the kinetic dynamics of plasmas. It can excite or damp waves, depending on the sign of the gradient of the distribution of particles in momentum direction. Indeed, when in the neighborhood of a resonance, the distribution of the faster particles is less than that of slower particles, there is a net transfer of energy from the wave to the plasma, leading to the well-known phenomenon of Landau damping 1 . In contrast, it is possible to excite a wave by taking energy from the resonating particles, in the case where there is a reversed gradient in the distribution function, such as when a beam of energetic particles is present in a cold plasma. This case is known as the bump-on-tail instability (referring to a bump of energetic particles on the tail of the distribution function of the thermal particles). This instability has been studied excluding 2 or including 3 in the model a finite dissipation for the wave. It was shown that without dissipation, the nonlinear dynamics would flatten the gradient around the position of the wave-particle resonance, thus leading to a saturation of the mode 2,3 . It was later shown that, with dissipative processes, that saturation would also occur, but there was the possibility of a change in frequency of the mode (called chirping) after its saturation in amplitude. This chirping was seen in particle simulations to be associated with the dynamics of energetic particles and the formation of coherent structures in the phasespace, moving in the momentum direction [3][4][5] . Such dynamics in phase space have been studied, not only for the bump-on-tail case, but also for various Energetic Particle Modes (EPMs), be it for specific modes such as the Alfvén waves 6 or in a more general way 7 .

In the context of Tokamak Physics, where an ultimate goal is to obtain a burning plasma, energetic particles play a key role in heating the thermal plasma. The way they deposit their energy onto the plasma is thus critical.

These energetic particles can be fusion produced alpha particles as well as energetic particles resulting from various heating mechanisms. An important concern for this energy deposition is that there are resonant instabilities that tend to take energy from energetic particles and, at the same time, to expel these particles from the plasma core, before they can transfer their energy to the plasma. One of these instabilities is the fishbone instability.

Fishbones were first observed in the Poloidal Divertor Experiment (PDX) tokamak 8 in 1983. This instability was observed to have a kink-like mode structure, with dominant mode numbers m = n = 1, and to occur in repetitive bursts, where the mode frequency decreases by a factor of about 2 during each burst. It was seen when the Neutral Beam Injection (NBI) was active, in a direction nearly perpendicular to the magnetic field. The mode frequency was comparable to the precession frequency of the deeply trapped particles generated by the NBI. Furthermore, these bursts were associated with losses of fast particles, i.e. an outward radial displacement of these particles.

Two different theoretical interpretations were given for two different regimes. The first regime 9 , of precessional fishbones, refers to the case where the mode frequency ω is much greater than the diamagnetic frequency ω * ,i associated with the radial pressure gradient of the thermal ions. In that regime, the trapped energetic particles destabilize a m = n = 1 mode with the frequency of their precessional motion ω D . This mode is purely kinetic and exists only when energetic particles are present, otherwise it is strongly damped by the Alvfén continuum 10,11 . The second regime 12 refers to a mode with frequency ω ≈ ω * ,i ≈ ω D . In that regime, the energetic particles excite a mode that already exists, a marginally stable internal kink mode that oscillates at the ion diamagnetic frequency.

It should be stressed that, in a torus, an m = n = 1 mode with an electric potential proportional to the minor radius within the q = 1 surface (kink-like mode) corresponds to a rigid displacement of the core region, without a bending of the field lines. Thus the mode, minimizing the work that should be done against the magnetic tension, is the easiest one to excite. Here q is, as usual, the safety factor.

One can make a few simplifications to analyze the essence of the precessional branch of the instability. As it was observed in the work by McGuire 8 , the mode frequency is close to the precessional frequency of deeply trapped energetic ions. It is thus possible to simplify the problem by assuming that all the particles are deeply trapped. With this simplification, it is possible to reduce the particle dynamics to a two-dimensional phase space where the particle-wave interaction occurs only through the electrostatic part of the wave field.

Again, with this simplification, and for a sufficiently large magnetic field, the toroidal angular momentum P ϕ , which is a constant of motion in the unperturbed field, is approximately the poloidal magnetic flux (to within a constant factor) and hence, the minor radius r is an approximate constant of motion. This is an important feature, as then a decrease in fast particle density along radius can be seen as a positive gradient of the distribution function in the toroidal momentum direction [START_REF] White | The theory of toroidally confined plasmas[END_REF] . A decreasing density of deeply trapped energetic particles play exacly the same role of a bump on the tail, in the context of electrostatic waves. Many different works have studied this problem, retaining more effects than we shall use. The more general effects can include: gyrokinetic ions [START_REF] Cole | [END_REF] , vertical plasma elongation 15 , or a different kind of energetic particles, such as barely trapped/circulating electrons 16,17 .

In this paper, we will study the linear phase of the precessional fishbone instability in the simplified model discussed above. This model, due to its simplicity and reduced dimensionality, has the great advantage of being analytically integrable in linear theory. It can thus give direct indication of the roles played by the different mechanisms that drive or suppress the instability. Furthermore it can serve as a benchmark for nonlinear codes, at least when the hypotheses are fulfilled. In section II, we will introduce the model for the bulk plasma as well as for the fast particles. This model is in the spirit of the one developed by Ödblom et al. 10 , with the notable difference that we have here a continuous radial distribution of fast particles but a single value for the magnetic moment, where they use a Heavyside function for the radial distribution and a Maxwellian for the perpendicular velocity distribution. In section III, the linear analysis of the model is presented, and an approximate analytic value is obtained for the instability threshold, the linear growth rate and the frequency. Section IV describes the numerical implementation of the model, and section V discusses the linear numerical benchmarking of the code. The summary and conclusions are found in section VI.

II. MODEL

The model will treat differently the dynamics of the thermal plasma and the dynamics of the energetic particles. In this model, we consider a small inverse aspect ratio tokamak where ε = r/R 1 (R is the major radius and r the minor radius). The thermal plasma will be described by the reduced MHD equations in a cylindrical geometry, meaning that in accord with the small inverse aspect ratio approximation, the toroidal curvature effects can be neglected, exept for the interaction with energetic particles. The fast particles will be described kinetically and the toroidal curvature, which is crucial for the dynamics of the energetic particles, will be included. The curvature will be used in the description of the phase space evolution and, together with the energetic particle pressure, in the vorticity equation.

A. MHD description of the thermal plasma

Let us define the magnetic field and the fluid velocity field as

B = B T + ∇ψ × e z (1)
v = ∇φ × e z = 1 B T ∇φ el × e z (2)
where φ el is the electric potential, φ the stream function, B T is the axial guide field in cylindrical geometry B T = B T e z and ψ = A where A is the axial component of the vector potential. The plasma density variations (ρ(r) = ρ 0 ) and the thermal pressure effects are neglected.

Keeping in the effects of energetic particle pressure, and noting that, as only deeply trapped particles are considered, the perpendicular component of the pressure is much larger than the parallel one P ⊥,h P ,h , the following two equations of the ideal reduced MHD are obtained

∂ψ ∂t = (B • ∇) φ (3)
∂ ∂t Ω + (v • ∇) Ω = 1 µ 0 ρ 0 (B • ∇) ∆ ⊥ ψ -(b × κ) • ∇ P ⊥,h B T ( 4
)
where b is the unit vector along the equilibrium magnetic field, κ mimics the toroidal curvature, Ω is the vorticity Ω = ∆ ⊥ φ and ∆ ⊥ is the cylindrical laplacian in the direction perpendicular to the magnetic guide field along e z .

The fields are expanded in an equilibrium part and a perturbation as ψ = ψ eq + ψ, φ = φ eq + φ. The perturbations are expressed in terms of the following Fourier series ( φ, ψ) = m,n ( φm,n (r), ψm,n (r))e -iωt+inϕ-imθ + c.c. (5) where θ is the poloidal angle and ϕ = 2πz Lz = z R is equivalent to the toroidal angle, implying a periodic axial direction z. For the linear study considered here, only the mode numbers m = n = 1 are relevant. As such, the perturbed fields become

( φ, ψ) = ( φ1,1 (r), ψ1,1 (r))e -iωt+iϕ-iθ + c.c. (6)
In the following, the "hat" notation, the indices, and the explicit radial dependence are dropped for simplicity. Equations ( 3) and ( 4) are now linearized assuming that there is no equilibrium flow (φ eq = 0)

-iωψ = i B T R 0 1 - 1 q φ (7) -iωΩ = i B T µ 0 ρ 0 R 0 1 - 1 q ∆ ⊥ ψ + B T µ 0 ρ 0 R 0 1 - 1 q φ r ∂ r J eq -(b × κ) • ∇ P ⊥,h B T (8)
where J eq is the equilibrium toroidal current, which can be expressed as J eq = ∆ ⊥ ψ eq = 1 r ∂ r r∂ r ψ eq . Combining equations ( 7) and ( 8) yields

1 r 2 ∂ r r 3 ω 2 A (r) -ω 2 ∂ r φ r = -iω B T ρ 0 (b×κ)•∇ P ⊥,h B T (9)
which, once integrated with respect to radius, gives

ω 2 A (r) -ω 2 ∂ r φ r = -iω ρ ( 10
)
where ρ = 1

r 3 r 0 r 2 B T ρ0 (b × κ) • ∇ P ⊥,h B T dr is the fast par- ticle contribution term and ω A (r) = v A,T R 1 -1 q
is the Alfvén frequency associated with the mode helicity.

B. Fast particle response

We now describe the kinetic evolution of the fast particles, which satisfies the Vlasov equation

∂ t f -{H, f } = 0 ( 11
)
where {g, h} = ∂ ϕ g∂ Pϕ h -∂ ϕ h∂ Pϕ g are the Poisson brackets, which operate on the toroidal angle ϕ and its conjugate momentum variable P ϕ . The fast particle distribution f = F eq + δf is decomposed into an equilibrium part F eq and a perturbation δf . The Hamiltonian H = H eq + H 1 is also decomposed into an equilibrium part H eq and a perturbation H 1 .

Then, it should be noted that we are considering a population of deeply trapped fast particles, for which v is negligible compared with v ⊥ and e M ψ eq , where parallel and perpendicular correspond to directions relative to the equilibrium magnetic field, e is the particle charge and M its mass.

These considerations give a few important results. The first one is that all fast particles lie close to the θ = 0 plane (the horizontal mid-plane of the tokamak, outside of the magnetic axis). Therefore, all interactions of the energetic particles with modes described as g(r, θ, ϕ) = Σ ∞ m,n=0 g m,n (r)e -iωt+inϕ-imθ + cc can be taken as g(r, θ = 0, ϕ) = Σ ∞ m,n=0 g m,n (r)e -iωt+inϕ + cc. Furthermore the perturbed Hamiltonian simplifies, from the form H 1 = e(φ -v δA ) to H 1 = eφ, where δA = A -A ,eq . Then the Hamiltonian is of the form

H = H eq (P ϕ , µ) + eφ el,n (P ϕ )e -iωt+inϕ (12)
Here P ϕ is the toroidal momentum, and it is equivalent to the radial variable. As it is defined P ϕ = eRψ, we find

∂P ϕ ∂r = - erB T q (13)
The equilibrium Hamiltonian is simply given by H eq = µB where µ is the magnetic moment.Thus the equilibrium precessional frequency of particles is simply

ω D = ∂H eq ∂P ϕ = µB T q ω C M Rr ( 14
)
where ω C is the cyclotron frequency of the particles considered. It can be seen here that, for a relatively flat q-profile, the precession frequency will decrease with radius.

III. LINEAR THEORY A. Linear particle response

Using this Hamiltonian in Eq. ( 11), one gets the linearized Vlasov equation

∂ t δf = ω D ∂ ϕ δf -∂ ϕ eφ el ∂ Pϕ F eq (15)
which, once solved, gives for the perturbation

δf = eφ el ∂ Pϕ F eq ω D -ω (16)
Now, the shape of the MHD mode must be taken into consideration. The dominant mode is kink-like, which means that the mode numbers are m = n = 1 and, well inside the q = 1 surface, the electric potential is, at the lowest order, proportional to the radius

φ el (r, θ, ϕ) = r φ el r r=r b e -iωt+iϕ-iθ (17)
Here r b is an arbitrary radius that is small enough compared to r * , the radial position of the q = 1 surface, so that φ el r has already reached the asymptotic value it has inside the core. Furthermore if we suppose that the energetic particles remain inside the volume defined by r b , the radial shape of the electric potential given by Eq.17, can be used for determining δf via Eq.16.

The contribution of fast particles to the MHD dynamics is given by

ρ = B T r 3 ρ 0 r 0 r 2 dr (b × κ) • ∇ P ⊥,h B T (18)
where

(b × κ) • ∇ ≈ e ω D µ ∂ ϕ (19)
which, it should be pointed out, is independent of µ. The perturbed pressure is given by

P ⊥,h = f µB T d 3 v = F eq (µ, v , P ϕ ) +δf (µ, v , P ϕ , ϕ) µB T d 3 v ( 20
)
where

d 3 v = 2πv ⊥ dv ⊥ dv = 2π B T M dµdv so that ρ = 2πeB 2 T M r 3 ρ 0 r 0 r 2 dr ω D µ ∂ ϕ µf dµdv = 2πieB 2 T M r 3 ρ 0 r 0 r 2 dr ω D µ µδf dµdv (21) ρ
= 2πieB 2 T M r 3 ρ 0 r 0 r 2 dr ω D µ × µ er φ el r r=r b ∂ r F eq ω D -ω ∂r ∂P ϕ dµdv (22)
Now let us define a general form for the equilibrium distribution function

F eq = Aδ µ µ * -1 δ v v ⊥ F r (r) ( 23
)
where δ is the Dirac delta function and we have chosen the single value of magnetic moment to be µ * . In this way Eq.22 reads ρ = - Starting back from Eq. ( 10), a streched variable can be defined as x = (r -r * )/ around the position where all variations of φ/r are located, i.e. around the q = 1 surface. These localized variations allow us to take the value of r in Eq. ( 25) to be r = r * i.e. y = 1. As a result, one can integrate in radius around the q = 1 surface, so that Linearizing ω A so that ω A = ω A x one finds

2 √ 2πieB 3/2 T A φ el r r=r b µ 3/2 * M 3/2 r 3 ρ 0 r 0 qr 2 dr ω D * ∂ r F r (r) ω D * -ω (
- φ r r=r b = π ω A ρ(r * ) (27)
which translates to the following dispersion relation

i = -K 1 0 y 2 q 2 ∂ y F y (y) q -ωT y dy (28)
where

K = 2 √ 2π 2 µ 3/2 * B 5/2 T eA M 3/2 r * ρ 0 ω A (29)
is a factor taking into account the energetic content of the fast particles as well as the thermal plasma, geometric and magnetic parameters.

C. Threshold, frequency and growth rate

In order to perform an analytic treatment, some additionnal assumptions need to be made about the distribution function. As it was mentioned, all fast particles are assumed to be well contained within the q = 1 surface, and in addition it is assumed that the distribution function's radial gradient goes to zero at the magnetic axis. Furthermore, in order to be able to get an analytic solution, it is assumed that q(y) 1 in the entire integration domain.

General distribution

Taking a F y in a form such that • it is monotonic • its gradient is centered around y = y 0 • its gradient is peaked enough around y = y 0 one can solve approximately the dispersion relation, getting a threshold condition

K 0 = - 1 πy 0 3 F y | y0 ; ω 0 = ω D * (y 0 ) (30)
This solution shows that the threshold value K 0 is inversely proportional to the maximum (negative) value of the distribution function's gradient, meaning that the mode will be more easily unstable if the gradient of the distribution function is more peaked. Furthermore, as expected, the mode frequency is consistent with the precessional frequency of particles at the location of the maximum gradient value.

Specific distribution

Let us now give an analytic profile for the radial distribution

F y (y) = 1 2 (1 -erf(β(y -y 0 )) (31)
with a β and a y 0 chosen so that F y goes to a constant at y = 0 and y = 1 (giving the "narrow gradient condition"

βy 0 1 ; (1 -y 0 )β 1)
. This distribution has a radial gradient

∂ y F y = - β √ π e -β 2 (y-y0) 2 (32)
which allows for an exact solution of the dispersion relation eq. ( 28) to be expressed as

i = K ωT 1 0 yF y (y)dy + K (ωT ) 2 1 0 F y (y)dy - K (ωT ) 2 1 0 F y (y) 1 -ωT y dy (33)
which can be simplified, noticing that [yF y ] 1 0 = 0 and that

1 0 F y (y)dy y 0 , to i = - K ωT y 0 - K (ωT ) 2 F y (0) - K (ωT ) 2 1 0 F y (y)dy 1 -ωT y (34)
The last integral term can be redefined as

1 0 F y (y)dy 1 -ωT y = - β √ π 1 0 e -β 2 (y-y0) 2 dy 1 -ωT y = - β ωT √ π β(1-y0) -βy0 e -z 2 dz X -z (35)
where z = β(y -y 0 ) and X = β ωT (1 -ωT y 0 ). Using again the narrow gradient condition, the integration domain can be stretched to infinity. Here the causal prescription is chosen, meaning that (ω) is taken positive (positive growth rate), giving a negative (X). In this case, the integral reads [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]

+∞

-∞ e -z 2 dz X -z = iπ(2e -X 2 -w(X)) (36)
where w(X) is the Faddeeva function. The exponential term is due to the pole coming from the causal prescription. The dispersion relation is then given as :

i = - Ky 0 ωT - K (ωT ) 2 + iKβ √ π (ωT ) 3 (2e -X 2 -w(X)) (37)
We can show that the frequency at the threshold is given by 1 y0T = ω D * (y 0 ) plus a small correction. Let us define δ = 1-ω 0 T y 0 1, so that X = β ω0T δ, and expand the dispersion relation while only keeping the first order terms in

X i = - K 0 y 2 0 1 -δ - K 0 y 2 0 (1 -δ) 2 + iK 0 βy 3 0 √ π (1 -δ) 3 (1 - 2iX √ π ) + O(X 2 )
(38) To the lowest order, this leads to

i = -2K 0 y 2 0 (1 -β 2 y 2 0 δ) + iK 0 βy 3 0 √ π ( 39
)
from where we get

K 0 = 1 √ πβy 3 0 ; ω 0 = 1 y 0 T (1-δ); δ = 1 β 2 y 2 0 (40)
These relations are the threshold values given for the fishbone mode in our model. Note that the correction to the frequency at threshold δ is indeed small. We can see that these relations are the same as the ones for a general distribution in Eq. 30, with a correction of order O(δ) related to the finite width of the gradient.

It is now possible to find values for the growth rate and the real frequency close to the threshold values. Let us expand the dispersion relation around the threshold values and write K and ω as K = K 0 + δK and ω = ω 0 + δω. The dispersion relation Eq. (37), written as i = D(K, ω) reads

i = D(K 0 , ω 0 )+δK ∂D ∂K K0,ω0 +δω ∂D ∂ω K0,ω0 +O(δK 2 , δω 2 ) (41) It is straightforward to show that ∂D ∂K = D K = i K . At the same time, to zeroth order in δ, ∂D ∂ω = -i √ πK 0 T βy 4 0 -2K 0 T β 2 y 5 0 ( 42
)
which leads to the relation

i δK K 0 -(i √ πK 0 T βy 4 0 + 2K 0 T β 2 y 5 0 )δω = 0 (43)
Taking δω = δω R + iγ and separating the dispersion relation into its imaginary part and its real part we obtain

γ δK K 0 √ π 2βy 0 ω 0 + O(δK 2 ) ; δω R 1 2 √ π βy 0 γ + O(δK 2 )
(44) that represent the growth rate and the real frequency correction close to the threshold.

IV. NUMERICAL IMPLEMENTATION OF THE MODEL

The numerical model is based on the same assumptions as the linear theory. In this case a domain decomposition is done, separating the area of interest in two regions. The first region is the core region, centered on the magnetic axis, and extending up to r b , a radius close to r * , that of the q = 1 surface. It is considered that all fast particles are contained in this region. The second domain is an annular region centered around the q = 1 surface, with r b < r < r b + 2∆ r , in which there are no fast particles. Here, a 2 dimensional code will be used instead of simply integrating numerically the linear equations, as it will serve for future nonlinear studies.

A. Core region

This region is the part of the θ = 0 plane that is located between the magnetic axis and the q = 1 surface. As it was said previously, the radius is parametrized by the toroidal momentum, and therefore this region can be described kinetically in angle-action variables, the action being the toroidal momentum P ϕ , and the angle the toroidal angle ϕ. As before the Hamiltonian is given by an equilibrium part H 0 = µ * B T and a time varying part, i.e. the electric potential of the mode:

H 1 (ϕ, P ϕ ) = er(P ϕ ) φ el (t) r(P ϕ ) r(Pϕ)=r b e iϕ (45)
The evolution of the distribution function f is given by Eq.11.

The Vlasov equation in this region is treated numerically using a semi-lagrangian scheme with the code VLAP, using a CIP-CSL method [START_REF] Nakamura | [END_REF]20 . It can be seen that this is a nonlinear equation, but for the study of interest here, we will look only at the linear phase of the evolution.

For consistency with the analytic prediction, the qprofile will be taken as q(r) = q 0 1 inside the core region. This allows us to get an explicit relation between P ϕ and r:

P ϕ (r) -P ϕ,0 = -e(r 2 -r 2 0 ) B T 2 (46)
where r 0 is the minimal radius of the core region, and P ϕ,0 is the minimal value for the toroidal momentum.

B. Annular region

This region is a narrow layer centered around r * . By considering a single helicity response (m = n = 1), and the fact that its width is small compared with its length (∆ r /r * 1), it is possible to simplify the geometry to a two dimensional slab. As there are no fast particles in this region, the dynamics is described by the usual reduced MHD model, without considering the fast particle contribution

∂ψ ∂t = [ψ, φ] + η∆ ⊥ ψ (47) ∂ ∂t Ω = [Ω, φ] + [ψ, ∆ ⊥ ψ] + ν∆ ⊥ Ω (48) where [f, g] = 1 r * ∂ r f ∂ θ g -1 r * ∂ θ f ∂ r g.
A small but finite resistivity η and viscosity ν have been added for numerical stability.

These equations have been cast in a dimensionless form by normalizing all length scales to the q = 1 radius r * (so that r = 1 corresponds to the q = 1 surface), all time scales to the Alfvén time, and the magnetic fields to the equilibrium magnetic field. Only the linearized equations for the mode m = 1 are retained in our simulations. Once linearized, they read

∂ψ ∂t = iφ∂ r ψ eq + η∆ ⊥ ψ (49) ∂ ∂t Ω = ∂ r ψ eq ∆ ⊥ ψ -iψ∂ rr ψ eq + ν∆ ⊥ Ω (50)
where 1 r * ∂ θ = ik θ = i. The MHD equations in this region are treated numerically using a semi-spectral finite differences scheme with the code AMON 21 . It should be noted here that, in accordance with the analytic theory, the equilibrium poloidal magnetic field profile is taken linear. It is of the form ∂ r ψ eq = B θ,eq = B T r * R r-r * L S , where L S is a characteristic shear length of the magnetic field. where r = 1 is the position of the q = 1 surface, during the linear growth. In solid blue, the real part, in dashed green, the imaginary part.

C. Coupling

The coupling between the two regions takes into account two effects. The first one is the influence of the core region on the layer region. It is calculated using Eq. (10) at the boundary between the two regions: the energetic particle term is computed numerically in the core region. In this way a time dependent boundary condition is obtained, fixing the value of the radial derivative of the electric potential at the boundary between the two domains. This value is used as an inner boundary condition to invert the Poisson equation ∆ ⊥ φ = Ω in the annular region. The outer boundary condition is simply given by φ = 0.

In the case of a slab geometry, the equation used is not actually Eq. ( 10), as it is considered that the radius is constant r = r * , the equation ends up being

ω 2 A (r) -ω 2 ∂ r φ| r b = -ir * ω ρ (51)
Noting here that the normalizations are such that µ 0 = ρ 0 = 1, we can see that the Alfvén frequency is ω A (r) = k θ v A,θ = k θ B θ,eq . Using that definition for ω A for the annular region one can express K as

K slab = 2 √ 2π 2 L S µ 3/2 * B 5/2 T eA M 3/2 r * ρ0k θ v A,θ
. The retroaction of the MHD dynamics on the particles' one is thus simply given by the amplitude and phase of φ/r at the inner boundary. Indeed the MHD evolution provides φ el r(Pϕ) r(Pϕ)=r b , which is used as the input for the perturbed Hamiltonian. These interactions between the two different domains are done at each time step.

V. LINEAR NUMERICAL BENCHMARKS

The numerical results for the linear analysis are presented in Fig. 123In the results presented below L S has been taken to be 2.5. The half-width of the annular layer is ∆ r = 0.15. That way, the Alfvén frequency continuum varies linearly from ω A (-∆ r ) = -0.06 to ω A (∆ r ) = 0.06. The distribution function parameters are set as β = 20 and y 0 = 0.35 (verifying the narrow gradient condition). For Fig. 1 and Fig. 3, µ * = 2 × 10 -2 , which leads to a theoretical linear frequency ω = 5.7 × 10 -3 . The small normalized values taken for the dissipation parameters are ν = η = 10 -9 .

It can be seen from Fig. 1 that the numerical results fit well the theory in the region linear in δK. The growth rate is seen to depart from the linear analytic value when K/K 0 grows far from 1. This is expected as the growth rate is calculated analytically for δK K 0 . It can be seen in Fig. 2 that the numerical mode frequency matches the analytic predictions as well for three cases near threshold, where the resonance frequency is changed via the parameter µ * . If, while keeping µ * constant, y 0 is changed, the mode frequency also follows the expected scaling.

It can be seen on the profile shown in Fig. 3 that, as expected, the strong variations of the mode φ(r) are only present around r = 1 (q = 1). The nearly constant profile on the left of the q = 1 surface describes, in a slab geometry, a rigid displacement of the plasma core, as expected. Around the q = 1 surface, one can observe the expected two MHD resonances located at ω A (r) = ±ω, whose width depends on the growth rate and on the resistivity and viscosity.

VI. CONCLUSION

A simplified model for the precessional fishbone instability, based on a separation between the thermal plasma, treated with a fluid approach, and an energetic particle population, treated kinetically, has been developed. A linear analytic treatment has been done, allowing us to find a threshold condition, a linear frequency and a linear growth rate for the instability within the scope of our model. It is shown that the threshold value depends on geometric and MHD parameters, as well as energetic particle parameters. A linear numerical study of this model has been performed and is in good agreement with the analytic theory, within the bounds of its validity.

It is shown that the internal kink mode can be destabilized by a negative radial gradient in a deeply trapped energetic particle population. This mode is easier to excite as this gradient gets steeper. The frequency of the mode is determined by the precessional frequency of the particles at the position of steepest gradient. It is seen that the kink mode radial structure, when it has a finite frequency, changes from a single step function centered at the q = 1 surface to a two-step structure, where the two steps are centered around the Alfvén resonances ω A = ±ω.

We calculated an instability threshold condition, showing that it is more easily exceeded with increase in the following:

• the radial density gradient of energetic particles

• the perpendicular velocity

• the radial position of the gradient.

The latter is likely due to the increase of the electric potential with the increase in minor radius. On the other hand, the mode will be stabilized by an increase in magnetic shear that causes a greater continuum damping, and by a higher bulk plasma density.

The longer term objectives of this work are to study the nonlinear evolution of the fishbone mode in a simplified model. The results obtained here have allowed us to benchmark the linear evolution of the mode in the numerical simulations for our model.

1 ωM 3 /2 y 3 r * ρ 0 y 0 y 2

 1302 24)Changing the integration variable to y = r r * where r * is the radius of the q = 1 surface, and expressing ω D * as ω D * = qr * rT where T is a characteristic period T = D * (r * ) , one gets ρ = -q 2 ∂ y F y (y ) q -ωT y dy(25)B. MHD response and dispersion relation

ω 2 A

 2 -ω 2 dx (26) Here ∆ r is then defined as ∆ r = r * -r b and φ goes to zero at r = r b + 2∆ r .

FIG. 1 .

 1 FIG.1. Growth rate of the mode as a function of K/K0. Solid line is the theoretical analytic result, stars linked with the dashed line are the numerical points.

FIG. 2 .FIG. 3 .

 23 FIG.2. Frequency of the mode as a function of µ * . Solid line is the theoretical result, stars linked with the dashed line are the numerical points.

Aknowledgements

The authors would like to thank Dr. Boris Breizman for fruitful discussions. This work has also benefited from discussions within the NLED Enabling Research project (the views and opinions expressed therein do not necessarily reflect those of the european commission). This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the Agence Nationale pour la Recherche.


Guillaume Ducoffe

The parallel complexity of coloring games

We wish to motivate the problem of finding decentralized lower-bounds on the complexity of computing a Nash equilibrium in graph games. While the centralized computation of an equilibrium in polynomial time is generally perceived as a positive result, this does not reflect well the reality of some applications where the game serves to implement distributed resource allocation algorithms, or to model the social choices of users with limited memory and computing power. As a case study, we investigate on the parallel complexity of a game-theoretic variation of graph coloring. These "coloring games" were shown to capture key properties of the more general welfare games and Hedonic games. On the positive side, it can be computed a Nash equilibrium in polynomial-time for any such game with a local search algorithm. However, the algorithm is time-consuming and it requires polynomial space. The latter questions the use of coloring games in the modeling of information-propagation in social networks. We prove that the problem of computing a Nash equilibrium in a given coloring game is PTIME-hard, and so, it is unlikely that one can be computed with an efficient distributed algorithm. The latter brings more insights on the complexity of these games.

Introduction

In algorithmic game theory, it is often the case that a problem is considered "tractable" when it can be solved in polynomial time, and "difficult" only when it is NP-hard or it is PLS-hard to find a solution. On the other hand, with the growing size of real networks, it has become a boiling topic in (non gametheoretic) algorithmic to study on the finer-grained complexity of polynomial problems [START_REF] Williams | Fine-Grained Algorithms and Complexity (Invited Talk)[END_REF]. In our opinion, the same should apply to graph games when they serve as a basis for new distributed algorithms. We propose to do so in some cases when it can be easily computed a Nash equilibrium in polynomial time.

The following case study will make use of well-established parallel and space complexity classes to better understand the hardness of a given graph game.

Precisely, we investigate on a "coloring game", first introduced in [START_REF] Panagopoulou | A game theoretic approach for efficient graph coloring[END_REF] in order to unify classical upper-bounds on the chromatic number. Since then it has been rediscovered many times, attracting attention on the way in the study of information propagation in wireless sensor networks [START_REF] Chatzigiannakis | Distributed game-theoretic vertex coloring[END_REF] and in social networks [START_REF] Kleinberg | Information-sharing in social networks[END_REF]. We choose to consider this game since it is a good representative of the separable welfare games -proposed in [START_REF] Marden | Distributed welfare games[END_REF] as a game-theoretic toolkit for distributed resource allocation algorithms -and the additively separable symmetric Hedonic games [START_REF] Burani | Coalition formation games with separable preferences[END_REF]. A coloring game is played on an undirected graph with each vertex being an agent (formal definitions will be given in the technical sections of the paper). Agents must choose a colour in order to construct a proper coloring of the graph. The individual goal of each agent is to maximize the number of agents with the same colour as hers. Furthermore, it can always be computed a Nash equilibrium in polynomial time with a simple local-search algorithm [START_REF] Ducoffe | The complexity of hedonic coalitions under bounded cooperation[END_REF][START_REF] Kleinberg | Information-sharing in social networks[END_REF][START_REF] Panagopoulou | A game theoretic approach for efficient graph coloring[END_REF].

However, for n-vertex m-edge graphs, the above-mentioned algorithm has O(m + n √ n)-time complexity and O(n + m)-space complexity. Therefore, when the graph gets larger, potential applications of coloring games as a computational mechanism design (e.g., in order to assign frequencies in sensor networks in a distributed fashion, or to model the behaviour of social network users with limited power and storage) can be questioned. In particular, the authors in [START_REF] Kearns | An experimental study of the coloring problem on human subject networks[END_REF] report on the limited abilities of human subject networks to solve a coloring problem. In this note, we will investigate on the belonging of our problemthe computation of a Nash equilibrium in coloring games -to some complexity classes that are related to parallel and space complexity. Our goal in doing so is to bring more insights on the complexity of the problem.

Related work.

Apart from lower-bounds in communication complexity [START_REF] Feigenbaum | Sharing the cost of multicast transmissions[END_REF], we are not aware of any analysis of decentralized complexity in game theory. Closest to our work are the studies on the sequential complexity of Hedonic games. Deciding whether a given Hedonic game admits a Nash equilibrium is NP-complete [START_REF] Ballester | NP-completeness in hedonic games[END_REF]. Every additively separable symmetric Hedonic games has a Nash equilibrium but it is PLS-complete to compute one [START_REF] Gairing | Computing stable outcomes in hedonic games[END_REF]. Coloring games are a strict subclass where the local-search algorithm terminates on a Nash equilibrium within a polynomial number of steps. We will go one step further by considering their parallel complexity, something we think we are the first to study.

In [START_REF] Chatzigiannakis | Distributed game-theoretic vertex coloring[END_REF], they introduced a distributed algorithm in order to compute the Nash equilibrium of a given coloring game. Their algorithm is a natural variation of the classical local-search algorithm for the problem, however, it does not speed up the computation of equilibria (at least theoretically). In addition, each agent needs to store locally the colouring of the graph at any given step, that implies quadratic space and communication complexity. Additional related work is [START_REF] Ducoffe | The complexity of hedonic coalitions under bounded cooperation[END_REF][START_REF] Kleinberg | Information-sharing in social networks[END_REF], where it is studied the number of steps of more elaborate local-search algorithms when up to k players are allowed to collude at each step. Informally, collusion means that the players can simultaneously change their colours for the same new colour provided they all benefit from the process (note that the classical local-search algorithm corresponds to the case k = 1).

Contributions. We prove that the problem of computing a Nash equilibrium in a given coloring game is PTIME-hard (Theorem 2). This is hint that the problem is inherently sequential, i.e., it is unlikely the computation of an equilibrium can be sped up significantly on a parallel machine with polynomially many processors. In particular, our negative result applies to the distributed setting since any distributed algorithm on graphs can be simulated on a parallel machine with one processor per edge and per vertex. By a well-known relationship between space and parallel complexity [START_REF] Papadimitriou | Computational complexity[END_REF], Theorem 2 also extends to show that no space efficient algorithm for the problem (say, within logarithmic workspace) can exist. Altogether, this may be hint that coloring games are a too powerful computational mechanism design for "lightweight" distributed applications.

Our reduction is from the standard Monotone Circuit Value problem. However, the gadgets needed are technically challenging, and we will need to leverage nontrivial properties of coloring games in order to prove its correctness. Definitions and useful background will be given in Section 2. We will detail our reduction in Section 3 before concluding this paper in Section 4.

Definitions and notations

We use the graph terminology from [START_REF] Bondy | Graph theory[END_REF]. Graphs in this study are finite, simple, and unweighted.

Coloring games. Let G = (V, E) be a graph. A coloring of G assigns a positive integer, taken in the range {1, . . . , n}, to each of the n vertices in V . For every i, let L i be the subset of vertices coloured i. We name L i a colour class in what follows. Nonempty colour classes partition the vertex set V . The partition is a proper coloring when no two adjacent vertices are assigned the same colour, i.e., for every 1 ≤ i ≤ n and for every u, v ∈ L i , {u, v} / ∈ E. Every graph G defines a coloring game whose n agents are the vertices in V . The strategy of an agent is her colour. Furthermore, every v ∈ L i receives payoff: -1 if there is u ∈ L i s.t. {u, v} ∈ E (in which case, the coloring is not proper), and |L i | -1 otherwise. We refer to Figure 1 for an illustration. Finally, a Nash equilibrium of the coloring game is any coloring of G where no agent can increase her payoff by changing her strategy. In particular, the proper coloring in Figure 1 is a Nash equilibrium. More generally, observe that a Nash equilibrium in this game is always a proper coloring of G. In what follows, we will focus on the computation of Nash equilibria in coloring games.

Theorem 1 ( [6, 12]). For any coloring game that is specified by an

n-vertex m-edge graph G = (V, E), a Nash equilibrium can be computed in O(m + n √ n)- time and O(n + m)-space.
Parallel complexity. Computations are performed on a parallel random-access machine (PRAM, see [START_REF] Greenlaw | Limits to parallel computation: P-completeness theory[END_REF]) with an unlimited amount of processors. However, as stated in the conclusion, our results also apply to more realistic parallel complexity classes. In what follows, we will use the fact that processors are numbered. We will handle with read/write conflicts between processors with the strategy CREW-PRAM (concurrent read, exclusive write). Let PTIME contain the decision problems that can be solved in sequential polynomial-time (that is, with a single processor). Problem A reduces to problem B if given an oracle to solve B, A can be solved in polylogarithmic-time with a polynomial number of processors. In particular, a problem B is PTIME-hard if every problem in PTIME reduces to B (this is formally defined as quasi-PTIME-hardness in [START_REF] Greenlaw | Limits to parallel computation: P-completeness theory[END_REF]). Such reductions are finer-grained than the more standard logspace reductions.

Main Result

Theorem 2. Computing a Nash equilibrium for coloring games is PTIME-hard.

In order to prove Theorem 2, we will reduce from a variation of the wellknown Monotone Circuit Value problem, defined as follows.

Problem 1 ( Monotone Circuit Value).

Input: A boolean circuit C with m gates and n entries, a word w ∈ {0, 1} n such that:

the gates are either AND-gates or OR-gates; every gate has exactly two entries (in-degree two); a topological ordering of the gates is given, with the m th gate being the output gate. Question: Does C output 1 when it takes w as input ?

Monotone Circuit Value is proved to be PTIME-complete in [START_REF] Greenlaw | Limits to parallel computation: P-completeness theory[END_REF].

The reduction

Let C, w be any instance of Monotone Circuit Value. We will reduce it to a coloring game as follows. Let G := (g 1 , g 2 , . . . , g m ) be the gates of the circuit, that are topologically ordered. Construction of the gate-gadgets. For every 1 ≤ j ≤ m, the j th gate will be simulated by a subgraph G j = (V j , E j ) with 12(n + j) -9 vertices. We refer to Figure 2 for an illustration. Let us give some intuition for the following construction of G j . We aim at simulating the computation of the (binary) output of all the gates in C when it takes w as input. To do that, given a supergraph G of G j (to be defined later), and a fixed Nash equilibrium for the coloring game that is defined on G, we aim at guessing the output of the j th gate from the subcoloring of G j . More precisely, the subcoloring will encode a "local certificate" that indicates which values on the two entries of g j cause the output.

Observe that to certify that an OR-gate outputs 1, it suffices to show that it receives 1 on any one of its two entries, whereas for an AND-gate it requires to show that it outputs 1 on its two entries. Since by de Morgan's laws, the negation of an AND-gate can be transformed into an OR-gate and vice-versa, therefore, we need to distinguish between three cases in order to certify the output of the gate. So, the vertices in V j are partitioned in three subsets of equal size 4(n + j) -3, denoted by V 1 j , V 2 j , V 3 j . Furthermore, for every 1 ≤ t ≤ 3, every vertex in V t j is adjacent to every vertex in V j \ V t j .

Fig. 2. Gadget subgraph Gj representing the j th gate. An edge between two subsets of vertices (delimited by an ellipse) denotes the existence of a complete bipartite subgraph.

Let us now describe the structure of the three (isomorphic) subgraphs G j [V t j ] = (V t j , E t j ) with 1 ≤ t ≤ 3. Informally, we will need this internal structure in order to ensure that every of the three subsets V t j will behave as a "truthful" certificate to decide on the output of the gate; i.e., only a few vertices of V j will be used to certify the output of the j th gate, while all others will be divided into artificial aggregates that we name "private groups" whose role is to ensure "truthfulness" of the certificate (this will be made clearer in the following). There are two nonadjacent vertices a t j , b t j ∈ V t j playing a special role. The other vertices in V t j \ {a t j , b t j } are partitioned in two subsets A t j , B t j of respective size 2(n + j) -3 and 2(n + j) -2. The sets A t j , B t j are called the private groups of a t j , b t j . Furthermore, every vertex in A t j is adjacent to every vertex in V t j \ (A t j ∪ {a t j }), similarly every vertex in B t j is adjacent to every vertex in V t j \ (B t j ∪ {b t j }). Since all edges are defined above independently the one from the other, the graph

G j [V 1 j ] = (V 1 j , E 1
j ) (encoded by its adjacency lists) can be constructed with |V 1 j | + |E 1 j | = 4(n + j) 2 -2(n + j) -2 processors simply by assigning the construction of each vertex and each edge to a different processor. Note that each processor can decide on the vertex, resp. the edge, it needs to compute from its number. Overall, it takes O(log(n + j))-time in order to construct G j [V 1 j ] in parallel. The latter can be easily generalized in order to construct G j in O(log(n+ j))-time with |V j |+|E j | processors. Therefore, the graphs

G 1 , G 2 , . . . , G m can be constructed in parallel in O(log(n + m))-time with m j=1 (|V j | + |E j |) processors, that is polynomial in n + m.
Construction of the graph. Let X = {x 1 , x 1 , . . . , x i , x i , . . . , x n , x n } contain 2n nonadjacent vertices, that are two vertices per letter in the binary word w. The graph G = (V, E) for the reduction has vertex-set V = X ∪ m j=1 V j . In particular, it has 2n -9m + 6m(m + 2n + 1) vertices. Furthermore, G[V j ] is isomorphic to G j for every 1 ≤ j ≤ m. In order to complete our reduction, let us now describe how our gadgets are connected the one with the other.

For technical reasons, we will need to make adjacent every vertex in the private group A t j (resp. B t j ), with 1 ≤ j ≤ m and 1 ≤ t ≤ 3, to every vertex in V \ V j . By doing so, note that every vertex in V \ (A t j ∪ {a t j }) is adjacent to every vertex in A t j (resp., every vertex in V \ (B t j ∪ {b t j }) is adjacent to every vertex in B t j ). Furthermore, each edge is defined independently the one from the other. Hence, similarly as above, Finally, we recall that for every j, there are three cases to distinguish in order to decide on the output of the j th gate, with each case being represented with some subset V t j . The union of subsets representing a positive certificate (output 1) is named Y j , while the union of those representing a negative certificate (output 0) is named N j . In particular, if the j th gate is an OR-gate, let Y j := {a 1 j , b 1 j , a 2 j , b 2 j } and N j := {a 3 j , b 3 j } (it suffices to receive 1 on one input). Else, the j th gate is an AND-gate, so, let Y j := {a 1 j , b 1 j } and N j := {a 2 j , b 2 j , a 3 j , b 3 j }. Suppose the j th gate is an OR-gate (the case when it is an AND-gate follows by symmetry, up to interverting Y j with N j , see also Figure 3). Let us consider the first entry of the gate. There are two cases. Suppose that it is the i th entry of the circuit, for some 1 ≤ i ≤ n. If w i = 0 then we make both x i , x i adjacent to both a 1 j , b 1 j ; else, w i = 1, we make both x i , x i adjacent to both a 3 j , b 3 j . Else, the entry is some other gate of the circuit, and so, since gates are topologically ordered, it is the k th gate for some k < j. We make every vertex in N k adjacent to both a 1 j , b 1 j , and we make every vertex in Y k adjacent to both a 3 j , b 3 j . The second entry of the gate is similarly considered, up to replacing above the two vertices a 1 j , b 1 j with a 2 j , b 2 j . We refer to Figure 3 for an illustration. In particular, observe that there is only a constant number of edges that are added at this step for each gate. Furthermore, the construction of these new edges only requires to read the two in-neighbours of the gate in the circuit C. As a result, the last step can be done in parallel in O(log(n + m))-time with m processors.

Structure of a Nash equilibrium

The graph G = (V, E) of our reduction (constructed in Section 3.1) defines a coloring game. Let us fix any Nash equilibrium for this game (that exists by Theorem 1). We will show that it is sufficient to know the colour of every vertex in Y m ∪ N m in order to decide on the output of the circuit C (recall that the m th gate is the output gate). To prove it, we will need the following technical claims in order to gain more insights on the structure of the equilibrium. More precisely, we will prove that there are exactly 6m + 1 colour classes, that are one colour class per private group A t j or B t j and an additional colour for the vertices in X. The intuition is that there are 2(n + m) vertices in one special colour class (including X) that simulates the computation of the output of C, whereas all other vertices are "trapped" with the vertices in their respective private group. We refer to Figure 4 for an illustration.

Claim 1. For every j, any colour class does not contain more than two vertices in every Y j ∪ N j . Furthermore, if it contains exactly two vertices in Y j ∪ N j then these are a t j , b t j for some 1 ≤ t ≤ 3. Proof. A Nash equilibrium is a proper coloring of G. Therefore, since any two vertices in different subsets among V 1 j , V 2 j , V 3 j are adjacent by construction, they cannot have the same colour. Since Y j ∪ N j = {a 1 j , b 1 j , a 2 j , b 2 j , a 3 j , b 3 j } and a t j , b t j ∈ V t j for every 1 ≤ t ≤ 3, the claim follows directly.

Claim 2. Any two vertices that are in a same private group have the same colour. Similarly, x i and x i have the same colour for every 1 ≤ i ≤ n.

Proof. Let S be either a private group (S = A t j or S = B t j for some 1 ≤ j ≤ m and 1 ≤ t ≤ 3), or a pair representing the same letter of word w (i.e., S = {x i , x i } for some 1 ≤ i ≤ n). Let v ∈ S maximize her payoff and let c be her colour. Note that v receives payoff |L c | -1 with L c being the colour class composed of all the vertices with colour c. Furthermore, every u ∈ S receives payoff lower than or equal to |L c | -1 by the choice of v. In such case, every u ∈ S must be coloured c, or else, since the adjacency and the nonadjacency relations are the same for u and v (they are twins), furthermore u, v are nonadjacent, the agent u would increase her payoff to |L c | by choosing c as her new colour, thus contradicting the hypothesis that we are in a Nash equilibrium.

The argument we use in Claim 2 is that twin vertices must have the same colour. In what follows, we will use the same argument under different disguises. Claim 3. Let 1 ≤ j ≤ m and 1 ≤ t ≤ 3. Either A t j or A t j ∪ {a t j } is a colour class, and in the same way either B t j or B t j ∪ {b t j } is a colour class. Furthermore, either B t j ∪ {b t j } is a colour class, or a t j and b t j have the same colour.

Proof. Recall that a Nash equilibrium is a proper coloring of G. Since a t j is the only vertex in V \A t j that is nonadjacent to A t j , furthermore every two vertices in A t j have the same colour by Claim 2, therefore, either A t j or A t j ∪ {a t j } is a colour class. Similarly, either B t j or B t j ∪ {b t j } is a colour class. In particular, suppose that b t j does not have the same colour as her private group. Then, she must receive payoff at least |B t j | = 2(n + j) -2 (else, she would increase her payoff by choosing the same colour as her private group, thus contradicting the hypothesis that we are in a Nash equilibrium). Furthermore, there can be only vertices in V \ (A t j ∪ B t j ) with the same colour c as b t j . Suppose for the sake of contradiction that a t j does not have colour c. There are two cases to be considered. Suppose that A t j ∪ {a t j } is a colour class. Then, a t j receives payoff |A t j | = 2(n + j) -3. In such case, since a t j and b t j are twin vertices in G \ (A t j ∪ B t j ), vertex a t j could increase her payoff to at least 2(n + j) -1 by choosing c as her new colour, thus contradicting the hypothesis that we are in a Nash equilibrium.

Else, a t j and b t j do not have the same colours as their respective private groups. In such case, A t j and B t j are colour classes, hence we can constrain ourselves to the subgraph G \ (A t j ∪ B t j ). In particular, the constriction of the Nash equilibrium to the subgraph must be a Nash equilibrium of the coloring game defined on G \ (A t j ∪ B t j ). Since a t j and b t j are twin vertices in G \ (A t j ∪ B t j ), they must have the same colour by a similar argument as for Claim 2.

As a result, a t j must have colour c in both cases, that proves the claim. We recall that we aim at simulating the computation of the output of all the gates in C. To do that, we will prove the existence of a special colour class containing X and some pair in Y j ∪ N j for every j. Intuitively, the two vertices of Y j ∪ N j are used to certify the output of the j th gate. However, this certificate is "local" in the sense that it assumes the output of the j -1 smaller gates to be already certified. Therefore, we need to prove that there can be no "missing gate", i.e., every gate is represented in the special colour class. Claim 4. Let c be a colour such that L c ⊆ X and L c does not intersect any private group (A t j or B t j for any 1 ≤ j ≤ m and 1 ≤ t ≤ 3). Then, X ⊆ L c and there exists an index j 0 such that the following holds true: |L c ∩ (Y j ∪ N j )| = 2 for every 1 ≤ j ≤ j 0 , and L c ∩ (Y j ∪ N j ) = ∅ for every j 0 + 1 ≤ j ≤ m.

Proof. By the hypothesis L c ⊆ X and L c does not intersect any private group, so, there is at least one vertex of m j=1 (Y j ∪ N j ) with colour c. Let j 0 be the largest index j such that there is a vertex in Y j ∪ N j with colour c. Since by Claim 1, there can be no more than two vertices of Y j ∪ N j that are in L c for every j, therefore, by maximality of j 0 we get |L c | ≤ |X| + 2j 0 = 2(n + j 0 ). In particular, observe that if |L c | = 2(n + j 0 ) then X ⊆ L c and for every 1 ≤ j ≤ j 0 there are exactly two vertices in Y j ∪N j with colour c. So, let us prove that |L c | = 2(n+j 0 ), that will prove the claim. By the choice of j 0 , there is some 1

≤ t ≤ 3 such that a t j0 ∈ L c or b t j0 ∈ L c . In particular, |L c | ≥ min{|A t j0 |, |B t j0 |} + 1 = 2(n + j 0 ) -2 or else, every vertex v t j0 ∈ L c ∩ {a t j0 , b t j0
} would increase her payoff by choosing the colour of the vertices in her private group (that is a colour class by Claim 3), thus contradicting the hypothesis that we are in a Nash equilibrium.

We prove as an intermediate subclaim that for any 1 ≤ j ≤ j 0 -1 such that L c ∩ (Y j ∪ N j ) = ∅, there is some 1 ≤ t ≤ 3 such that a t j , b t j ∈ L c . Indeed, in this situation, there is some t such that a t j ∈ L c or b t j ∈ L c . If b t j ∈ L c then we are done as by Claim 3, a t j ∈ L c . Otherwise, b t j / ∈ L c and we prove this case cannot happen. First observe that a t j ∈ L c in this case. Furthermore, since a t j and b t j do not have the same colour we have by Claim 3 that B t j ∪ {b t j } is a colour class. In this situation, b t j receives payoff 2(n + j) -2 ≤ 2(n + j 0 -1) -2 < |L c |. Since in addition a t j and b t j are twins in G \ (A t j ∪ B t j ), vertex b t j could increase her payoff by choosing colour c, thus contradicting that we are in a Nash equilibrium. This proves a t j , b t j ∈ L c , and so, the subclaim. By the subclaim, there is an even number 2k of vertices in j0-1 j=1 (Y j ∪ N j ) with colour c, for some k ≤ j 0 -1. Similarly, since by Claim 2 the vertices x i , x i have the same colour for every 1 ≤ i ≤ n, |X ∩ L c | = 2n for some n ≤ n. Now there are two cases to be considered.

Suppose that b t j0 ∈ L c . Then, by Claim 3 a t j0 ∈ L c . Furthermore |L c | ≥ 2(n + j 0 ) -1 or else, vertex b t j0 would increase her payoff by choosing the colour of the vertices in B t j0 (that is a colour class by Claim 3), thus contradicting the hypothesis that we are in a Nash equilibrium. As a result, Altogether, |L c | ≥ 2(n + j 0 ), that proves the claim.

|L c | = 2(n + k + 1) ≥ 2(n + j 0 ) -1, that implies n + k ≥ n + j 0 -1,
We point out that by combining Claim 1 with Claim 4, one obtains that for every 1 ≤ j ≤ m, there are either zero or two vertices in Y j ∪ N j in each colour class not containing a private group, and in case there are two vertices then these are a t j , b t j for some 1 ≤ t ≤ 3.

Claim 5. Any two vertices in X have the same colour. Furthermore, for every 1 ≤ j ≤ m, every vertex in Y j ∪ N j either has the same colour as vertices in X or as vertices in her private group.

Proof. Let L c be any colour class with at least one vertex in m j=1 (Y j ∪ N j ). Let j 0 be the largest index j such that there is a vertex in Y j ∪ N j with colour c. In order to prove the claim, there are two cases to be considered. Suppose that We suggest through this case study a more in-depth analysis of the complexity of computational mechanism designs. We would find it interesting to pursue similar investigations for other games. Experiments in the spirit of [START_REF] Kearns | An experimental study of the coloring problem on human subject networks[END_REF] could be helpful for our purposes. Further, we note that PRAM is seen by some as a too unrealistic model for parallel computation. Thus, one may argue that proving our reduction in this model casts a doubt on its reach. However, we can leverage on the stronger statement that Monotone Circuit Value is strictly PTIMEhard [START_REF] Condon | A theory of strict P-completeness[END_REF]. It implies roughly that the sequential time and the parallel time to solve this problem cannot differ by more than a moderate polynomial-factor (unless the solving of all problems in PTIME can be sped up on a parallel machine by at least a polynomial-factor). Our reduction directly shows the same holds true for the problem of computing a Nash equilibrium in a given coloring game, that generalizes our hardness result to more recent parallel complexity classes (e.g., [START_REF] Karloff | A model of computation for MapReduce[END_REF]).

Fig. 1 .

 1 Fig. 1. Proper coloring of a graph G. Each colour class is represented by an ellipse. Every agent receives unit payoff.

m j=1 3 t=1

 3 (|A t j | + |B t j |)|V \ V j | processors are sufficient in order to construct these edges in O(log(n + m))-time, that is polynomial in n + m.

Fig. 3 .

 3 Fig. 3. Edges in G to simulate the two connections of an AND-gate in the circuit.

Fig. 4 .

 4 Fig. 4. A boolean circuit (left) with a Nash equilibrium of the coloring game from our reduction (right). Each colour class is represented with an ellipse. Intuitively, vertices in the central colour class simulate the computation of the output. Other colour classes contain a private group and they are "inactive".

  and so, |L c | ≥ 2(n + j 0 ), as desired. However, since a t j0 and b t j0 do not have the same colour, B t j0 ∪ {b t j0 } is a colour class by Claim 3. In particular, b t j0 receives payoff 2(n + j 0 ) -2 < |L c |. Since a t j0 , b t j0 are twins in G \ (A t j0 ∪ B t j0 ), vertex b t j0 could increase her payoff by choosing colour c, thus contradicting that we are in a Nash equilibrium.

	Else, b t j0 / ∈ L c and we prove this case cannot happen. First observe that j0 ∈ L a t

c . Furthermore, |L c | = 2(n + k) + 1 ≥ 2(n + j 0 ) -2, that implies n + k ≥ n + j 0 -1,

and so, |L c | ≥ 2(n + j 0 ) -1.

This work is partially supported by ANR project Stint under reference ANR-13-BS02-0007 and ANR program "Investments for the Future" under reference ANR-11-LABX-0031-01.

L c = A t j0 ∪ {a t j0 } and L c = B t j0 ∪ {b t j0 } for any 1 ≤ t ≤ 3. We will prove that X ⊆ L c , that will imply that L c is unique in such a case, and so, will prove the claim. By the choice of colour c, L c ⊆ X. Further, observe that there can be no private group with a vertex in L c . As a result, this case follows directly from Claim 4.

Else, either L c = A t j0 ∪ {a t j0 } or L c = B t j0 ∪ {b t j0 } for some 1 ≤ t ≤ 3, and we may assume that it is the case for any colour class L c that contains at least one vertex in m j=1 (Y j ∪ N j ) (or else, we are back to the previous case). So, let us constrain ourselves to the subgraph G[X]. In particular, the constriction of the Nash equilibrium to the subgraph must be a Nash equilibrium of the coloring game defined on G[X]. Since the vertices in X are pairwise nonadjacent, they must form a unique colour class in such case, that proves the claim.

We will need a "truthfulness" property to prove correctness of our reduction. Namely, the value of the output of any gate in the circuit must be correctly guessed from the agents with the same colour as vertices in X. Claim 6. Let 1 ≤ j 0 ≤ m such that for every 1 ≤ j ≤ j 0 , there is at least one vertex in Y j ∪ N j with the same colour c 0 as all vertices in X. Then for every 1 ≤ j ≤ j 0 , L c0 ∩ Y j = ∅ if and only if the output of the j th gate is 1.

Proof. In order to prove the claim by contradiction, let 1 ≤ j 1 ≤ j 0 be the smallest index j such that either Y j ∩ L c0 = ∅ and the output of the j th gate is 1 (false negative) or Y j ∩ L c0 = ∅ and the output of the j th gate is 0 (false positive). We will show that in such case, there is an edge with two endpoints of colour c 0 , hence the coloring is not proper, thus contradicting the hypothesis that we are in a Nash equilibrium. Note that since by de Morgan's laws, the negation of an AND-gate can be transformed into an OR-gate and vice-versa, both cases are symmetrical, and so, we can assume w.l.o.g. that the j th 1 gate is an OR-gate. There are two subcases to be considered.

Suppose that the output of the j th 1 gate is 0 (false positive). In such case, Y j1 ∩ L c0 = ∅. Let us consider the first entry of the gate. If it is the i th entry of the circuit for some 1 ≤ i ≤ n then w i = 0 (because the output of the j th 1 gate is 0) and so, by construction, x i , x i ∈ L c0 are adjacent to a 1 j1 , b 1 j1 . Else, it is the k th gate of the circuit for some k < j 1 . By minimality of j 1 , since the output of the k th gate must be 0 (because the output of the j th

j1 have a neighbour in L c0 in this subcase. We can prove similarly (by considering the second entry of the gate) that a 2 j1 , b 2 j1 have a neighbour in L c0 in this subcase. The latter implies the existence of an edge with both endpoints in

Else, the output of the j th 1 gate is 1 (false negative). In such case, Y j1 ∩L c0 = ∅, hence N j1 ∩ L c0 = ∅. Since the output of the gate is 1, there must be an entry of the gate such that: either it is the i th entry of the circuit for some 1 ≤ i ≤ n, and w i = 1 (in which case, the two vertices x i , x i ∈ Lc 0 are adjacent to both a 3

j1 , b 3 j1 by construction); or it is the k th gate of the circuit for some k < j 1 and this gate outputs 1. In the latter case, by minimality of j 1 , Y k ∩ L c0 = ∅. By construction, every vertex in Y k is adjacent to a 3 j1 , b 3 j1 . As a result, a 3 j1 , b 3 j1 have a neighbour in L c0 in this subcase. The latter implies the existence of an edge with both endpoints in L c0 since N j1 = {a 3 j1 , b 3 j1 }.

Proof of Theorem 2

Proof of Theorem 2. Let C, w be any instance of Monotone Circuit Value.

Let G = (V, E) be the graph obtained with our reduction from Section 3.1, which can be constructed in polylogarithmic-time with a polynomial number of processors. The graph G defines a coloring game. We fix any Nash equilibrium for this game, that exists by Theorem 1. By Claim 5, any two vertices in X have the same colour c 0 . We will prove that there is at least one vertex in Y m with colour c 0 if and only if the circuit C outputs 1 when it takes w as input. Since Monotone Circuit Value is PTIME-complete [START_REF] Greenlaw | Limits to parallel computation: P-completeness theory[END_REF], the latter will prove that computing a Nash equilibrium for coloring games is PTIME-hard. By Claim 6, we only need to prove that for every 1 ≤ j ≤ m, there is at least one vertex in Y j ∪ N j with colour c 0 . To prove it by contradiction, let j 0 be the smallest index j such that no vertex in Y j ∪ N j has colour c 0 . By Claim 5, every vertex in Y j0 ∪ N j0 has the same colour as her private group. In particular, the three of a 1 j0 , a 2 j0 , a 3 j0 receive payoff 2(n + j 0 ) -3. We will prove that one of these three agents could increase her payoff by choosing c 0 as her new colour, thus contradicting that we are in a Nash equilibrium. Indeed, by the minimality of j 0 , it follows by Claim 4 that for any 1 ≤ j ≤ j 0 -1, there are exactly two vertices of Y j ∪ N j with colour c 0 , while for every j 0 ≤ j ≤ m there is no vertex in Y j ∪ N j with colour c 0 . As a result, |L c0 | = 2(n + j 0 ) -2. In particular, any agent among a 1 j0 , a 2 j0 , a 3 j0 could increase her payoff by choosing c 0 as her new colour -provided she is nonadjacent to every vertex in L c0 . We will show it is the case for at least one of the three vertices, that will conclude the proof of the theorem. Assume w.l.o.g. that the j th 0 gate is an OR-gate (indeed, since by de Morgan's laws, the negation of an AND-gate can be transformed into an OR-gate and vice-versa, both cases are symmetrical). There are two cases.

Suppose that the output of the j th 0 gate is 1. In such case, there must be an entry of the gate such that: it is the i th entry of the circuit, for some 1 ≤ i ≤ n, and w i = 1; or it is the k th gate of the circuit for some k < j 0 and the output of that gate is 1. In the latter case, we have by Claim 6 that the two vertices of Y k ∪ N k with colour c 0 are in the set Y k . Assume w.l.o.g. that the abovementioned entry is the first entry of the gate. By construction, the two vertices a 1 j0 , b 1 j0 are nonadjacent to every vertex in L c0 . Else, the output of the j th 0 gate is 0. Therefore, for every entry of the gate: either it is the i th entry of the circuit, for some 1 ≤ i ≤ n, and w i = 0; or it is the k th gate of the circuit for some k < j 0 and the output of that gate is 0. In the latter case, we have by Claim 6 that the two vertices of Y k ∪ N k with colour c 0 are in the set N k . By construction, the two vertices a 3 j0 , b 3 j0 are nonadjacent to every vertex in L c0 . In both cases, it contradicts that we are in a Nash equilibrium.


Hongjiang Chen
email: hongjiangchen@163.com

Analysis of the Rigidity of Arc Flexible Hinge with Different Geometrical Parameters

Keywords: arc flexure hinge, finite element analysis, rigidity

Flexible hinges are widely used in micro robotic. Its rigidity directly influences an organization's terminal localization. Its actual structure geometry size cannot satisfy the theoretical analysis completely in a theoretical supposition condition. In this paper, we analyzed the rotation rigidity of arc flexible hinges in different parameters using finite element software ANSYS. The errors are discovered and compared with theoretical result. An analysis of changes of parameters on the performance of arc flexible hinge was carried out, though the graph of the flexible hinge's parameters and its performance. The key manufacture parameters that affect the performance of arc flexible hinge the most and rules of design are given, which can provide directions of design precision for the flexible hinge.

Introduction

Flexible hinge have some characteristics, such as small volume, without rubs, ceaseless, good rigidity and high sensitivity. With microcomputer electrical system series (MEMS) technical rapidly expanding, flexible hinge are widely applied in the displacement which requests small angular and high-precision rotation, such as gyroscopes, accelerometers, precision instruments and so on. It has broad application prospects in the micron level domain.

The common flexible hinge is in two kinds: beam-shape flexible hinge and arcshaped flexible hinge. The beam-shaped flexible hinge has a big slewing area, but the movement precision is bad. The arc-shaped flexible hinge's movement precision is high, but the slewing range is relatively small [START_REF] Zuo | Calculation and analysis of rotational stiffness for three types of flexure hinges[END_REF]. In order to take into account the movement precision and scope, the following several rotation flexible hinges have been generated: parabolic flexure hinge, an arc flexure hinge and a hyperbola-shaped flexure hinge, etc [START_REF] Lobontiu | Analytical model of displacement amplification and stiffness optimization for a class of flexure-based compliant mechanisms[END_REF].The properties of flexible hinges are rigidity, precision and stress characteristic, etc. The rigidity performance reflects the stress ability and also manifests movement to a vice-flexible degree. In 1965, Paros et al [START_REF] Paros | How to design flexure hinge[END_REF] announced his design development of the circular flexible hinge for the first time, and gave the rigidity formula. Smith et al used the similar method to obtain an elliptic flexible hinge mechanics expression [START_REF] Smith | Arc flexure hinges[END_REF]. Nicolae Lodonitu inferred the parabola and the hyperbolic flexible hinge's rigidity formula [START_REF] Lobontiu | Corner-filleted flexure hinges[END_REF].Wei Xu and Tim King analyzed the rectangular and ellipse flexible hinge's rigidity and rotation precision using the finite element method [START_REF] Xu | Flexure hinges for piezo-actuator displacement amplifiers: flexibility, accuracy and stress considerations[END_REF].

In this paper the arc flexible hinge stiffness to different geometrical parameters is analyzed with software ANSYS10.0. Compared with results of theoretical analysis and finite element analysis (FEA), the errors are analyzed. Through the graph of the flexible hinge parameters and its performance, an analysis of changes of parameters on the performance of the arc flexible hinge was carried out. The key manufacture parameters that affect the performances of an arc flexible hinge the most and rules of design are given, which can give directions of design precision for the flexible hinge.

Rigidity formula of the arc flexible hinge

An arc flexible hinge, as shown in Figure 1, is a particular type of flexure that consists of a necked down section. Parameters t, h, b are flexible hinge's smallest thickness, height and width, respectively, parameter R is the radius of arc, and m  is the centre angle of arc.

Fig. 1. Model of arc flexible hinge

As shown in Figure 1, the infinitesimal is intercepted in the abscissa axis that is

) sin (     x
. To begin, the infinitesimal section is vertical to the abscissa axis. The flexible hinge's angular deformation a z is generated under torque M Z as given in Eq.1 [START_REF] Zuo | Calculation and analysis of rotational stiffness for three types of flexure hinges[END_REF].

        d R t EbR M dx x EI x M m m m m z R R z z z         3 2 sin sin ) cos 2 2 ( cos 12 ) ( ) ( (1)
) tan 1 4 arctan( ) 1 4 ( ) 1 2 ( 12 ] tan ) 1 4 ( 1 [ ) 1 4 ( tan ) 1 3 6 ( 4 ] tan ) 1 4 ( 1 [ ) 1 4 ( tan ) 1 2 ( 8 ) cos 2 2 ( cos 2 2 / 5 4 2 2 2 2 3 3 2 2 2 2 2 4 3 1 m m m m m m m d R t f                                            (2)
Where the parameters r can be expressed by Eq.3.

t R   (3)
The rotation rigidity formula of arc flexible hinge is given by Eq.4.

1 2 1 12 f EbR M K Z Z    (4)

FEA model of arc flexible hinge

ANSYS has some characteristics that the general finite element analysis technology, powerful computing, and reliable result. The arc flexible hinge's basic structure size is b=10mm, t=1mm, R=2.5mm, θm.=90 o The material is Spring steel, E=200GPa, The theoretical calculation and FEA rotational stiffness are obtained through changing arc flexure hinge parameters E、b、t and R as well asθm, as shown in Figure 3 to7. The theoretical value of arc flexible hinge stiffness is simplified theory. FEA is FEA value of arc flexible hinge stiffness. (1)From Figures 3 and4 it can be observed that the rotation stiffness is a linearly increasing relation with material young's modulus E and width b. The theoretical value is bigger than the FEA value. When E and b are smaller, the FEA value and the theoretical value are closer.

(2) From Figure 5 it can observed that the rotation of the stiffness FEA value and the theoretical value is a curve increasing with thickness t, end the speed-up is getting quicker and quicker. When t is bigger, the difference of the theoretical value and the FEA value is bigger. When t<2mm, the FEA value and the theoretical value are close.

(3) From Figure 6 it can be observed that the relation of rotation rigidity and radius R is a decreasing curve, and the rate of reduced scope is gradually decreased. When R>2mm, the FEA value and the theoretical value are closer.

(4) From Figure7 it can be observed that the relation of rotation rigidity and centre angleθm is a decreasing curve, but the decreasing scope is small. The theoretical value is bigger than the FEA value. Whenθm>20 o the FEA value and the theoretical value are closer.

From Figures 3 to 7 it can be observed that the influence of flexure hinge's parameter to its rotation stiffness is: the influence of thickness t is biggest, followed by radius R, centre angleθm, width b and E.

The theoretical value and the FEA value of an arc flexible hinge rotation rigidity is not equal, even if has a big differential value. The reasons are:

(1) The flexible hinge theoretical model that is established using materials mechanics' bending strain theory is built on the basis of certain assumptions.

(2) From Figure 2 it can be observed that the FEA model not only has the displacement in the y axis direction, but also has the displacement change on the z axis direction, when the torque Z M exerted on the z axis for the model. In other words, the flexible hinge not only has the bending strain, but also will have the shearing force to cause upward deformation. Under certain design parameters, the flexible hinge's theory and FEA solution achieve a good match.

Conclusions

The different design parameters to the flexible hinge rotation rigidity influence and the linear relationship are obtained by comparing the rotation rigidity theory solution and the finite element analysis. It can be observed that the influence of flexible hinge's parameter to its rotation stiffness is: the influence of thickness t is biggest, followed by radius R, centre angleθm, width b and E. The reasons that the theoretical value and the FEA value of an arc flexible hinge's rotation rigidity is not equal is that the flexible hinge not only has the bending strain, but also will have the shearing force to cause upward deformation. It is helpful to further analyze the movement of the mechanical deformation using the finite element technology to simulate the flexible hinge performance.

Fig. 2 .

 2 Fig. 2. FEA model of arc flexure hinge The arc flexible hinge FEA model is shown in Figure 2.The model left end surface is restrained completely; the right end surface exerts bending moment m N M . 1 . 0  . The special node 1 of the grid model right end surface represents the output displacement. The unit type chooses 3-D the entity SOLID92 unit model, the entire model uses smart size to free mesh. The FEA model and the output displacement have been obtained with changing arc flexible hinge's parameters E, b, t and R as well as m  .

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .

 34567 Fig. 3. Comparison of FEA value and the theoretical value of flexure hinge with changing E

Please note that the LNCS Editorial assumes that all authors have used the western naming convention, with given names preceding surnames. This determines the structure of the names in the running heads and the author index.

Acknowledgments. This research is financially supported by the Hangzhou Binteng Technology Co.,LTD, Hangzhou Province, China.


Claudiu Tiberiu
email: claudiu.albulescu@upt.ro

Albulescu 1

Dominique Pépin

The loss of interest for the euro in Romania

Keywords: money demand, open economy model, currency substitution, Romania JEL codes: E41, E52, F41

We generalize a money demand micro-founded model to explain Romanians' recent loss of interest for the euro. We show that the reason behind this loss of interest is a severe decline in the relative degree of the euro liquidity against that of the Romanian leu.

Introduction

Romania joined the European Union (EU) in 2007 being now one of the Euro area candidate countries. Long before Romania's entrance to the EU, the leu and the euro went hand in hand as the main transactions and savings currencies. However, since September 2001, the euro holding has considerably diminished as compared to that of the domestic currency. Let t M denote the Romanian domestic money holding, while * t M is the euro holding.

Assuming that t S is the exchange rate, then * t t M S represents the euro holding denominated in lei. Figure 1 shows the drop of the * t t M S / t M ratio over the period 2001:M9-2015:M11. The literature on money demand provides no explanation for such a trend. The money demand in the CEE countries is empirically investigated by [START_REF] Van Aarle | Currency substitution and seignorage in Eastern Europe[END_REF], [START_REF] Mulligan | Shortage and currency substitution in transition economies: Bulgaria, Hungary, Poland, and Romania[END_REF], [START_REF] Dreger | Long-run money demand in the new EU member states with exchange rate effects[END_REF], [START_REF] Fidrmuc | Money demand and disinflation in selected CEECs during the accession to the EU[END_REF], or [START_REF] Dritsaki | A panel data approach to the demand for money in Bulgaria and Romania[END_REF]. Single-country analyses of money demand are conducted inter-alia by [START_REF] Komárek | Money Demand in an Open Transition Economy: The Czech Republic, 1993-2001[END_REF] for the Czech Republic, by [START_REF] Siliverstovs | Dynamic Modelling of the Demand for Money in Latvia[END_REF] for Latvia, or by [START_REF] Hsieh | Tests of currency substitution, capital mobility and nonlinearity of Hungary's money demand function[END_REF] for Hungary. The money demand in Romania was investigated by [START_REF] Andronescu | Long-run estimates of money demand in Romania[END_REF] and [START_REF] Ruxanda | Household money demand in Romania. Evidence from cointegrated VAR[END_REF]. None of these papers provides a micro-founded theoretical model to justify the specification of their empirical money demand functions.

Albulescu and Pépin (2016) represent an exception.

Our contribution to the existing literature is twofold. First, we generalize the micro-founded model of [START_REF] Albulescu | The money demand in an open economy model with microeconomic foundations: An application to the CEE countries[END_REF] by assuming that the relative liquidity degree of the euro against that of the leu is changing. Second, we apply the new model on the Romanian case and explain the loss of interest for the euro during the last period.

A money demand model in an open economy

Generalizing [START_REF] Albulescu | The money demand in an open economy model with microeconomic foundations: An application to the CEE countries[END_REF], we suppose that the lifetime utility function of the domestic agent is:

                             i t , ; P M S , P M , P X U E V i t i t * i t i t i t i t i t i t 0 i i t t , ( 1
)
where t X is the monetary consumption spending denominated in lei, t P is the price index, t  is a stationary stochastic process and t is a deterministic trend.

 .

E t is the expectation conditional upon the information available at time t and the presence of t  and of trend t in the utility function indicates that its properties are subject to changes. This utility specification is based on the assumption that the representative agent holds foreign and domestic money in relation with his total consumption, with no distinction between his consumption of foreign and domestic goods. Now suppose that the utility function takes the form:

                                                              t *
, P X U ,      1 , (2) where
) 1 /( 1    
is the elasticity of substitution between the leu and the euro and   t , t   is a function of t and t  (the share parameter).

If the elasticity  is high, it is easier to replace one currency by another, which represents a proof of monetary integration [START_REF] Fidrmuc | Money demand and disinflation in selected CEECs during the accession to the EU[END_REF]. Therefore, if 1   we have substitutability between currencies, while a value 1   indicates their complementarity. In their simplified model, where  

t , t t     
, [START_REF] Albulescu | The money demand in an open economy model with microeconomic foundations: An application to the CEE countries[END_REF] find that the elasticity of substitution between the leu and the euro is weak, ranging between 0.3 and 0.5 under different estimations, and reject thus the hypothesis of monetary integration with the Euro area.

The expression

                                         1 t * t t t t t t P M S , t 1 P M , t
is the liquidity production function and the term

      t t , t / , t 1     
measures the liquidity degree of the euro against the leu in the eyes of the Romanian representative agent. It indicates the amount of euros needed to produce the same liquidity service as one leu. The liquidity service provided by a currency is likely to be influenced by any institutional change, by financial innovations or technological developments in payment systems, or, in this particular case, by blockages in the real estate market, hence the hypothesis of a changing degree of liquidity.

The real consumption and the CES liquidity production function are next combined according to a Cobb-Douglas utility function with a consumption elasticity of   1 and a liquidity elasticity of  , while the parameters are restricted so that

  1 , t 0 , 1 0 t        and     0 .
The representative agent faces the budget constraint: The parameter  represents the proportional cost the agent faces for holding money (for simplification this cost is considered the same for cash and deposits and is fixed as in [START_REF] Lucas | On the stability of money demand[END_REF] at 1% per year, that is 0.082953% on a monthly basis). It stands for the charges related to the use of money (cost of a bank card, cash theft or loss…).

        * t t t * t t t t t * t * 1 t t t 1 t * 1 t t 1 t B S B M S M X Z i 1 B S i 1 B 1 M S 1 M                    , (
                                                            1 t t t 1 t t t 1 t * t t t t t t P M ) , t ( oc P X P M S ) , t ( 1 P M ) , t ( 1 , (4)
where

1 t 1 t 1 t i 1 i oc       
is the opportunity cost of domestic money holding, while the remaining two first-order conditions,

    0 P M S V E t *
) δ(t, θ 1 θ                                                    , (5)
where

* 1 t * 1 t * 1 t i 1 i oc       
is the opportunity cost of euro holding.

As the left-hand terms of equations ( 4) and ( 5) are the same, it comes:

      * 1 t 1 t t t t * t t oc ln oc ln , t , t 1 ln M M S ln                             . (6)
Equation ( 6) describes the euro holding relative to the leu holding. If the opportunity cost of holding euros increases more than the opportunity cost of holding lei, the Romanians abandon the euro in favor of the leu, and this is what really happened. This substitution effect is even stronger if the elasticity of substitution between currencies is high. Figure 2 shows that the opportunity cost spread

   if   t , t t     
. However, this is not the case because the correlation coefficient is considerably smaller than 1. Consequently, other factors might explain this strong loss of interest for the euro, like a changing liquidity degree of the euro.

Indeed, equation ( 6) also shows that the euro demand depends on

      t t , t / , t 1      ratio. A
decrease in this ratio is associated with a decrease in the liquidity service provided by the euro and thus a reduction in the demand for euro. It remains to assess the relevance of this explanation by estimating the function  

    t 2 2 1 0 t t t t , t , t 1 ln                       . (7)
Inserting ( 7) in ( 6), we obtain the cointegration equation which allows the estimation of the elasticity of substitution  and the ratio

      t t , t / , t 1      :   t * 1 t 1 t 2 2 1 0 t * t t oc ln oc ln t t M M S ln                      . ( 8
)

Empirical evidence

We use monthly statistics from 2001M9 to 2015M11 to estimate equation ( 8). The money market rate is used to compute the discounted interest rate (data come from IMF International Financial Statistics). The structure of bank deposits is used as a proxy of the money demand structure (data are extracted from the National Bank of Romania monthly bulletins). apply the [START_REF] Hansen | Tests for Parameter Instability in Regressions with I(1) Processes[END_REF] test to check the existence of a long-run relationship and we estimate the cointegration equation ( 8) by the Fully-Method Ordinary Least Squares (FMOLS) method.

Before estimating equation ( 8), we verify if our series are I(1). The application of the ADF and PP unit root tests prove the existence of a unit root for both processes (Table 1). [START_REF] Hansen | Tests for Parameter Instability in Regressions with I(1) Processes[END_REF] and Probability Lc test is the p-value of test (the null hypothesis is cointegration).

All the parameters are significant and the  elasticity of substitution is positive, as expected.

However, similar to [START_REF] Albulescu | The money demand in an open economy model with microeconomic foundations: An application to the CEE countries[END_REF], we report a low value for  (0.20), which

shows that the leu and the euro are rather complements than substitutes.

The parameters 1  and 2  are significant, underlining that the euro liquidity service compared to that of the leu       A simple examination of this figure highlights Romanians' loss of interest for the euro. If at the end of 2001 the degree of liquidity of the euro was about 0.25 (indicating that four euros produce the same liquidity service as one leu), it has continuously declined until 2004, with a stabilization around a level of about 0.02. The complementary role of the euro became marginal nowadays.

Conclusions

We build a generalized money demand micro-founded model which allows a change in the relative liquidity degree of the euro against the Romanian leu. The long-run money demand is influenced by the interest rate spread, a drop in this spread generating the loss of interest for the euro.

Figure 1 .

 1 Figure 1. Romanians' euro holding to domestic money holding ratio

Table 1 .

 1 Unit root testsThe results of the FMOLS estimation for equation (8) are highlighted in Table2. The Hansen's test accepts the existence of a cointegrating relationship at 10%, result that validates the theoretical specification of our model.

	Variables				Test ADF	Test ADF	Test PP	Test PP
											(with intercept)	(with trend and	(with intercept)	(with trend and
												intercept)	intercept)
	ln	 S	t	M	* t	M	t			-2.89*	-2.26***	-3.01*	-2.22***
	ln	oc	t	1  	ln	* t oc 	1	-1.74***	-2.00***	-1.74***	-1.99***
	Notes: (i) the null hypothesis is the presence of unit root; (ii) *, **, *** means a p-value for the t-statistic >1%,
	>5% and >10% respectively. Otherwise said, ** implies that null hypothesis of unit root cannot be rejected at
	5% significance level.

Notes: (i) ***, **, * means significance at 1%, 5% et 10% significance level; (ii) Probability means p-value of the test of nullity ; (iii) Lc statistic is the statistic of the Hansen's cointegration test

Acknowledgements

This work was supported by a grant of the Romanian National Authority for Scientific Research and Innovation, CNCS -UEFISCDI, project number PN-II-RU-TE-2014-4-1760.