J Barré

D Métivier

Y Y Yamaguchi

Trapping scaling for bifurcations in Vlasov systems

Keywords: numbers: 05, 45, -a Nonlinear dynamics and chaos 02, 30, Oz Bifurcation theory

We study non oscillating bifurcations of non homogeneous steady states of the Vlasov equation, a situation occurring in galactic models, or for Bernstein-Greene-Kruskal modes in plasma physics. Through an unstable manifold expansion, we show that in one spatial dimension the dynamics is very sensitive to the initial perturbation: the instability may saturate at small amplitude -generalizing the "trapping scaling" of plasma physics-or may grow to produce a large scale modification of the system. Furthermore, resonances are strongly suppressed, leading to different phenomena with respect to the homogeneous case. These analytical findings are illustrated and extended by direct numerical simulations with a cosine interaction potential.

I. INTRODUCTION

The Vlasov equation models the dynamics of a large number of interacting particles, when the force acting on them is dominated by the mean-field. It is a fundamental equation of plasma physics and galactic dynamics, but also describes a number of other physical systems, such as wave-particles interactions [START_REF] Elskens | Microscopic dynamics of plasmas and chaos IoP[END_REF], free electron lasers [START_REF] Bonifacio | Collective instabilities and high-gain regime in a free electron laser[END_REF], certain regimes of nonlinear optics [START_REF] Picozzi | Optical wave turbulence: Towards a unified nonequilibrium thermodynamic formulation of statistical nonlinear optics[END_REF], and wave propagation in bubbly fluids [START_REF] Smereka | A Vlasov equation for pressure wave propagation in bubbly fluids[END_REF]. Understanding its qualitative behavior has proved to be a formidable challenge, from the first linear computations of Landau [START_REF] Landau | On the vibrations of the electronic plasma[END_REF] to the most recent mathematical breakthroughs [6]. We address here a problem of bifurcations in these systems.

The study of bifurcation of homogeneous stationary solutions to the Vlasov equation (when particles are homogeneously distributed in space, with a certain velocity distribution) has a long and interesting history, mainly related to plasma physics; a detailed account can be found in [START_REF] Balmforth | Pattern formation in Hamiltonian systems with continuous spectra; a normal-form single-wave model[END_REF]. One of the main question is to determine the asymptotic behavior of the system, when started close to an unstable stationary state. The basic mechanism, ultimately responsible for the instability saturation, is the resonance phenomenon: particles with a velocity close to that of the growing mode are strongly affected by the small perturbation. This was qualitatively recognized in the 1960's [START_REF] Baldwin | Perturbation Method for Waves in a Slowly Varying Plasma[END_REF][START_REF] Neil | Collisionless Damping of Nonlinear Plasma Oscillations[END_REF], leading these authors to predict the famous "trapping scaling", that is the asymptotic electric field, after nonlinear saturation, would be of order λ 2 , with λ the instability rate. However, formal computations relying on standard nonlinear expansions often led to the very different λ 1/2 scaling (see discussion and refs in [START_REF] Balmforth | Pattern formation in Hamiltonian systems with continuous spectra; a normal-form single-wave model[END_REF]). In the 1990's, Crawford's careful unstable manifold computations finally recovered unambiguously the trapping scaling [START_REF] Crawford | Amplitude equations for electrostatic waves: Universal singular behavior in the limit of weak instability[END_REF]. Yet, his expansion is plagued by strong divergences in the λ → 0 + limit. Soon after, Del-Castillo-Negrete derived an infinite dimensional reduced model for the instability development, using different techniques [START_REF] Del-Castillo-Negrete | Nonlinear evolution of perturbations in marginally stable plasmas[END_REF]. This is currently the best understanding we have, and shows a striking degree of universality.

Understanding the non homogeneous case is further challenging, and has a practical importance for various physical systems. A first key example comes from astrophysics. Radial Orbit Instability destabilizes a spherical self-gravitating system when the number of low angular momentum particles, usually stars, grows. The fate of the initially unstable system has been qualitatively and numerically studied in [START_REF] Palmer | Neighbouring equilibria to radially anisotropic spheres: possible end-states for violently relaxed stellar systems[END_REF], providing a rare example of such a bifurcation study. In the plasma physics context, Bernstein-Greene-Kruskal exhibited a large family of exact solutions of the Vlasov-Poisson equation [START_REF] Bernstein | Exact Nonlinear Plasma Oscillations[END_REF], now commonly called BGK modes. Some of these solutions are stable, others unstable [START_REF] Schwarzmeier | Stability of Bernstein-Greene-Kruskal equilibria[END_REF][START_REF] Manfredi | Stability of Bernstein-Greene-Kruskal modes[END_REF][START_REF] Lin | Nonlinear instability of periodic BGK waves for Vlasov-Poisson system[END_REF][START_REF] Pankavich | Instability conditions for some periodic BGK waves in the Vlasov-Poisson system[END_REF]. The non linear behavior in case of instability is an open question.

To tackle the bifurcation in the non homogeneous case we will consider in this article as a first step the case of a real unstable eigenvalue, for a one dimensional system. We will show that the number of particles around zero frequency is negligible in the non homogeneous case compared to the homogeneous case. This is a crucial remark, since it implies that for non oscillating instabilities, the resonance phenomenon will be absent, or strongly suppressed. Nevertheless, similarly to the homogeneous case, the unstable eigenvalue still bifurcates from a marginally stable continuous spectrum. Thus, one should expect a new type of bifurcation. Our strategy is to attack this problem by revisiting Crawford's expansion, in order to describe the dynamics on the unstable manifold. The basic questions are: Does the instability saturate at small amplitude? What is the analog of the "trapping scaling"? Can we expect the same kind of universality as in the homogeneous context? Our main results include: i) A striking asymmetry on the unstable manifold: one branch leads to a small amplitude saturation, the other escapes far away from the original stationary point. ii) The λ 2 "trappping scaling", well known in the homogeneous case, generalizes for the small amplitude saturation; however, the reason for this scaling is very different. This result can also be reached by a different theoretical approach, following [START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF]. iii) Since these features depend on some generic physical characteristic of the problem, we expect our results to be generic, at least for 1D systems.

The outline of the article is as follows: Section II deals with a general one-dimensional Vlasov equation. We first explain formally how to obtain the unstable manifold expansion and emphasize the main physical points (Section II A), then present a more detailed computation (Section II B). In Section III, we restrict to a particular cosine interaction potential. We are then able to perform more explicit computations, and provide precise predictions, as well as push the computations at higher orders, which is crucial to understand the divergence pattern of the expansion (Section III B). We also show how some of these predictions can also be obtained via the "rearrangement formula" idea, recently introduced [START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF]; this provides an independent theoretical approach (Section III C). In Section IV, we compare all previous predictions with direct numerical simulations. The choice of the cosine potential allows us to perform very precise computations. Section V is devoted to some final comments and open questions.

II. UNSTABLE MANIFOLD REDUCTION (SPATIALLY PERIODIC 1D SYSTEMS)

A. General picture

The Vlasov equation associated with the one-particle Hamiltonian H[f ] is

∂ t f = {f, H[f ]}, H[f ](q, p, t) = p 2 /2 + V [f ](q, t) (1)
where the one-particle potential V [f ] is defined by using the two-body interaction v(q) as

V [f ](q, t) = v(q -q ′ )f (q ′ , p ′ , t)dq ′ dp ′ (2) and the Poisson bracket {•, •} is {u, v} = ∂ p u∂ q v -∂ q u∂ p v. (3)
Denote by f µ a family of unstable stationary states depending on the parameter µ, and by g a perturbation around f µ . The equation for g consists of a linear and a nonlinear parts:

∂ t g = L µ g + N [g], (4)
where

L µ g = {g, H[f µ ]} + {f µ , V [g]} , N [g] = {g, V [g]}. ( 5
)
Let L µ have an unstable real nondegenerate eigenvalue λ depending on µ, and Ψ λ be the corresponding eigenfunction. We assume that the family {f µ } has a critical point µ c , where the eigenvalue is λ = 0, and we will study the system in the limit µ → µ c which implies λ → 0 + . Notice that the Hamiltonian structure implies that L µ also has the eigenvalue -λ. In addition, by space translation symmetry when f µ is homogeneous in space, λ is in this case twice degenerate. When f µ is non homogeneous, the symmetry is broken, and a neutral Goldstone mode Ψ N appears.

We use the standard L 2 Hermitian product, and introduce L † µ , the adjoint operator of L µ , and Ψ λ the adjoint eigenfunction with eigenvalue λ * = λ, where λ * is the complex conjugate of λ. We normalize it as Ψ λ , Ψ λ = 1, where the scalar product is defined by

f, g = f * g dqdp, (6)
and note Π µ (•) = Ψ λ , • Ψ λ the projection onto the unstable eigenspace Span{Ψ λ }. Depending on the behavior of Ψ ±λ and Ψ N when λ → 0 + , Π µ can be singular in this limit. As will become clear later on, this is the origin of an important difference between homogeneous and non homogeneous cases. We are interested in the dynamics on the unstable manifold of f µ , and expand g as g(q, p, t) = A(t)Ψ λ (q, p) + S(q, p, A(t)).

The term S parameterizes the unstable manifold, and we have the estimation S = O(A 2 ) for small A by assuming that the unstable manifold is tangent to the unstable eigenspace. Note that this parametrization may be only local.

Performing the projection, we obtain for A(t) the following dynamical equation, expanded in powers of A:

Ȧ = λA + Ψ λ , N [f ] = λA + a 2 (λ)A 2 + a 3 (λ)A 3 + • • • , (8) where the coefficient a 2 is a 2 (λ) = Ψ λ , {Ψ λ , V [Ψ λ ]} . (9)
In the homogeneous case, the unstable eigenspace is actually two-dimensional, and the reduced equation involves A, which is complex, and A * ; since an equation similar to (8) can be written for |A|, we omit this slight difference here. In this case we have by symmetry a 2l = 0 for any l ∈ N * and a 2l+1 ∝ λ 1-4l [START_REF] Crawford | Amplitude equations for electrostatic waves: Universal singular behavior in the limit of weak instability[END_REF], this strong divergence preventing the truncation of the series [START_REF] Baldwin | Perturbation Method for Waves in a Slowly Varying Plasma[END_REF]. By contrast, in the non homogeneous case, a 2 = 0; we show however below that a similar divergence pattern when λ → 0 + is expected. Nevertheless, keeping term up to A 2 , (8) looks like the normal form of a transcritical bifurcation; it is not a standard one however, because λ is always positive, and a 2 diverges when λ → 0 + . From the truncated equation, we can conjecture the general phenomenology:

• on the unstable manifold, in one direction the dynamics is bounded and attracted by a stationary state close to f µ , characterized by A ∞ = -λ/a 2 (λ) (notice that this state is stationary and stable for the truncated unstable manifold dynamics; it may not be so for the unconstrained dynamics);

• in the other direction, the dynamics leaves the perturbative regime, and thus the range of validity of our analysis.

Estimating a 2 (λ) in the λ → 0 + limit gives access to A ∞ ; this requires more technical work performed in the next subsection. Note that the formal computations leading to ( 8) and ( 9) are essentially independent of the dimension; the following estimate for a 2 (λ) is valid in 1D.

B. Technical computation

We present now the explicit form of a 2 (λ), the order A 2 coefficient in the reduced dynamics, for the non homogeneous case of a generic 1-dimensional system. The computation is somewhat technical, but the physical message is rather simple, and summarized here:

• a 2 ∝ 1/λ when λ → 0 + • the normalization condition requires Ψ λ = O(1/λ),
and this divergence of the amplitude in the neighborhood of λ = 0 is responsible for the divergence of a 2 .

• the integrals appearing in the scalar product in [START_REF] Neil | Collisionless Damping of Nonlinear Plasma Oscillations[END_REF] do not diverge when λ → 0 + .

The last two points are major differences with the homogeneous case in which the amplitude of Ψλ does not diverge, and the resonances translate into pinching singularities for the integrals analogous to (9) [START_REF] Crawford | Amplitude equations for electrostatic waves: Universal singular behavior in the limit of weak instability[END_REF]. These pinching singularities are responsible for the divergences in the unstable manifold expansion. From these results, we can now conclude that the stationary state close to f µ on the unstable manifold is characterized by

A ∞ ∝ λ 2 ,
which is a central result of this article.

The spectrum matrix

Let us turn to the explicit computation. We first need to introduce quantities regarding the linearized problem; we follow a classical method [START_REF] Kalnajs | Dynamics of flat galaxies. IV -The integral equation for normal modes in matrix form[END_REF][START_REF] Lewis | Linearized analysis of inhomogeneous plasma equilibria: General theory[END_REF] and use here the notations of [START_REF] Barré | Algebraic damping in the one-dimensional Vlasov equation[END_REF]. For preparation, we introduce the biorthogonal functions {d i } i∈I and {u k } k∈K , on which we expand the density and potential respectively:

f (q, p, t)dp = i∈I a i (t)d i (q), (10)
V [f ](q, t) = k∈K a k (t)u k (q). (11)
The functions {d i } and {u k } satisfy

v(q -q ′ )d k (q ′ )dq ′ = u k (q), (12)
and

(d i , u k ) = d i (q) * u k (q)dq = ν k δ ik , (13)
Since all the functions expanded are real, we can choose the d i 's and u k 's to be real functions; the ν k 's are then also real. In a generic system the spectrum function is a matrix denoted by Λ µ (λ) [START_REF] Kalnajs | Dynamics of flat galaxies. IV -The integral equation for normal modes in matrix form[END_REF][START_REF] Lewis | Linearized analysis of inhomogeneous plasma equilibria: General theory[END_REF], whose (l, k) element is

(Λ µ ) lk (λ) = ν k δ lk -2π m∈Z im∇ J F µ λ + imΩ c k,m c * l,m dJ, ( 14
)
where we have introduced the action-angle variables (J, θ) for an integrable Hamiltonian system H[f µ ] with 1 degree of freedom, and assumed that f µ depends only on J; this is denoted by

f µ (q, p) = F µ (J). The function c k,m (J) is defined by c k,m (J) = 1 2π u k (q)e -imθ dθ, (15)
and the frequency Ω(J) by Ω = dH[f µ ]/dJ. We remark that det Λ µ (λ) = 0 implies λ to be the eigenvalue of the linear problem. The matrix's size is #K, the number of biorthogonal functions in the expansion of the potential; it is usually infinite, but only 2×2 for the cosine potential introduced later.

Properties of the spectrum matrix

To analyze Λ, we need to specify a little bit more the frequency Ω(J). We assume that Ω(J) does not vanish at finite J, or does so only logarithmically. This includes the cases where the stationary potential V [f µ ](q) has a single minimum and is infinite for |q| infinite (such as for 1D gravity), and the generic situation with periodic boundary conditions; indeed, in the latter situation, local minima of the stationary potential give rise to separatrices, on which the action is constant. At these specific values of the action Ω vanishes, but generically it does so only logarithmically (an example of this situation is given in the next section). In particular, under the above assumption for Ω, the functions (Λ µ ) lk (λ) defined in [START_REF] Schwarzmeier | Stability of Bernstein-Greene-Kruskal equilibria[END_REF] are well defined for λ = 0 (the integral over J converges), and continuous for λ ∈ R. Furthermore, the same reasoning applies to the derivatives of (Λ µ ) lk (λ): these functions are continuous for λ ∈ R, even for λ = 0. The fact that integrals over J converge is the technical counterpart of a physical phenomenon: the absence, or weakness, of resonances for λ = 0. Furthermore, starting from ( 14), changing the summation variable m to -m, and using c k,-m (J) = c * k,m (J), it is not difficult to see that Λ µ is a real matrix, and that (Λ µ ) lk (λ) = (Λ µ ) kl (-λ). Hence Λ µ (0 + ) is symmetric. This will be useful below.

Comparison with the homogeneous case

The properties of the spectrum matrix detailed in the above paragraph are in sharp contrast with the better known situation for a homogeneous stationary state. It is useful to highlight the comparison. We assume here periodic boundary conditions. For a homogeneous stationary state, a Fourier transform with respect to the space variable diagonalizes the spectrum matrix. The diagonal elements are [START_REF] Crawford | Amplitude equations for electrostatic waves: Universal singular behavior in the limit of weak instability[END_REF]:

Λ (hom) k (λ) = 1 -2iπkv k F ′ 0 (p) λ + ikp dp, (16)
where F 0 (p) is the momentum distribution of the considered stationary state, and vk the Fourier of the two-body potential. Note that since 1/p is not integrable close to p = 0, (16) may not be well defined for λ = 0, and in any case is not differentiable along the imaginary axis.

In particular, the function defined by ( 16) for λ ∈ R is not differentiable for λ = 0. This is the technical counterpart of the strong resonance phenomenon between a non oscillating perturbation and particles with vanishing velocities. We come back now to the non homogeneous case.

Eigenfunctions

The eigenfunctions with eigenvalue λ of L µ and L † µ (recall that λ is real) are respectively expressed as

Ψ λ = k∈K m∈Z im∇ J F µ λ + imΩ c k,m e imθ b k (17)
and

Ψ λ = - k∈K m∈Z c k,m λ -imΩ e imθ bk , (18)
where the real vectors b = (b k ) and b = ( bk ) depending on λ must satisfy the conditions

Λ µ (λ)b(λ) = 0, T b(λ)Λ µ (λ) = 0. ( 19
)
where the upper T represents the transposition. The assumption that λ is a nondegenerate eigenvalue implies that dim KerΛ µ (λ) = 1. Otherwise, there would be two linearly independent eigenfunctions for λ.

Analysis of the a2 coefficient

From the above eigenfunctions the a 2 (λ) term is written as

a 2 (λ) = -2π j∈K k∈K l∈K bj b k b l × m∈Z n∈Z c * j,m+n ϕ m,n λ + i(m + n)Ω dJ (20)
where

ϕ m,n = in∇ J im∇ J F µ c k,m λ + imΩ c l,n - im∇ J F µ c k,m λ + imΩ im∇ J c l,n . (21)
We study the λ dependence of the a 2 (λ) term in the limit λ → 0 + . First, we estimate the order of magnitude of b. The amplitude of b can be chosen arbitrarily but b cannot, since Ψλ must satisfy the normalization condition, which is expressed as

1 = - k∈K l∈K bl (Λ µ ) ′ lk (λ)b k = -b, Λ ′ µ (λ)b , (22)
with •, • the standard scalar product. Making use of the facts that the (Λ µ ) kl are regular around λ = 0 + , the normalization condition reads

-1 = b, Λ ′ µ (0)b 0 + λ b, Λ ′′ µ (0)b 0 + Λ ′ µ (0)b 1 + O(λ 2 ), (23)
where we assumed that b can be expanded in λ:

b(λ) = b 0 + λb 1 + O(λ 2 ). ( 24
)
We will show that the amplitude of b(λ) diverges in the limit λ → 0 + , and accordingly, we assume that b(λ) can be expanded as

b(λ) = λ -y Ỹ (λ), Ỹ (λ) = b0 + λ b1 + O(λ 2 ). ( 25
)
The normalization condition [START_REF] Campa | Statistical mechanics and dynamics of solvable models with long-range interactions[END_REF] together with (24) imposes y ≥ 0. First, note that the definition of b (19) reads, at order λ:

Λ µ (0)b 1 = -Λ ′ µ (0)b 0 . (26)
Since Λ µ (0) is symmetrical, this implies that Λ ′ µ (0)b 0 is orthogonal to the kernel of Λ µ (0). Now, writing at leading order the definition of b [START_REF] Kalnajs | Dynamics of flat galaxies. IV -The integral equation for normal modes in matrix form[END_REF], we obtain Λ µ (0) b0 = 0. Hence, we conclude that b0 , Λ ′ µ (0)b 0 = 0, and the first term of ( 23) is estimated as

b(λ), Λ ′ µ (0)b 0 = O(λ 1-y ). ( 27
)
The second term of ( 23) is also estimated in the same ordering as

λ b, Λ ′′ µ (0)b 0 + Λ ′ µ (0)b 1 = O(λ 1-y ), (28)
and the remaining terms are higher than O(λ 1-y ). Consequently, to satisfy the normalization condition [START_REF] Campa | Statistical mechanics and dynamics of solvable models with long-range interactions[END_REF], the exponent y should be unity and b = O(1/λ).

In order to show that b is the unique source of divergence in the a 2 (λ) term, we have to check two cases which may give an extra diverging factor 1/λ in the integrand of a 2 (λ) (see [START_REF] Lewis | Linearized analysis of inhomogeneous plasma equilibria: General theory[END_REF]): m = 0 and m + n = 0. It is easy to find that the contribution from the former vanishes. Expanding a 2 (λ), we find that the contribution from the latter is canceled out at leading order between ±m, and the integral gives a finite value. Therefore, we conclude that a 2 (λ) is of order 1/λ in the limit λ → 0 + .

III. AN EXPLICIT EXAMPLE: COSINE POTENTIAL (HMF MODEL)

The previous computations, for a generic 1D potential, are a bit abstract, and complicated enough already at order A 2 . In order to perform higher order computations, and to prepare for precise numerical tests, we specialize in this section to a particular interaction potential v(q) =cos q, with q ∈ [0, 2π[ and periodic boundary conditions; this model is sometimes called the Hamiltonian mean field (HMF) model [START_REF] Antoni | Clustering and relaxation in Hamiltonian long-range dynamics[END_REF][START_REF] Campa | Statistical mechanics and dynamics of solvable models with long-range interactions[END_REF]. All computations then become more explicit, which will be useful to analyze higher orders, see Sec.III B. In addition, we can also make use in this context of the "rearrangement" approach described in [START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF]: this will provide an alternative analytical approach to test the results, see Sec.III C.

A. Cosine potential: computation at order A 2

The mean-field potential created by a phase space distribution f is in this case V MF [f ](q) = -M x cos q -M y sin q with

M = M x + iM y = f (q, p)e iq dpdq.
M is usually called the "magnetization". We choose M µ , the magnetization of f µ , to be real, and assume a symmetric initial condition, such that M y (t) = 0 for all t (test simulations with M y (0) = 0 -not shown-showed a similar behavior). We also assume M µ = 0, that is f µ non homogeneous. The one particle dynamics in V MF is a simple pendulum dynamics, obviously integrable; we introduce the associated action-angle variables (J, θ). The biorthogonal functions in the expansion of the potential are {cos q, sin q} and the size of matrix Λ µ (λ) is Frequency Ω as a function of the action J in the HMF model. The action is Js at the separatrix (see the upper-left inset representing the phase space (q, p)), where Js is different if the separatrix is approached from below or from above [START_REF] Barré | Dynamics of perturbations around inhomogeneous backgrounds in the HMF model[END_REF]. Ω goes to Ω0 in the limit J → 0 + . The lower-right inset shows the logarithmic divergence of 1/Ω around the separatrix action Js for both inside (red solid upper line) and outside (blue broken lower line).

2 accordingly. The matrix Λ(λ) is diagonal: Λ µ (λ) = diag(Λ c µ (λ), Λ s µ (λ)) [START_REF] Barré | Dynamics of perturbations around inhomogeneous backgrounds in the HMF model[END_REF], where

Λ c µ (λ) = 1 + 2π m =0 imF ′ µ (J)|c m (J)| 2 λ + imΩ(J) dJ (29)
Λ s µ (λ) = 1 + 2π m =0 imF ′ µ (J)|s m (J)| 2 λ + imΩ(J) dJ (30)
and

c m (J) = 1 2π 2π 0 cos q(θ, J)e -imθ dθ, (31)
s m (J) = 1 2π 2π 0 sin q(θ, J)e -imθ dθ. ( 32
)
The spectrum function satisfies Λ s µ (0) = 0 [START_REF] Campa | A dynamical stability criterion for inhomogeneous quasi-stationary states in longrange systems[END_REF], which corresponds to the Goldstone mode. The critical point µ c is, therefore, determined by the other part and the root of the equation Λ c µ (λ) = 0 tends to 0 + in the limit µ → µ c . Hereafter we consider the cosine part Λ c µ only, and denote it by Λ µ for simplicity.

In the above setting the eigenfunctions are

Ψ λ = m∈Z imF ′ µ (J)c m (J) λ + imΩ(J) e imθ (33)
Ψ λ = - 1 Λ ′ µ (λ) * m∈Z c m (J) λ -imΩ(J) e imθ . ( 34
)
It is easy to see that Λ ′ µ (λ) ∼ Cλ for µ → µ c , inducing a divergence in Ψ λ . The possibility of resonances is readily seen on Eqs.( 33) and (34): they correspond to actions J such that λ ± imΩ(J) = 0. Since we have assumed λ ∈ R, this can happen only at the separatrix J = J s , where Ω(J s ) = 0 (see Fig. 1). The divergence is very weak, Ω(J) -1 = O(ln |J -J s |), thus resonances are strongly suppressed since the number of particles around J = J s is negligible. On the other hand, for a complex λ strong resonances are possible, a priori similar to those occurring in the homogeneous case.

For the cosine potential (HMF model), the above expression (20) for a 2 (λ) simplifies a little bit:

a 2 (λ) = -2π Λ ′ µ (λ) m,n∈Z c * m+n ϕ m,n (J) λ + i(m + n)Ω dJ ( 35
)
where Λ µ (λ) is given by ( 29), and

ϕ m,n (J) = m n ∂ ∂J F ′ µ c m λ + imΩ c n -m F ′ µ c m λ + imΩ c ′ n . ( 36
)
It is straightforward to see that a 2 (λ) is of order 1/λ in the limit λ → 0 + due to the divergent factor 1/Λ ′ µ (λ). Some details on how to compute numerically with a good accuracy the function Λ µ and its roots are given in the appendix.

B. Computations at higher orders

We explain here how to compute the higher order terms. While the structure is standard, the effective computations are intricate. The first paragraph is valid for a general potential, but in order to explicitly work out the order of magnitude of a 3 , we restrict to the cosine potential. We then show that a 3 ∼ C/λ 3 , as λ → 0 + .

Structure of the computation for a general potential

Let us parameterize the unstable manifold of f µ as

g = k≥1 S k A k , S 1 = Ψ λ . ( 37
)
Since the nonlinear part N [g] of the Vlasov equation,

∂ t g = L µ g + N [g]
, is bilinear, we write it as

N [g] = B(g, g) , with B(g, h) = ∂ J g∂ θ V [g] -∂ θ g∂ J V [g].
Recalling that the reduced dynamics on the unstable manifold is

Ȧ = k≥1 a k A k , a 1 = λ, ( 38
)
our goal is to provide a formal expression for a k . This requires computing at the same time the S k 's. We write the time evolution equation for g in two ways:

∂ t g = k≥1 l≥1 A k+l-1 ka l S k = k≥1 A k L µ S k + k≥1 l≥1 A k+l B(S k , S l ), (39)
where we have used (38). Picking up the terms order by order, we have for any k

(ka 1 -L µ )S k = k-1 l=1 [B(S k-l , S l ) -(k -l)a l+1 S k-l ] .
(40) This equation for k = 1 is simply L µ Ψ λ = λΨ λ , and we focus on k ≥ 2.

We now project these equations onto span{Ψ λ } and span{Ψ λ } ⊥ ; we note the corresponding projection operators Π and Π ⊥ = I -Π, where I is the identity. The projection operators work as

(Π • L µ )Ψ λ = λΨ λ , (Π ⊥ • L µ )Ψ λ = 0, (Π • L µ )S k = 0, (Π ⊥ • L µ )S k = L µ S k , k ≥ 2. ( 41
)
The projection operator Π induces

a k Ψ λ = k-1 l=1 Π • B(S k-l , S l ), (42) and Π
⊥ yields (kλ -L µ )S k = k-1 l=1 Π ⊥ • B(S k-l , S l ) - k-2 l=1 (k -l)a l+1 S k-l .
(43) Equation (42) determines a k from S l , and (43) S k from S l and a l with l ∈ {1, • • • , k -1}.

Cosine potential -Order of magnitude of a3

We now specialize to the cosine potential and focus on the third order coefficient a 3 .

Calling G k the r.h.s. of Eq.( 43), we have

S k = R(kλ) G k , (44)
where

R(z) = (z -L µ ) -1 is the resolvent of L µ . Our first task is to determine the λ dependence of R(kλ) G k .
Let us then consider the equation

(z -L µ )X = G ,
and solve for X(θ, J). Denoting the m-th Fourier component of X in θ as Xm , we find after some computations, for all m ∈ Z:

Xm = Ĝm z + imΩ(J) -C(z) imc m F ′ µ (J) z + imΩ(J) (45)
IV. DIRECT NUMERICAL SIMULATIONS

Our analytical description of the bifurcation can be accurately tested in the HMF case. The time-evolved distribution function is obtained via a GPU parallel implementation of a semi-Lagrangian scheme for the Vlasov HMF equation with periodic boundary conditions [START_REF] Rocha Filho | Solving the Vlasov equation for onedimensional models with long range interactions on a GPU[END_REF]. We use a 2 n × 2 n grid in position momentum phase space truncated at |p| = 2 with n up to 12; the time step is usually 10 -2 . We use two families of reference stationary states

F 0 µ (E) = N -1 F /[1 + e β(E-µ) ], G 0 µ (E) = N -1 G E 2 e -µE , (55)
with the one-particle energy E(q, p) = p 2 /2 + M µ (1cos(q)), where the magnetization M µ has to be computed self-consistently, and N F and N G are the normalization factors. Fig. 2 presents M µ and its bifurcations, for the family F 0 µ . For both families, a real positive eigenvalue appears at a critical value of µ: for F 0 µ , this is at point a, see Fig. 2. The initial perturbation is ǫT (q, p) = ǫ cos(q) exp(-βp 2 /2). Note that this is not proportional to the unstable eigenvector Ψ λ : this allows us to test the robustness of our unstable manifold analysis with respect to the initial condition. For all simulation results presented in Fig. 3, the size of the perturbation ǫ was chosen small enough such that the saturated solution reached for ǫ > 0 does not depend on ǫ. On the other hand, the smaller ǫ, the more accurate computations are required to avoid numerical errors. In particular, we have observed that numerical errors may drive the system far away from the reference stationary solution, following a dynamics similar to the one with ǫ < 0. In such cases, we have used a finer phase space grid: GPU computational power was crucial to reach very fine grids. For example the initial Fermi distribution F 0 µ was very sensitive to numerical errors and to ǫ: we took ǫ = 1.8 • 10 -6 with a 4096 × 4096 grid; for G 0 µ distribution, which is much smoother, ǫ = 10 -5 and a 1024 × 1024 grid was enough (except for the point corresponding to λ = 0.032 where more precision was needed, and we took ǫ = 3 • 10 -6 ).

Typical evolutions for the order parameter M (t) are shown in Fig. 3 (insets), for positive and negative ǫ. The asymmetry is clear: for one perturbation the change in δM = M (t) -M µ remains small, for the other it is O(1). In the case where δM remains small, we compute its saturated value by averaging the small oscillations; the result is plotted as a function of λ on Fig. 3: the δM ∝ λ 2 behavior is clear, for both families. The fact that numerical simulations are able to reach this stationary state suggests that it is a genuine stationary state, indeed stable with respect to the whole dynamics, and not only on the unstable manifold. However, longer simulations, or with smaller grid sizes (not shown), indicate it is also easily destabilized by numerical noise. We conclude that the O(λ 2 ) state is thus probably close to the instability µ with β = 40. Dotted lines correspond to unstable stationary solutions and solid lines to stable ones. We are interested in this article in the neighborhood of point a, where a branch of stable inhomogeneous Mµ = 0 stationary solutions becomes unstable. At point l, a branch of stable homogeneous stationary solutions (Mµ = 0, red solid line) becomes unstable; we do not consider this bifurcation in this article.

threshold. Notice that for the initial unstable reference state |Λ µ (0)| = O(λ 2 ) is very small; hence the very close nearby stationary state with δM ∝ λ 2 may have a different stability, although (in)stability of non homogeneous stationary states is rather robust comparing with the homogeneous case (see [START_REF] Barré | On the neighborhood of an inhomogeneous stable stationary solution of the Vlasov equation -Case of an attractive cosine potential[END_REF] for a discussion).

We provide as Supplementary Material 3 videos, in order to better illustrate the phase space dynamics. The video "Fermi eps +.mp4" shows the time evolution of the distribution function in phase space with initial condition F 0 µ (E(q, p)) + ǫT (q, p), for ǫ = +1.8 • 10 -6 , β = 40, (M µ = 0.328, µ = 0.658) (it corresponds to the dashed blue curve in the upper inset of Fig. 3) from t = 0 to t = 600. Since the system reaches a new stationary state close to the original one, we observe almost no change in the distribution function. It is important to notice however that this picture is very different from that of the saturation of an instability over an homogeneous background: in that case, resonances would create small "cat's eyes" structures, which do not appear here. To better appreciate the dynamics in this case, we also provide the video "Fermi eps + diff.mp4", which is the same as the previous one, except that the reference state F 0 µ has been substracted; hence the evolution of the perturbation is more clearly shown.

The video "Fermi eps -.mp4" shows the time evolution of the distribution function with the same parameter values except that ǫ = -1.8 • 10 -6 (it corresponds to the green curve in the upper inset of Fig. 3). This time the distribution changes completely its shape and seems to approach a periodic solution, far away from the original stationary distribution. In all simulations the relative error between the total energy of the system at a given . Note that the reference stationary states F 0 µ follow the curve from point a towards point l on Fig. 2. Along this curve, λ starts from 0 at point a, then grows and reaches an upper limit, about 0.15, before decreasing and reaching 0 again at point l. This is why two values of µ may correspond to the same value of λ (but different values of δM ), as can be seen around λ = 0.15.

time and the total initial energy is at most of the order of 10 -7 .

V. CONCLUSION

To summarize, we have elucidated for one-dimensional Vlasov equations the generic dynamics of a non oscillating instability of a non homogeneous steady state. We have shown in particular: i) due to the absence of resonances, the physical picture is very different from the homogeneous case; ii) a striking asymmetry: in one direction on the unstable manifold, the instability quickly saturates, and the system reaches a nearby stationary state, in the other direction the instability does not saturate, and the dynamics leaves the perturbative regime; iii) for an instability rate λ, the nearby stationary state is at distance O(λ 2 ) from the original stationary state. Direct numerical simulations of the Vlasov equation with a cosine potential have confirmed these findings, and furthermore suggested that this phenomenology may have relevance for initial conditions that are not on the unstable manifold.

These results are generic for 1D Vlasov equations. It will be of course important to assess if they can be extended to two-and three-dimensional systems. The structure of the reduced equation on the unstable manifold, which underlies the asymmetric behavior, should be valid in any dimension. However, the crucial physical ingredient underlying our detailed computations leading to the O(λ 2 ) scaling is the absence (or the weak effect) of resonance between the growing mode and the particles' dynamics. This has to be carefully analyzed in higher dimensions, and is left for a future study. Let us note that the findings of [START_REF] Palmer | Neighbouring equilibria to radially anisotropic spheres: possible end-states for violently relaxed stellar systems[END_REF] regarding the nonlinear evolution of the Radial Orbit Instability are in agreement with this picture: they start from a weakly unstable spherical state; close to this reference state, they find an axisymmetric weakly oblate stationary state; at finite distance from the reference state (outside the perturbative regime), they find a stable prolate stationary state. This is consistent with the asymmetric behavior we predict on the unstable manifold. An open question regards the stability of the new weakly oblate stationary state: in our numerical examples, it seems stable, but very close to threshold. In [START_REF] Palmer | Neighbouring equilibria to radially anisotropic spheres: possible end-states for violently relaxed stellar systems[END_REF] they claim it is unstable, but since their numerical simulations are much less precise, this could be a numerical artefact.

Finally, the possibility of an infinite dimensional reduced dynamics, as in [START_REF] Del-Castillo-Negrete | Nonlinear evolution of perturbations in marginally stable plasmas[END_REF], remains open. We have also restricted ourselves in this article to a real unstable eigenvalue. It is obviously worthwhile to investigate the bifurcation with a complex pair of eigenvalues. Resonances should appear in this case, similar to those in the homogeneous case, but with a different symmetry: new features are thus expected. function) are known [START_REF] Milne | Infinite families of exact sums of squares formulas, Jacobi elliptic functions, continued fractions, and Schur functions[END_REF], and can be used with [START_REF] Ogawa | Spectral and formal stability criteria of spatially inhomogeneous stationary solutions to the Vlasov equation for the Hamiltonian mean-field model[END_REF] to obtain the following expressions:

c 0 (k) =        2E(k) K(k) -1, k < 1 2k 2 E(1/k) K(1/k) + 1 -2k 2 , k > 1 (57)
and for m = 0

           c 2m (k) = 2π 2 K(k) 2 mq(k) m 1 -q(k) 2m , k < 1 c 2m+1 (k) = 0, k < 1 c m (k) = 2π 2 k 2 K(1/k) 2 mq(1/k) m 1 -q(1/k) 2m , k > 1 (58)
where K(k) and E(k) are the elliptic functions of first and second kind respectively and q(k) = exp -πK( √ 1k 2 )/K(k) . Note that one may also use the sn × dn and sn × cn Fourier expansions to obtain similar expressions for the Fourier transform of sin q with respect to θ. These explicit expressions allow efficient computations of the c m and help in choosing a right truncation for the m summation. A standard numerical solver then provides precise values of λ.

FIG. 1 :

 1 FIG. 1: (Color online)Frequency Ω as a function of the action J in the HMF model. The action is Js at the separatrix (see the upper-left inset representing the phase space (q, p)), where Js is different if the separatrix is approached from below or from above[START_REF] Barré | Dynamics of perturbations around inhomogeneous backgrounds in the HMF model[END_REF]. Ω goes to Ω0 in the limit J → 0 + . The lower-right inset shows the logarithmic divergence of 1/Ω around the separatrix action Js for both inside (red solid upper line) and outside (blue broken lower line).

FIG. 2 :

 2 FIG. 2: (Color online) Phase diagram of the Fermi distribution F 0µ with β = 40. Dotted lines correspond to unstable stationary solutions and solid lines to stable ones. We are interested in this article in the neighborhood of point a, where a branch of stable inhomogeneous Mµ = 0 stationary solutions becomes unstable. At point l, a branch of stable homogeneous stationary solutions (Mµ = 0, red solid line) becomes unstable; we do not consider this bifurcation in this article.

FIG. 3 :

 3 FIG. 3: (Color online) δM (λ) for F 0 µ (circles) and G 0 µ (crosses) with associated quadratic fit; β = 40. For each function we show two runs of Mx(t) = M (t) (here My(t) = 0) with positive and negative ǫ. δM in the main diagram is computed as a long time average for ǫ > 0. For F 0 µ , |ǫ| = 1.8•10 -6 ; for G 0 µ , |ǫ| = 10 -5 , except for λ = 0.032 where |ǫ| = 3 • 10 -6. Note that the reference stationary states F 0 µ follow the curve from point a towards point l on Fig.2. Along this curve, λ starts from 0 at point a, then grows and reaches an upper limit, about 0.15, before decreasing and reaching 0 again at point l. This is why two values of µ may correspond to the same value of λ (but different values of δM ), as can be seen around λ = 0.15.

Acknowledgments D.M. thanks T. Rocha for providing his GPU Vlasov-HMF solver. Y.Y.Y. acknowledges the support of JSPS KAKENHI Grant Number 23560069. J.B. acknowledges the hospitality of Imperial College, where part of this work was performed.

where

Taking z = 2λ in (45), we see that Λ µ (2λ) -1 introduces a 1/λ 2 divergence, and that the l = 0 term in the sum yields an extra 1/λ divergence, unless Ĝl=0 = 0. Thus, the C factor (46) gives the leading singularity. Now, recalling (43), we have to apply the resolvent to G = Π ⊥ B(Ψ λ , Ψ λ ). Using the definition of a 2 , we have

We first note that B(Ψ λ , Ψ λ ) m=0 = 0 (B contains two terms, each containing a derivative with respect to θ; hence the zeroth Fourier mode vanishes) and that Ψλ,m=0 = 0. Hence the possible divergence related to l = 0 in (45) does not exist, and the resolvent introduces only a 1/λ 2 divergence. We conclude that in the r.h.s. of (47), the first term is O(λ -2 ). The second is a priori O(λ -3 ), since a 2 ∝ 1/λ, but we show now it is actually also O(λ -2 ). The C factor (46) for (R(2λ

where the last line is for λ → 0 + . Hence, the second term in the r.h.s. of (47) is O(λ -2 ), and so is S 2 . From (42)

where the r.h.s. is Π applied to a O(λ -2 ). Now, the projection Π contains a diverging 1/λ factor, coming from the normalization factor 1/Λ ′ µ (λ), needed in Ψ λ to ensure that Ψ λ , Ψ λ = 1. Hence, except for a restricted set of functions ϕ such that Ψ λ , ϕ = O(1), we have (for ϕ independent of λ): Πϕ ∝ 1 λ . The exceptional ϕ such that the projection Π does not introduce a diverging 1/λ factor lie close to the kernels of Π and Π ⊥ (the latter is just Span(Ψ λ )). With this in mind, it is not difficult to conclude that a 3 ∝ 1/λ 3 . This is the same divergence strength as the one which appears in the homogeneous case [START_REF] Crawford | Universal Trapping Scaling on the Unstable Manifold for a Collisionless Electrostatic Mode[END_REF], although the mechanism inducing the divergence is different.

C. Rearrangement formula for the cosine potential

The rearrangement formula is a powerful tool to predict the asymptotic stationary state f A µ from a given ini-tial state f µ . We show here that it allows to recover some of the above results, following a very different route.

The main result from the rearrangement formula is, roughly speaking,

where the symbol • A represents the average over angle variable at fixed action, angle and action being associated with the asymptotic Hamiltonian H A µ = H[f A µ ] whose potential part is determined self-consistently [START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF]. We show that this formula predicts

, where M µ and M A µ are magnetizations in f µ and f A µ respectively: this is consistent with the unstable manifold analysis.

We assume that f µ = F µ (H µ ). The average • A can be removed when applied to any function of energy:

) by the definition [START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF]. Thus, we can expand f µ A with respect to small δH = H A µ -H µ , and have

) Adding and subtracting F ′ (H µ ) δH I where the average • I is taken over angle variable associated with H µ , multiplying by cos q and integrating over (q, p), we obtain a self-consistent equation for δM

where N (δM ) is of higher order and we used the Parseval's equality [START_REF] Ogawa | Spectral and formal stability criteria of spatially inhomogeneous stationary solutions to the Vlasov equation for the Hamiltonian mean-field model[END_REF] to derive the coefficient Λ µ (0). Straightforward computations give

where the homogeneous case is rather singular [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF], and the leading nonzero solution to (51) is

The relation between Λ µ (0) and the instability rate λ is obtained by recalling that λ satisfies Λ µ (λ) = 0 and expanding Λ µ (λ) with respect to λ. This expansion implies

and the former and the latter correspond to the homogeneous and non homogeneous cases respectively. Summarizing, the rearrangement formula gives the scaling δM = O(λ 2 ) for both homogeneous and non homogeneous cases, although the mechanisms for these scalings are different in each case, as already found in the unstable manifold reduction. We remark that, in the non homogeneous case, the leading three terms of (51) gives two nonzero solutions of order O(λ 2 ) and O(1), whose stability is not clear by the rearrangement formula only. The latter solution is a priori out of range of the method; nevertheless, later we will numerically observe that one direction of perturbation drives the system to the O(λ 2 ) state, and the other direction to O(1).

VI. APPENDIX: COMPUTING THE FUNCTION Λµ(λ) AND ITS ROOTS

Testing the predicted scaling A ∞ ∝ λ 2 , as done in Sec.IV, requires computing λ with a good accuracy. We present our method in this appendix. In order to compute the eigenvalue λ associated with a given initial distribution f µ (q, p) = F µ (J), one has to find a positive root of the dispersion function

The two main numerical obstacles are to compute efficiently the functions c m (J), and the infinite sum in m. Fourier expansions of sn 2 (sn is the sine Jacobi elliptic