Guilherme Mazanti

Islam Boussaada

Silviu-Iulian Niculescu

Yacine Chitour

Effects of Roots

Keywords: Time-delay equations, stability analysis, spectral methods, root assignment

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Effects of Roots of Maximal Multiplicity on

the Stability of Some Classes of Delay Differential-Algebraic Systems: The Lossless Propagation Case Guilherme Mazanti * * , * Islam Boussaada * * , * Silviu-Iulian Niculescu * Yacine Chitour *

INTRODUCTION

Time delays are useful modeling tools in a wide range of scientific and technological domains such as biology, chemistry, economics, physics, or engineering. They may represent, for instance, the reaction time of an engineering system, the transfer time of material, energy, or information between parts of a system, the duration of a chemical reaction, or the duration of maturation processes in biology. Due to these applications and the challenging mathematical problems arising in their analysis, systems with time delays have been the subject of much attention by researchers in several fields, in particular since the 1950s and 1960s, such as, for instance, in [START_REF] Bellman | Differentialdifference equations[END_REF]; [START_REF] Halanay | Differential equations: Stability, oscillations, time lags[END_REF]; [START_REF] Pinney | Ordinary difference-differential equations[END_REF]. We refer to [START_REF] Diekmann | Delay equations: Functional-, complex-, and nonlinear analysis[END_REF]; [START_REF] Gu | Stability of time-delay systems[END_REF]; [START_REF] Hale | Introduction to functional-differential equations[END_REF]; [START_REF] Insperger | Semi-discretization for time-delay systems[END_REF]; [START_REF] Michiels | Stability and stabilization of time-delay systems: An eigenvalue-based approach[END_REF]; [START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF] for details on time-delay systems and their applications.

This paper is interested in the analysis of stability and asymptotic behavior of systems with delays of the general form

             x (t) = Ax(t) + N k=1 B k y(t -τ k ), y(t) = Cx(t) + N k=1 D k y(t -τ k ), (1)
Corresponding author: Guilherme Mazanti (guilherme.mazanti@l2s.centralesupelec.fr). G.M., I.B., and S.I.N. are also with Inria, DISCO Team, Saclay -Île-de-France research center.

where x(t) ∈ R dx , y(t) ∈ R dy , N , d x , and d y are positive integers, τ 1 , . . . , τ N are positive delays, and, for k ∈ {1, . . . , N }, A, B k , C, and D k are matrices of appropriate dimensions. Systems such as (1) are delay differentialalgebraic systems since they are written as a system of delay differential equations coupled with a system of algebraic equations with delays. They correspond to a particular class of delay differential-algebraic systems in which delays appear only in the y variable, which is the case in particular in models of lossless propagation phenomena (see Section 2 for an example).

Differential-algebraic systems have been extensively studied in the delay-free setting (see, e.g., [START_REF] Brenan | Numerical solution of initial value problems in differential-algebraic equations[END_REF]; Coleman (1998); [START_REF] Griepentrog | Differential-algebraic equations and their numerical treatment[END_REF]; [START_REF] Kumar | Control of nonlinear differential algebraic equation systems[END_REF]; [START_REF] Kunkel | Differentialalgebraic equations[END_REF]). This kind of system arises naturally in several situations, such as in some electronic circuit models, in some control problems with constraints, or in the limiting behavior of singularly perturbed systems. Delay differential-algebraic systems have also been considered in the literature, arising in general from systems of transport partial differential equations representing some propagation phenomenon and coupled with static and dynamic boundary conditions (see, e.g., [START_REF] Niculescu | Delay effects on stability: A robust control approach[END_REF]; [START_REF] Halanay | Stability radii for some propagation models[END_REF]; [START_REF] Hale | Introduction to functional-differential equations[END_REF]; [START_REF] Brayton | Small-signal stability criterion for electrical networks containing lossless transmission lines[END_REF]). More specific motivating examples are described in Section 2.

The stability analysis of time-delay systems has attracted much research effort and is an active field (see, e.g., [START_REF] Cooke | On zeroes of some transcendental equations[END_REF]; [START_REF] Gu | Stability of time-delay systems[END_REF]; [START_REF] Michiels | Stability and stabilization of time-delay systems: An eigenvalue-based approach[END_REF]; [START_REF] Olgac | An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems[END_REF]; [START_REF] Sipahi | Stability and stabilization of systems with time delay: limitations and opportunities[END_REF]). Similarly to the delay-free situation, one may address the asymptotic behavior of a linear time-invariant time-delay system through spectral methods by considering the corresponding characteristic function, whose complex roots determine the asymptotic behavior of solutions of the system (see, e.g., [START_REF] Hale | Introduction to functional-differential equations[END_REF]; [START_REF] Michiels | Stability and stabilization of time-delay systems: An eigenvalue-based approach[END_REF]; [START_REF] Mori | On an estimate of the decay rate for stable linear delay systems[END_REF]). The characteristic function of (1) is the function ∆

: C → C defined for s ∈ C by ∆(s) = det sE -A - N k=1 e -sτ k B k , (2)
where E, A, B 1 , . . . , B N are (d x + d y ) × (d x + d y ) matrices defined by blocks as

E = Id dx 0 0 0 , A = A 0 C -Id dy , B k = 0 B k 0 D k ,
for k ∈ {1, . . . , N }. Notice that (1) can be rewritten in terms of these matrices as

Ez (t) = Az(t) + N k=1 B k z(t -τ k ),
where z(t) = (x(t) T , y(t) T ) T . Similarly to the delay-free case, the exponential behavior of (1) is determined by the spectral abscissa γ of ∆, defined by γ = sup{Re s | s ∈ C and ∆(s) = 0}, and all solutions of (1) converge exponentially to 0 if and only if γ < 0.

The spectral abscissa of ∆ is related to the notion of dominant root, defined as follows.

Definition 1. Let Q : C → C and s 0 ∈ C be such that Q(s 0 ) = 0. We say that s 0 is a dominant (respectively, strictly dominant) root of Q if, for every s ∈ C \ {s 0 } such that Q(s) = 0, one has Re s ≤ Re s 0 (respectively, Re s < Re s 0 ).

It follows immediately from the above definition that, if Q admits a dominant root s 0 , then γ = Re s 0 , but dominant roots may not exist in general.

Functions of the form (2) are particular instances of quasipolynomials, whose definition is the following. Definition 2. A quasipolynomial is an entire function Q which can be written under the form

Q(s) = k=0 P k (s)e λ k s ,
where is a positive integer, λ 0 , . . . , λ are pairwise distinct real numbers, and, for k ∈ {0, . . . , }, P k is a non-zero polynomial of degree d k . The integer

D = + k=0 d k is called the degree of Q.
The above definition of the degree of a quasipolynomial is motivated by a classical property, provided in (Pólya and Szegő, 1998, Problem 206.2) and known as the Pólya-Szegő bound, which implies that, given a quasipolynomial Q of degree D ≥ 0, the multiplicity of any root of Q does not exceed D. Recent works such as Boussaada and Niculescu (2016a,b) have provided characterizations of multiple roots of quasipolynomials using approaches based on Birkhoff and Vandermonde matrices.

It has been recently remarked (see, e.g., [START_REF] Boussaada | Tracking the algebraic multiplicity of crossing imaginary roots for generic quasipolynomials: a Vandermonde-based approach[END_REF]; [START_REF] Boussaada | Further remarks on the effect of multiple spectral values on the dynamics of time-delay systems. Application to the control of a mechanical system[END_REF][START_REF] Boussaada | Multiplicity-induceddominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF]; Mazanti et al. (2020a,b)) that, for quasipolynomials coming from some systems with time-delays, real roots of maximal mul-tiplicity are often dominant, a property usually referred to as multiplicity-induced-dominancy (MID for short). This property has been shown to hold, in particular, for scalar single-delay differential equations of retarded type (see [START_REF] Mazanti | On qualitative properties of single-delay linear retarded differential equations: Characteristic roots of maximal multiplicity are necessarily dominant[END_REF]), a result we present and explain below in Section 3, and also extended for some systems to the case of complex roots of maximal multiplicity in [START_REF] Mazanti | Spectral dominance of complex roots for single-delay linear equations[END_REF].

One of the applications of the MID property is in the design of stabilizing feedback controllers for control systems with time delays, as in [START_REF] Boussaada | Multiplicity-induceddominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF]. A major difficulty when addressing this question is that, except in degenerate situations, quasipolynomials have infinitely many roots, and one usually disposes only of finitely many parameters that can be chosen in the controller design. If, however, the MID property holds and these parameters are chosen in such a way as to guarantee the existence of a root of maximal multiplicity, then this root is dominant, and hence determines the asymptotic behavior of the system, allowing for stabilization if one chooses this root with negative real part.

The aim of this paper is to start the investigation of whether the MID property holds for delay differentialalgebraic systems of the form (1). For this purpose, we restrict our attention to the first non-trivial situation, corresponding to the case d x = d y = N = 1 in which both equations in (1) are scalar and the system contains a single delay. This simple-looking, low-dimensional case illustrates many of the subtleties in the analysis of the MID property for delay differential-algebraic systems. Our main result for this system, Theorem 4, shows that the MID property does hold in this setting, providing necessary and sufficient conditions on the system parameters for having a real root of maximal multiplicity and characterizing further the other roots of the characteristic quasipolynomial when these conditions are satisfied.

The paper is organized as follows. Section 2 presents some examples of systems which can be put under the form (1). We then present, in Section 3, a previous result from [START_REF] Mazanti | On qualitative properties of single-delay linear retarded differential equations: Characteristic roots of maximal multiplicity are necessarily dominant[END_REF] on the MID property for retarded delay differential equations, which can be seen as a particular case of (1) in which N = 1 and D 1 = 0. We briefly recall the strategy of its proof, which serves as inspiration for the analysis of the MID property for a delay differential-algebraic system with scalar unknowns and a single delay in Section 4. An application to one of the examples from Section 2 is provided in Section 5.

Notation. For a given complex number s, we denote by s, Re s, and Im s its complex conjugate, real part, and imaginary part, respectively. Given nonnegative integers n, k with 0 ≤ k ≤ n, the notation n k represents the usual binomial coefficient

n! k!(n-k)! .

MOTIVATING EXAMPLES

In this section, we briefly describe two systems that can be modeled by delay differential-algebraic equations.

E(t) R 0 v(t, x) i(t, x) C 1 R 1 L, C
Fig. 1. Electrical circuit with a transmission line

Electrical circuit with a transmission line

Consider the electrical circuit from Figure 1, in which a voltage source E(t) with internal resistance R 0 is connected, through a lossless transmission line of normalized length 1 and normalized characteristic impedance L and capacitance C, to a parallel association between a capacitor of capacitance C 1 and a resistor of resistance R 1 . This electrical circuit can be described by the system

           ∂ x v(t, x) + L∂ t i(t, x) = 0, t ≥ 0, x ∈ (0, 1), ∂ x i(t, x) + C∂ t v(t, x) = 0, t ≥ 0, x ∈ (0, 1), v(t, 0) = E(t) -R 0 i(t, 0), t ≥ 0, i(t, 1) = C 1 ∂ t v(t, 1) + 1 R 1 v(t, 1), t ≥ 0,
(3) with suitable initial conditions. Performing the classical change of variables into Riemann invariants

u 1 (t, x) = 1 2 v(t, x) + L C i(t, x) , u 2 (t, x) = 1 2 v(t, x) - L C i(t, x) ,
and setting y 1 (t) = v(t, 1) and y 2 (t) = u 2 (t, 1), one verifies that (3) is rewritten as

                         y 1 (t) = - 1 C 1 1 R 1 + C L y 1 (t) - 2 C 1 C L ρy 2 (t -τ ) + 2 C 1 C L E(t -τ /2) 1 + R 0 L C , y 2 (t) = y 1 (t) + ρy 2 (t -τ ) - E(t -τ /2) 1 + R 0 L C , (4)
where ρ =

1-R0 √ L C 1+R0 √ L C
is the reflection coefficient at the voltage source and τ = 2 √ LC represents the roundtrip travel time along the transmission line. If one is interested in the stability analysis of (4) under no input, i.e., when E(t) = 0 for every t, then (4) reduces to a delay differential-algebraic system of the form (1) with d x = d y = 1 and N = 1. We refer the interested reader to [START_REF] Brayton | Small-signal stability criterion for electrical networks containing lossless transmission lines[END_REF]; [START_REF] Halanay | Stability radii for some propagation models[END_REF]; [START_REF] Niculescu | Delay effects on stability: A robust control approach[END_REF] and references therein for more details on these kind of circuits and their applications, as well as for further examples of engineering systems which can be put under the form of a delay differential-algebraic system (1).

Delayed output feedback for proper control systems

Consider the linear control system ẋ(t) = Ax(t) + Bu(t), y(t) = Cx(t) + Du(t),

(5) where x(t) ∈ R dx is the state, u(t) ∈ R du is the control, y(t) ∈ R dy is the output, d x , d u , and d y are positive integers, and A, B, C, and D are matrices of appropriate dimensions. We consider the problem of stabilizing ( 5) by a delayed output feedback of the form u(t) = Ky(t -τ ), where τ > 0 is the delay and K is a matrix of appropriate dimension. The closed-loop system is then ẋ(t) = Ax(t) + BKy(t -τ ),

y(t) = Cx(t) + DKy(t -τ ), (6)
which is under the form (1) with N = 1 delay.

PRELIMINARY RESULTS ON RETARDED DELAY DIFFERENTIAL EQUATIONS

Consider the delay differential-algebraic system (1) in the particular case N = 1 and D 1 = 0. Letting τ = τ 1 and B = B 1 C, this system can be written under the form x (t) = Ax(t) + Bx(t -τ ). ( 7) System (7) is a system of delay differential equations said to be of retarded type since the derivative of highest order only appears in the non-delayed term x (t).

An important particular case which can be put into the form ( 7) is that of a scalar retarded delay differential equation of order n,

x (n) (t) + n-1 k=0 a k x (k) (t) + n-1 k=0 α k x (k) (t -τ ) = 0, ( 8
)
where n is a positive integer and the coefficients a k and α k are real numbers for k ∈ {0, . . . , n -1}. The corresponding characteristic quasipolynomial is the function ∆ : C → C given by

∆(s) = s n + n-1 k=0 a k s k + e -sτ n-1 k=0 α k s k . ( 9
)
Notice that, according to Definition 2, ∆ is of degree 2n, and hence any root of ∆ has multiplicity at most 2n.

The MID property has been extensively studied for some retarded differential equations under the form (8), such as in [START_REF] Boussaada | Further remarks on the effect of multiple spectral values on the dynamics of time-delay systems. Application to the control of a mechanical system[END_REF][START_REF] Boussaada | Multiplicity-induceddominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF]; [START_REF] Boussaada | Tracking the algebraic multiplicity of crossing imaginary roots for generic quasipolynomials: a Vandermonde-based approach[END_REF]; Mazanti et al. (2020a,b). In particular, [START_REF] Mazanti | On qualitative properties of single-delay linear retarded differential equations: Characteristic roots of maximal multiplicity are necessarily dominant[END_REF] presents the following result. Theorem 3. Consider the quasipolynomial ∆ from (9) and let s 0 ∈ R.

(a) The number s 0 is a root of multiplicity 2n of ∆ if and only if, for every k ∈ {0, . . . , n -1}, 10) is satisfied, then s 0 is a strictly dominant root of ∆.

                     a k = n k (-s 0 ) n-k + (-1) n-k n! n-1 j=k j k 2n -j -1 n -1 s j-k 0 j!τ n-j , α k = (-1) n-1 e s0τ n-1 j=k (-1) j-k (2n -j -1)! k!(j -k)!(n -j -1)! s j-k 0 τ n-j . ( 10
) (b) If (
The proof of this result is detailed in the case n = 2 in [START_REF] Mazanti | On qualitative properties of single-delay linear retarded differential equations: Characteristic roots of maximal multiplicity are necessarily dominant[END_REF], the full proof being provided in the extended version of that reference. A first step of the proof is to remark that it suffices to consider the case s 0 = 0 and τ = 1, since the general case can be reduced to this setting by performing the translation and scaling of the spectrum represented by the change of variables z = τ (s -s 0 ). Part (a) of Theorem 3 can then be obtained by straightforward computations, imposing that ∆ (k) (0) = 0 for every k ∈ {0, . . . , 2n -1}. The dominance proof for establishing (b) is then carried by providing first a priori bounds on the imaginary part of non-roots with nonnegative real part and then using a suitable factorization of ∆ to show, using this bound, that such roots cannot exist.

MID FOR A DELAY DIFFERENTIAL-ALGEBRAIC EQUATION

We consider in this section the delay differential-algebraic equation

x (t) = ax(t) + by(t -τ ), y(t) = cx(t) + dy(t -τ ), (11)
where x(t) ∈ R, y(t) ∈ R, and a, b, c, d are real coefficients. System (11) corresponds to (1) with d x = d y = N = 1 and the explicit computation of its characteristic quasipolynomial ∆ from (2) yields ∆(s) = s -a -e -sτ (sd -ad + bc) .

(12) Note that ∆ is a quasipolynomial of degree 3. The main result we prove in this paper is the following counterpart of Theorem 3. Theorem 4. Consider the quasipolynomial ∆ from (12) and let s 0 ∈ R.

(a) The number s 0 is a root of multiplicity 3 of ∆ if and only if the coefficients a, b, c, d, the root s 0 , and the delay τ satisfy the relations 13) is satisfied, then s 0 is a dominant root of ∆. Moreover, for every other complex root s of ∆, one has Re s = s 0 .

a = s 0 + 2 τ , d = -e s0τ , bc = - 4 τ e s0τ . ( 13
) (b) If (
(c) Let Ξ = {ξ ∈ R | tan ξ = ξ}. If (13) is satisfied, then the set of roots of ∆ is {s 0 + i 2 τ ξ | ξ ∈ Ξ}. Remark 5.
With respect to Theorem 3, Theorem 4 provides, in its part (c), additional information on the location of the other roots of ∆. On the other hand, s 0 is not strictly dominant in this case. Remark 6. In the particular case s 0 = 0 and τ = 1, (13) yield a = 2, d = -1, and bc = -4. The corresponding quasipolynomial ( 12) is then given by ∆(z) = z -2 + e -z (z + 2).

(14) For general s 0 ∈ R and τ > 0, one may reduce to the above setting by performing the translation and scaling of the spectrum represented by the change of variables z = τ (s -s 0 ).

The proof of Theorem 4 follows the same general line of that of Theorem 3, but extra properties should be proved in order to obtain the additional conclusions of Theorem 4. The main properties we need for the proof are provided in Appendix A.

Proof of Theorem 4. Let ∆ be the quasipolynomial obtained from ∆ by setting

∆(z) = τ ∆ z τ + s 0 (15) for z ∈ C. Then ∆(z) = z + b 0 + e -z (β 1 z + β 0 ) with b 0 = τ (s 0 -a), β 1 = -de -s0τ , β 0 = τ e -s0τ (ad -bc -ds 0 ). ( 16
)
It follows immediately from relation ( 15) that s 0 is a root of multiplicity 3 of ∆ if and only if 0 is a root of multiplicity 3 of ∆. Since ∆ is a quasipolynomial of degree 3, 0 is a root of multiplicity 3 of ∆ if and only if ∆(0) = ∆ (0) = ∆ (0) = 0. We compute

∆ (z) = 1 + e -z (-β 1 z -β 0 + β 1 ), ∆ (z) = e -z (β 1 z + β 0 -2β 1 ),
and thus 0 is a root of multiplicity 3 of ∆ if and only if b 0 + β 0 = 0, 1 -β 0 + β 1 = 0, β 0 -2β 1 = 0. One immediately verifies that the above linear system of equations on (b 0 , β 1 , β 0 ) admits a unique solution, given by (b 0 , β 1 , β 0 ) = (-2, 1, 2). Using ( 16), one concludes that s 0 is a root of multiplicity 3 of ∆ if and only if (13) holds, concluding the proof of (a). Notice moreover that, under (13), one has ∆ = ∆, where ∆ is the quasipolynomial defined in (14).

To prove (b), it suffices to show that every root of ∆ lies on the imaginary axis. Note first that

∆(z) = z 3 1 0 t(1 -t)e -zt dt, (17)
as one immediately verifies by integrating by parts. Assume, to obtain a contradiction, that there exists a root z 0 ∈ C of ∆ such that Re z 0 = 0. Writing z 0 = σ + iω for σ, ω ∈ R with σ = 0, one may assume, with no loss of generality thanks to Corollary 10 in Appendix A, that σ > 0 and ω > 0. By Lemma 11 in Appendix A, one has 0 < ω < 2.

Using the fact that z 0 is a non-zero root of ∆, one obtains from ( 17) and taking the imaginary part that 1 0 t(1 -t)e -σt sin(ωt) dt = 0.

Since 0 < ω < 2, the function t → t(1 -t)e -σt sin(ωt) is strictly positive in (0, 1), which contradicts the above equality. Hence (b) is proved.

Finally, (c) follows immediately from the relation between ∆ and ∆ under (13), the fact that all roots of ∆ lie on the imaginary axis, and Lemma 12 in Appendix A.

Remark 7. With respect to other results on multiplicityinduced-dominancy for delay differential equations such as Theorem 3, Theorem 4 provides, in its part (c), the additional information of the location of all roots of ∆.

The set Ξ of the real roots of the equation tan ξ = ξ is infinite, discrete, and can be written as Ξ = {ξ k | k ∈ Z}, where (ξ k ) k∈Z is the increasing sequence of the roots of tan ξ = ξ with the convention that ξ 0 = 0. In particular, for every k ∈ Z, one has ξ k ∈ -π 2 + kπ, π 2 + kπ and ξ -k = -ξ k . We also recall that ξ k = kπ + π 2 + o(1) as k → +∞.

APPLICATION TO THE DELAYED OUTPUT FEEDBACK FOR PROPER CONTROL SYSTEMS

Consider system (6) with

d x = 1, d y = 2, d u = 1, which we write as    ẋ(t) = ax(t) + bk 1 y 1 (t -τ ) + bk 2 y 2 (t -τ ), y 1 (t) = c 1 x(t) + d 1 k 1 y 1 (t -τ ) + d 1 k 2 y 2 (t -τ ), y 2 (t) = c 2 x(t) + d 2 k 1 y 1 (t -τ ) + d 2 k 2 y 2 (t -τ ). (18)
From (2), we compute its characteristic quasipolynomial 18) is not under the form (11), its characteristic quasipolynomial ( 19) is of the same form of that of (12), and thus Theorem 4 can be applied to (19). We wish to design the parameters k 1 and k 2 and obtain conditions on the delay τ in order to achieve maximal multiplicity of some root s 0 < 0. Then Theorem 4 will ensure the dominance of this root, implying the exponential stability of the system. Conditions (13) can be rewritten in the context of ( 19) as

∆(s) = s -a -e -sτ (d 1 k 1 + d 2 k 2 )s -k 1 (ad 1 -bc 1 ) -k 2 (ad 2 -bc 2 ) . (19) Even though (
a = s 0 + 2 τ , d 1 k 1 + d 2 k 2 = -e s0τ , (ad 1 -bc 1 )k 1 + (ad 2 -bc 2 )k 2 = 2 τ -s 0 e s0τ . (20)
Note that the first equation in (20) can be rewritten as s 0 = a -2 τ , hence one may stabilize the system by a root of maximal multiplicity only if a ≤ 0 or τ < 2 a . Let us denote δ i = ad i -bc i for i ∈ {1, 2}. The second and third equations of (20) can be rewritten as

d 1 d 2 δ 1 δ 2 k 1 k 2 = -1 4 τ -a e aτ -2 . ( 21
)
This system admits at least one solution if and only if d 1 δ 2 = d 2 δ 1 or δ i = a -4 τ d i for i ∈ {1, 2}, with exactly one solution in the first case and infinitely many solutions in the second case. This discussion can be concentrated in the following result. Proposition 8. Consider system (18) for given real parameters a, b, c 1 , c 2 , d 1 , and d 2 and a positive delay τ , and assume that either a ≤ 0 or τ < 2 a . Let δ i = ad i -bc i for i ∈ {1, 2} and assume moreover that either d 1 δ 2 = d 2 δ 1 or δ i = a -4 τ d i for i ∈ {1, 2}. Then there exist real parameters k 1 , k 2 such that ( 18) is exponentially stable, with exponential decay rate a -2 τ . Moreover, k 1 and k 2 are solutions of the linear system (21). Consider, as an example, the case a = 1, b = 1, c 1 = 2, c 2 = 1, d 1 = 1, d 2 = 2, and τ = 3 2 . Since a = 1, the corresponding open-loop system is unstable. Note that the inequality τ < 2 a is indeed satisfied and the exponential decay rate one may obtain with Proposition 8 is a -2 τ = -1 3 . One computes, using Proposition 8, the feedback parameters k 1 ≈ -0.87610 and k 2 ≈ 0.13478. Figure 2 presents a numerical simulation of the solutions of (20) with these parameters and with initial conditions x(0) = 1 and y 1 (t) = y 2 (t) = 0 for t < 0. One observes, as expected, that solutions converge to 0 exponentially.

Fig. 2 .

 2 Fig. 2. Trajectories of (18) with the designed delayed output feedback.

ACKNOWLEDGEMENTS

This work is supported by a public grant overseen by the French National Research Agency (ANR) as part of the "Investissement d'Avenir" program, through the iCODE project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. The authors also acknowledge the support of Institut Polytechnique des Sciences Avancées (IPSA).

Appendix A. SOME TECHNICAL RESULTS

We present in this appendix technical results used in the proof of Theorem 4. We start by providing some properties of the quasipolynomial ∆ from ( 14). The first one is the following identity, whose proof is straightforward. Lemma 9. Let ∆ be given by ( 14). Then, for every z ∈ C, one has ∆(-z) = -e z ∆(z).

As a consequence of the previous identity, one immediately obtains the following symmetry property of the roots of ∆.

Corollary 10. Let ∆ be given by ( 14) and assume that

We also need the following a priori bound on the imaginary part of the roots of ∆ outside of the imaginary axis. Lemma 11. Let ∆ be given by ( 14) and assume that

Proof. Let z 0 ∈ C be as in the statement. We write z 0 = σ + iω with σ, ω ∈ R and σ = 0. Thanks to Corollary 10, we assume, with no loss of generality, that σ > 0. Since z 0 is a root of ∆, one has e -z0 (z 0 +2) = 2-z 0 and thus, in particular,

Since f is differentiable and σ > 0, the mean value theorem yields the existence of

and thus

Letting g : R → R be the polynomial g(x) = x 2 -4x + ω 2 , since lim x→±∞ g(x) = +∞, the above inequality implies that g must admit two distinct real roots, and thus its discriminant is positive, i.e., 16 -4ω 2 > 0, which is equivalent to ω 2 < 4. Thus |ω| < 2, as required.

As a final technical result, we provide the following characterization of the roots of ∆ on the imaginary axis.

Lemma 12. Let ∆ be given by ( 14), Ξ be as in the statement of Theorem 4(c), and

Proof. Note that iζ is a root of ∆ if and only if iζ -2 + e -iζ (iζ + 2) = 0, which is the case if and only if

where, for θ ∈ R, R θ is the rotation matrix in R 2 , defined by

.

One then immediately verifies that the above system is equivalent to


D Grasso

D Borgogno

E Tassi

A Perona

Asymmetry effects driving secondary instabilities in two-dimensional collisionless magnetic reconnection

In the framework of the studies on magnetic reconnection, much interest has been recently devoted to asymmetric magnetic configurations, which can naturally be found in solar and astrophysical environments as well as in laboratory plasmas. Several aspects of this problem have been investigated, mainly in a two-dimensional geometry and by means of Particle-in-Cell simulations.

Still, there are open questions concerning the onset and the effects of secondary instabilities in the nonlinear phase of an asymmetric reconnection process. In this work we focus on the conditions that lead to the appearance of the Kelvin-Helmholtz instability following an asymmetric reconnection event in a collisionless plasma. The investigation is carried out by means of two-dimensional numerical simulations based on a reduced fluid model assuming a strong guide field. We show that, unlike the symmetric case, in the presence of asymmetry, a Kelvin-Helmholtz-like instability can develop also for a finite equilibrium electron temperature. In particular, simulations indicate the formation of steep velocity gradients, which drive the instability, when the resonant surface of the equilibrium magnetic field is located sufficiently far from the peak of the equilibrium current density. Moreover, a qualitative analysis of the vorticity dynamics shows that the turbulent behavior induced by the secondary instability is not confined inside the island but can also affect the plasma outside the separatrices. The comparison between simulations carried out with an adiabatic closure and a Landau-fluid closure for the electron fluid, indicates that the latter inhibits the secondary instability by smoothing velocity gradients.

I. INTRODUCTION

Magnetic reconnection (MR) is a change in the topology of the magnetic field, usually accompanied by a fast conversion of magnetic energy into plasma kinetic energy and heat,as well as by a release of accelerated particles following the formation of intense current and vorticity layers (also referred to as sheets, in the literature and in the present paper). Due to its relevance in a wide variety of physical phenomena ranging from space to laboratory plasmas, MR has benefit of a vast amount of studies [START_REF] Biskamp | Magnetic Reconnection in Plasmas[END_REF][START_REF] Priest | Magnetic Reconnection[END_REF] since the formulation of the first model by Sweet and Parker [START_REF] Parker | Sweet's mechanism for merging magnetic fields in conducting fluids[END_REF][START_REF] Sweet | The production of high-energy particles in solar flares[END_REF]. In the past three decades much attention has been paid to the problem of fast reconnection. Starting in the 90's many studies [START_REF] Birn | Geospace environmental modeling (gem) magnetic reconnection challenge[END_REF][START_REF] Biskamp | Collisional and collisionless magnetic reconnection[END_REF][START_REF] Ottaviani | Nonlinera collisionless magnetic reconnection[END_REF][START_REF] Zakharov | Twofluid magnetohydrodynamic description of the internal kink mode in tokamaks[END_REF] focused on the physical mechanisms that can drive fast reconnection in high temperature and low collisional plasmas. Electron inertia became the first candidate to explain the occurrence of MR in this context, where, in the absence of dissipation, the topological constraint on the magnetic field is transferred to generalized fields associated with infinite families of Casimirs invariants, while the topological transition in the magnetic field leads to a continuous cascade of energy towards small scales, thus causing the narrowing of the current and vorticity layers.

It is namely MR caused by electron inertia that we will consider in our analysis.

When the spectrum of the magnetic field is dominated by a single mode, the study of MR can be carried out under the two-dimensional (2D) approximation. In such 2D framework it is possible to distinguish between two different configurations, depending on the presence of finite gradients of the equilibrium current density at the resonant surface, where the MR instability develops: a symmetric one, in which the peak of the current density is located at the resonant surface and an asymmetric configuration where the peak is not located at the resonant surface.

The symmetric case has been extensively studied, analytically and numerically, in the linear and nonlinear regimes [START_REF] Biskamp | Magnetic Reconnection in Plasmas[END_REF]. At the beginning of the 2000s a number of works based on reduced fluid models, which are suitable when a strong component of the magnetic field in one direction is present, were devoted to the nonlinear evolution of the current and vorticity sheets in symmetric configurations in the cold ion limit and for large values of the standard instability parameter ∆ [START_REF] Furth | Finite resistivity instabilities of a sheet pinch[END_REF]. In particular, these works focused on the onset of secondary instabilities of the fluid type, like the Kelvin Helmholtz (KH) one, occurring nonlinearly on top of a MR instability and leading to the generation of turbulence in the current and vorticity sheets. In the low-β e limit, with β e indicating the ratio between the equilibrium kinetic pressure of the electrons and the magnetic pressure of the guide field, past studies carried out in this context [START_REF] Sarto | Secondary instabilities and vortex formation in collisionless-fluid magnetic reconnection[END_REF][10][START_REF] Sarto | Electron parallel compressibility in the nonlinear development of two-dimensional collisionless magnetohydrodynamic reconnection[END_REF][START_REF] Grasso | Secondary instabilities in two-and three-dimensional magnetic reconnection in fusion relevant plasmas[END_REF] agreed on two points. In the first place, it is possible to develop secondary fluid-like instabilities only when the electron temperature is negligible (or, more in general, when ρ s d e , with ρ s indicating the ion sound Larmor radius and d e the electron skin depth). In the second place, the fluid-like instability develops inside the magnetic island generated by the MR process. Still in the context of symmetric reconnection, further studies, based on reduced fluid models, showed that secondary KH-like instabilities develop in the plasma vorticity also in the presence of finite electron temperature, provided that β e is sufficiently high [START_REF] Grasso | Coupling between reconnection and Kelvin-Helmholtz instabilities in collisionless plasmas[END_REF][START_REF] Tassi | Hamiltonian four-field model for magnetic reconnection: nonlinear dynamics and extension to three dimensions with externally applied fields[END_REF].

How the current and vorticity layers are affected by such secondary instabilities in an asymmetric configuration is still an open question, in spite of the relevance of MR in asymmetric systems. Indeed, the asymmetric case is important for instance in strongly driven, laser-generated plasmas [START_REF] Rosenberg | A laboratory study of asymmetric magnetic reconnection in strongly driven plasmas[END_REF] and in the framework of Tokamak physics, due to the fact that gradients of the current density at the rational surface naturally occur in toroidal and cylindrical geometry [START_REF] Bertin | Effect of local current gradients on magnetic reconnection[END_REF][START_REF] Militello | Asymmetric tearing mode in the presence of viscosity[END_REF][START_REF] Militello | Effects of local features of the equilibrium current density profile on linear tearing modes[END_REF]. The asymmetric problem has recently been studied also in endogenous reconnection processes relevant to the electron temperature heating in a fusion burning plasma [START_REF] Coppi | Endogenous and asymmetric magnetic reconnection with assoicated processes of relevance to fusion burning plasmas[END_REF]. A renewed interest in the asymmetric configurations was also generated, after 2015 by the observations in the subsolar magnetosheath region [START_REF] Phan | Electron magnetic reconnection without ion coupling in Earths turbulent magnetosheath[END_REF] and from the Magnetosphere Multiscale Spacecraft (MMS), of reconnection occurring at the Earth's dayside between the solar wind and the Earth magnetic field, where the asymmetry arises both from the different plasma densities on the two sides of the reconnecting layers and from the configuration of the magnetic field itself [START_REF] Burch | Electron-scale measurements of magnetic reconnection in space[END_REF]. Moreover, the KH instability has been observed in conjunction with MR events in several astrophysical and space contexts such as in the magnetosphere, in the boundaries of the coronal mass ejections, in the solar corona and in the solar blowout jets (see Ref. [START_REF] Li | Observing Kelvin-Helmholtz instability in solar blowout jet[END_REF] and references therein). Many numerical simulations have been devoted to this particular 2D asymmetric setting, aimed at highlighting the difference with a fully 3D system and explaining the data reported from the MMS mission [START_REF] Le | Enhanced electron mixing and heating in 3d asymmetric reconnection at the earths magnetopause[END_REF][START_REF] Price | The effects of turbulence on threedimensional magnetic reconnection at the magnetopause[END_REF]. In these studies, mainly conducted using PIC codes, a lot of attention has been paid to the turbulence developing along the current density layers inside and around the magnetic island that forms during the reconnection event. The main result of these investigations is that secondary instabilities leading to turbulence may develop mainly when a fully 3D setting is considered, while the 2D configurations remain laminar.

The evolution of the narrow vorticity sheets following MR in asymmetric configurations is the subject of the present paper. We describe asymmetric MR by adopting the 2D limit of the reduced fluid model of Ref. [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF], which accounts for the parallel temperature evolution and for Landau damping effects as well. The asymmetry is introduced in the equilibrium magnetic field, while no equilibrium density gradients are assumed. The basic idea of our approach is to reduce the problem to a few essential ingredients in order to evaluate the impact of the asymmetry on the stability of the current density and vorticity sheets. We will show that, contrary to what happens in the symmetric case, turbulence driven by KHlike instabilities may develop also when d e ≤ ρ s . Moreover the turbulence, that in the symmetric case is confined inside the magnetic island, in the asymmetric case may cross the separatrix, opening new transport scenarios. We investigate the onset of the KH instability as a function of various parameters characterizing the equilibrium magnetic field. We find that the crucial parameter for the onset of the instability appears to be the ratio x s /λ, between the distance x s of the of the rational surface from the location of the peak of the equilibrium current density, and the equilibrium scale length λ. Finally, the effect of an electron Landau closure on the development of secondary instabilities is also studied.

The paper is organized as follows. In Sec. II the 2D reduced Landau-gyrofluid model is reviewed in its cold ion limit. In Sec. III, after introducing the numerical setup, we describe the results of our analysis of secondary instabilities following asymmetric MR. We conclude in Sec. IV.

II. MODEL EQUATIONS

Our analysis is based on the following reduced fluid model for a plasma consisting of electrons and of a single ionized particle species:

∂∇ 2 ⊥ ϕ ∂t + [ϕ, ∇ 2 ⊥ ϕ] - 2 β e [A , ∇ 2 ⊥ A ] = 0, (1)
∂ ∂t A - 2 β e δ 2 ∇ 2 ⊥ A + [ϕ, A - 2 β e δ 2 ∇ 2 ⊥ A ] + [A , ∇ 2 ⊥ ϕ] + [A , t e ] = 0, (2)
∂t e ∂t + [ϕ, t e ] -[A , q e ] - 4 β e [A , ∇ 2 ⊥ A ] = 0. (3)
In Eqs. ( 1)-( 3), the fields ϕ, A , t e and q e are all functions of the time coordinate t and of the Cartesian spatial coordinates (x, y)

∈ D = {(x, y) : -L x ≤ x ≤ L x , -L y ≤ y ≤ L y }, with
L x and L y positive constants. Translational invariance of the system along the z coordinate of the Cartesian system is assumed.

The model equations ( 1)-(3) are formulated in terms of the following normalized variables:

t = Ω i t, x = x ρ s , y = ŷ ρ s , ϕ = e φ T e , A = Â B 0 ρ s , t e = t e T e , q e = q e n 0 T e ρ s Ω i , (4)
where the caret symbol denotes dimensional quantities. We indicate with ϕ the electrostatic potential whereas A is the magnetic flux function, related to the normalized magnetic field B by B(x, y, t) = ∇A (x, y, t) × ẑ + ẑ, with ẑ indicating the unit vector along the z direction. The model assumes indeed the presence of a constant and uniform component of the magnetic field (referred to as guide field) along the z direction. The amplitude of the guide field is assumed to be much larger than that of the x and y components of the magnetic field, i.e. |∇A | 1, for all x, y ∈ D.

The fields t e and q e indicate the fluctuations of the electron temperature and heat flux, respectively, in the direction parallel to the guide field. In Eq. ( 4) we indicate with B 0 the amplitude of the (dimensional) guide field, with T e the constant and uniform equilibrium electron temperature, with e the proton charge and with n 0 the constant and uniform equilibrium particle density (equal for ions and electrons). The symbols Ω i and ρ s , on the other hand, refer to the ion cyclotron frequency based on the guide field, and to the ion sound Larmor radius, defined by ρ s = T e /m i /Ω i , with m i indicating the ion mass. The system (1)-( 3) is characterized by the two parameters β e and δ. Indicating with m e the electron mass, these two parameters are defined as β e = 8πn 0 T e /B 2 0 and δ = m e /m i . They correspond to the ratio between the equilibrium internal pressure and the magnetic pressure exerted by the guide field, and to the mass ratio, respectively. The symbol ∇ 2 ⊥ = ∂ xx + ∂ yy refers to the Laplacian operator with respect to the coordinates x and y, whereas the bracket

[ , ] is defined by [f, g] = ∂ x f ∂ y g -∂ y f ∂ x g
for two functions f and g.

The system (1)-( 3) is evidently not closed, as long as one does not provide an expression for the parallel heat flux q e in terms of lower order moments. For the present analysis, we will consider two possible closures. The first one corresponds to the adiabatic closure

q e = 0. (5)
In this case the heat transport along the magnetic field lines is neglected and the resulting system possesses a Hamiltonian structure [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF].

The second closure is given by

q e = - 2 δ 2 π Lt e . (6)
with L an operator accounting for the Landau damping. More precisely, we follow the Landau closure previously implemented in Ref. [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF] and derived from Refs. [START_REF] Kobayashi | Three-dimensional simulations and spacecraft observations of sub-ion scale turbulence in the solar wind: influence of Landau damping[END_REF][START_REF] Passot | A fluid approach for ion scale plasmas with weakly distorted magnetic fields[END_REF][START_REF] Sulem | Influence of the Nonlinearity Parameter on the Solar Wind Sub-ion Magnetic Energy Spectrum: FLR-Landau Fluid Simulations[END_REF]. When periodic boundary conditions are imposed along the x and y directions we can indicate with

f (x, y, t) = +∞ k=-∞ f k (t) exp(i(k x x + k y y))
the Fourier series for a generic time-dependent function f satisfying periodic boundary conditions over D. The symbol k is defined by k = (k x , k y ) where k x = πp/L x and k y = πm/L y are the wave numbers along the x and y direction, respectively, with p, m ∈ Z. In the present two-dimensional setting, the action of the operator L on a function f is given by Lf (x, y, t)

= +∞ k=-∞ (Lf ) k (t) exp(i(k x x + k y y)) where (Lf ) k (t) = -[(k 0 • τ • k 0 ) -1/2 ][A , f ] k (t). (7)
In Eq. ( 7), τ = b ⊗ b with indicating the spatial average over the domain D and b = B/|B|. We indicate with k 0 the vector (k x , k y , 0) and with [A , f ] k (t) the Fourier coefficient, associated with the wave number k, of the Fourier series of the function [A , f ](x, y, t). The expression [START_REF] Burch | Electron-scale measurements of magnetic reconnection in space[END_REF] descends from the original formulation of the Landau closure of Ref. [START_REF] Scott | Derivation via free energy conservation constraints of gyrofluid equations with finite-gyroradius electromagnetic nonlinearities[END_REF] and applies, as mentioned above, to the case where periodic boundary conditions are imposed along both the x and y directions. In the present analysis, however, periodic boundary conditions will be imposed only along the y direction. Also, the ordering on which the model ( 1)-( 3) is based, admits only quadratic nonlinearities. Therefore, as already applied in Ref. [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF], a variant of the expression [START_REF] Burch | Electron-scale measurements of magnetic reconnection in space[END_REF], specified in Sec. III A, will be adopted.

The model ( 1)-( 3) applies to collisionless plasmas with β e ∼ δ 2 1. Ions are assumed to be cold and the characteristic frequency of the phenomena under consideration has to be much smaller than the ion cyclotron frequency Ω i . For the derivation of the model we refer to Ref. [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF], where the model is obtained from the more general gyrofluid model of Ref. [START_REF] Scott | Derivation via free energy conservation constraints of gyrofluid equations with finite-gyroradius electromagnetic nonlinearities[END_REF] by making use of an asymptotic ordering. Eq. ( 1) can be seen as the evolution equation for the vorticity ∇ 2 ⊥ ϕ associated with the E×B velocity v E×B = ẑ ×∇ϕ. From the quasi-neutrality equation, it follows that for this model the relation n e = ∇ 2 ⊥ ϕ holds [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF],

with n e indicating the fluctuations of the electron density with respect to the equilibrium value n 0 . Therefore, Eq. ( 1) can also be seen as a continuity equation for the electron species. Eq. ( 2) corresponds to the parallel electron momentum equation or, equivalently, to the parallel Ohm's law, where the terms associated with the electron inertia, i.e. those multiplied by the coefficient δ 2 , break the frozen-in condition and allow for MR. Finally, Eq. ( 3) governs the evolution of the parallel electron temperature. Perpendicular electron temperature fluctuations as well a fluctuations of the ion gyrofluid moments are neglected.

III. SIMULATION RESULTS AND ANALYSIS A. Numerical setup

We solve numerically the system (1)-( 3) (adopting the closure ( 5) or ( 6)) over the domain D with L x = 46.26 and L y = 12π. The initial condition on the fields ϕ, A and t

e corresponds to ϕ(x, y, 0) = 0, A (x, y, 0) = A eq (x) + A 0 x -x s λ exp(ik y y) + c.c., t e (x, y, 0) = 0. (8)
In Eq. ( 8), A eq (x) = -A ln cosh (x/λ) + αx, where A, λ > 0 and α are three constant parameters. The function A eq is the equilibrium magnetic flux function, which yields the following asymmetric sheared component of the equilibrium magnetic field:

B eq (x) = ∇A eq (x) × ẑ = A λ tanh x λ -α ŷ, (9)
where ŷ is the unit vector along the y direction. The parameter A is the amplitude of the equilibrium magnetic field, while λ and α control its scale of variation and asymmetry, respectively. In the numerical simulations, the mode number m, associated with the wave number k y of the initial perturbation, has been fixed by the condition m = 1. We indicate with x = x s the position of the resonant surface which, in this 2D setting, is identified by the condition k y B eq (x s ) = 0, i.e. by

x s = λ arctanh λα A . (10)
In particular, for α = 0 one recovers a symmetric equilibrium, with the resonant surface at x = 0, equidistant from the boundaries x = L x and x = -L x . For α = 0 one obtains an asymmetric equilibrium, in which the equilibrium parallel current density (whose amplitude is given by j eq (x) = -(d 2 A eq /dx 2 )(x) = -A/(λ 2 cosh 2 (x/λ))) does not have an extremum value at the resonant surface. We recall that a 2D problem for fields depending on x, y and t with an asymmetric equilibrium such as that of Eq. ( 9), can also be mapped to a single-helicity problem with fields depending on x, y + αz and t, where α = k z /k y , and k z indicates the wave number of the perturbation along the z direction [START_REF] Borgogno | Aspects of three-dimensional magnetic reconnection[END_REF] . We recall, however, that the three-dimensional (3D) extension of the model ( 1)-( 3) requires |α| 1.

In Eq. ( 8) we also indicated with A 0 a function, peaked around (x -x s )/λ = 0, which determines the profile of the amplitude of the perturbation.

At a given time t, the values of the fields correspond to ϕ(x, y, t) = φ(x, y, t), A (x, y, t) = A eq (x) + à (x, y, t) and t e (x, y, t) = t e (x, y, t), where φ, à and t e indicate the perturbations of the fields, which satisfy the relation φ(x, y, 0) = t e (x, y, 0) = 0 and à (x, y, 0) = A 0 ((x -x s )/λ) exp(ik y y) + c.c.. We impose that such perturbations vanish at x = -L x and x = L x and that they are periodic along the y direction.

As anticipated in Sec. II, a variant of the original Landau closure [START_REF] Borgogno | Test-electron analysis of the magnetic reconnection topology[END_REF], adapted to this model, is used for non-adiabatic simulations. This form respects the requirement on quadratic nonlinearities and permits to apply the operator L to functions which are periodic only along the y direction. The explicit expression of the operator can be provided in terms of the coefficients of the Fourier series of the function Lf (x, y, t). The expression for the generic Fourier coefficient associated with the wave number k y is given by

(Lf ) ky (x, t) = = 1 (A 2 /λ 2 -(A 2 /(λL x )) tanh(L x /λ) + α 2 )k 2 y ((A/λ) tanh(x/λ) -α)ik y ) f ky (x, t). ( 11
)
The expression [START_REF] Sarto | Electron parallel compressibility in the nonlinear development of two-dimensional collisionless magnetohydrodynamic reconnection[END_REF] is analogous to the one adopted in Ref. [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF] and is obtained from the original expression (7) by removing nonlinearities in the perturbations, by adopting the expression (9) for the magnetic equilibrium and by using differential operators in the real space along the x coordinate.

Simulations have been carried out with the numerical code described in Ref. [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF]. This code is 3D, explicit in time, assuming a third order Adams-Bashfort scheme, and parallelized along both the periodic directions y and z, while a compact finite difference scheme is implemented on a non-equispaced grid along the magnetic field shear direction [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF]. For the purpose of this paper the code has been run in the 2D limit, avoiding dependence on the zdirection, while a resolution of n y = 960 grid points has been adopted along the y-direction.

For the x-direction n x = 1200 grid points have been adopted on the non equispaced grid so as to guarantee a resolution of dx = 0.02 around x = 0, where the secondary instabilitiy takes place.

B. Secondary instability induced by asymmetry in the adiabatic case

We consider first the model ( 1)-( 3) with the adiabatic closure [START_REF] Borgogno | Aspects of three-dimensional magnetic reconnection[END_REF].

We fix the values of the parameters δ and β e by setting δ = 0.1 and β e = 0.05. The values of these two parameters were not changed throughout the whole paper because we focused on the dependence on parameters characterizing the equilibrium magnetic field and, in particular, its asymmetry.

We begin our analysis illustrating the results from a simulation of an asymmetric case with the following parameters: A = 1, α = 1/4, λ = 3. In this Subsection, these results will be presented in comparison with a symmetric case possessing the same values of the parameters except for α = 0. In Fig. 1 a comparison of the growth of the perturbed magnetic flux at the X-point, A X , in the two cases, is shown, up to the nonlinear stage when the saturation of the magnetic island is reached. Note that the symmetric and asymmetric case differ in the initial x coordinate of the X-points, which are at x = 0 and x = 2.92 respectively.

The simulation times have been rescaled in order to start with the same amplitude. The perturbed magnetic flux A X provides also a measure of the linear growth rate. We see then that the asymmetry causes a decrease of the growth rate. Incidentally we mention (although we do not show it here) that for the m = 1 mode considered here, which is the most unstable one for our equilibrium configuration, we observe in general a decrease of the growth rate as the value of the asymmetry parameter increases. We point out that linearly this is not true for the less unstable modes, as shown in Ref. [START_REF] Borgogno | Test-electron analysis of the magnetic reconnection topology[END_REF]. In the same plot we show also the drift, present in the asymmetric case, of the X-point location from the initial position at x = 2.92 towards x = 0.

The drift along the x-direction can be clarified by equating the time evolution, in the Lagrangian and Eulerian frames of reference, of the generalized magnetic flux function

F = A -2 βe δ 2 ∇ 2 ⊥ A , evaluated at the magnetic island X-point of coordinates (x X (t), y X (t)). This yields dF (x X (t), y X (t), t) dt = ∂F (x, y, t) ∂t | X + v X • ∇F | X ( 12
)
where the subscript X indicates the value of the field at the istantaneous position of the X-point and v X is the drift velocity of the X-point.

At each time, the eq. ( 2), evaluated at the X-point, reads as

∂F (x, y, t) ∂t | X = -[ϕ, F ]| X (13)
where

[A , ∇ 2 ⊥ ϕ]| X = [A , t e ]
| X = 0 because of the cancellation of the reconnecting component of the magnetic field at the X-point, i.e. ∇A (x, y, t) × ẑ| X = 0. Since the magnetic flux A is an even function of y it follows that the first derivative of F along y vanishes at the magnetic field null point. Hence the X-point position varies only with respect to the x-coordinate, and its velocity along the x direction, denoted as v X x , becomes

v X x = dF (x X (t), y X (t), t) dt + [ϕ, F ]| X ∂F ∂x | X -1 (14)
Note that, in the symmetric case, the denominator ∂F/∂x| X in Eq. ( 14) vanishes, which means that also the numerator must be zero for the relation not to become singular. This condition is fulfilled only for a constant value of the X-point position.

The blue circle in Fig. 1 shows the time evolution of the X-point position according to the theoretical prediction obtained at each time step by the relation x X (t + ∆t) = x X (t) + v X x (t)∆t. The agreement with the the curve obtained from the simulation is quite relevant as long as the fields exhibit a rather smooth nonlinear behavior, while the post-processing evaluation of the field gradients (which appear at the denominator of Eq. ( 14) across the X-line are not reliable any more in the advanced nonlinear phase.

The reader might be surprised by the abrupt jump in the position of the X-point visible at around t = 330 in Fig. 1. We claim that this displacement is a consequence of the nonlinear evolution, described later on in this section, in which the turbulence propagating along the separatrix has reached the X-point.

In Fig. 2 we show a comparison between the contour plots of the vorticity for the asymmetric and symmetric cases at a comparable time in terms of the island width. The current densities, not shown here, have analogous distribution. It is worth mentioning that in the low β e regime considered here the current density and vorticity layers are coupled together and undergo a similar nonlinear evolution [START_REF] Sarto | Secondary instabilities and vortex formation in collisionless-fluid magnetic reconnection[END_REF][START_REF] Grasso | Coupling between reconnection and Kelvin-Helmholtz instabilities in collisionless plasmas[END_REF]. Contour plots show that, in both cases, vorticity tends to concentrate along the separatrices of the magnetic island, although the asymmetry leads to an island with a D-shaped structure. In the symmetric case the phase mixing process of the Lagrangian invariants associated with the infinite families of Casimirs, characterizing the adiabatic and isothermal limits of this model [START_REF] Grasso | Phase mixing and satu-ration in Hamiltonian reconnection[END_REF], generates laminar vorticity layers also inside the magnetic island, as the nonlinear evolution proceeds. This is no longer true in the asymmetric case, where the layers remain concentrated mainly along the separatrix. Moreover, as the color bar reveals, the asymmetric vorticity layers are more intense. This latter feature can be better appreciated when analyzing the velocity fields shown in Fig. 3. In the top right panel of Fig. 3 the profiles, for the asymmetric case, of the x and y components of the velocity field v = v E×B = ẑ × ∇ϕ are plotted as a function of x at y = 10 in comparison with the symmetric case (left top panel). The dashed lines represent the location of the two branches of the separatrix for the chosen y value. We first highlight the difference between the peaks in the v y profiles. In the symmetric case they are larger and contained within the island, while in the asymmetric case the profile of v y is more peaked on the lower branch of the separatrix (the one located at a position with x slightly less than x = 0), where the thinnest velocity layers appear. Second, focusing on this layer close to x = 0, we observe that the v y profile is steeper in the inner side of the magnetic island with respect to the outer side. The shape of this peak resembles a kind of asymmetric Bickley jet which, according to the shear flow instability theory [START_REF] Biskamp | Magnetic Reconnection in Plasmas[END_REF], may be prone to the KH instability. We recall that the Bickley jet, defined as v y (x) = 1/ cosh(x) 2 , is the standard shear flow adopted in fluid mechanics for studying the onset of the KH instability.

According to the linear theory, in absence of magnetic field this jet is unstable for any wave vector below a marginal value. However, in the case of a magnetized plasma a stabilizing role can be exerted by a magnetic field parallel to the flow. In this case in fact the plasma motion associated with the instability is counteracted by the tension that is developed by the distortion of magnetic field lines trapped into the plasma. As a consequence for an ideal plasma when the jet is embedded in a sheared magnetic field like B y (x) = B 0 tanh(x) the linear growth rate of the KH instability progressively decreases with increasing B 0 , up to a critical value, as well as the number of linearly unstable modes, both at small and high k y .

In the nonlinear regime, even for field intensities significantly below the critical threshold derived in the linear limit, the stabilizing effect is dramatic.

In the bottom panel of Fig. 3 the profiles of the x and y components of the magnetic field are also plotted. We observe that, as expected, the B x component is almost negligible inside the magnetic island in both the symmetric and asymmetric case. On the contrary the absolute value of the B y component has a different behavior. In the symmetric case it has comparable values on the two sides of the lower branch of the separatrix, while in the asymmetric case it shows a steep gradient through the lower branch itself dropping to almost a constant value inside the island. This balance between the velocity and magnetic fields justifies the onset of the KH instability [START_REF] Biskamp | Magnetic Reconnection in Plasmas[END_REF], that we observe later on in the nonlinear evolution, on the side of the layer facing the inner part of the magnetic island where the stabilizing effect of B y is weaker.

The development of the KH instability is illustrated by three snapshots taken at subsequent times in Fig. 4 and by the movie provided in the supplemental material. From the dynamical point of view the velocity layers can be seen as counter streaming plasma jets starting from the X-point and propagating at the beginning towards y = 0, where they collide forming vortices. Subsequently, jets propagate backward towards the X-point undergoing KH instability, which increases the amplitude of the vortices. When the first vortex has reached the X-point we observe an acceleration of the process of vortex formation and propagation.

Moreover the vortices start propagating also on the other branch of the separatrix, where they are stretched by the more intense magnetic field. The current density layers, not shown here, exhibit a similar behavior.

This indicates a crucial difference with respect to the symmetric case. We point out that the values of the parameters correspond to a regime where ρ s /d e = β e /(2δ 2 ) ≈ 1.6. According to previous analyses [START_REF] Sarto | Secondary instabilities and vortex formation in collisionless-fluid magnetic reconnection[END_REF][10][START_REF] Sarto | Electron parallel compressibility in the nonlinear development of two-dimensional collisionless magnetohydrodynamic reconnection[END_REF][START_REF] Grasso | Secondary instabilities in two-and three-dimensional magnetic reconnection in fusion relevant plasmas[END_REF] in such regime with ρ s ≥ d e , even though the current density and vorticity fields exhibit small scale structures, in the symmetric case their behavior remains laminar, as in Fig. 2(a), throughout the whole nonlinear evolution and no secondary instabilities develop. Also, we remark that here we are adopting an adiabatic closure, unlike

Refs. [START_REF] Sarto | Secondary instabilities and vortex formation in collisionless-fluid magnetic reconnection[END_REF][10][START_REF] Sarto | Electron parallel compressibility in the nonlinear development of two-dimensional collisionless magnetohydrodynamic reconnection[END_REF][START_REF] Grasso | Secondary instabilities in two-and three-dimensional magnetic reconnection in fusion relevant plasmas[END_REF], where an isothermal closure was assumed for the electron fluid. When the same initial conditions, as is the case here, are considered for ∇ 2 ⊥ ϕ and t e , the solution t e = 2∇ 2 ⊥ ϕ can be taken for the system (1)-( 3) with q e = 0. Replacing this solution in the last term of Eq. (2) casts the system in the form

∂∇ 2 ⊥ ϕ ∂t + [ϕ, ∇ 2 ⊥ ϕ] - 2 β e [A , ∇ 2 ⊥ A ] = 0, (15)
∂ ∂t A - 2 β e δ 2 ∇ 2 ⊥ A + [ϕ, A - 2 β e δ 2 ∇ 2 ⊥ A ] + 3[A , ∇ 2 ⊥ ϕ] = 0. (16)
This form differs from that of the isothermal system (which descends from replacing Eq.

(3) with t e = 0) by the factor 3 that multiplies the last term on the left-hand side of Eq. ( 16). This term is denoted as the electron compressibility term and is the one responsible for suppressing the secondary instability in the cases studied in Refs. [START_REF] Sarto | Secondary instabilities and vortex formation in collisionless-fluid magnetic reconnection[END_REF][10][START_REF] Sarto | Electron parallel compressibility in the nonlinear development of two-dimensional collisionless magnetohydrodynamic reconnection[END_REF][START_REF] Grasso | Secondary instabilities in two-and three-dimensional magnetic reconnection in fusion relevant plasmas[END_REF]. With the adiabatic assumption, the relative strength of this term is thus even greater than in the isothermal case. In spite of this, the asymmetry appears to be able to overcome the damping effect associated with this term.

We remark that a secondary KH instability was indeed observed also in the symmetric case but in a different regime of the physical parameters, more precisely in the cold electron limit, where ρ s /d e << 1. In the latter case, however, the head of the colliding jets have a typical symmetric (mushroom) shape (see for instance Ref. [START_REF] Sarto | Secondary instabilities and vortex formation in collisionless-fluid magnetic reconnection[END_REF]). These jets collide at the O-point and then propagate along a direction perpendicular to the initial one. This behavior is inhibited in the asymmetric configuration, where only one side of the jets exhibits vortex formation and cannot propagate outward the separatrix, but is forced to move along the separatrix itself.

In Fig. 7 a snapshot of the advanced nonlinear phase shows that the turbulent behavior spreads out of the magnetic island when the deformation of the magnetic island separatrix at the O-point latitude is faster than the dynamics of the vorticity layers, which temporarily remain outside the island before being reabsorbed at later times. This vortices motion through the magnetic island separatrix opens a new transport scenario as suggested also in Refs. [START_REF] Le | Enhanced electron mixing and heating in 3d asymmetric reconnection at the earths magnetopause[END_REF][START_REF] Lee | Electron debye scale kelvin-helmholtz instability: Electrostatic particle-in-cell simulations[END_REF][START_REF] Price | The effects of turbulence on threedimensional magnetic reconnection at the magnetopause[END_REF].

In Ref. [START_REF] Lee | Electron debye scale kelvin-helmholtz instability: Electrostatic particle-in-cell simulations[END_REF] the electron KH instability was the mechanism allowing for the transport of particles, while in Refs. [START_REF] Le | Enhanced electron mixing and heating in 3d asymmetric reconnection at the earths magnetopause[END_REF][START_REF] Price | The effects of turbulence on threedimensional magnetic reconnection at the magnetopause[END_REF] the authors suggested that the transport of particles from the solar wind across the magnetopause into the magnetosphere is favoured by the turbulence generated by the lower hybrid drift instability driven by the steep density gradients at the magnetopause. We recall that in the cold ion limit adopted here the electron density coincides with the vorticity [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF], which also exhibits steep gradients, associated with the above described velocity profiles that generate the KH instability leading to turbulence in our model. We also note that previous 2D simulations described in the literature report a very weak turbulence [START_REF] Le | Enhanced electron mixing and heating in 3d asymmetric reconnection at the earths magnetopause[END_REF][START_REF] Price | The effects of turbulence on threedimensional magnetic reconnection at the magnetopause[END_REF] around the separatrix. Here we have shown that also in 2D configurations a significant level of turbulence may arise in the presence of a strongly asymmetric configuration.

C. Role of the location of the rational surface

We stress that the secondary KH instability previously observed in the asymmetric case is purely due to the geometry of the equilibrium magnetic field, α being the only parameter that has been changed between the two cases presented so far. Nonetheless we observed that, in general, a finite value of α is a necessary but not sufficient condition to trigger this instability. Indeed we performed different simulations keeping the same values of δ = 0.1 and β e = 0.05, but varying α, λ and the equilibrium magnetic field amplitude A. The results are summarized in Table I and discussed hereafter. By keeping A = 1, λ = 3 but with α = 1/8, i.e. by decreasing the value of the asymmetry parameter α with respect to the previously investigated asymmetric case, no instability has been found. This shows that, as above anticipated, a finite value of α is not sufficient to trigger the instability.

Assuming A = 1.35, but keeping α = 1/4 and λ = 3 we do not find any instability either, which indicates that the amplitude of the equilibrium magnetic field can suppress the instability. A simulation with A = 1, λ = 3 and α = 6/25 shows the occurrence of the KH instability, confirming the crucial role of the α parameter. The simulations with A = 0.67, λ = 2, α = 1/4, and A = 0.33, λ = 1, α = 1/4 lead again to the occurrence of the KH instability. We note that all the parameters α, A and λ contribute to the definition of the location of the resonant surface given in Eq. ( 10). However, the data in the fourth column of Table I show that the value of the distance of the rational surface from the location of the current density peak alone is not sufficient to determine whether the instability takes place or not. Indeed, the last case in the table exhibits the KH instability but is characterized by a value x s = 0.97, smaller than the values of x s of the second and third cases, which did not lead to a KH instability. It emerges that, among the seven cases discussed here, only those that have a value of x s /λ sufficiently large, exhibit the turbulent behavior induced by the instability, while the other cases persist in a laminar behavior throughout all the advanced nonlinear phase. To understand the physical reason of the observed difference, it is interesting to compare the profiles of the y-component of the velocity field shown in Fig. 5 for different values of the resonant surface of the corresponding equilibrium. In particular, we selected one stable case and two cases leading to instability. The profiles have been normalized to their maximum value and have been shifted along the x axis in such a way that their maxima are located at the same value of the x coordinate. We note that, as the resonant surface goes further away from x = 0, the width of the peak narrows and the gradient on the right side becomes steeper. In particular, this makes the part of the layer inside the magnetic island more prone to the KH instability as the value of x s /λ exceeds a certain threshold, in agreement with the subsequent development of the instability shown by the simulations. and on the ratio x s /λ. One observes that the instability takes place if x s /λ is sufficiently large.

α λ A x s x s /λ KH inst. 0 3 1 0 0 No 1 8 3
All these considerations permit to formulate a further condition that must be satisfied for an asymmetric case to develop the KH instability. The resonant surface must be located far enough from the x = 0 surface. How to define what is "enough" is a subtle question, which requires a more detailed parametric analysis, which is beyond the scope of this paper. However we note that, for the simulation campaign carried out here, the x s threshold value is of the order of the equilibrium scale length, i.e.

x s /λ ∼ 1 appears to be an estimate for the instability threshold in the adiabatic case.

Although the results reported here have been obtained keeping constant the box length 2L y = 24π, other simulations carried out with 2L y = 12π have shown an analogous behavior.

D. Effect of Landau dissipation on the secondary instabilities

In this Subsection we analyze the effect of the dissipation by switching on the Landau damping term in a simulation with the same parameters of the case with α = 1/4 reported in Subsection III B. We recall that α = 1/4 corresponds to k z = 1/48 ≈ 0.021, according to the mapping with the single helicity problem mentioned in Sec. III A. We confirm that for such a small value of k z the effect on the linear growth rate is almost negligible, as already shown in Ref. [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF].

The major effects of the Landau damping term can be appreciated when looking into the nonlinear phase of the reconnection process, where the presence of dissipation inhibits the secondary instabilities observed in the adiabatic case. Indeed the evolution of the vorticity field remains laminar until the end of the simulation. In Fig. 6 a comparison of the y component of the velocity profiles is shown for the two cases under consideration, at a time just before the onset of the KH instability occurring in the adiabatic case. The profiles have been normalized to their maximum. It is clear that the Landau effect prevents the velocity peaks from narrowing and smooths their gradients across the separatrix (dashed line), thus eliminating the conditions that make the colliding jet unstable to the KH instability.

We find it appropriate also to comment our results with respect to the role played by the invariants that the model possesses in the adiabatic limit and whose conservation is violated by the presence of the Landau term. Indeed, the model ( 1)-(3), in the adiabatic limit q e = 0 can be cast in the form [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF] ∂ν

+ ∂t + [ϕ + , ν + ] = 0, (17)
∂ν - ∂t + [ϕ -, ν -] = 0, (18)
∂s ∂t + [ϕ, s] = 0, (19)
where

ν ± = A - 2 β e δ 2 ∇ 2 ⊥ A ± δ √ 3 (∇ 2 ⊥ ϕ + t e ), s = ∇ 2 ⊥ ϕ - t e 2 , (20)
ϕ ± = ϕ ± √ 3 δ A . (21)
The fields ν ± and s correspond to Lagrangian invariants advected by generalized incompressible velocity fields. The presence of such Lagrangian invariants is associated with the existence of infinite families of Casimir invariants given by

C ± = d 2 x f ± (ν ± ), C = d 2 x f (s), (22)
where f ± and f are arbitrary functions. In particular, for the initial conditions considered in the present analysis, as already mentioned, the relation t e = 2∇ 2 ⊥ ϕ holds. In this case one has

ν ± = A - 2 β e δ 2 ∇ 2 ⊥ A ± √ 3δ∇ 2 ⊥ ϕ, s = 0, (23)
ϕ ± = ϕ ± √ 3 δ A . (24)
Thus, only the two Lagrangian invariants ν ± lead to a non trivial dynamics. In previous works, such as those of Refs. [START_REF] Sarto | Electron parallel compressibility in the nonlinear development of two-dimensional collisionless magnetohydrodynamic reconnection[END_REF][START_REF] Grasso | Phase mixing and satu-ration in Hamiltonian reconnection[END_REF], the isothermal two-field model reduction of our threefield model was used. Such two-field model also admits a reformulation in terms of two Lagrangian invariants analogous to ν ± . In these References it was shown that the system tends to form small scale structures concentrated in particular along the separatrices of the magnetic island, caused by the phase mixing process of these invariants. This prevents the formation of strong gradients representing the necessary condition for a secondary KHlike instability. However, when the system cannot be completely reformulated in terms of advection equations for Lagrangian invariants, for instance due to the assumption of cold electrons [START_REF] Sarto | Secondary instabilities and vortex formation in collisionless-fluid magnetic reconnection[END_REF] or of higher β e [START_REF] Grasso | Coupling between reconnection and Kelvin-Helmholtz instabilities in collisionless plasmas[END_REF][START_REF] Tassi | Hamiltonian four-field model for magnetic reconnection: nonlinear dynamics and extension to three dimensions with externally applied fields[END_REF], the secondary instability can take place. Our conclusion in Sec. III B can then be re-expressed in this context, stating that a secondary KH-like instability can take place even if the system can be recast in the form of advection equations for Lagrangian invariants. This can happen if a sufficiently asymmetric equilibrium (in the sense described in Sec. III B) is adopted.

On the other hand, in the presence of the Landau damping term, a reformulation of the model in the form ( 17)-( 19) is no longer possible. For instance, one can easily see that the Landau term tends to damp temperature fluctuations, without affecting in the same manner the vorticity [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF]. As a consequence, the difference between t e and ∇ 2 ⊥ ϕ is no longer constant along the flow, which expresses the fact that s is no longer a Lagrangian invariant.

As mentioned above, the Landau damping term appears to inhibit the KH instability.

Clearly, this inhibition is of different nature with respect to that due to the phase mixing of the Lagrangian invariants observed in the symmetric case as it is due to the smoothing of the velocity gradients that form in the asymmetric adiabatic case.

In order to have an estimate of the impact of Landau damping on the non-conservation of Casimir invariants, we considered the temporal evolution of the quantities d 2 x ν+ and d 2 x ν 2 + both in the adiabatic and Landau case. Because these two integral quantities are of the form C + indicated in Eq. ( 22), they are Casimir invariants in the adiabatic limit.

It turns out that, for our choice of the equilibrium parameters, the conservation of the quantities d 2 x ν + and d 2 x ν 2 + does not differ significantly between the adiabatic and the Landau cases. As shown in Fig. 8, the curves are almost coincident. This behavior reflects also in the feeble differences between the energy distributions (not shown here) and mainly depends on the strength of the sheared equilibrium magnetic field (for a higher equilibrium magnetic field strength the Landau damping effect on the energy is more evident, see [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF]).

However, while the Landau damping seems not to play a significant role in the conservation of integral quantities, it significantly affects the plasma behavior locally by preventing the Kelvin-Helmoltz instability onset.

Finally, we would like to point out that the linearized Landau-fluid model exhibits a very good agreement, in terms of the estimation of the linear growth rates for reconnecting modes, with the analytical theory based on a hybrid fluid-kinetic model [START_REF] Tassi | A reduced Landaugyrofluid model for magnetic reconnection driven by electron inertia[END_REF]. This could justify the choice of the Landau-fluid closure as an effort in the direction of providing fluid models retaining features of the kinetic theory. On the other hand, the quality of the agreement between the Landau-fluid and the kinetic dynamics in the nonlinear phase is a more delicate issue. Also, when considering the present Landau-fluid model as a step toward kinetic models, it should be pointed out that, as mentioned above, the Landau damping term in the fluid model breaks the Hamiltonian structure, which, on the other hand, is present in the drift-kinetic model from which the Landau-fluid model can be derived. In particular, as shown in Ref. [START_REF] Liseikina | Foliation and mixing of the electron drift-kinetic distribution function in nonlinear two-dimensional magnetic reconnection[END_REF], according to a Hamiltonian drift-kinetic model, the phase space of gyrocenters can be foliated into leaves, each hosting a dynamics analogous to that of the Lagrangian invariants of the fluid model takes place. Thus, although the Landau-fluid closure provides in general a better matching with the linear kinetic theory with respect to the adiabatic closure, one has to be aware of possible discrepancies emerging in the advanced nonlinear phase.

FIG. 1: Time evolution of the perturbed magnetic flux at the X-point for the symmetric (dashed line) and the asymmetric case (solid line). The presence of asymmetry decreases the linear growth rate. The drift along the x coordinate (represented on the right axis) of the X-point location is plotted in blue. This drift, typical of the asymmetric case, shows how the X-point tends to migrate toward the x = 0 line. The blue circles represent the approximation of the X-point displacement according to the velocity field in Eq. ( 14). The fit between the curves is restricted to a limited section of the total simulated time, because the post-processing treatement of the data does not allow a reliable evaluation of the F field gradient across the X-line (which appear at the denominator of the Eq. ( 14)) in the advanced nonlinear phase.

FIG. 2: Contour plots of the vorticity ∇ 2 ⊥ ϕ, at comparable times in terms of the island width, for different values of the asymmetry parameters. The left panel corresponds to a symmetric case, α = 0, while the right panel corresponds to an asymmetric case with α = 1/4. In both cases vorticity concentrates along the separatrices but in the symmetric case an elliptic structure is visible also inside the island. Note that the asymmetric layers are more intense.

IV. CONCLUSIONS

In this paper we analyzed how the intense vorticity layers, that form when a MR event takes place, evolve in an asymmetric low β e 2D configuration. In the finite electron temperature regime with ρ s /d e > 1 that we considered here, these layers are aligned along the separatrices of the magnetic island and have a laminar behavior in the symmetric configuration [START_REF] Sarto | Electron parallel compressibility in the nonlinear development of two-dimensional collisionless magnetohydrodynamic reconnection[END_REF]. We have shown that this paradigm breaks down when an asymmetric equilibrium magnetic field is taken into account. In the asymmetric configuration the vorticity may undergo a secondary instability of the KH type, which drives a turbulent plasma behavior around the magnetic island. In this case the velocity gradients, which develop firstly along the branch of the separatrix close to the x = 0 axis, tend to peak like an asymmetric Bickley jet that is well known to be prone to secondary instabilities. However, the presence of an asymmetry in the equilibrium configuration turns out to be a necessary but not sufficient condition for the occurrence of secondary instabilities. We have shown that the actual driving parameter for the onset of such an instability depends on the radial position of the resonant surface in the linear phase, i.e. in our 2D setting, where the sheared component of the equilibrium magnetic field goes to zero. In particular, for the instability to take place, the value of the rational surface, x s , must be of the order of the equilibrium scale length. Moreover in the asymmetric configuration the turbulence is not confined into the magnetic island but is allowed to spread out of the separatrix, opening transport channels otherwise closed. Our results highlight that a fully developed turbulent behavior can be achieved in 2D asymmetric configurations, differently from previous results mainly obtained with PIC simulations, where a strong turbulent behavior was observed only in 3D settings [START_REF] Le | Enhanced electron mixing and heating in 3d asymmetric reconnection at the earths magnetopause[END_REF][START_REF] Price | The effects of turbulence on threedimensional magnetic reconnection at the magnetopause[END_REF].

We also analyzed the stabilizing effect of a non-adiabatic closure on these secondary insta-FIG. 5: Profiles of v y for different values of the resonant surface of the equilibrium magnetic field.

The profiles have been normalized to their maximum value and have been shifted along the x axis so that the location of their maxima coincides at the same x coordinate. Note the narrowing of the peak as the value of x s increases. The cases with x s = 1.95 and x s = 2.92 correspond to values of x s /λ sufficiently large to yield the KH instability.

bilities, by including a Landau damping term in our model. This term prevents the velocity peaks from narrowing, so inhibiting their destabilization.

In the future we intend to address this problem in a full 3D setting, where several modes can be nonlinearly excited driving a secondary instability on their own. Moreover, the presence of several modes may influence this scenario because of the magnetic chaos that can be generated by their interaction.

Another interesting aspect worth of a deeper attention concerns the particle density, which, in our model, is proportional to the vorticity. In the minimal asymmetric framework where we have described the onset and the behavior of the KH instability, the equilibrium density is uniform and, considering the magnetospheric context, the equilibrium magnetic asymmetry depends only on the direction of the solar wind. In other works with a different equilibrium configuration, equilibrium density gradients are assumed and secondary Rayleigh-Taylor modes associated with the density inhomogeneities are induced by the KH instability in the nonlinear phase ( [START_REF] Faganello | Magnetized KelvinHelmholtz instability: theory and simulations in the Earths magnetosphere context[END_REF][START_REF] Matsumoto | Onset of turbulence induced by a KelvinHelmholtz vortex[END_REF] and references therein). In Ref. [START_REF] Faganello | Magnetized KelvinHelmholtz instability: theory and simulations in the Earths magnetosphere context[END_REF], the role of the small-scale KH structures in the mass transfer between the solar wind and the magnetosphere is also discussed. In this perspective, a further development of our investigation is represented by the particle transport in the turbulent region produced by the KH instability and, in particular, across the island's boundaries, which could be carried out through a test particle approach ( [START_REF] Perona | A test electron model for the study of three dimensional magnetic reconnection effects[END_REF] and references therein). In conclusion, the use of the presently adopted reduced fluid model implies evident limitations, in terms of the range of applicability, when compared to kinetic models. In order to reduce the gap with kinetic descriptions and in view of possible applications to asymmetric reconnection events in the magnetosphere, we also intend to adopt more refined fluid models valid for larger β e values and including finite ion Larmor radius effects as well as equilibrium temperature anisotropies.

FIG. 3 :

 3 FIG. 3: Profiles of v x , v y (top panels), B x and B y (bottom panels) plotted as a function of x at y = 10. The symmetric case, α = 0, is on the left and the asymmetric case, α = 1/4, is on the right. The dashed lines represent the location of the two branches of the separatrix for the chosen y value.

FIG. 4 :

 4 FIG. 4: Contour plots of the vorticity at three subsequent times: t = 300 (top frame), t = 306 (middle frame), t = 318 (bottom frame). The arrows indicate the direction of propagation of the instability front. The dashed line in each figure shows the magnetic island separatrix at the corresponding time.

FIG. 6 :

 6 FIG.6: Profiles of v y , normalized to their maximum, comparing the effect of the Landau damping term with the adiabatic case. The values of the parameters are α = 1/4, λ = 3, A = 1 for both cases. The presence of the Landau damping leads to weaker velocity gradients.

FIG. 7 :

 7 FIG.7: Contour plots of the vorticity at three subsequent times: t = 354 (top frame), t = 360 (middle frame), t = 375 (bottom frame). The white dashed lines identify the corresponding magnetic island separatrices.

FIG. 8 :

 8 FIG. 8: Time evolution of the quantities C 1 = d 2 x ν + (black lines) and C 2 = d 2 x ν 2 + (blue lines) for the adiabatic (solid lines) and the Landau (dashed lines) cases. The times are normalized with respect to the linear growth rate of the reconnection process. The values of the parameters are α = 1/4, λ = 3, A = 1 for both cases.

TABLE I :

 I Table summarizing whether the secondary KH instability occurs or not, depending on the values of the equilibrium parameters α, λ and A, on the location of the resonant surface x s

	1 1.18 0.39	No
	1 4 3 1.35 1.88 0.63	No
	6 25 3 1 2.72 0.91	Yes
	1 4 2 0.67 1.95 0.97	Yes
	1 4 3 1 2.92 0.97	Yes
	1 4 1 0.33 0.97 0.97	Yes

Acknowledgments

Computations have been performed on the Mesocentre SIGAMM machine, hosted by Observatoire de la Côte d'Azur. The authors thank P.L. Sulem and T. Passot for providing useful hints and for discussing the results. DG and DB thank prof. B. Coppi and the Ignitor theory group at CNR for stimulating discussions.

Supplementary material

See the movie in the supplementary material for a complete description of the evolution of the KH instability. The movie shows the time evolution of the vorticity spatial distribution U (x, y, t) = ∇ 2 ⊥ ϕ(x, y, t) for the case α = 1/4, A = 1 and λ = 3 in the left frame. The superposed dashed lines correspond to the magnetic island separatrices. The time evolution of A X is plotted on the right frame, in order to highlight the link between KH and MR instability.