Julien Barré
email: jbarre@unice.fr

Alain Olivetti

Yoshiyuki Y Yamaguchi

Algebraic damping in the one-dimensional Vlasov equation

Keywords: numbers: 05.20.Dd, 45.50.-j, 52.25.Dg, 98.10.+z

We investigate the asymptotic behavior of a perturbation around a spatially non homogeneous stable stationary state of a one-dimensional Vlasov equation. Under general hypotheses, after transient exponential Landau damping, a perturbation evolving according to the linearized Vlasov equation decays algebraically with the exponent -2 and a well defined frequency. The theoretical results are successfully tested against numerical N -body simulations, corresponding to the full Vlasov dynamics in the large N limit, in the case of the Hamiltonian mean-field model. For this purpose, we use a weighted particles code, which allows us to reduce finite size fluctuations and to observe the asymptotic decay in the N -body simulations.

Introduction

Systems of particles interacting through long-range forces are usually described over a certain time scale by Vlasov equations. This situation is encountered in various fields of physics: plasma physics, self-gravitating systems, wave-particles interactions for instance. One may add here two dimensional fluid dynamics, since the 2D Euler equation shares many properties with the Vlasov equation.

A Vlasov equation usually admits a continuous infinity of stationary states. Investigating the stability of these states is a natural question. In a celebrated paper [1], Landau considered stationary states of a plasma which are homogeneous in space, and addressed the issue of the asymptotic behavior of a small perturbation through a Laplace transform analysis of the Vlasov equation linearized around the stationary state. This was the starting point of an extremely abundant research on Landau damping in plasma physics, and more generally on the fate of perturbations around stationary states of the Vlasov equation, usually homogeneous in space. We focus now on the known mathematically rigorous results, all of them obtained in the context of homogeneous stationary solutions of the Vlasov equation. A rigorous linear treatment "à la Landau" is provided for instance in [2,3]: under very strong regularity hypothesis for the stationary state and the perturbation (both should be analytic functions of the velocity), it proves the exponential asymptotic decay of a solution of the linearized Vlasov equation in a bounded spatial domain, as predicted by Landau. However, it is known that such a linear solution may decay at a much slower rate, and even not decay at all, when the analyticity hypothesis for the perturbation (see [START_REF] Weitzner | Magneto-Fluid and Plasma Dynamics[END_REF][START_REF] Crownfield | [END_REF] for references in the physics literature) or the bounded spatial domain hypothesis [6] is not satisfied. Mouhot and Villani proved recently the asymptotic exponential decay of a perturbation evolving according to the full non linear equation in a bounded domain, using analytic norms (which implies the analyticity of the stationary state and its perturbation) [7,8]. However, it is shown in [9] (see also [10] for a previous non rigorous treatment) that such a non linear damping may fail if weaker norms are used.

In many cases, one would like to ask the same question in the case of non homogeneous stationary solutions. In particular, this is the case in the astrophysical context, where the concept of Landau damping is often used [11]. [2] contains a discussion of the linearized Vlasov equation around a non homogeneous stationary solution. However, the analysis of this situation faces some technical difficulties. At a formal level, these difficulties may be partially overcome using the "matrix" formulation introduced in the works of Kalnajs [12] and Polyachenko and Schukhman [13]. Since then, many papers have computed linear instability rates in an astrophysical context using this method (see [14,15], to mention just a few). Recently, some new methods to investigate the stability and compute such growing rates for unstable non homogeneous stationary states have been introduced and tested on toy models [16,17,18,19]. Beyond these unstable states, stable oscillating modes in the context of 1D self gravitating models have also been investigated in some details [20,21]. However, there seems to be very few analytical studies of what would be the strict analog of Landau damping: the decay of a perturbation close to a stable non homogeneous solution of the Vlasov equation. The reason is probably that such a study is technically more difficult, since it requires the further step of an analytical continuation. Weinberg has performed this step using a numerical approximation and computed the analog of a Landau damping rate in some cases [22] (see also [23]). However, at variance with the homogeneous case, we shall see that this Landau damping rate never controls the asymptotic decay of the perturbation.

It is well known that the 2D Euler equation, as well as other related conservative 2D fluid equations, share a lot of similarities with the Vlasov equation. There have been many linear investigations of perturbations around stationary base flows such as a vortex or a shear flow, starting with Rayleigh [24]. Transposing these studies in the context of the Vlasov equation, linearized around a non homogeneous stationary state, one would expect the following generic picture: the dispersion relation has branch point singularities on the real axis, which make the analytical continuation procedure used in Landau's work trickier; if the stationary state is stable, the asymptotic decay of the perturbation is controlled by these branch point singularities, and is in general algebraic. An exponential decay of the perturbation "à la Landau" may be visible, but it is restricted to an intermediate time window, before the algebraic decay kicks in. In particular, it seems reasonable to expect that the exponential damping studied by Weinberg in the gravitational case [22] is actually a transient effect (this statement does not preclude of course its potential physical relevance). We remark that such an exponential decay followed by an algebraic decay has also been studied in a dissipative system of coupled oscillators [25]; the decay mechanism in this case is similar to the one in Vlasov systems.

Smereka [26] investigated the asymptotic behavior of a linearized 1D Vlasov equation around non homogeneous stationary states, and reached the conclusions outlined above; his analysis is however limited to a particular class of interactions, and does not contain direct numerical simulations. Moreover, his conclusion is affected by a non generic choice of the perturbations, as we will discuss later on. Besides [26], we are not aware of other studies tackling this problem.

To be more specific, our goals in this paper are (i) Study as generally as possible the linear asymptotic decay of perturbations around non homogeneous stationary states of a 1D Vlasov equation.

(ii) Perform explicitly the computations in the simple case of the Hamiltonian meanfield (HMF) model.

(iii) Compare the analytical results with detailed numerical computations; this is made possible by the use of a simple toy model such as HMF.

This paper is organized as follows. All theoretical results are described in Sec.2. From the linearized Vlasov equation, we derive two equations for the Fourier-Laplace components of perturbation and of potential in Sec.2.1. In Sec.2.2 these two equations are formally solved in the Laplace space with the help of biorthogonal functions; an example of such functions is given in Sec.2.3 for spatially periodic systems. The asymptotic dynamics of a perturbation is determined by its singularities in the Laplace space; we classify the singularities in Sec.2.4. Focusing on one type of singularities, which is in many cases the relevant one, we show that the perturbation asymptotically decays algebraically with the exponent -2 in Sec.2.5. These general results are applied to the HMF model in Sec.2.6; in this case, we show that the two components of the magnetization vector decay with exponents -3 and -2 respectively. We note that the exponent -3 comes from a special cancellation due to a symmetry of the HMF model. The theoretical results for the HMF model are numerically examined in Sec. 3. In order to test the two exponents separately, we introduce two types of perturbation in Sec.3.1 and in Sec.3.2 respectively. The last section 4 is devoted to conclusion.

Theory

We consider the Vlasov equation in one dimension for the one-particle distribution function f (x, p, t),

∂ t f + ∂ p H∂ x f -∂ x H∂ p f = 0, (1)
where x ∈ D ⊂ R is the position variable, p ∈ R the conjugate momentum variable, and H is the one-particle Hamiltonian defined by

H[f ](x, p, t) = p 2 2 + φ[f ](x, t) + φ ext (x). (2)
The potential φ[f ](x, t) is defined from the two-body interaction potential v and the distribution f as

φ[f ](x, t) = D R v(x -y)f (y, p, t)dpdy, (3)
and φ ext (x) derives from an external force F ext (x) as

F ext (x) = -∂ x φ ext (x). ( 4
)
Although it would be interesting to consider a time dependent and/or non-potential external force, we focus in this article on a static and potential force such as [START_REF] Weitzner | Magneto-Fluid and Plasma Dynamics[END_REF]. The domain D is typically R or [0, 2π]:

(i) Case 1, D = R:
The model is defined on the whole real line. A typical example is a 1D self-gravitating system, a very much studied caricature of the more realistic 3D self-gravitating systems.

(ii) Case 2, D = [0, 2π]: The system has periodic boundary conditions. Such boundary conditions are sometimes used in plasma physics; it is also the setting of the HMF model, a paradigmatic toy model for long range interacting systems (see Sec. 2.6).

We give in the following subsections a general analysis of the linearized Vlasov equation, which applies to the both cases.

Two equations to be solved

Let us call f 0 a stationary solution of the Vlasov equation ( 1). The one-particle Hamiltonian H[f 0 ](x, p) is integrable, and its trajectories are level curves of H[f 0 ]. If D = R, one could imagine trajectories that are unbounded and not periodic. However, the stationarity of f 0 imposes that f 0 (x, p) is constant along the trajectories; this imposes that such unbounded non periodic trajectories either do not exist, or are not populated by the density f 0 . As a consequence, we can always define the angle θ and the action J from the original coordinate (x, p). Strictly speaking, this change of variables is not always one-to-one: whenever there exists a separatrix in the one-particle Hamiltonian H[f 0 ](x, p), there are distinct trajectories with the same value of the action. Thus, although we will formally use this change of variables, a careful treatment may be needed for specific potentials (see for instance [28]). We add a perturbation f 1 to the stationary solution, and start from an initial condition f 0 + f 1 (t = 0). The Hamiltonian H is linear with respect to f , and we have

H[f 0 + f 1 ] = H[f 0 ] + φ[f 1 ], where φ[f 1 ](x, t) = φ 1 (x, t) = D R v(x -y)f 1 (y, p, t)dpdy. (5)
The stationary solution f 0 must be constant along trajectories of the Hamiltonian H[f 0 ]; a sufficient condition for this is to take f 0 a function of the action alone. It is not a necessary condition if there is a separatrix, since two disjoint trajectories correspond to the same action. Thus, with a slight loss of generality, we will assume in the following that f 0 may be written as f 0 (J). The unperturbed part of one-particle Hamiltonian is also a function of the action alone, and written as H[f 0 ](J). Using the angle-action variables (θ, J), we have the linearized Vlasov equation

∂ t f 1 + Ω(J)∂ θ f 1 -f ′ 0 (J)∂ θ φ 1 = 0 , (6)
where we have defined the frequency Ω(J) = dH[f 0 ]/dJ. Notice that the external potential does not enter in this linear equation; it appears implicitly of course through the definition of the angle-action variables.

To analyze the linearized Vlasov equation ( 6), we introduce the Fourier-Laplace transform û(m, J, ω) of a function u(θ, J, t) as

û(m, J, ω) = π -π dθ e -imθ +∞ 0 dt e iωt u(θ, J, t) (7)
where m is an integer and Im(ω) large enough to ensure convergence. The inverse transform is then

u(θ, J, t) = 1 (2π) 2 ∞ m=-∞ Γ dω û(m, J, ω)e -iωt e imθ ( 8
)
where Γ is a Bromwich contour running from -∞ + iσ to +∞ + iσ, and the real value σ is larger than the imaginary part of any singularity of û(m, J, ω) in the complex ω-plane.

Performing a Fourier transform with respect to θ and a Laplace transform with respect to time on (6), we obtain, after simple algebraic manipulations f1 (m, J, ω) = A(m, J, ω) φ1 (m, J, ω) + B(m, J, ω),

where

A(m, J, ω) = mf ′ 0 (J) mΩ(J) -ω , (10)
B(m, J, ω) = g(m, J) mΩ(J) -ω , (11)
and ig(m, J) is the Fourier transform of the initial perturbation f 1 (θ, J, t = 0) with respect to θ. We assume

f 1 (θ, J, t = 0)dθdJ = 0 ( 12
)
and hence g(0, J) = 0. The two equations ( 5) and ( 9) relate f 1 to φ 1 . The strategy is now to combine these two equations to compute f 1 and φ 1 . One sees however that f 1 is easily obtained in angle-action variables whereas φ 1 is more easily expressed in the original (x, p) variables. To overcome the difficulty of the two natural coordinate basis (x, p) and (θ, J), we follow the standard procedure and introduce two families of biorthogonal functions [27,12,13].

Biorthogonal functions

We introduce the linear mapping L v by

L v : d → u (13)
where

u = L v • d is defined by u(x) = D v(x -y)d(y)dy. (14)
We assume that there exist two index sets I and I ′ which satisfy I ′ ⊂ I ⊂ Z, and two families {d j (x)} j∈I and {u k (x)} k∈I ′ which satisfy the following conditions:

(i) {d j (x)} j∈I is linearly independent and any density function ρ(x) may be expanded as

ρ(x) = j∈I a j d j (x), (15)
(ii) {u k (x)} k∈I ′ is linearly independent and spans Im(L v ); any function g(x) ∈ Im(L v ) may be expanded as

g(x) = k∈I ′ b j u j (x), (16)
(iii) the two families are orthogonal to each other:

(d j , u k ) = D d j (x)ū k (x)dx = λ k δ jk , (j ∈ I, k ∈ I ′ ) (17)
with λ k = 0, and where δ jk is the Kronecker δ.

(iv) For all k ∈ I ′ , d k and u k satisfy the property

u k = L v • d k . For all k ∈ I \ I ′ , L v • d k = 0.
Let us start from the definition of the potential φ 1 (x, t) (5). This definition may be rewritten as

φ 1 (x, t) = D v(x -y)ρ 1 (y, t)dy (18)
by using the perturbation density

ρ 1 (x, t) = R f 1 (x, p, t)dp. ( 19
)
From the assumption (i), the perturbation density is expanded in the form

ρ 1 (x, t) = j∈I a j (t)d j (x). ( 20
)
Substituting ( 20) into (18) and using (iv), we obtain an expansion for φ 1 in the form

φ 1 (x, t) = k∈I ′ a k (t)u k (x). ( 21
)
The Fourier-Laplace transform of ( 21) is expressed by

φ1 (m, J, ω) = k∈I ′ ãk (ω)c km (J), (22)
where ãk (ω) is the Laplace transform of a k (t),

ãk (ω) = ∞ 0 a k (t)e iωt dt, (23)
and c km (J) is the Fourier transform of u k (x),

c km (J) = π -π u k (x)e -imθ dθ. (24)
Substituting ( 22) into (9), we obtain

f1 (m, J, ω) = A(m, J, ω) k∈I ′ ãk (ω)c km (J) + B(m, J, ω). (25)
The inverse Fourier transform of (25) gives f1 (θ, J, ω), the Laplace transform of f 1 (θ, J, t):

f1 (θ, J, ω) = 1 2π m∈Z f1 (m, J, ω)e imθ (26) = 1 2π m∈Z A(m, J, ω) k∈I ′ ãk (ω)c km (J) + B(m, J, ω) e imθ . ( 27
)
We observe that f1 (θ, J, ω), and hence f 1 (θ, J, t), are determined by the {ã k (ω)} k∈I ′ , so that the subset {ã j (ω)} j∈I\I ′ is not necessary. We therefore seek a solution for the subfamily {ã k (ω)} k∈I ′ instead of the whole family {ã j (ω)} j∈I . For this purpose, we multiply (26) by ūl (x) (l ∈ I ′ ) and integrate it over θ and J. Using the biorthogonality relation (17) and noting that the change of variable (x, p) → (θ, J) is symplectic dx ∧ dp = dθ ∧ dJ, the left-hand-side of (26) becomes

f1 (θ, J, ω)ū l (x)dθdJ = R D f1 (x, p, ω)ū l (x)dxdp (28) = D ρ1 (x, ω)ū l (x)dx = ãl (ω)λ l . (29)
To derive the last equality, we have used the fact that the Laplace transform of

ρ 1 (x, t) (20) is ρ1 (x, ω) = j∈I ãj (ω)d j (x). ( 30
)
On the other hand, the right-hand-side of ( 26), submitted to the same operations, becomes 1 2π m∈Z A(m, J, ω)

k∈I ′ ãk (ω)c km (J) + B(m, J, ω) ūl (x)e imθ dθdJ = 1 2π m∈Z k∈I ′ ãk (ω) A(m, J, ω)c km (J)c lm (J)dJ + 1 2π m∈Z B(m, J, ω)c lm (J)dJ. (31)
Remembering the definitions of A(m, J, ω) and B(m, J, ω), we introduce the functions

F lk (ω) = 1 2π m∈Z mf ′ 0 (J) mΩ(J) -ω c km (J)c lm (J)dJ, l, k ∈ I ′ (32) G l (ω) = 1 2π m∈Z g(m, J) mΩ(J) -ω clm (J)dJ, l ∈ I ′ (33)
where contributions from m = 0 vanish not only for F lk but also for G l , thanks to assumption (12). We further define the (♯I ′ )×(♯I ′ ) matrices Λ and

F (ω) = (F lk (ω)) l,k∈I ′ ,
where Λ is diagonal with elements {λ k } k∈I ′ , and the (♯I

′ )-dimensional vectors G(ω) = (G l (ω)) l∈I ′ and ã(ω) = (ã k (ω)) k∈I ′ .
Using the above matrices and vectors, the equation for ã(ω) reads in matrix form:

[Λ -F (ω)]ã(ω) = G(ω). (34)
The equation ( 34) is formally solved as

ã(ω) = [Λ -F (ω)] -1 G(ω), (35)
and the temporal evolution of the {a k (t)} k∈I ′ , and of f 1 , is obtained from the inverse Laplace transform of ã(ω). The inverse matrix [Λ -F (ω)] -1 does not always exist since the determinant det(Λ -F (ω)) is not always non-zero. This determinant is sometimes called the dispersion function, and its roots are poles of ã(ω). We will discuss the singularities of ã(ω) in Sec.2.4 after giving an example of the two families {d j } j∈I and

{u k } k∈I ′ in the domain D = [0, 2π].

Example: D = [0, 2π]

If we consider the domain D = [0, 2π] with periodic boundary condition, the interaction potential v(x) must be also 2π-periodic and is expanded in Fourier series as

v(x) = 1 2π m∈Z v m e imx , (36)
where the coefficients v m are determined by

v m = π -π v(x)e -imx dx. ( 37
)
We can choose the two families {d j (x)} j∈I and {u k (x)} k∈I ′ as

d j (x) = 1 2π e ijx , j ∈ I = Z, (38)
and

u k (x) = 2π 0 v(x -y)d k (y)dy = 1 2π v k e ikx , k ∈ I ′ (39)
where the index set I ′ is defined by

I ′ = {k ∈ Z | v k = 0}. ( 40
)
From the definition of u k (x) and I ′ the assumption (iv) is satisfied. Expansions on the two families {d j (x)} j∈Z and {u k (x)} k∈I ′ are essentially Fourier expansions, and hence the assumptions (i)-(iii) are also satisfied. The factors

λ k are λ k = v k /2π for k ∈ I ′ .
Generically, all v k are non zero, and I = I ′ = Z. However, for the HMF model, which we will introduce in section 2.6, v(x) =cos x and Im(L v ) is 2-dimensional. Accordingly the index set I ′ is I ′ = {1, -1} and the matrix Λ -F (ω) in ( 34) is a 2 × 2 matrix.

Singularities of ã(ω)

From the knowledge of the functions ã(ω), one may easily compute the time evolution of the potential and density perturbations, through an inverse Laplace transform. The asymptotic in time behavior of this inverse Laplace transform will be determined by the singularities of the functions ãk (ω) (k ∈ I ′ ) in the complex plane. We turn now to the study of these singularities. The matrix coefficients F lk (ω) and the vector coefficients G l (ω) are defined through integrals over the real variable J, see [START_REF] Buyl | [END_REF] and (33). These integrals are naturally defined in the whole half plane Im(ω) > 0. Note however, that expressions [START_REF] Buyl | [END_REF] and (33) are in general not properly defined for ω ∈ R, as in this case mΩ(J)ω may vanish. For Im(ω) ≤ 0, we will actually have to consider rather the analytical continuations of the expressions [START_REF] Buyl | [END_REF] and (33). We may have two kinds of singularities, described in the following.

The first kind of singularities is poles, coming from roots of the dispersion function det(Λ -F (ω)). If such a root exists in the half plane Im(ω) > 0, it corresponds to an eigenvalue of the linearized Vlasov operator, and it yields an exponential growth of the perturbation. In the upper half plane Im(ω) > 0, the only possible singularities for ã(ω) are the poles. We assume in the following that the reference stationary state f 0 is linearly stable, so that the determinant of Λ -F (ω) does not have any roots with Im(ω) > 0. If det(Λ -F (ω)) has a root on the real axis Im(ω) = 0, it corresponds to a purely oscillating mode. This is compatible with a linearly stable f 0 , but corresponds to a non decaying perturbation; we also assume in the following that this does not happen. The analytical continuation of det(Λ -F (ω)) in the lower half plane Im(ω) < 0 may also have roots. They correspond to "Landau poles" for ã(ω), and give rise to an exponential damping of the perturbation. This damping behavior is known as Landau damping, and has been studied in the gravitational case in [22,23], and in the HMF case in [28]. As anticipated in the introduction, we will see that this exponential damping, if it exists, is subdominant in the large time region.

The second kind of singularities comes from the integral in ( 32) and ( 33), and appears on the real axis of Im(ω) = 0. We have to study the singularities of functions

ϕ(z) = b a ψ(J) mΩ(J) -z dJ (41)
properly defined for Im(z) > 0, where ψ and Ω are real functions, and [a, b] ⊂ R∪{+∞}. J has to be thought of as the action coordinate, and Ω as the associated frequency. We assume that ψ is analytic. We may set m = 0 since the contributions from m = 0 in (32) and ( 33) vanish. We now show that ϕ(z) is regular for z ∈ R except for special points, and will classify the special points into three types. Notice that both functions F kl and G k fit in this framework.

If the equation mΩ(J) = z has no solutions in [a, b] for any z in a neighborhood of real z 0 , then ϕ is analytic in a neighborhood of z 0 . Assume now that the equation mΩ(J) = z has one or several branches of solutions

J * i (z) ∈ [a, b] in a neighborhood of z 0 , where all J * i (z) ∈ [a, b]
are regular as functions of z. In this case, ϕ can be analytically continued from the open half plane Im(z) > 0 to the neighborhood of z 0 by taking into account the possible residue contributions of the roots J * i (z), in a straightforward generalization of the "Landau prescription". Thus, generically, no singularity of ϕ appears at z = z 0 .

The singularities of ϕ are hence associated with special points z 0 , such that the branches of solutions J * i (z) ∈ [a, b] of the equation mΩ(J) = z undergo a bifurcation or are singular. This may happen in the following three types, illustrated on Fig. 1: (i) z 0 = mΩ(a) or z 0 = mΩ(b), when one or both are finite. This is a common situation, generically encountered in 1D self-gravitating systems [20] as well as in the HMF model [28]. For instance, this mechanism creates a singularity around the frequency mω 0 = mΩ(J = 0), where J = 0 corresponds to the minimum of the effective potential φ[f 0 ] + φ ext , see Fig. 1. We will show in Sec.2.5 that this singularity is logarithmic.

(ii) z 0 corresponds to a J such that Ω(J) is singular: this may be an action J s corresponding to a separatrix; an illustration is given on Fig. 1, for J = J s , ω = 0.

ω 2 ω 1 ω 0 ω =0 Ω J J s Figure 1.
An example of function Ω(J), with the special points at the origin of singularities for the associated function ϕ defined as in Eq. (41). This is not a common situation for the 1D self-gravitating models we have in mind. This generically happens however for periodic systems, where the trajectories in the one-particle Hamiltonian H[f 0 ] may be oscillating or librating, and the two regions are delimited by a separatrix. We will see such a situation in the HMF case.

(iii) z 0 corresponds to a J such that Ω ′ (J) = 0: this corresponds generically to a local maximum or minimum of the frequency (see frequencies ω 1 and ω 2 on Fig. 1). We are not aware of any model studied in the literature where this phenomenon happens. This is certainly an interesting case to study, especially in view of the results obtained in an analogous situation for the diocotron instability of a magnetically confined electron column [29], and recently for the 2D Euler equation [30]. We will not be concerned with this type of singularities in the following.

In the next subsection we estimate the asymptotic relaxation of a k (t) by considering contributions from singularities of the first type, since they should be in many cases the relevant singularities. We will confirm whether the estimation is valid by performing direct N-body simulations in Sec.3.

Contribution from singularities of the first type

We concentrate now on the singularities of the first type at z = mΩ(a) = mω a . Notice that if b = +∞, and ω ∞ = Ω(+∞) is finite, z = mω ∞ is a singular point; indeed, the number of solutions of the equation mΩ(J) = z changes from z > mω ∞ to z < mω ∞ . We expect this singularity to be irrelevant if f 0 (J) decreases rapidly enough as J → ∞, and we neglect it in the following.

To analyze the function ϕ(z) (41) around z = mω a , we expand Ω(J) and ψ(J) in power series around J = a as

Ω(J) = ω a + c(J -a) + O (J -a) 2 (42)
and

ψ(J) = d(J -a) ν + o ((J -a) ν ) ( 43
)
for some ν ∈ N. The exponent ν will be determined for a given m later. Substituting these assumptions into ϕ and changing Ja to J, this yields the following singular part for ϕ:

ϕ(z) = C b-a 0 J ν J -ζ dJ + (regular part) (44)
with the constant C = d/mc and ζ = (zmω a )/mc. Using the equality

J ν J -ζ = J ν-1 + ζ J ν-1 J -ζ (45)
recursively, the singular part of ϕ(z), which comes from the lower bound of the integral (44), behaves as a logarithm times a power around ζ = 0:

ϕ(z) = cste × (z -mω a ) ν ln(z -mω a ) + (regular part). ( 46
)
Thus, at ω = mω a , the matrix elements F kl (ω) and the vector elements G k (ω) have singularities of the type (ωmω a ) ν ln(ωmω a ). From (35), ãk (ω) is expressed as a sum, product and ratio of F kl (ω) and G k (ω) functions. Thus, the leading singularity of ãk (ω) at ω = mω 0 is also of the type (ωmω a ) ν ln(ωmω a ) for some ν. a k (t) is the inverse Laplace transform of ãk (ω). Since ãk has only logarithmic singularities on the real axis, we may deform the contour Γ of the inverse Laplace transform down to the real axis. The inverse Laplace transform then becomes an inverse Fourier transform. The asymptotic decay of a k (t) is then determined by the strongest singularity of ãk (ω) on the real axis (see [START_REF] Lighthill | Introduction to Fourier Analysis and Generalized Functions[END_REF] p.52). A singularity such as (46) yields an asymptotic decay as (see [START_REF] Lighthill | Introduction to Fourier Analysis and Generalized Functions[END_REF] p.42)

a k (t) ∼ cste e -imωat t ν+1 . (47)
See Appendix A for a heuristic explanation on how to obtain estimates such as (47).

To obtain ν for a given m, we need to introduce some assumptions about the system we consider. We assume that the potential created by the stationary state (interaction potential + external potential φ[f 0 ](x) + φ ext (x)) has a single minimum at J = a = 0, and is quadratic with respect to x around its minimum at leading order, except for an irrelevant constant term. This is the case in many situations of interests, such as stationary states for a 1D self-gravitating system. We also assume that f 0 (J) is analytic, and decays fast enough at infinity (for instance exponentially). This assumption excludes for instance truncated f 0 , or compactly supported stationary states. As an example, the thermal equilibria of a 1D self-gravitating system or of the HMF model satisfy all assumptions. Finally, we assume that the perturbation is also analytic.

Under the above assumptions, we now estimate the exponent ν. From ( 32) and ( 33), we see that the function ψ(J) reads ψ(J) = f ′ 0 (J)c km (J)c lm (J) for F kl (ω) functions and

ψ(J) = g(m, J)c lm (J)
for G(ω) functions; we now expand these functions with respect to J. Let us start from the function c km (J). In the limit of J → 0, motion is harmonic with the frequency ω 0 , and the position x is written in polar coordinates using angle-action variables as x ∝ J 1/2 sin θ. Expanding u k (x) with respect to x, and substituting the above expression of x into the expansion, the function c km (J) reads

c km (J) ∼ 2π 0 ∞ n=0 u (n) k (0) n! J n/2 sin n θ e -imθ dθ,
where u

(n) k
is the n-th derivative of u k . The first non-vanishing term corresponds to n = |m|, and hence the leading order for c km (J) in a small J expansion is

c km (J) ∼ J |m|/2 . ( 48
)
c km (J) is the Fourier transform of u k (x); in a similar way, the leading order in a small J expansion for g(m, J), the Fourier transform of the initial perturbation, is

g(m, J) ∼ J |m|/2 , ( 49
)
since the perturbation is regular. The function f ′ 0 (J) is regular and hence the leading order is constant. Consequently, the function ψ(J) is, at leading order

ψ(J) ∼ J |m| (50)
both for F and G functions. Hence, for each term in the infinite series defining the coefficient F kl (ω) [START_REF] Buyl | [END_REF] and G l (ω) (33), we have ν = |m|.

Going back to (32), (33) and making use of section 2.4, we see that the strongest singularities for F kl (ω) and G l (ω) come from the m = ±1 terms in the sum over m, so that the exponent is ν = 1. We conclude using (47) that the functions a k (t), under the hypothesis of this section, decay as e -iω 0 t /t 2 .

This result has to be compared with [26], which finds a decay exponent 3/2. Since we have performed the same kind of analysis as Smereka does in [26], this discrepancy is surprising even if the class of Hamiltonians studied is different. It may be traced back to the fact that this author uses the following hypothesis for the perturbation: lim J→0 g(m, J) = 0 even for m = 0; this would correspond to a singular perturbation, since the initial perturbation f 1 is not well-defined in the limit J → 0, because lim J→0 e imθ g(m, J) = lim J→0 e imθ ′ g(m, J) as soon as e im(θ-θ ′ ) = 1. This singular perturbation implies that g(m, J) ∼ J 0 in the limit J → 0 instead of (49), and hence ψ(J) = g(m, J)c lm (J) ∼ J |m|/2 . Accordingly, the strongest singularities for G l (ω) (33), coming from the m = ±1 terms, corresponds to ν = 1/2. Using considerations similar to the ones described in Appendix A Smereka showed that this gives a decay exponent 3/2. Considering in Smereka's setting a regular perturbation, as is more natural (and as actually assumed in Eq. ( 21) of [26]), would yield the -2 exponent also for the class of Hamiltonians studied in [26].

We may also compare this result to the asymptotic decay of perturbations around a stationary 2D shear flow or a vortex. In these cases, the longitudinal (resp. transverse) velocity perturbation asymptotically decays as 1/t (resp. as 1/t 2 ), without temporal oscillations.

Example of the HMF model

In this section, we analyze in more details a specific example, the HMF model whose Hamiltonian is

H(x, p) = N j=1 p 2 j 2 + 1 2N N j=1 N k=1 [1 -cos(x j -x k )]. (51)
The canonical equation of motion of the HMF model is described through the magnetization M = (M x , M y ) defined by

M x = 1 N N j=1 cos x j , M y = 1 N N j=1 sin x j , (52)
and hence the computational cost is O(N) for each time step, although the number of interactions between the N particles is O(N 2 ). This advantage allows precise numerical tests of the predictions. The associated Vlasov equation reads:

∂f ∂t + p ∂f ∂x - ∂φ[f ] ∂x ∂f ∂p = 0 , (53)
with

φ[f ](x, t) = -M x [f ] cos x -M y [f ] sin x (54)
and

M x [f ](t) = ∞ -∞ 2π 0 cos xf (x, p, t)dxdp, (55)
M y [f ](t) = ∞ -∞ 2π 0 sin xf (x, p, t)dxdp. ( 56
)
Note that this is a pendulum potential, so that the dynamics admits a separatrix. The action-angle variables are explicitly written in terms of elliptic integrals [28]. Without loss of generality, we consider a stationary solution with M y = 0, and write M

x (t) and M

(1) y (t) the magnetization perturbations.

For convenience, we choose real functions for the family {u k } k∈I ′ : u c (x) = cos x and u s (x) = sin x and I ′ = {c, s}. The coefficients a c (t) and a s (t) of the potential φ 1 (x, t),

φ 1 (x, t) = -a c (t) cos x -a s (t) sin x, correspond to M (1)
x (t) and M

(1) y (t) respectively. The expansions of u c and u s in the Fourier series of the angle variable θ define the coefficients {c cm (J)} m∈Z and {c sm (J)} m∈Z as in (24); to simplify the notations, we rename these coefficients c m (J) and s m (J) respectively.

This choice for the family {u k } makes the matrix Λ -F (ω) diagonal [28], with

F cc (ω) = 1 2π m∈Z mf ′ 0 (J) mΩ(J) -ω |c m (J)| 2 dJ, (57)
F ss (ω) = 1 2π m∈Z mf ′ 0 (J) mΩ(J) -ω |s m (J)| 2 dJ, (58)
F cs (ω) = F sc (ω) = 0. ( 59
)
We also have the following expressions for G:

G c (ω) = 1 2π m∈Z g(m, J) mΩ(J) -ω cm (J)dJ, (60)
G s (ω) = 1 2π m∈Z g(m, J) mΩ(J) -ω sm (J)dJ. ( 61
)
A priori, according to the general discussion in the previous section, the leading singularity of F cc and F ss is located at ω = ω 0 = √ M 0 , and has an index ν = 1, since ω 0 here plays the role of ω a in the previous section 2.5. However, due to the symmetries of the system, a further cancellation occurs: the function c m (J) identically vanishes for all m odd and all J < J s , with J s the action at the separatrix. Thus, the strongest singularity for F cc and G c actually comes from the m = 2 term in ( 57) and (60), and is located at ω = 2ω 0 ; its index is ν = 2. This has an interesting consequence on the asymptotic behavior of M x (t): since it is now governed by a singularity with index ν = 2, we expect M (1) x (t) ∼ e -2iω 0 t t 3 for t → ∞ (62)

at variance with M

(1) y (t), which is still governed by a ν = 1 singularity

M (1) y (t) ∼ e -iω 0 t t 2 for t → ∞. (63)
This feature makes the HMF model particularly suitable for a numerical test of the theory developed in this section, as we should be able to probe two different asymptotic behaviors for M

(1)

x (t) and M

(1) y (t).

Numerical simulations

In this section we test numerically the linear predictions of the previous section, on the example of the HMF model, by solving the whole (non linear) Vlasov equation. Solving the Vlasov equation over long times may be a very heavy numerical task, or even impossible with current computers. In this case, several features help: the model is one dimensional, it is particularly simple, and we only need to solve the Vlasov equation close to a stationary state.

A natural strategy could be to solve directly the Hamiltonian N-body dynamics, with N large enough; we know that this provides an approximation to the continuous Vlasov evolution. We have found that the finite-N fluctuations were too big to allow a test of the asymptotic in time regime. Another strategy would be to use a standard Vlasov solver; for instance, a semi-Lagrangian method has already been used for HMF [START_REF] Buyl | [END_REF]. This resulted in very heavy computations. Finally, we have chosen to introduce an algorithm relying on a Hamiltonian simulation of appropriately weighted particles [33], it provides a very convenient tool to test the theoretical predictions. The algorithm is described and discussed in Appendix B. We have tested the weighted particles algorithm against (i) a semi-lagrangian code (ii) a simple unweighted N-body code (results not reported). The temporal evolution of the magnetization from the three codes are in good agreement up to a certain time. The unweighted N-body simulation becomes dominated by finite size fluctuations much earlier than the weighted particles' one.

All simulations discussed in the following were performed using the weighted particles algorithm, close to a thermal equilibrium stationary state, parametrized by the temperature T = 1/β:

f 0 (x, p) = N e -β(p 2 /2-M 0 cos x) ( 64
)
where N is the normalization and the magnetization M 0 (T ) is solution of a consistency equation [34]. The thermodynamical equilibrium state of the HMF is non homogeneous (that is M 0 = 0) as soon as T < 0.5. In the following we only use T = 0.1; larger temperatures resulted in increased fluctuations and made it more difficult to reach the asymptotic in time regime. The magnetization is M 0 = 0.946 and the harmonic frequency is ω 0 = √ M 0 = 0.972 for T = 0.1.

Cosine perturbation

We consider first a cosine perturbation of the thermal equilibrium:

f 0 (x, p) + f 1 (x, p) = N a e -β(p 2 /2-M 0 cos x) (1 + a cos x) , (65)
where |a| is small enough. This perturbation is compatible with the (x, p) → (-x, -p) symmetry of the canonical equation of motion, and M

y (t) is then identically equal to zero. We may restrict the initial points such that p i (0) > 0, and obtain the temporal evolutions of particles which are initially in the lower half of µ space without direct computations, using this symmetry.

To estimate the perturbed magnetization M

x (t), we subtract from M x (t) its long time average. A typical temporal evolution of M

x (t) is shown in Fig. 2. Finite size effects, which are visible on the curve for N = 10 7 points, may prevent the study of the asymptotic behavior, so that it is usually necessary to use a very large number of points (see the curve for N = 10 8 points). We fit the envelop of the decaying M x (t) curve by a power-law, using the least square method (see figure (3)). We find an exponent -2.95, which is in very good agreement with the prediction (62) .

Sine perturbation

We consider now a sine perturbation of the thermal equilibrium:

f 0 (x, p) + f 1 (x, p) = N a e -β(p 2 /2-M 0 cos x) (1 + a sin x) . (66)
In this case, the (x, p) → (-x, -p) symmetry is broken, and we have to compute the evolution of the whole µ space. A slow rotating motion of the magnetization appears, which makes it difficult to define an asymptotic average value of M

x (t) and M

y (t). Rather, we used a running average to eliminate the rotation effect (see Appendix C).

x (t) for an initial condition as in (65). We used the weighted particles code with T = 0.1, a = 0.1, N = 10 7 (black curve) and N = 10 8 (gray curve). The initial weighted points are equally distributed in ]π, π] × [-3, 3], and we used the (x, p) → (-x, -p) symmetry in order to reduce the computations.

x (t) envelop for an initial condition perturbed by a cosine, as in (65). Red points: numerical simulation using the weighted particles code with T = 0.1, a = 0.1 and N = 10 9 . The initial weighted points are equally distributed in ]π, π] × [-3, 3]. Black dashed line: power law fitting from t = 600 to 6000 (∝ t -2.95 ).

Finally, we obtain curves similar to the case of a cosine perturbation to confirm the prediction (62) and (63).

We then fit the envelops of M

x (t) and M

y (t) with power laws. The exponents are -2.94 and -1.77, to be compared with the predicted -3 and -2, see (62) and (63). The relative error for M x (t) and M

(1) y (t), with frequency respectively ω 0 and 2ω 0 . On Fig. 5, we plot the power

x (t) (upper panel) and M

(1) y (t) (lower panel) envelops for an initial condition perturbed by a sine, as in (66). Red points: numerical simulation using the weighted particles code with T = 0.1, a = 0.1 and N = 10 9 . The initial weighted points are equally distributed in ]π, π] × [-3, 3]. Black dashed line: power law fit from t = 600 to 6000. We obtain an exponent equal to -2.94 for M (1) x (t) and -1.77 for M

spectra of M

(1)

x (t) and M

(1) y (t). We observe that indeed in both cases a single frequency is selected in the long time regime, with numerical values respectively 1.944 and 0.973. This is in almost perfect agreement with the theoretical prediction 2ω 0 = 1.945 and ω 0 = 0.972.

Figure 5 also allows to observe the cross-over between two different dynamics explored by the system: the short time evolution is driven by the Landau pole contribution, and the asymptotic behavior by (62) and (63). Indeed, as shown in [28], it is possible to compute the dominant Landau pole for the parameters of Fig. 5; one finds a frequency Re(ω L ) ≃ 1.8. This is in good agreement with the short time power spectrum of M x (t), which is maximal around ω = 1.85. We therefore conclude that Landau damping occurs in the short time regime, before the algebraic decay dominates. The peak position at short time in the M y 's spectrum may be a signature of the growing ω 0 peak in a long time regime, since it is close to ω 0 . It might also be related to a Landau pole associated to M y : since the matrix Λ -F (ω) is diagonal, Landau poles for M x and M y may be different. Finally notice that the power spectrum divergence close to ω = 0 is only due to the rotation of the magnetization.

x (t) (upper panel) and M

(1) y (t) (lower panel) for an initial condition perturbed by a sine, as in (66). The numerical simulation was done using the weighted particles code with T = 0.1, a = 0.1 and N = 10 9 . The initial weighted points are equally distributed in ]π, π] × [-3, 3]. Each inset includes a magnification around the maxima and the short time power spectrum. The maxima of the long time power spectra of M

Conclusion

We have investigated the asymptotic dynamics of perturbations around stable non homogeneous backgrounds in spatially one-dimensional Vlasov equations.

The dispersion relation of the linearized Vlasov equation has poles in the lower half of the complex plane and logarithmic branch points on the real axis. The poles yield exponentially decaying contributions: this is a form of Landau damping. The branch points yield algebraically decaying contributions. If the dominant branch point is related to a minimum of the potential created by the stationary state, the leading order of the potential is quadratic with respect to position around the minimum, and unless some special cancellation occurs, the perturbation potential behaves asymptotically as e iω 0 t /t 2 , where ω 0 is the harmonic frequency of the potential well. We expect this situation to happen in many cases of interest.

We have tested the theory on the HMF model by performing N-body simulations which correspond to the full Vlasov equation in the large N limit; these simulations used the weighted particles code. The exponent and frequency of the decay have been confirmed by these simulations, including in one case where a special symmetry imposes a decay as e i2ω 0 t /t 3 . This summarizes in the following scenario for decaying perturbations around stable non homogeneous background in one-dimensional systems: i) The perturbation potential first roughly behaves as an oscillating decaying exponential, with the frequency and decay rate related to the poles of the dispersion relation, as usual Landau damping on homogeneous backgrounds. ii) After this transient, the algebraic decay sets in, and the frequency changes; the decay exponent and the new frequency are now governed by the dominant branch point singularity on the real axis of the dispersion relation.

These results prompt several questions, which this work does not answer. First, when the stationary distribution has a compact support in action variable, the edge of the support may create a singularity stronger than the bottom of the potential well. In some cases, one would then expect a 1/t asymptotic decay, with a frequency corresponding to the action at the edge of the support. However, a stationary distribution with compact support may also sustain purely oscillatory modes [20], which do not decay at all. We have not been able to find a stationary state with compact support and no oscillatory mode to test the possibility of a 1/t decay. Second, the local extrema of the function Ω(J) may also create a different type of singularity, and thus modify the asymptotic behavior of a perturbation. Recently, Bouchet and Morita have studied a similar situation in the context of the 2D Euler equation, unveiling the phenomenon of "vorticity depletion" close to these local extrema [30]. The possibility of a similar behavior for the Vlasov equation is an open question. Third, we picked up the slowest decaying contribution among Fourier modes with respect to the angle variable, and we have not investigated if the sum over the infinite number of Fourier modes affects the asymptotic decay.

It is also of primary interest to understand the asymptotic behavior of a perturbation in a three dimensional setting. When the potential created by the stationary state is integrable, angle-action variables can be defined, and the method used in this article is viable: one would have to study the singularities of the dispersion relation in this case. In a situation where the potential created by the stationary state is not integrable, the strategy would fail.

Finally, our numerical study required to compute a Vlasov evolution with good precision for a long time, which raises difficulties. To overcome them, we have introduced a particle method, the weighted particles approach. It proved particularly well-suited for our purpose, and it would be interesting to investigate the reasons for this.

to each lattice point a weight proportional to the initial distribution we want to sample. The concrete algorithm is as follows.

We span the µ space with a regular lattice, having N points. The upper and lower boundaries ±p max for the lattice in the p-direction must be set such that the initial distribution f (x, p) for |p| > p max is negligible. Let i = 1, • • • , N denote the lattice points, and (x i , p i ) their coordinates. We assign the weight w i = Cf (x i , p i ) to each lattice point, where the constant C is defined by the normalization

N i=1 w i = 1.
We put one particle on each lattice point, and the particles move on µ space following by the canonical equations of motion for the HMF model ẋi (t) = p i (t), ṗi (t) = -Mx (t) sin x i (t) + My (t) cos x i (t), where the suffix i runs from 1 to N, and Mx (t) and My (t) are defined by

Mx (t) = N i=1 w i cos x i (t), My (t) = N i=1 w i sin x i (t). (B.1)
We note that the lattice is used only to define the weight w i , and to set the initial condition for the i-th particle as (x i (0), p i (0)) = (x i , p i ). It is worth stressing that particles are not fixed at lattice points, but move in the whole µ space, keeping their initially assigned weights w i . The magnetization (M x , M y ) is computed using the evolving positions and the initially fixed weight as (B.1). In a semi-Lagrangian code [START_REF] Buyl | [END_REF], the distribution is defined at fixed lattice points. To evolve the distribution over a time step ∆t, one computes the inverse temporal evolution of a particle during ∆t with initial condition given by the lattice (x i (0), p i (0)) = (x i , p i ). Integrating the Vlasov equation, one gets f (x i , p i , ∆t) = f (x i (-∆t), p i (-∆t), 0). The point (x i (-∆t), p i (-∆t)) does not coincide with a lattice point generally, so that the value of the distribution at this point is obtained by interpolation. The semi-Lagrangian code thus requires three discretization parameters: the time step ∆t, and two spatial and velocity discretization parameters ∆x, ∆p. In order to insure numerical stability, the time step must become small as spatial and velocity discretizations become small. As a result, we need a small time step for computations with high spatial resolution, which results in heavy computations to reach the long time regime.

In the weighted particles code, we may set the ∆x and ∆p discretizations independently of the time step ∆t, and hence the computational burden may be reduced by taking a rather large ∆t, still giving a good enough accuracy.

The weighted particles code has further advantages against the semi-Lagrangian code: (1) The Vlasov equation has an infinite number of conserved quantities which are (f (x, p)) l dxdp for l ∈ N, but it is known that the semi-Lagrangian code cannot preserve them for l ≥ 2. In the weighted particles code, the weight w i is a fixed value and hence all such quantities, approximated by i (w i ) l , are preserved exactly. (2) The semi-Lagrangian code uses interpolation. It is a delicate step to obtain the temporal evolution of the distribution function, and the accuracy of the code depends on the algorithm of interpolation. The weighted particles code does not require any interpolation.

Let us remark that the weighted particles algorithm is highly parallelizable and its convenient structure makes it possible to use a lot of parallelization methods. In particular, it allows to take advantage of the available computer architecture, be it a cluster with distributed memory or shared memory. In the case of a distributed memory, the mean field property allows to restrict the communication between node to the magnetization, which can be computed piece by piece on each node. It is thus possible to compute the long time evolution of the system for a very large number of weighted particles.

One of the disadvantages of the weighted particles codes is that this code cannot compute the temporal evolution of the distribution directly. We can obtain a coarsegrained distribution, but its resolution is lower than the initially given lattice. Another disadvantage may be the limitation of objects for which the weighted particles code works well. The weight w i on a lattice point corresponds to set several particles with the same initial condition. If these particles were given slightly different initial conditions, they would eventually separate as time goes by. In the weighted particles code, they remain together. Consequently, weighted particles code might have to be improved if it is to be used in order to observe more drastic changes of the distribution function, such as violent relaxation from a waterbag initial state to a Lynden-Bell quasi-stationary state. Further investigations are needed to understand why the weighted particles code seems to work so well in our case. The numerical tests and theoretical arguments given in [33] may be a first step in this direction.

Appendix C. Extraction of a rotating part

When we use a asymmetric perturbation such as the one given in (66), we observe a small rotation of the x and y magnetizations M x (t) and M y (t) (see Fig. C1). In this case, it is not possible to directly define M x (t) and M

(1) y (t) by substracting the long time average. Our method is then to use a running average. We define M x (t) and M

(1) y (t) such that M (1) x (t) = M x (t) - In order to compute the right exponent of a power law fit, we have to choose the parameter ∆t. However our different tests (see Fig. C2) show that modifying ∆t does not change much the result. We have taken ∆t = 5 for Fig. 4. This is roughly the time needed to observe one oscillation of M y (t), and two oscillations of M x (t). x (t) (upper panel) and M

(1) y (t) (lower panel) for an initial condition perturbed by a sine, as in (66). The numerical simulation was done using the weighted particles code with T = 0.1, a = 0.1 and N = 10

Figure 2 .

 2 Figure 2. Temporal evolution of M

Figure 3 .

 3 Figure 3. (color online) Temporal evolution of the M

( 1 )

 1 y (t) exponent is close to 10%. Different explanations are possible: the fit range does not completely lie in the asymptotic regime and/or there are numerical errors. (62) and (63) predict an oscillating decay of the perturbations M

Figure 4 .

 4 Figure 4. (color online) Temporal evolutions of the M

Figure 5 .

 5 Figure 5. (color online) Power spectra of the M

  M x (t ′ )dt ′ , (C.1)andM (1) y (t) = M y (t) -1 2∆t t+∆t t-∆t M y (t ′ )dt ′ . (C.2)

Figure C1 .

 C1 Figure C1. Slow rotation of M

  9 . The initial weighted points w i are equally distributed in ]π, π] × [-3, 3].

Figure C2 .

 C2 Figure C2. (color online) Temporal evolution of the M

( 1 )

 1 y (t) envelop for an initial condition perturbed by a sine, as in (66), using different window sizes for the running average. The numerical simulation was done using the weighted particles code with T = 0.1, a = 0.1 and N = 10 9 . The initial weighted points are equally distributed in ]π, π] × [-3, 3].

Acknowledgments

The authors thank Freddy Bouchet and Bruno Marcos for helpful discussions. This work is supported by the Ministry of Educations, Science, Sports and Culture, Grantin-Aid for Young Scientists (B), 19760052 and by the ANR-09-JCJC-009401 INTERLOP project. Computations have been partially performed on the "Mésocentre SIGAMM" machine, hosted by the Observatoire de la Côte d'Azur.

We assume that A(z) rapidly decreases for |z| → ∞. We take a branch cut in the lower half of the imaginary axis as shown in Fig. A1. To compute the integral over the Bromwich contour Γ, we add the paths Γ R , Γ B and Γ L to the Bromwich contour Γ, and make a closed path. From the assumption no singularity is enclosed by the closed path, and hence

Thanks to the factor e -izt of the integrand, contributions from Γ R and Γ R vanish, and hence

) ν ln(re -iπ/2 )e -ire -iπ/2 t e -iπ/2 dr (A.5)

Due to the branch cut, we have to distinguish A(re i3π/2 ) from A(re -iπ/2 ), and we denote them by A L (r) and A R (r) respectively. The integral is hence written as

The two functions A L (r) and A R (r) coincide in the limit r → 0. Around the singularity r = 0, we can therefore estimate

Using a scaling of the variable as y = rt, the integral is expressed by

y ν e -y dy.

(A.9)

The first term of the right-hand-side yields t -ν-2 decay and t -ν-2 ln t decay, and the second term t -ν-1 decay. Returning to (A.2), we obtain e -imωat /t ν+1 as the slowest decay. From (A.1), we now see that the asymptotic decay of a k (t) is governed by the smallest ν m .

Appendix B. Weighted particles code

The N-body simulations are performed by the weighted particles code. In a standard N-body simulations, we would prepare N initial positions and momenta by drawing random numbers according to a given initial distribution. An initial condition prepared in this way has fluctuations of order O(1/ √ N ). In the weighted particles code, we prepare N initial positions and momenta as lattice points of a square lattice, and give