Shun Ogawa
email: shun.ogawa@riken.jp

Stability criterion of spatially inhomogeneous solutions to Vlasov equation

Stability of spatially inhomogeneous stationary solutions to the Vlasov equation for multi-dimensional systems is investigated. The formal stability criterion taking into account all of the Casimir invariants is derived in an explicit form without use of the angle-action variables of a single particle's Hamiltonian associated with a stationary solutions to the Vlasov equation, when biorthogonal basis can be defined.

Introduction

Temporal evolution of systems with long range interaction (LRI) is approximately described by the Vlasov equation (also called the collisionless Boltzmann equation) which is an evolution equation of a single particle density function, when the number of particles are large enough [START_REF] Balescu | Statistical Dynamics Matter Out of Equilibrium[END_REF][START_REF] Campa | Physics of Long-Range Interacting Systems[END_REF][START_REF] Campa | [END_REF][START_REF] Ichimaru | Basic principles of plasmas physics a statistical approach[END_REF][START_REF] Binney | Galactic dynamics[END_REF]. This approximation gets to be exact when one takes a large population limit (mean-field limit) [START_REF] Braun | [END_REF][7][START_REF] Spohn | Large scale dynamics of interacting particles[END_REF]. This equation is hence a basic tool for analyzing mesoscopic dynamics of Hamiltonian systems with LRI, for example, plasmas [START_REF] Ichimaru | Basic principles of plasmas physics a statistical approach[END_REF] and self-gravitating systems [START_REF] Binney | Galactic dynamics[END_REF]. One of remarkable phenomena observed in the LRI systems is existence of long-lasting out-of-equilibrium stationary states called quasistationary states (QSSs) and slowly equilibration in the relaxation process. The QSSs are regarded as stable stationary solutions to the Vlasov equation [START_REF] Yamaguchi | [END_REF][10][11][12] . The LRI system is often trapped in a QSS and its duration diverges as the number N of elements consisting the system. After that, the finite N effect appears and the temporal evolution is no longer described by the Vlasov equation and is described by the kinetic equation such as Balescu-Lenard equation [START_REF] Balescu | Statistical Dynamics Matter Out of Equilibrium[END_REF][START_REF] Campa | Physics of Long-Range Interacting Systems[END_REF][START_REF] Campa | [END_REF][START_REF] Ichimaru | Basic principles of plasmas physics a statistical approach[END_REF]. Then, the system leaves the QSS and relaxes to the thermal equilibrium state. Such a relaxation process is directly observed in the N body simulations for the Hamiltonian mean field (HMF) model (globally coupled XY model with inertial terms) [13] in Refs. [START_REF] Yamaguchi | [END_REF]14,15] and for the α-HMF model with α ≤ 1 (the coupling constant decays as r -α where r denotes distance between lattice points) [16] in Refs. [17,18]. The duration τ QSS of QSS diverges as, for instance, τ QSS ∼ N ν with ν > 0 for the HMF model and the α-HMF model with α < 1 [START_REF] Yamaguchi | [END_REF]10,14,15,18], τ QSS ∼ ln N for the α-HMF model with α = 1 [18], and τ QSS ∼ N / ln N for self-gravitating systems [START_REF] Binney | Galactic dynamics[END_REF].

The stability analysis for the Vlasov equation is considered to be a first step of study on QSSs because it clarifies which kind of solutions can be QSSs. For the spatially homogeneous stationary solutions, the simple criteria of formal and spectral stability have been already derived in Refs. [START_REF] Yamaguchi | [END_REF][19][START_REF] Inagaki | gravitating systems[END_REF][21][22]. It is shown rigorously that some stable homogeneous states are valid for the finite N system within a time scale τ st ∼ N 1/8 [12]. The asymptotic stability in the nonlinear regime is rigorously shown for the Vlasov-Poisson system in the study on nonlinear Landau damping [23].

The aim of this article is to derive a precise formal stability criterion for spatially inhomogeneous stationary solutions to the Vlasov equation in an explicit form. The stability of inhomogeneous states has been investigated in the context of astrophysics for a long time [START_REF] Binney | Galactic dynamics[END_REF][24][25][26][27], and it has been investigated for simple systems with long-range interaction [28,29]. The nonlinear (orbital) stability is rigorously investigated recently for the spherical gravitating systems [30] and for the HMF model [31].

The spectral stability has been investigated by use of the Nyquist methods [19][START_REF] Inagaki | gravitating systems[END_REF][21][22]. In these works the back ground stationary state is assumed to be spatially homogeneous, and an explicit form of the plasmas dispersion function is required. Thus it is difficult to look into stability of spatially inhomogeneous states by using the Nyquist method. There exists another way to derive a stability criterion for spatially inhomogeneous states of the HMF model [28], and the explicit form of dispersion function is necessary in this procedure too, so that it cannot be applied for stationary states whose effective Hamiltonian is non-integrable.

In this article, the formal stability [START_REF] Yamaguchi | [END_REF]28,29,32,[START_REF] Arnold | Mathematical methods of classical mechanics[END_REF] is investigated. Let f 0 be a stationary solution to the Vlasov equation and also be a critical point of some invariant functional (The detail is exhibited in Sec. 3.). A solution f 0 is said to be formally stable when a second variation ("Hessian") of the invariant functional around f 0 is positive or negative definite. In this case, L 2 -like-norm of perturbation given by this "Hessian" is an invariant of the linearized Vlasov equation, so that the formally stable solution is linearly stable.

The Vlasov equation has infinitely many constraint conditions coming from the Casimir invariant functionals. Then a physically relevant perturbation satisfies these constraint conditions within a linear order. The problem of dealing these constraints appeared in a stability analysis of perfect fluid flows [START_REF] Arnold | Mathematical methods of classical mechanics[END_REF][START_REF] Ellis | [END_REF]. More recently, this is a main problem to derive the asymptotic form of the linear response of the Vlasov systems around spatially inhomogeneous states [14,35,36]. The QSSs are looked as stable fixed points on a constraint "surface" and these constraints prevent systems relaxing to the equilibrium state [14], and the solution to the linearized Vlasov equation around f 0 should be on the tangent "plane" at f 0 in a constraint "surface." When f 0 is spatially homogeneous [START_REF] Yamaguchi | [END_REF] or is inhomogeneous whose effective single particle Hamiltonian is integrable [28], the formal stability criterion taking into account all the Casimir constraints is derived. Stationary states with nonintegrable effective Hamiltonians are out of scope in these previous studies. However, for a system in two or three dimensional spaces, the effective Hamiltonians are nonintegrable in general. In Ref. [29], the linear stability is investigated without use of the angle-action variables, and less-refined formal stability criteria taking into account a finite number of invariants have been derived. But the formal stability criterion taking into account all of the Casimir constraints has not been derived, because it requires to determine an infinite number of Lagrangian multipliers associated with the Casimir constraints [29], and it is impossible in general. Some stationary states are stable against the perturbation violating some Casimir constraints. For instance, if an only normalization condition is taken into account, the inhomogeneous thermal equilibria are judged to be stable in the HMF model [29]. On the other hand, some perturbation violating the constraint conditions can bring about instability. For instance, if one does not take into account any constraint condition, all equilibria in the ordered phase are not judged to be stable in the HMF model [29]. Further, as shown in Ref. [28], some inhomogeneous solutions are not judged to be stable when one use the stability criterion taking into account only normalization of density functions, even though this is formally stable. Thus, it is important in practice to derive the stability criterion taking into account infinitely many Casimir invariants without use of the angle-action variables.

We then derive the formal stability criterion without use of the angle-action variables by expressing the constraint conditions in a simpler manner as in Ref. [36] (a similar formula is found in Ref. [25,30]). The formal stability criterion is explicitly written in terms of the positive or negative definiteness of a matrix in a similar manner with the previous result [28]. Everything is given a priori in this matrix, so that we can judge the stability of f 0 by substituting it into the criterion, in principle. When this formula is applied for the one-dimensional (1d) systems, this matrix coincides with the dispersion matrix [START_REF] Binney | Galactic dynamics[END_REF]37] with null frequency. For the multi-dimensional systems, the inverse matrix appears in the linear response and the zero-field susceptibility diverges when the matrix has a null eigenvalue [36]. The marginally formally stable solutions are the marginally spectrally stable ones not only in the 1d systems but also in the multidimensional systems.

The rest part of this article is organized as follows: We introduce the Vlasov equation and the constraint conditions of Vlasov dynamics in Sec. 2. In Sec. 3, we exhibit a definition of formal stability and derive its criterion. We check the present result includes the previous one for the 1d HMF model, and derive the stability criterion for the 2d HMF model. The article is summarized in Sec. 4.

Vlasov equation 2.1 Model

A model dealt in this article is an N -body system with Hamiltonian, in the d dimensional real space,

H N = N ∑ i =1 ∥p i ∥ 2 2 + 1 2N N ∑ i , j =1 V ( q i -q j ) , (1)
where

q i ∈ D d ⊂ R d and its conjugate momentum p i ∈ R d , for i = 1, 2, • • • , N .
We assume that the domain D d is bounded and the potential V is smooth. In the limit N → ∞, the temporal evolution is well described by the Vlasov equation [START_REF] Braun | [END_REF][7][START_REF] Spohn | Large scale dynamics of interacting particles[END_REF],

∂ f ∂t + { H [ f ], f } = 0, H [ f ] = ∥p∥ 2 2 + V [ f ](q), V [ f ] = V (q -q ′ ) f (q ′ , p ′ , t )dq ′ dp ′ , (2)
where f (q, p, t ) is a single particle density function, and the Poisson bracket is given by {a, b}

= ∂a ∂p • ∂b ∂q - ∂a ∂q • ∂b ∂p .
The Vlasov equation has infinitely many stationary solutions, and we here deal mainly with a stationary solution f 0 depending only on the one particle effective Hamiltonian H 0 ,

f 0 (q, p) = F 0 ( H 0 (q, p) ) , H 0 ( q, p ) = H [ f 0 ], (3)
where F 0 is a monotonically decreasing function. In the similar manner, we can deal with a stationary solution f 0 depending on two integrals as follows

f 0 (q, p) = F 0 ( H 0 (q, p) + νL(q, p) ) , H 0 ( q, p ) = H [ f 0 ], (4)
where L is an additional invariant of the Hamiltonian system with H 0 , and where F 0 is monotonic with respect to H 0 + νL where ν is a constant. The stationary solution f 0 is assumed to depend on a linear combination of two integrals because we focus on the stationary solution f 0 which is also a solution to the optimization problem with some constraint conditions in the present article. We additionally assume that biorthogonal basis can be defined in this model, and this is introduced in a forthcoming section.

Casimir invariants

The Vlasov equation has infinitely many constraint conditions coming from Casimir functionals,

C [ f ] = c ( f (q, p, t ) ) dqdp, (5)
which are conserved in the temporal evolution for any smooth function c. Then solutions to the linearized Vlasov equation around f 0 should be on the tangent "plane" at f 0 on the iso-Casimirs "surface," so that the physically relevant perturbation δ f satisfies, within a linear order,

c ′ ( f 0 (q, p) ) δ f (q, p)dqdp = 0, (6)
for any c, where c ′ (x) = dc/dx. Since f 0 = F 0 (H 0 ) is assumed to be monotonic with respect to H 0 , it is shown that this constraint condition is equivalent to the condition:

R ( H 0 (q, p) ) δ f (q, p)dqdp = 0 (7)
for any function R [28]. This is equivalent to the condition:

R ( H 0 (q, p) ) 〈δ f (q, p)〉 H 0 dqdp = 0 (8)
(see A for the derivation). The constraint condition are expressed as follows:

〈δ f 〉 H 0 = 0, (9)
for almost every H 0 , where 〈•〉 H 0 denotes the average over the iso-H 0 set. Thus, the accessible perturbation can be written as

δ f (q, p) = g (q, p) -〈g 〉 H 0 (q ,p ) . ( 10
)
It should be remarked that the condition ( 6)-( 9) is not for the initial state of the perturbation δ f (t = 0) but for the relevant dynamical part of the perturbation whose constrained on a constraint surface. Even if the δ f (t = 0) does not satisfies the constraint condition ( 9), the linearized dynamics should be constrained on a tangent plane of the constraint surface and the directions that are not parallel to the tangent plane are not accessible. Thus, it is worthwhile to check the formal stability to show the linear stability of the stationary solution f 0 .

When the system has an additional invariant L, we should find a condition likes Eq. ( 9) for each steady state. If the stationary solution f 0 depends on only H 0 , the physically relevant perturbation satisfies Eq. ( 9). When f 0 is written as Eq. ( 4), the accessible perturbation satisfies

〈δ f 〉 H 0 +νL = 0. ( 11
)
Then the stability criterion for f 0 in Eq. ( 4) is derived as done for the stationary solution depending only H 0 by replacing H 0 with H 0 + νL.

3 Formal stability

Definition of formal stability

The concept of the formal stability [32] is introduced here. Let us consider the stationary solution f 0 to the Vlasov equation, which is a solution to an optimization problem with constraints, normalization and energy conservation;

maximizing S [ f ] = s( f )dqdp, subject to N [ f ] = f dqdp = 1, E [ f ] = ∥p∥ 2 2 f dqdp + 1 2 V [ f ] f dqdp = E , ( 12
)
where s is a strictly concave function. This is similar to the maximizing entropy principle. A solution to the optimization problem (12) is

f 0 = (s ′ ) -1 (βH 0 + α), ( 13
)
where s ′ (x) = ds/dx, and α and β are Lagrangian multipliers corresponding to the normalization and the energy conservation respectively. Since s is concave, f 0 satisfies the condition (3) when β > 0. Taking the second variation of the functional F = S -αN -βE around f 0 , we have

δ 2 F [ f 0 ][δ f , δ f * ] = β ( |δ f | 2 γ(H 0 ) -δ f * V [δ f ] ) dqdp, ( 14
)
where z * denotes complex conjugate of z ∈ C and γ(H 0 ) = β/s ′′ ( f 0 ) = d f 0 /dH 0 . The functional ( 14) is an invariant of the linearized Vlasov dynamics. When the quadratic form ( 14) is negative for any perturbations δ f satisfying Eq. ( 9), the stationary solution f 0 is said to be formally stable. The formally stable solution is linearly stable [32].

Let us consider the case that there is an additional invariant,

L [ f ] = L(q, p) f dqdp = L. ( 15
)
The solution to the optimization problem is f 0 = (s ′ ) -1 (βΨ + α), where Ψ = H 0 + νL/β and ν is a Lagrangian multiplier associated with the condition (15). The second variation of the functional ( 15) is zero, so that it does not appear in Eq. ( 14). Then, we just put γ(H 0 , L) = d f 0 /dΨ in Eq. ( 14) instead of γ(H 0 ).

Bi-orthogonal basis

A useful tool for analyzing perturbations of the Vlasov systems are biorthogonal basis, {d i (q)} i ∈I and {u i (q)} i ∈I ′ , where I ′ ⊂ I ⊂ Z, and we assume a number of elements of I ′ , ♯I ′ < ∞ [START_REF] Binney | Galactic dynamics[END_REF]37,42,43]. Any spatial density perturbations δρ(q) = ∫ δ f (q, p)dp is spanned with {d i (q)} i ∈I , and {u i (q)} i ∈I ′ are given by

u i (q) = (V * d i )(q) ≡ ∫ V (q -q ′ )d j (q ′ )dq ′ . ( 16
)
They are supposed to satisfy

∫ u * i (q)d j (q)dq = λ i δ i j , (17)
where λ i ̸ = 0 for i ∈ I ′ and λ i = 0 for i ∈ I \ I ′ . By use of the biorthogonal basis, we can express the second term of right hand side in Eq. ( 14) as

δ f * V [δ f ]dqdp = ∑ i ∈I ′ 1 λ i u * i δ f dqdp 2 . ( 18
)
The derivation is exhibited in B.

Formal stability criterion

Let us derive the formal stability criterion explicitly. We suppose that the single particle Hamiltonian H 0 has only one integral, H 0 per se. Substituting Eq. ( 10) into Eq. ( 14), we express the formal stability condition as

G [g , g * ] = δ 2 F [ f 0 ] [ g -〈g 〉 H 0 , g * -〈g * 〉 H 0 ] /β < 0, (19)
for any g . Our strategy is to find a g max maximizing G under a constraint,

M [g , g * ] = ∑ i ∈I ′ u * i (g -〈g 〉 H 0 )dqdp 2 = 1, (20)
and to derive the condition so that G [g max , g * max ] < 0 which is equivalent to Eq. ( 19). The value 1 itself in the right hand side of Eq. ( 20) is not important. The condition M [g , g * ] > 0 is really important. This is because the rescaling g → c g with a constant c just changes the values of functionals as G → |c| 2 G and M → |c| 2 M and it does not affect the sign of G . When M [g , g * ] = 0, then G [g , g * ] < 0 because γ(H 0 ) < 0. Thus such a perturbation does not bring about instability. By using a vector ξ = (ξ i ) i ∈I ′ given by

ξ i [g ] = u * i ( g -〈g 〉 H 0 ) dqdp, (21)
the constraint ( 20) is expressed as ∥ξ∥ = 1, where ∥•∥ denotes Euclidean norm. We write ξ i [g max ] = ξi and ξ[g max ] = ξ. Let  † denote the Hermitian conjugate of a matrix  and a † = (a * ) T where the superscript T denotes transposition. By use of the equality (see A for the derivation),

a〈b〉 H 0 dqdp = 〈a〉 H 0 bdqdp, (22)
for any functions a and b for which the integrals are defined, we have g max holding δG /δg -ηδM /δg = 0 as

g * max -〈g * max 〉 H 0 = γ(H 0 ) ∑ i ∈I ′ ( 1 λ i + η ) ξ * i ( u * i -〈u * i 〉 H 0 ) , ( 23
)
where η is a Lagrangian multiplier with respect to the constraint condition [START_REF] Inagaki | gravitating systems[END_REF]. Let 23) into Eq. ( 19), we have the condition of the formal stability of f 0 as

ζ i = (λ -1 i + η) ξi and ζ = (ζ i ) i ∈I ′ . Substituting Eq. (
∑ i , j ∈I ′ ζ * i ζ j γ(H 0 ) ( u * i u j -〈u * i 〉 H 0 〈u j 〉 H 0 ) dqdp - ∑ i ∈I ′ λ -1 i | ξi | 2 < 0. ( 24
)
This inequality is similar to Eq. (1.26) in Ref. [30]. For simplicity, we introduce several notations;

F i j = γ(H 0 ) ( u * i u j -〈u * i 〉 H 0 〈u j 〉 H 0 ) dqdp, (25)
F = ( F i j ) (i , j )∈I ′ ×I ′ , and Λ = diag(λ 1 , λ 2 , • • • , λ ♯I ′ )
. By use of them, a relation ξ = F ζ is derived from Eq. ( 23). We thus have the formal stability condition as

ζ † F ( 1 -Λ-1 F ) ζ < 0, where 1 ∈ C ♯I ′ ×♯I ′
is the unit matrix and we have used the fact F = F † . The matrix F is negative definite, because

F i j = γ(H 0 ) ( 〈u * i u j 〉 H 0 -〈u * i 〉 H 0 〈u j 〉 H 0 ) dqdp, (26)
the matrix

( 〈u * i u j 〉 H 0 -〈u * i 〉 H 0 〈u j 〉 H 0 ) (i , j )∈I ′ ×I ′ is
[ f 0 ] = 1 -Λ-1 F is positive definite. ■
It should be remarked that the matrix D[ f 0 ] is the matrix form of dispersion functions with null frequency when the spatial dimension d = 1 [START_REF] Binney | Galactic dynamics[END_REF]37,42]. Further, the inverse matrix of D[ f 0 ] appears in a linear response formula for the 1d and multi-dimensional systems [36,38]. Thus, the zero minimum eigenvalue of D[ f 0 ] might also correspond to the marginal spectral stability, because the zero field susceptibility diverges when det D[ f 0 ] = 0.

We mention the relation between the present main claim and the classical results on the spherical stellar system. The Antonov's variational principle [START_REF] Binney | Galactic dynamics[END_REF] says: For the stellar system in the state

f 0 = F 0 (H 0 ) (F ′ 0 < 0) is stable if and only if W A [ f 0 ][h] ≡ -δ 2 F [ {h, f 0 }, {h, f 0 } ] β > 0, ( 27
)
for any h(q, p). Furthermore, it is shown that F 0 (H 0 ) with

F ′ 0 < 0 is stable by showing W A [ f 0 ][h] > 0 [5]
. We call this statement the classical statement. Let us recall that f 0 = F 0 (H 0 ) (F ′ 0 < 0) is formally stable if and only if (Eq. ( 19))

-G [g , g * ] = -δ 2 F [ g -〈g 〉 H 0 , g * -〈g * 〉 H 0 ] /β > 0.
This is similar to Eq. ( 27). Let W and G be sets of δ f such that

∃h(q, p), s.t. δ f = {h, f 0 }, ( 28
)
∃g (q, p), s.t. δ f = g -〈g 〉 H 0 , (29)
respectively. Since 〈{h, f 0 }〉 H 0 = 0, we have W ⊂ G. It is hence shown that W A [ f 0 ][h] > 0 for any h and Antonov's variational principle holds true, when Eq. ( 19) holds, or equivalently, D[ f 0 ] is positive definite. If there exists h satisfying

F ′ 0 (H 0 ){h, H 0 } = g -〈g 〉 H 0 , (30)
for each g , we have W = G. Thus, the positive definiteness of D[ f 0 ] is equivalent to the condition ( 27) and we have: For the spherical stellar system, if F 0 (H 0 ) satisfies F ′ 0 (H 0 ) < 0 for all H 0 , then D[ f 0 ] is positive definite. On the other hand, F ′ 0 (H 0 ) < 0 is assumed in our main claim, so that the present criterion is not so powerful for the spherical stellar systems, when W = G.

The set I ′ is assumed to be a finite set in the present section. One problem is that the matrix D[ f 0 ] becomes infinite matrix, i. e. I ′ is not a finite set in general. Thus, some truncation or approximation of the model

V (q -q ′ ) = ∑ i ∈I ′ u i (q)u * i (q ′
) is necessary. When the biorthogonal basis are the Fourier modes, and V is smooth enough, it may be possible to truncate the matrix elements D[ f 0 ] i j for large i or j when one computes the determinant. We do not deal with this problem in this article. In the precedence subsections, two examples are exhibited, in which ♯I ′ < ∞, the 1d and two dimensional (2d) HMF models [13,39,40]. Both of them are out of scope of the classical statement. The 2d one is out of scope of the previous paper [28], because the angle-action variables are not constructed.

One-dimensional Hamiltonian mean-field model

Let us examine the obtained criterion includes the one obtained in the previous article [28], the formal stability criterion for the HMF model whose Hamiltonian is

H = N ∑ i =1 p 2 i 2 - 1 2N ∑ i ̸ = j cos(q i -q j ) -h N ∑ i =1 cos q i , (31)
where h is an external field,

p i ∈ R, q i ∈ [-π, π) for i = 1, 2, • • • , N .
The effective Hamiltonian H 0 of the stationary state f 0 = F 0 (H 0 ) is

H 0 = p 2 2 - cos(q -q ′ ) f 0 (q ′ , p ′ )dq ′ dp ′ -h cos q. ( 32
)
In this case, I = N, I ′ = {1, 2} and d i and u i are given respectively by d 2k = π -1 cos(kq), d 2k-1 = π -1 sin(kq), u 1 =sin q, and u 2 =cos q, and λ 1 = λ 2 = -1. Substituting them into the general form (26) of the matrix D = (D i j ) (i , j )∈{1,2} 2 , we have D 12 = D 21 = 0 [41], and

D 11 = 1 + d F 0 (H 0 ) d H 0 ( sin 2 q -〈sin q〉 2 H 0 (q,p)
) dqdp,

D 22 = 1 + d F 0 (H 0 ) d H 0 ( cos 2 q -〈cos q〉 2 H 0 (q,p) ) dqdp. ( 33
)
The stationary solution f 0 is formally stable if and only if D 11 > 0 and D 22 > 0. Using the angle-action variables, one shows that these coincide with the dispersion function (also called a dielectric function in the context of plasma physics) with null frequency. This criterion is hence the same with the one obtained in the previous study [28], and with the criterion of the orbital stability [31].

Two-dimensional Hamiltonian mean-field model

The 2d HMF model is out of scope of the previous study [28], because it is impossible to construct the angle-action variables for the effective Hamiltonian of the inhomogeneous stationary states. The N -body Hamiltonian of this model is

H = N ∑ i =1 ∥p i ∥ 2 2 -h x N ∑ i =1 cos x i -h y N ∑ i =1 cos y i - 1 2N ∑ i ̸ = j [ cos(x i -x j ) + cos(y i -y j ) + cos(x i -x j ) cos(y i -y j ) ] , (34)
where [39,40]. The effective Hamiltonian for the stationary state with the external field h x = h y = h is

p i = (p xi , p yi ) ∈ R 2 and q i = (x i , y i ) ∈ [-π, π) 2 for i = 1, 2, • • • , N
H 0 = p 2 x + p 2 y 2 + V (x, y), V (x, y) = -(M + h) cos x -(M + h) cos y -P cc cos x cos y, (35)
where we assume that the stationary state f 0 (q, p) = F 0 (H 0 ) is symmetric with respect to both x and y, and we define the order parameters in the stationary state f 0 as M = cos x f 0 dqdp = cos y f 0 dqdp and P cc = cos x cos y f 0 dqdp. We chose {u i } i ∈{1,2,••• ,8} as (

u 1 = cos x, u 2 =
) 36
For the symmetric potential V , we have 〈u i 〉 H 0 = 0, for i = 4, 5, • • • , 8, and F i j = 0 when u i u j is odd with respect to x or y. Let 0 m,n be an m × n null matrix. In the similar manner in Ref. [36], we explicitly write the matrix D = (D i j

) (i , j )∈{1,2,••• ,8} 2 as D =      De 0 3,4 0 3,1 0 4,3 Do 0 2,2 0 2,2 Do 0 4,1 0 1,3 0 1,4 D 88      = diag ( De , Do , Do , D 88 ) , (37)
where ) ,

D 88 = 1 + d F 0 (H 0 ) d H 0 sin 2 x sin 2 ydqdp, (38)
and where

G e 1 = d F 0 (H 0 ) d H 0 ( cos 2 x -〈cos x〉 2 H 0 ) dqdp, G e 2 = d F 0 (H 0 ) d H 0 ( cos x cos y -〈cos x〉 2 H 0 ) dqdp, G e 3 = d F 0 (H 0 ) d H 0 ( cos 2 x cos y -〈cos x〉 H 0 〈cos x cos y〉 H 0 ) dqdp, G e 4 = d F 0 (H 0 ) d H 0 ( cos 2 x cos 2 y -〈cos x cos y〉 2 H 0 ) dqdp, (39)
(see Appendix of Ref. [36] for a way to compute numerically them.) and

G o 1 = d F 0 (H 0 ) d H 0 sin 2 xdqdp, G o 2 = d F 0 (H 0 ) d H 0 sin 2 x cos ydqdp, G o 3 = d F 0 (H 0 ) d H 0 sin 2 x cos 2 ydqdp, (40)
Thus, for the 2d HMF model, F 0 is formally stable if and only if De and Do is positive definite and D 88 > 0.

The matrix De appears in the linear response formula [36] and if one of its eigenvalues is 0, the response diverges. Thus, the marginally formally stable solution f 0 is the marginally spectrally stable solution.

Summary

The formal stability criterion of stationary solutions to the Vlasov equation for systems in multi-dimensional space is derived in an explicit form, where the biorthogonal basis can be defined. This criterion takes into account all of the Casimir invariants.

The practical usefulness of the refined criterion has been already mentioned in the previous paper [28]. The less refined criteria taking into account on the finite numbers of invariants cannot be used when the stationary state is close to the marginal stable one with spatial inhomogeneity. Furthermore, we overcome a weak point of the result obtained in Ref. [28]; the criteria has been derived only for the HMF model or its generalized ones such as α-HMF model [16,17]. We have checked that our general result includes what is obtained in the previous study for the HMF model [28], and the matrix D[ f 0 ] coincides with the dispersion matrix when the system is in 1d real space.

For the 1d systems the minimum eigenvalue of D[ f 0 ] is 0 if and only if the operator L[ f 0 ] in the linearized Vlasov equation

∂ t δ f = L[ f 0 ][δ f ]
around f 0 has an (embedded) eigenvalue 0. Thus, the marginally formally stable solution f 0 is also the marginally spectrally stable solution. In the multi-dimensional systems, although the dispersion relation for the non-zero frequency has not been derived, we can obtain the same statement as follows: The matrix D[ f 0 ] -1 appears in the linear response formula [36]. Thus, the linear response diverges when D[ f 0 ] has a null eigenvalue, and the marginally formally stable f 0 is the marginally spectrally stable. It should be noted that the spectral stability criterion is not obtained and this is an open problem. It should be noted that when the matrix D[ f 0 ] is not positive definite, and has both positive and negative eigenvalues, it is possible that a growing (unstable) mode exists around f 0 .

A Derivation of Eq. (22) We exhibit a derivation of Eq. (22). 〈X 〉 H 0 is explicitly written as

〈X 〉 H 0 = 1 S(H 0 ) X (q, p)δ ( H 0 -H 0 (q, p) ) dqdp, (41)
where S denotes a measure of iso-H 0 in the µ-space, the phase space of single particle,

S(E ) = δ ( E -H 0 (q, p) ) dqdp. ( 42
)
We assume that ∥∇H 0 ∥ ̸ = 0, so that the iso-H 0 sets are null sets in 2d -dimensional µ space for each level. Hereafter, for any function X (q, p), X , X ′ , and X ′′ denote X (q, p), X (q ′ , p ′ ), and X (q ′′ , p ′′ ) respectively. The primes do not mean derivatives hereafter. The equation ( 22) is shown as follows: The right hand side in Eq. ( 22) is equal to

a [ b ′ δ ( H 0 -H ′ 0 ) S (H 0 ) dq ′ dp ′ ] dqdp = b ′ [ aδ ( H 0 -H ′ 0 ) S ( H ′ 0 ) dqdp ] dq ′ dp ′ , ( 43
)
and this is the same with the left hand side in Eq. (22). We have used the fact, for functions f ∈ R and G,

δ( f -f ′ )G( f ′ )dq ′ dp ′ = G( f ) δ( f -f ′ )dq ′ dp ′ , ( 44
)
where the range of f is in R and δ( ff ′ )dq ′ dp ′ < ∞.

B Derivation of Eq. (18) The Eq. ( 18) is shown as follows: The spatial density perturbation is

δρ(q) = ∫ δ f (q, p)dp = ∑ i ∈I a i d i (q). ( 45
)
We firstly show that

δ f * V [δ f ]dqdp = ∫ δρ * (V * δρ)dq = ∑ i ∈I ′ a i ∫ δρ * u i dq = ∑ i ∈I ′ λ i |a i | 2 . ( 46
)
It is also shown that

∑ i ∈I ′ 1 λ i ∫ u * i δ f dqdp 2 = ∑ i ∈I ′ ∫ u * i δρdq 2 = ∑ i ∈I ′ λ i |a i | 2 . ( 47
)
We therefore have Eq. ( 18).

  positive definite, and γ(H 0 ) < 0. It is therefore shown the main claim of this article.Main ClaimLet f 0 = F 0 (H 0 ) (F ′ 0 (H 0 ) < 0)be a stationary solution to the Vlasov equation and H 0 have only one integral of motion, H 0 itself. A stationary solution f 0 to the Vlasov equation is formally stable if and only if the matrix D

cos y, u 3 =

 3 cos x cos y, u 4 = sin x, u 5 = sin x cos y, u 6 = sin y, u 7 = cos x sin y, u 8 = sin x sin y.

  De ∈ R 3×3 , Do ∈ R 2×2 , and D 88 are respectively De =

Acknowledgments

The author acknowledges the financial support of the A * MIDEX project (n • ANR-11-IDEX-0001-02) funded by the "investissements d'Avenir" French Government program, managed by the French National Research Agency (ANR).


Shun Ogawa

Linear and nonlinear response of Vlasov system with non-integrable Hamiltonian

Linear and nonlinear response of the Vlasov system with nonintegrable Hamiltonian

I. INTRODUCTION

Long range interaction systems show several phenomena which are out of scope of the equilibrium statistical mechanics [START_REF] Campa | Physics of Long-Range Interacting Systems[END_REF][START_REF] Campa | Statistical mechanics and dynamics of solvable models with long-range interactions[END_REF]. One of them is that such a system is often trapped in out-of-equilibrium quasi-stationary states (QSSs) whose duration gets to be longer as the number N of elements in a system and diverges when the large population limit N → ∞ is taken [START_REF] Campa | Physics of Long-Range Interacting Systems[END_REF][START_REF] Campa | Statistical mechanics and dynamics of solvable models with long-range interactions[END_REF][START_REF] Binney | Galactic dynamics[END_REF][START_REF] Yamaguchi | Stability criteria of the Vlasov equation and quasi-stationary states of the HMF model[END_REF][START_REF] Barré | The Vlasov equation and the Hamiltonian mean-field model[END_REF]. Then, if the system of interest is huge enough, the relaxation time is so long that one cannot see the thermal equilibrium state. It is hence interesting to investigate the non-equilibrium statistical mechanics or thermodynamics of QSSs. In particular, the topic of this article is effect of external forces in the QSSs.

When N is huge enough, temporal evolution of the longrange interaction system is well described by the Vlasov equation [START_REF] Braun | The Vlasov dynamics and its fluctuations in the 1/N limit of interacting classical particles[END_REF][START_REF] Dobrushin | Vlasov equations[END_REF][START_REF] Spohn | Large scale dynamics of interacting particles[END_REF] (also called the collisionless Boltzmann equation [START_REF] Binney | Galactic dynamics[END_REF]) which describes evolution of a density function defined on a µ-space, a single particle phase space. The QSSs are interpreted as stable stationary solutions to the Vlasov equation [START_REF] Yamaguchi | Stability criteria of the Vlasov equation and quasi-stationary states of the HMF model[END_REF]. The Vlasov equation has unique solutions for each given initial state [START_REF] Dobrushin | Vlasov equations[END_REF][START_REF] Spohn | Large scale dynamics of interacting particles[END_REF], and thus the QSSs depend not only on macroscopic variables such as temperature and energy but also on mesoscopic things, details of the single particle density function. Thus the study on responses to external forces in QSSs should be based on the Vlasov equation.

The linear response theory for the Vlasov systems has been developed for stability analysis in self-gravitating systems [START_REF] Binney | Galactic dynamics[END_REF], for looking into plasmas responses in magnetically confined plasmas [START_REF] Boozer | Physics of magnetically confined plasmas[END_REF], for computing time-asymptotic response to the external forces of long-range interaction systems in both spatially homogeneous [START_REF] Patelli | Linear response theory for long-range interacting systems in quasistationary states[END_REF] and inhomogeneous [START_REF] Ogawa | Linear response theory in the Vlasov equation for homogeneous and for inhomogeneous quasistationary states[END_REF] QSSs and of a fluid systems [START_REF] Chavanis | Linear response theory for hydrodynamic and kinetic equations with long-range interactions[END_REF]. By use of this theory, critical phenomena in QSSs [START_REF] Ogawa | Non-mean-field critical exponent in a mean field model: Dynamics versus statistical mechanics[END_REF] are investigated, and some informations of unforced systems are extracted by observing responses to oscillating external forces [START_REF] Ogawa | Linear response theory in the Vlasov equation for homogeneous and for inhomogeneous quasistationary states[END_REF]. Further, the nonlinear response theory has been developed * shun.ogawa@riken.jp, Author's current address: RIKEN BSI, 2-1 Hirosawa, Wako, 351-0198 Japan to investigate the response to the finite size external forces in QSSs near or on the critical point in which the linear response theory does not work [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF][START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF]. These response formulae have been derived only when the single particle effective Hamiltonian is integrable and the angle-action variables are used for solving test particle dynamics.

Analyzing the linearized equation with the nonintegrable effective Hamiltonian is practically important, because systems in multi-dimensional spaces are more realistic (for example, self-gravitating systems in the three dimensional (3D) space, magnetically confined hot plasmas) and their effective Hamiltonians are non-integrable in general. To tackle this problem, one method to take into account constraints a posteriori is proposed to obtain linear response formulae approximately [START_REF] Patelli | General linear response formula for non integrable systems obeying the Vlasov equation[END_REF]. This method provides canonical (taking into account the normalization) and micro-canonical (taking into account the normalization and the energy concervation) linear responses and other kinds of linear responses with a finite number of constraints systematically. However one cannot obtain isolated one with this method because it is practically impossible to take into account infinitely many Casimirs constraints with this method. The same problem lies in the stability analysis of the Vlasov equation [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]. If the effective Hamiltonian is not integrable, it was impossible to obtain the precise stability criterion with a finite time step procedure in general, since one should obtain an infinite number of Lagrangian multipliers associated with the Casimir's constraints [START_REF] Campa | A dynamical stability criterion for inhomogeneous quasi-stationary states in long-range systems[END_REF].

The above problems in the linearized equation should be solved as the first step to understand the dynamics around the QSSs with non-integrable effective Hamiltonian. After that, we will be able to continue tackling more difficult problems on nonlinear Landau damping, nonlinear stability, nonlinear response, critical phenomena and their universality, and finite N effects.

In this article, we firstly obtain the linear response of the Vlasov system in the multi-dimensional space without solving the linearized Vlasov equation and without using the angle-action variables. Let an initial state without external field be f 0 (q, p) and a final state f h (q, p) after exerting external field h. The linear response is obtained by restricting the form of accessible perturbation by assuming smoothness of f h with respect to h, and by taking into account the constraint conditions that the perturbation should be on a tangent "plane" of a constraint surface at f 0 .

Further, we shall mention the nonlinear response formula [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF] derived via the transient (T-) linearization method developed by Lecellotti and Dorning [START_REF] Lancellotti | Critical Initial States in Collisionless Plasmas[END_REF][START_REF] Lancellotti | Time-asymptotic wave propagation in collisionless plasmas[END_REF][START_REF] Lancellotti | Nonlinear Landau damping[END_REF] to analyse the plasmas oscillation and nonlinear Landau damping. In this theory, the Vlasov equation is linearized around an "unknown" asymptotic stationary state. Solving this equation, we obtain a self-consistent equation determining the asymptotic state. The asymptotic solution is obtained by redistributing the initial density function along iso-asymptotic effective Hamiltonian sets, and this formula is called rearrangement formula [START_REF] Yamaguchi | Conditions for predicting quasistationary states by rearrangement formula[END_REF]. The same formula is also derived for predicting the QSSs in a 1D system [START_REF] Ribeiro-Teixeira | Ergodicity breaking and quasistationary states in systems with long-range interactions[END_REF] and a 3D self-gravitating system [START_REF] Benetti | Nonequilibrium Stationary States of 3D Self-Gravitating Systems[END_REF] via the another consideration when non-stationary initial states satisfy a (generalized) viral condition and there is no parametric resonance. It is derived the rearrangement formula keeps Casimir invariants at the order of T-linearization. Then, the nonlinear response formula is derived in a similar manner to derive the linear response formula in this article.

This article is organized as follows: The model and the dynamics in a mean-field limit N → ∞ are firstly introduced in Sec. II, and the explicit form of constraint condition coming from Casimir invariants is derived in Sec. III. Based on this constraint condition, the linear response formula is derived in Sec. IV and several examples are exhibited in Sec. V. We derive the nonlinear response formula and make a brief comment on problems in the T-linearized Vlasov equation for spatially multi-dimensional systems in Sec. VI, and summarize this article in Sec. VII.

II. MODEL AND ITS DYNAMICS

A. Model and Vlasov equation

Let us consider a system with long-range interaction whose Hamiltonian is

H N = N ∑ i =1 ∥p i ∥ 2 2 + 1 2N N ∑ i , j =1 V ( q i -q j ) + h(t ) N ∑ i =1 Φ(q i ), (1)
where q i denotes configuration of the i -th particle, p i its conjugate momentum, V the inter-particles (-cites) potential, and h(t )Φ(q i ) interaction between the external field h and the i -the particle. Taking the mean-field limit N → ∞ [START_REF] Braun | The Vlasov dynamics and its fluctuations in the 1/N limit of interacting classical particles[END_REF][START_REF] Dobrushin | Vlasov equations[END_REF][START_REF] Spohn | Large scale dynamics of interacting particles[END_REF], the temporal evolution of this system can be described in terms of the single-particle density function f (q, p, t ) which is a solution to the Vlasov equation,

∂ f ∂t + { H [ f ], f } = 0, (2)
where H [ f ] is an effective single-particle Hamiltonian,

H [ f ] = ∥p∥ 2 2 + V [ f ](q) + h(t )Φ(q), (3)
and {a, b} is the Poisson bracket given by {a, b} = ∂a ∂p

• ∂b ∂q - ∂a ∂q • ∂b ∂p . ( 4
)
The system is initially in a QSS, f 0 , and the effective Hamiltonian

H 0 (q, p) = H [ f 0 ](q, p) = ∥p∥ 2 2 + V [ f 0 ](q) (5)
has only one integral of motion, H 0 itself. The f 0 is expressed as f 0 (q, p) = F 0 ( H 0 (q, p) ) by use of a monotonically decreasing function F 0 . This assumption is reasonable when we are interested in the asymptotic behavior of perturbations around the formally stable solutions to the Vlasov equation [START_REF] Yamaguchi | Stability criteria of the Vlasov equation and quasi-stationary states of the HMF model[END_REF][START_REF] Barré | The Vlasov equation and the Hamiltonian mean-field model[END_REF][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][START_REF] Holm | Nonlinear stability of fluid and plasma equilibria[END_REF]. The formal stability is defined in terms of positive or negative definiteness of second variation of an invariant functional around f 0 which is a solution to the optimization problem:

maximizing S [ f ] = s( f )d qd p, subject to 1 = N [ f ] = f d qd p, E = E [ f ] = ∥p∥ 2 2 f d qd p + 1 2 V [ f ] f d qd p ( 6
)
where s is a convex function. A formally stable solution is linearly stable [START_REF] Holm | Nonlinear stability of fluid and plasma equilibria[END_REF]. By solving the optimization problem, we have a solution,

f 0 (q, p) = (s ′ ) -1 (βH 0 + α) ≡ F 0 (H 0 ), (7)
where s ′ denotes d s(x)/d x, and α and β are Lagrangian multipliers with respect to the normalization and the energy conservation. Since s is convex, then the inverse of its first derivative (s ′ ) -1 is a strictly decreasing function. The parameter β must be positive.

B. Linear response

The external field h(t ) is turned on and it converges to a constant h(t ) → h as t → ∞. In the previous studies [START_REF] Patelli | Linear response theory for long-range interacting systems in quasistationary states[END_REF][START_REF] Ogawa | Linear response theory in the Vlasov equation for homogeneous and for inhomogeneous quasistationary states[END_REF][START_REF] Chavanis | Linear response theory for hydrodynamic and kinetic equations with long-range interactions[END_REF] the asymptotic linear response δ f is obtained by solving the linearized Vlasov equation around f 0 ,

∂g p ∂t + {H 0 , g p } + {V [g p ] + h(t ))Φ, f 0 } = 0,
where g p (t ) ∼ O(h) is a perturbation around f 0 , and by taking the limit, δ f = lim t →∞ g p (t ). The angle-action variables of the Hamiltonian H 0 is necessary to solve the linearized Vlasov equation, but it is impossible in general for multidimensional systems. To avoid this problem, we focus on constraint conditions restricting a form of perturbations, and we obtain the linear response δ f without solving the linearized Vlasov equation.

III. CASIMIR INVARIANTS

We assume that f and its derivatives converge to 0 rapidly enough as ∥p∥ → ∞. Further, f and its derivatives are assumed to vanish on the boundary of spatial domain or the system has periodic boundary condition with respect to q. Under these assumptions, it is shown (see Appendix. A) that the Vlasov equation keeps values of Casimir functionals,

C [ f ] = c ( f (q, p, t ) ) d qd p, (8)
for any smooth function c. The linearized Casimir conservation condition is expressed as that the accessible perturbation δ f satisfies

c ′ ( f 0 (q, p) ) δ f (q, p)d qd p = 0, (9)
where c ′ (x) = d c/d x for any smooth function c. Since f 0 (q, p) = F 0 ( H 0 (q, p) ) , F 0 is a monotonically decreasing function, and c is chosen arbitrarily, the constraint condition ( 9) is equivalent to the condition,

0 = R ( H 0 (q, p) ) δ f (q, p)d qd p = R ( H 0 (q, p) ) 〈δ f 〉 H 0 (q ,p ) d qd p, (10)
for any function R [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]. The second equality is shown as follows:

R ( H 0 (q, p) ) 〈δ f 〉 H 0 (q ,p ) d qd p = R ( H 0 (q, p) ) δ f (q ′ , p ′ )δ ( H 0 (q, p) -H 0 (q ′ , p ′ ) ) d q ′ d p ′ δ ( H 0 (q, p) -H 0 (q ′′ , p ′′ ) ) d q ′′ d p ′′ d qd p = R ( H 0 (q, p) ) [ δ f (q ′ , p ′ )δ ( H 0 (q, p) -H 0 (q ′ , p ′ ) ) S ( H 0 (q, p) ) d q ′ d p ′ ] d qd p = δ f (q ′ , p ′ ) [ R ( H 0 (q, p) ) δ ( H 0 (q, p) -H 0 (q ′ , p ′ ) ) S ( H 0 (q, p) ) d qd p ] d q ′ d p ′ = δ f (q ′ , p ′ ) R ( H 0 (q ′ , p ′ ) ) S ( H 0 (q ′ , p ′ ) ) S ( H 0 (q ′ , p ′ ) ) d q ′ d p ′ = R ( H 0 (q, p) ) δ f (q, p)d qd p, ( 11
)
where S denotes a volume of iso-H 0 set,

S(E ) = δ ( E -H 0 (q ′ , p ′ ) ) d q ′ d p ′ . ( 12
)
More generally it is possible to show as in the 1D case [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF] that 〈a〉 H 0 (q ,p ) b(q, p)d qd p = a(q, p)〈b〉 H 0 (q ,p ) d qd p, [START_REF] Ogawa | Non-mean-field critical exponent in a mean field model: Dynamics versus statistical mechanics[END_REF] for the functions a and b when 〈a〉 H 0 b and a〈b〉 H 0 are integrable. Thus it has been shown that the condition ( 10) is equivalent to

〈δ f 〉 H 0 = 0, for almost every (q, p). ( 14
)

IV. LINEAR RESPONSE FORMULA

A. Implicit form of linear response

After that the external field is exerted and the limit t → ∞ is taken, the effective Hamiltonian becomes

H h (q, p) = H 0 (q, p) + δV (q) + hΦ(q) + O(h 2 ), ( 15
)
where the linear response δV ≡ V [δ f ]. Let the initial and final states be respectively as:

f 0 (q, p) = F 0 ( H 0 (q, p) ) = G 0 ( H 0 (q, p) ) 〈G 0 〉 µ , f h (q, p) = F h ( H h (q, p) ) = G h ( H h (q, p) ) 〈G h 〉 µ , (16)
where 〈a〉 µ = ad qd p and

G h (H h ) = G 0 (H h ) + hG 1 (H h ) + O(h 2 ) = G 0 (H 0 ) +G ′ 0 (H 0 )(H h -H 0 ) + hG 1 (H 0 ) + O(h 2 ). ( 17
)
Expanding f h around f 0 , we have

f h = f 0 + G ′ 0 (H 0 )(δV + hΦ) + G 1 (H 0 ) 〈G 0 (H 0 )〉 µ - ⟨ G ′ 0 (H 0 )(δV + hΦ) +G 1 (H 0 ) ⟩ µ G 0 (H 0 ) 〈G 0 (H 0 )〉 2 µ + O(h 2 ). ( 18
)
We then obtain the linear response,

δ f ≡ F ′ 0 (H 0 )(δV + hΦ) + G 1 (H 0 ) 〈G 0 (H 0 )〉 µ - ⟨ F ′ 0 (H 0 )(δV + hΦ) ⟩ µ F 0 (H 0 ) - 〈G 1 (H 0 )〉 µ F 0 (H 0 ) 〈G 0 (H 0 )〉 µ (19)
by taking the linear order. The function G 1 is determined so that 〈δ f 〉 H 0 = 0, that is,

G 1 〈G 0 〉 µ = 〈G 1 〉 µ 〈G 0 〉 µ F 0 -F ′ 0 (H 0 )〈δV + hΦ〉 H 0 + F 0 (H 0 ) ⟨ F ′ 0 (H 0 ) (δV + hΦ) ⟩ µ , ( 20
)
where

F ′ 0 (H 0 ) = G ′ 0 (H 0 )/〈G 0 (H 0 )〉 µ .
The response is therefore implicitly given by

δ f = F ′ 0 (H 0 ) ( δV (q) - ⟨ δV (q) ⟩ H 0 ) + hF ′ 0 (H 0 ) ( Φ(q) - ⟨ Φ(q) ⟩ H 0 ) . ( 21
)
Solving the implicit linear response formula ( 21) by using bi-orthogonal basis, we obtain explicitly the linear response taking into account the constraint conditions. We make a comment on the case that there exist two integrals and f 0 depends on the both of them. Let L = L(q, p) f d qd p be an additional integral (the angular momentum density for example). We consider the optimization problem (6) and we add the additional constraint L = Const. to Eq. ( 6). A solution f 0 depends on H 0 and L as

f 0 = F 0 ( βH 0 + νL )
where β and ν are Lagrangian multipliers. Thus the accessible perturbation satisfies 〈δ f (q, p)〉 (βH 0 (q ,p )+νL(q ,p )) = 0, [START_REF] Yamaguchi | Conditions for predicting quasistationary states by rearrangement formula[END_REF] where the bracket 〈•〉 (βH 0 (q ,p )+νL(q ,p )) means the average taken over iso-(βH 0 +νL) set. It should be noted that a form of constrained perturbation depends on how f 0 depends on H 0 and L. We should find ways to restrict the form of perturbations for each stationary state.

If H 0 has three independent integrals of motion, we can use angle-action variables and can solve the linearized Vlasov equation.

B. Explicit form of linear response

We introduce the bi-orthogonal basis [START_REF] Binney | Galactic dynamics[END_REF][START_REF] Barré | Algebraic damping in the one-dimensional Vlasov equation[END_REF][START_REF] Barré | On algebraic damping close to inhomogeneous Vlasov equilibria in multi-dimensional spaces[END_REF], {d i (q)} i ∈I and {u i (q)} i ∈I ′ , where the sets I ′ and I satisfy I ′ ⊂ I ⊂ Z. A perturbation of spatial density is spanned by the base {d i } i ,

δρ(q) = ∫ δ f (q, p)d p = ∑ i ∈I a i (t )d i (q). ( 23
)
The base {u i } i is introduced as

u i (q) ≡ (V * d i )(q) ≡ ∫ V (q -q ′ )d i (q ′ )d q ′ , (24)
and the orthogonal relation ∫ d i (q) ū j (q)d q = λ j δ i j [START_REF] Holm | Nonlinear stability of fluid and plasma equilibria[END_REF] holds, where λ i ̸ = 0 when i ∈ I ′ and it vanishes otherwise, and δ i j is the Kronecker delta. The upper bar denotes complex conjugate. Integrating the terms including δ f or δV in Eq. ( 21) with respect to p, we have

∫ δ f d p - ∫ F ′ 0 (H 0 ) ( δV (q) - ⟨ δV (q) ⟩ H 0 ) d p = ∑ i ∈I a i [ d i - ∫ F ′ 0 (H 0 ) ( u i (q) - ⟨ u i (q) ⟩ H 0 ) d p ] . (26)
Multiplying ū j to both sides and integrating them with respect to q, we have,

∑ i ∈I ′ a i [ λ j δ i j - ∫ F ′ 0 (H 0 ) ( u i ū j -〈u i 〉 H 0 ⟨ ū j ⟩ H 0 ) d pd q ] = ∑ i ∈I ′ [ λ j δ j i - ∫ F ′ 0 (H 0 ) ( ū j u i - ⟨ ū j ⟩ H 0 〈u i 〉 H 0 ) d pd q ] a i . ( 27
)
Let F = (F j i ) (i , j )∈I ′ ×I ′ be a matrix whose elements are given by

F j i = ∫ F ′ 0 (H 0 ) ( ū j u i - ⟨ ū j ⟩ H 0 〈u i 〉 H 0 ) d pd q. ( 28
)
We further assume that the term coupling with external force can be expanded as

Φ(q) = ∑ i ∈I ′ b i u i (q). ( 29
)
We then have

h ∫ F ′ 0 (H 0 ) ( Φ(q) - ⟨ Φ(q) ⟩ H 0 ) d p = h ∑ i ∈I ′ b i ∫ F ′ 0 (H 0 ) ( u i (q) - ⟨ u i (q) ⟩ H 0 ) d p, (30)
by integrating the term coming from external force in Eq. ( 21) with respect to p. Multiplying ū j to it and integrating with respect to q, we have, (as we have already done)

h ∑ i ∈I ′ F j i b i . ( 31
)
Combining Eqs. ( 27) and ( 31), we get the linear equation de-

termining {a i } i ∈I ′ , ∑ j ∈I ′ ( λ j δ i j -F i j ) a j = h ∑ j ∈I ′ F i j b j . ( 32
)
Introducing symbols x = (x i ) i ∈I ′ for x = a, b and Λ = diag(λ i ) i ∈I ′ , we can simplify the equation as follows:

(

1 -Λ -1 F ) a = hΛ -1 Fb, ( 33
)
and it is solved as

a = h(1 -F) -1 Fb = hD -1 (1 -D)b, (34)
where 1 denotes the unit matrix and D = 1 -Λ -1 F.

The maximal-eigenvalue of D is zero when f 0 might be marginally stable, and corresponds to the critical point.

When we apply Eq. ( 34) to 1D systems, This explicit response formula formally coincides with what is derived in Ref. [START_REF] Ogawa | Linear response theory in the Vlasov equation for homogeneous and for inhomogeneous quasistationary states[END_REF].

V. EXAMPLES: HAMILTONIAN MEAN-FIELD MODELS A. One dimensional case

Let us examine the proposed theory by use of the Hamiltonian mean-field (HMF) model whose Hamiltonian is

H = N ∑ i =1 p 2 i 2 - 1 2N ∑ i ̸ = j cos(q i -q j ) -h N ∑ i =1 cos q i . ( 35
)
where h is an external field,

p i ∈ R, q i ∈ [-π, π) for i = 1, 2, • • • , N .
In the equilibrium state, this model shows second order phase transition at the temperature T = 0.5, where the Boltzmann's constant k B = 1. By use of the linear and nonlinear response formulae, the critical phenomena in the equilibrium state and QSSs for the isolated HMF model are investigated and it is shown that the Casimir constraints bring about the non-classical critical exponents [START_REF] Ogawa | Non-mean-field critical exponent in a mean field model: Dynamics versus statistical mechanics[END_REF][START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF][START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF]. We here check that the present method yields the linear response formula same with the previously obtained one for the 1D HMF model. The effective Hamiltonian is

H = p 2 2 - cos(q -q ′ ) f (q ′ , p ′ , t )d q ′ d p ′ -h cos q. ( 36
)
Applying Eq. ( 21) for the HMF model, we can derive the linear response formula

δ f = (-δM -h)F ′ 0 (H 0 ) ( cos q -〈cos q〉 H 0 ) , (37)
where δM = cos qδ f d qd p. Multiplying cos q to both sides and integrating over the µ-space, we obtain the linear response as

δM = 1 -D D h, ( 38
)
where

D = 1 + F ′ (H 0 ) ( cos 2 q -〈cos q〉 2 H 0 ) d qd p. ( 39
)
This is equivalent to the linear response formula obtained in Ref. [START_REF] Ogawa | Non-mean-field critical exponent in a mean field model: Dynamics versus statistical mechanics[END_REF].

For more general 1D systems it is obvious that Eq. ( 21) is equivalent to the linear order of the nonlinear response formula derived in Ref. [START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF].

V(x,y)

x y FIG. 1. (Color online) The effective potential of the Hamiltonian (41). A minimum point is at (0, 0) and min V = -2M -P . Saddle points are at (±π, 0) and (0, ±π) and V = P . A maximum point are on (±π, ±π) and max V = 2M -P , where M > P > 0.

B. Two dimensional case

We next examine our theory in the 2D system whose Hamiltonian is

H = N ∑ i =1 ∥p i ∥ 2 2 -h x N ∑ i =1 cos x i -h y N ∑ i =1 cos y i - 1 2N ∑ i ̸ = j [ cos(x i -x j ) + cos(y i -y j ) + cos(x i -x j ) cos(y i -y j ) ] , (40)
where [START_REF] Antoni | Anomalous diffusion as a signature of a collapsing phase in two-dimensional self-gravitating systems[END_REF][START_REF] Torcini | Equilibrium and dynamical properties of two-dimensional N -body systems with long-range attractive interactions[END_REF]. Let us assume the initial QSS f 0 with h x = h y = 0 is even with respect to both x and y. Then, the effective Hamiltonian is

q i = (x i , y i ) ∈ [-π, π) 2 for i = 1, 2, • • • , N
H 0 = p 2 x + p 2 y 2 + V (x, y), V (x, y) = -M x cos x -M y cos y -P cc cos x cos y, (41)
where

M x = cos xρ(q)d q, M y = cos yρ(q)d q, P cc = cos x cos yρ(q)d q, ( 42
)
and where ρ(q) = f (q, p)d p. The effective potential is shown in Fig. 1.

To compute the linear response of the macroscopic observables M x , M y , and P cc to the external field h x = h y = h, it is necessary to compute 〈cos x〉 H 0 , 〈cos y〉 H 0 , and 〈cos x cos y〉 H 0 . We here set M x = M y = M and P cc = P .

For any smooth function g depending only on q, 〈g 〉 E is expressed as follows (see derivation for Appendix B.)

〈g 〉 E = ∫ R 2 d p ∫ [-π,π) 2 g (q)δ ( H (q, p) -E ) d q = 2π ∫ [-π,π) 2 g (q)Θ ( E -V (x, y) ) d q. ( 43
)
Since H 0 is even with respect to both x and y, we have The following notations are introduced for simplicity;

〈sin x〉 H 0 = 〈sin y〉 H 0 = 〈sin x cos y〉 H 0 = 〈sin x sin y〉 H 0 = 〈cos x sin y〉 H 0 = 0. (44) When H 0 > 2M -P = max V (x,
G 1 = - F ′ 0 (H 0 ) ( cos 2 x -〈cos x〉 2 H 0 ) d qd p = - F ′ 0 (H 0 ) ( cos 2 y -〈cos y〉 2 H 0 ) d qd p, G 2 = - F ′ 0 (H 0 ) ( cos x cos y -〈cos x〉 H 0 〈cos y〉 H 0 ) d qd p = - F ′ 0 (H 0 ) ( cos x cos y -〈cos x〉 2 H 0 ) d qd p, G 3 = - F ′ 0 (H 0 ) ( cos 2 x cos y -〈cos x cos y〉 H 0 〈cos x〉 H 0 ) d qd p = - F ′ 0 (H 0 ) ( cos x cos 2 y -〈cos x cos y〉 H 0 〈cos y〉 H 0 ) d qd p, G 4 = - F ′ 0 (H 0 ) ( cos 2 x cos 2 y -〈cos x cos y〉 2 H 0 ) d qd p. ( 47
)
By use of them and Eq. ( 21), we have

  1 -G 1 -G 2 -G 3 -G 2 1 -G 1 -G 3 -G 3 -G 3 1 -G 4     δM x δM y δP cc   =   h x G 1 + h y G 2 h x G 2 + h y G 1 (h x + h y )G 3   (48)
We therefore obtain the explicit linear response formula as follows:

δM x = χ 1 h x + χ 2 h y , δM y = χ 2 h x + χ 1 h y , δP cc = χ 3 (h x + h y ), ( 49
)
where explicit expressions of χ 1 , χ 2 , and χ 3 are

χ 1 = 1 det G (G 1 -G 2 1 -G 1 G 4 + G 2 2 + G 2 3 + G 2 1 G 4 -2G 1 G 2 3 + 2G 2 G 2 3 -G 2 2 G 4 ), ( 50
)
χ 2 = 1 det G ( G 2 +G 2 3 -G 2 G 4 ) , ( 51
)
χ 3 = 1 det G G 3 (1 -G 1 +G 2 ) = G 3 1 -G 1 -G 2 -G 4 +G 1 G 4 + G 2 G 4 -2G 2 3 , ( 52
)
respectively, where the determinant of G, the matrix in the left hand side of Eq. ( 48), is

det G = (1 -G 1 + G 2 ) × ( 1 -G 1 -G 2 -G 4 + G 1 G 4 +G 2 G 4 -2G 2 3 ) . ( 53
)
A way to compute terms including 〈cos x〉 H 0 , 〈cos y〉 H 0 , and

〈cos x cos y〉 H 0 is exhibited in Appendix D. When F 0 (H 0 ) is spatially homogeneous, that is, M = P = 0, we have G 2 = G 3 = 0 and G 1 and G 4 (1 -G 1 < 1 -G 4 when M = P = 0) do not vanish. Thus the susceptibilities are χ 1 = G 1 1 -G 1 , χ 2 = χ 3 = 0 (54)
in the disordered phase. We numerically confirm the linear response formula. The initial state is the Maxwell-Boltzmann type:

f MB (q, p) = exp(-H 0 /T ) 〈exp(-H 0 /T )〉 µ . ( 55
)
This system shows the first order phase transition [START_REF] Antoni | Anomalous diffusion as a signature of a collapsing phase in two-dimensional self-gravitating systems[END_REF] and there is no (meta-)stable homogeneous state with T < 0.5. The initial values of order parameters for T = 0.3 and 0.4 are exhibited in Table I. The external field h x = h y = h is exerted. Theoretically obtained susceptibilities are exhibited in Table I when the temperature T = 0.3 and 0.4, so that the initial equilibria are spatially inhomogeneous. We integrate an equation of motion derived from the Hamiltonian (40) by using a fourth order symplectic integrator [START_REF] Yoshida | Recent progress in the theory and application of symplectic integrators[END_REF], and compute the order parameters of N body systems, given respec-

M N xc (t , h) = 1 N N ∑ i =1 cos x i (t , h), M N yc (t , h) = 1 N N ∑ i =1 cos y i (t , h), M N xs (t , h) = 1 N N ∑ i =1 sin x i (t , h), M N ys (t , h) = 1 N N ∑ i =1
sin y i (t , h), (56)

P N cc (t , h) = 1 N N ∑ i =1 cos x i (t , h) cos y i (t , h), P N cs (t , h) = 1 N N ∑ i =1 cos x i (t , h) sin y i (t , h), P N sc (t , h) = 1 N N ∑ i =1 sin x i (t , h) cos y i (t , h), P N ss (t , h) = 1 N N ∑ i =1 sin x i (t , h) sin y i (t , h), ( 57
)
for null amplitude h = 0 and non-zero h. We compare the theoretically obtained linear response δM x , δM y , and δP cc with the numerically obtained responses given respectively by

δM N x (h) = M N x (h) -M N x (0), δM N y (h) = M N y (h) -M N y (0), δP N x y (h) = P N x y (h) -P N x y (0), (58)
where M N x , M N y , and P N x y are given by

M N x = √ M N xc 2 + M N xs 2 , M N y = √ M N yc 2 + M N ys 2 , P N x y = √ P N cc 2 + P N cs 2 + P N sc 2 + P N ss 2 , (59)
and where upper bars in Eq. (58) denote the time average

M N xc (h) = 1 τ ∫ t 0 +τ t 0 M N xc (t , h)d t . ( 60
)
For the 2D HMF model there is error between M N x (0) and M which is a solution to the self-consistent equation. We then focus on the difference M N x (h) -M N x (0) rather than M N x (h) -M . We set t 0 = 200, τ = 200, N = 4 × 10 6 , and the time step δt = 0.05. Figure 2 shows that these t 0 and τ are appropriate, and Fig. 3 shows that the numerically obtained results confirm the theory.

VI. NONLINEAR RESPONSE FORMULA

The nonlinear response formula [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF] which is called the rearrangement formula in Ref. [START_REF] Yamaguchi | Conditions for predicting quasistationary states by rearrangement formula[END_REF] keeps the Casimir invariants within an order of the T-linearization, the linearization around an asymptotic (A-) state f A (q, p) = lim t →∞ f (q, p, t ) assumed to be stationary. We derive the nonlinear response formula via the same strategy for deriving the linear response formula in the present article. We assume that the asymptotic effective Hamiltonian

H A = ∥p∥ 2 /2 + V A + hΦ, V A = V [ f A ], (61)
have only one integral of a single particle motion and f A depends only on H A . Further f A = F A (H A ) is assumed to be monotonically decreasing with respect to H A . Under these assumptions, expanding Eq. ( 8) around f A as done in Sec. III, the constraint condition coming from Casimir invariants within an order of T-linearization can be expressed as

R ( H A (q, p) ) ( f 0 (q, p) -f A (q, p) ) d qd p = 0, ( 62
)
for any smooth function R on R. By use of Eq. ( 13), it is shown that Eq. (62) holds true if and only if

f A = 〈 f 0 〉 H A (63)
for almost every (q, p) in the µ-space. Deriving a selfconsistent equation

V A (q) = V (q -q ′ )〈 f 0 〉 H A (q ′ ,p ′ ) d q ′ d p ′ (64)
from Eq. ( 63) and solving it, one can obtain H A and the nonlinear response δ f = f Af 0 . Is it possible to derive Eq. ( 63) for the multi-dimensional systems as done in Refs. [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF][START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF]? There is some difficulty to derive the same formula from the T-linearization method for the multi-dimensional systems. To see this, let us exhibit a sketch of this T-linearization method (see Refs. [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF][START_REF] Lancellotti | Time-asymptotic wave propagation in collisionless plasmas[END_REF][START_REF] Lancellotti | Nonlinear Landau damping[END_REF][START_REF] Yamaguchi | Conditions for predicting quasistationary states by rearrangement formula[END_REF] for details.) We firstly divide f (q, p, t ) in two ways: One is the naive perturbation decomposition,

f = f 0 + g p , ( 65
)
where g p is the perturbation around f 0 , and the other one is the asymptotic-transient (AT) decomposition

f = f A + g T , (66)
where g T is the T-term satisfying lim t →∞ g T = 0. According to the AT decomposition, the potential is also decomposed as

V [ f ] + h(t )Ψ = V A + V T + hΦ, V T = V [g T ] + (h(t ) -h)Φ. ( 67
)
Substituting Eqs. ( 65) and (67) into the Vlasov equation, and omitting the nonlinear term coupling with the T-field V T , we have the T-linearized Vlasov equation

∂ f ∂t + {H A , f } + {V T , f 0 } = 0. ( 68
)
It should be noted that the nonlinearity still remains in the term {H A , f }. A solution f TL to the T-linearized equation is implicitly given by f TL (q, p, t ) = f ON (q, p, t ) + f LA (q, p, t ), f ON (q, p, t ) = e -t {H A ,•} f 0 (q, p),

f LA (q, p, t ) = - ∫ t 0 e -(t -s){H A ,•} ( F T • ∂ f 0 ∂p ) d s, ( 69
)
where we introduce the operator {H A , •}a = {H A , a} for any function a(q, p), and F T = -∂V T /∂q. The terms f ON and f LA are called O'Neil term and Landau term respectively [START_REF] Lancellotti | Time-asymptotic wave propagation in collisionless plasmas[END_REF][START_REF] Lancellotti | Nonlinear Landau damping[END_REF].

When the asymptotic stationary state f A exists, it can be picked by taking the long-time average of f TL , and we have

f A (q, p) = lim τ→∞ 1 τ ∫ τ 0 f TL (q, p, t )d t (70)
within an order of the T-lineatization method. Then, our next job is to compute the long-time average of f ON and f LA , but there are several difficulties to this for the multidimensional systems.

In the 1D systems, it is shown that

lim τ→∞ 1 τ ∫ τ 0 e -t {H A ,•} a(q, p)d t = 〈a〉 H A (q ,p ) , (71)
by use of the angle-action variables of H A [START_REF] Lancellotti | Time-asymptotic wave propagation in collisionless plasmas[END_REF][START_REF] Lancellotti | Nonlinear Landau damping[END_REF]. However, in our case, the angle-action variables cannot be constructed, so that it is unclear that this ergodic like formula holds true or not.

There is another problem, slowly algebraic damping of Tforce-field F T . In Ref. [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF] we use the fact that F T damps rapidly (∼ t -ν with ν ≥ 2) for the 1D systems [START_REF] Barré | Algebraic damping in the one-dimensional Vlasov equation[END_REF] and

lim t →∞ ∫ ∞ t F T (t )d t = 0 (72)
when we compute lim t →0 f LA . Meanwhile, in the multidimensional Vlasov systems [START_REF] Barré | On algebraic damping close to inhomogeneous Vlasov equilibria in multi-dimensional spaces[END_REF] and the 2D Euler equations [START_REF] Bouchet | Large time behavior and asymptotic stability of the 2D Euler and linearized Euler equations[END_REF], the T-force-field F T damps as or slower than t -1 , so that the integral ∫ ∞ 0 F T (t )d t is not defined in the L 1 meaning apparently. Sometimes, the transient part is asymptotically F T ≍ e -i Ωt t -γ (0 < γ ≤ 1) with Ω ̸ = 0, and the integral ∫ ∞ 0 F T (t )d t exists in the Riemannian meaning. In this case, one should be more careful when one computes the integrals and takes the limit. It should be remarked that there exists a case that Ω = 0 [START_REF] Barré | On algebraic damping close to inhomogeneous Vlasov equilibria in multi-dimensional spaces[END_REF], so that it should be checked for each system. The relation between the nonlinear response formula obtained by considering the constraint conditions and a solution to the T-linearized Vlasov equation might be an interesting future problem.

VII. SUMMARY AND PERSPECTIVE

The linear response formula has been derived without use of the analytic solution of the single particle orbit or the angle-action variables of the effective Hamiltonian. The present method improves the generalized linear response formula obtained in Ref. [START_REF] Patelli | General linear response formula for non integrable systems obeying the Vlasov equation[END_REF] when the back ground density function is a monotonically decreasing function of the effective Hamiltonian H 0 . The response formula (21) results in the one obtained in the previous studies for 1D systems [START_REF] Patelli | Linear response theory for long-range interacting systems in quasistationary states[END_REF][START_REF] Ogawa | Linear response theory in the Vlasov equation for homogeneous and for inhomogeneous quasistationary states[END_REF][START_REF] Ogawa | Non-mean-field critical exponent in a mean field model: Dynamics versus statistical mechanics[END_REF][START_REF] Ogawa | Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems[END_REF], and is numerically confirmed by use of the 2D HMF model. Further the nonlinear response formula [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF] has been derived via the same strategy, when the asymptotic solution f A to the T-linearized Vlasov equation is monotonically decreasing function of the effective Hamiltonian H A .

The nonlinear response theory based on the Tlinearization method deals with the nonlinearity of order O(h ν ) with 1 < ν < 2 [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF]. It should be noted that it is difficult to obtain the nonlinear response of order higher than O(h 2 ) successively via the proposed method, so that the error O(h 2 ) is unavoidable up to now. Let δ f n be a response of order O(h n ) for n ∈ Z. The condition Eq. [START_REF] Patelli | Linear response theory for long-range interacting systems in quasistationary states[END_REF] in the nonlinear regime O(h 2 ) is written as

0 = ( R 1 (H 0 )δ f 2 + R 2 (H 0 )δ f 2 1 ) d qd p = ( R 1 (H 0 )〈δ f 2 〉 H 0 + R 2 (H 0 )〈δ f 2 1 〉 H 0 ) d qd p, (73)
where R 1 = c ′ ( f 0 ) and R 2 = c ′′ ( f 0 )/2. It is quite difficult to obtain explicitly δ f 2 satisfying this equation for any c unlike the linear regime. Then, the error O(h 2 ) is unavoidable in both naive perturbation and T-liearization methods.

In the present article, the form of perturbation is restricted so as to subject to the constraint conditions coming from Casimir invariants at the linear order. By use of the form of constraint conditions [START_REF] Ogawa | Nonlinear response for external field and perturbation in the Vlasov system[END_REF], it is possible to take into account on the Casimir constraints when we derive the formal stability criterion without use of angle-action variables and this is a topic of forthcoming paper [START_REF] Ogawa | Stability criterion of spatially inhomogeneous solutions to Vlasov equation[END_REF].

In this article, we exert the uniform external force to the systems without integrability. We may also consider the case that unperturbed system is integrable but an external force breaks its integrability. It might be an interesting future work, how the local chaos induced by the static external field affects meso-or macro-scopic properties of systems. Such a phenomenon is found in a toy-model with one charged particle confined in cylindrical or toroidal magnetic fields [START_REF] Cambon | Chaotic motion of charged particles in toroidal magnetic configurations[END_REF][START_REF] Ogawa | Full particle orbit effects in regular and stochastic magnetic fields[END_REF].

as follows respectively: When -2M -P < H 0 < P , we have The integral F ′ 0 (H 0 )〈cos x〉

FIG. 2 .

 2 FIG. 2. (Color online) Time series of order parameters: The panel (a) is the time series for M x and the panel (b) for P x y . The temperature T = 0.3, the number of particles N = 4 × 10 6 , and the time step δt = 0.05. For each panel, the upper (red) curve is for h = 0.025 and the lower (blue) one for h = 0.

FIG. 3 .

 3 FIG. 3. (Color online) δM x , δM y , and δP x y as functions of h. The lines are the linear responses obtained theoretically and the crosses are responses obtained numerically. We set temperature of the initial states as T = 0.3 (Left column, panels (a, c, e)) and T = 0.4 (Right column, panels (b, d, f)), and a number of particles N = 4 × 10 6 and the time step δt = 0.05.

  〈cos x〉 H 0 = 〈cos y〉 H 0 range of the arccosine function is [0, π]. When P < H 0 < 2M -P , we have〈cos x〉 H 0 = 〈cos y〉 H 0 in G n (n = 1, 2,[START_REF] Binney | Galactic dynamics[END_REF][START_REF] Yamaguchi | Stability criteria of the Vlasov equation and quasi-stationary states of the HMF model[END_REF]

  y), we have〈cos x〉 H 0 = 〈cos y〉 H 0 = 〈cos x cos y〉 H 0 = 0. (45)Then we have to compute 〈cos x〉 H 0 , 〈cos y〉 H 0 , and 〈cos x cos y〉 H 0 when H 0 < 2M -P , and these are exhibited in Appendix C. We next derive an explicit form of linear responses,

	δM x =	cos xδ f d qd p, δM y =	cos yδ f d qd p,
				(46)
		δP cc =	cos x cos yδ f d qd p.

TABLE I .

 I Initial equilibria and zero-field susceptibilities

	tively by
	T	M	P d δM x/y /d h| h=0 d δP cc /d h| h=0
	0.3 0.90223 0.81556	0.034428	0.059089
	0.4 0.84269 0.71910	0.071298	0.099709

  + V (x, y) )and σ(H 0 ) is defined in Eqs. (C6) and (C9). The similar terms in G 2 and G 4 are computed in the same manner.

						2 H 0	d qd p included in G 1 and
	G 2 is computed as follows;
						∫	∫
	F ′ 0 (H 0 )〈cos x〉 2 H 0 d qd p = 2π	pd p	F ′ 0 (H 0 )〈cos x〉 2 H 0 d q
	= 2π = 2π	∫ 2M -P -2M -P ∫ 2M -P -2M -P	d H 0 F ′ 0 (H 0 )〈cos x〉 2 H 0 F ′ 0 (H 0 )〈cos x〉 2 H 0 σ(H 0 )d H 0 , ∫ Θ ( H 0 + V (x, y)	)	d q
	where p = ∥p∥ =	√ 2	(	H 0

ACKNOWLEDGMENTS

The author is grateful to Yoshiyuki Y. Yamaguchi and Xavier Leoncini for valuable discussions. He acknowledges the financial support of the A * MIDEX project (n • ANR-11-IDEX-0001-02) funded by the "investissements d'Avenir" French Government program, managed by the French National Research Agency (ANR).

Appendix A: Conservation of the Casimir functionals [START_REF] Spohn | Large scale dynamics of interacting particles[END_REF] It is shown that the Casimir functional [START_REF] Spohn | Large scale dynamics of interacting particles[END_REF] is conserved in the Vlasov dynamics. Taking the time derivative of C [ f ], we have

Under the conditions asserted above Eq. ( 8), the boundary terms vanish and the left hand side of Eq. (A1) is

Appendix B: Derivation of Eq. (43)

We derive Eq. ( 43). The all we have to do is to perform integration with respect to p in the left hand side of Eq. ( 43

where Θ(x) = 0 (resp.1) when x < 0 (resp. x ≥ 0) is the Heaviside step function, and we have used the relation,

Thus we have

Appendix C: Computation of 〈cos x〉 H 0 and 〈cos x cos y〉 H 0 On the iso-H 0 curve, x and y satisfy

for H 0 ∈ [-2M -P, P ] and