Eric Cancès
email: eric.cances@enpc.fr

Gaspard Kemlin
email: gaspard.kemlin@u-picardie.fr

Antoine Levitt
email: antoine.levitt@universite-paris-saclay.fr

A PRIORI ERROR ANALYSIS OF LINEAR AND NONLINEAR PERIODIC SCHRÖDINGER EQUATIONS WITH ANALYTIC POTENTIALS

This paper is concerned with the numerical analysis of linear and nonlinear Schrödinger equations with periodic analytic potentials. We prove that, for linear equations, when the potential is analytic in a strip of width A of the complex plane, the solution is analytic in the same strip, ensuring an exponential convergence of the planewave discretization of the equation with rate A. On the other hand, for nonlinear equations, we find that the solution may be analytic only in a strip of width smaller than A. This behavior is illustrated by two examples using a combination of numerical and analytical arguments.

Introduction

In this paper, we study models inspired by Kohn-Sham Density Functional Theory (KS-DFT). KS-DFT is currently the most popular model in quantum chemistry and materials science as it offers a good compromise between accuracy and computational efficiency. It aims at computing, for a given configuration of the nuclei of the molecular system or material of interest, the electronic ground-state energy and density. From the latter, it is possible to compute the effective forces acting on the nuclei in this configuration, and thus to identify the (meta)stable equilibrium configurations of the system, or to simulate the dynamics of the molecular system in various thermodynamic conditions. In materials science applications, computations are commonly done in a periodic simulation cell, which can be either the unit cell of a crystal (for the special case of perfect crystals), or a supercell (for all the other cases: crystals with defects, disordered alloys, glassy materials, liquids...).

We denote by L = Za 1 + Za 2 + Za 3 the periodic lattice, where (a 1 , a 2 , a 3 ) is a non-necessarily orthonormal basis of R 3 , and by Ω = [0, 1)a 1 + [0, 1)a 2 + [0, 1)a 3 the simulation cell. Let us denote by L 2 per,L := {u ∈ L 2 loc (R 3 , C) | u is L-periodic} the Hilbert space of complex-valued L-periodic locally square integrable functions on R 3 , endowed with its usual inner product. The KS-DFT equations read

H ρ ϕ i = λ i ϕ i , ϕ i , ϕ j L 2 per,L = δ ij , ρ(x) = 2 Np i=1 |ϕ i (x)| 2 , (1)
where H ρ is the Kohn-Sham Hamiltonian, a self-adjoint operator on L 2 per,L bounded below and with compact resolvent. The ϕ i 's are the Kohn-Sham orbitals, and the λ i 's their energies. Since H ρ depends on ρ, which in turn depends on the eigenfunctions ϕ i , 1 i N p , (1) is a nonlinear eigenproblem. The parameter N p represents physically the number of valence electron pairs per simulation cell and ρ the ground-state electronic density. We assume here, as is the case for most physical systems, that λ 1 λ 2 • • • λ Np are the lowest N p eigenvalues of H ρ (Aufbau principle). The Kohn-Sham Hamiltonian with pseudopotentials reads

H ρ = - 1 2 ∆ + V nl + V Hxc loc,ρ
where V nl is a finite-rank self-adjoint operator (the nonlocal part of the pseudopotential), and V Hxc loc,ρ (x) = V loc (x) + V H,ρ (x) + V xc,ρ (x) is a periodic real-valued function depending (nonlocally) on ρ. The function V loc is the local component of the pseudopotential, the Hartree potential V H,ρ is the unique solution with zero mean to the periodic Poisson equation

-∆V H,ρ (x) = 4π ρ(x) - 1 |Ω| ˆΩ ρ , ˆΩ V H,ρ = 0,
and the function V xc,ρ , called the exchange-correlation potential, depends on the chosen approximation of the exchange-correlation energy functional. In the simple Xα model [START_REF] Slater | A Simplification of the Hartree-Fock Method[END_REF], V xc,ρ (x) = -C D ρ(x) 1/3 , where C D > 0 is the Dirac constant.

It is not mandatory to use pseudopotentials in KS-DFT calculations. Some softwares allow for all-electrons calculations in which the total pseudopotential operator V nl +V loc is replaced with a local potential with Coulomb singularities at the positions of the nuclei. However, most calculations are done with pseudopotentials, or use the formally similar Projector Augmented Wave (PAW) method [START_REF] Blöchl | Projector augmented-wave method[END_REF], for three reasons: (i) core electrons are barely affected by the chemical environment and can usually be considered to occupy "frozen states", (ii) in heavy atoms, core electrons must be dealt with relativistic quantum models which makes the simulation more expensive from a computational viewpoint, (iii) due to the Coulomb singularities, all-electron Kohn-Sham orbitals have cusps at the positions of the nuclei and are therefore only Lipschitz continuous, while the Kohn-Sham orbitals computed with pseudopotentials are much more regular and can be well approximated with Fourier spectral methods (usually called planewave discretization methods in the field).

Several methods for constructing pseudopotentials have been proposed in the literature, leading to local and nonlocal functions of different regularity. As expected, the rate of convergence of the planewave discretization method is directly linked to the regularity of these functions. The a priori error analysis of this problem was performed in [START_REF] Cancès | Numerical analysis of the planewave discretization of some orbital-free and Kohn-Sham models[END_REF] for pseudopotential with Sobolev regularity. It was proved in particular, for the simple Xα exchange-correlation functional, but also for the much more popular local density approximation (LDA) exchange-correlation functional, that if the local and nonlocal parts of the pseudopotential are in the periodic Sobolev space of order s > 3/2, then the Kohn-Sham orbitals ϕ i and the density ρ are in the periodic Sobolev space of order s + 2, and (optimal) polynomial convergence rates were obtained in any Sobolev spaces of order r with -s < r < s + 2. In addition, as for linear second-order elliptic eigenproblems, the error on the eigenvalues converges to zero as the square of the error on the eigenfunctions evaluated in H 1 -norm. The analysis in [START_REF] Cancès | Numerical analysis of the planewave discretization of some orbital-free and Kohn-Sham models[END_REF] covers for example the case of Troullier-Martins pseudopotentials [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF], for which s = 7 2ε for all ε > 0. On the other hand, these estimates are not sharp in the case of Goedecker-Teter-Hutter (GTH) pseudopotentials [START_REF] Goedecker | Separable dual-space Gaussian pseudopotentials[END_REF][START_REF] Hartwigsen | Relativistic separable dual-space Gaussian pseudopotentials from H to Rn[END_REF], for which the local and nonlocal contributions are periodic sums of Gaussian-polynomial functions, and therefore have entire continuations to the whole complex plane. Such pseudopotentials are implemented in different DFT softwares, such as BigDFT [START_REF] Ratcliff | Flexibilities of wavelets as a computational basis set for large-scale electronic structure calculations[END_REF], Quantum Espresso [START_REF] Giannozzi | QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials[END_REF] or Abinit [START_REF] Gonze | The Abinit project: Impact, environment and recent developments[END_REF][START_REF] Romero | ABINIT: Overview and focus on selected capabilities[END_REF], as well as DFTK, a recent electronic structure package in the Julia language [START_REF] Herbst | DFTK: A Julian approach for simulating electrons in solids[END_REF].

The purpose of this paper is to investigate this case. While it has been known for a long time (see e.g. [START_REF] Bernstein | Sur la nature analytique des solutions des équations aux dérivées partielles du second ordre[END_REF][START_REF] Friedman | On the Regularity of the Solutions of Non-Linear Elliptic and Parabolic Systems of Partial Differential Equations[END_REF][START_REF] Morrey | On the Analyticity of the Solutions of Analytic Non-Linear Elliptic Systems of Partial Differential Equations: Part I. Analyticity in the Interior[END_REF][START_REF] Morrey | On the Analyticity of the Solutions of Analytic Non-Linear Elliptic Systems of Partial Differential Equations: Part II. Analyticity at the Boundary[END_REF][START_REF] Petrovskii | Sur l'analyticité des solutions des systèmes d'équations différentielles[END_REF] and references therein for historical insight or [START_REF] Blatt | On the analyticity of solutions to non-linear elliptic partial differential systems[END_REF][START_REF] Hashimoto | A Remark on the Analyticity of the Solutions for Non-Linear Elliptic Partial Differential Equations[END_REF] for more recent developments) that the solutions to elliptic equations on R d with real-analytic data have an analytic continuation in a complex neighborhood of R d , the size of this neighborhood is a priori unknown. In the periodic case we are considering, the latter directly impact the decay rate of the Fourier coefficients of the solution, hence the convergence rate of the planewave discretization method. For pedagogical reasons, we will work most of the time with one dimensional linear or nonlinear Schrödinger equations, because (i) it is easier to visualize analytic or entire continuations of functions originally defined on the real space R d when d = 1, and (ii) exponential convergence rates of planewave discretization methods are easier to spot in 1D. However, most of our arguments extend to the multidimensional case. In Section 2, we introduce a hierarchy of spaces (H A ) A>0 of complex-valued 2πperiodic functions on the real line having analytic continuations to the strip R + i(-A, A) of the complex plane. We then pick a real-valued function V ∈ H B for some B > 0 and consider the one-dimensional Schrödinger operator H = -∆ + V . A low vs high-frequency decomposition of the periodic L 2 space allows to prove that for all 0 < A < B, the solution u to the linear equation Hu = f is in H A whenever f ∈ H A (see Section 3.1), and that the eigenfunctions of H are in H A (see Section 3.2). We rely on this result to prove in Section 3.3 that the planewave discretization method converges exponentially fast in this case and we provide a numerical illustration of these results. We turn in Section 4 to the nonlinear setting, where we present two examples for which we show, using a combination of numerical and analytical tools, that the analyticity strip of the solution may be much smaller than the one of the data. Finally, we consider in the Appendix the multidimensional case, which is an immediate extension, and its application to Kohn-Sham models.

Spaces of analytic functions

Let us first introduce some notation. We denote by L 2 per (R, C) the space of locally square integrable complexvalued 2π-periodic functions on R, endowed with its natural inner product

(u, v) L 2 per := ˆ2π 0 u(x) v(x) dx,
and by S per (R, C) the space of tempered complex-valued 2π-periodic distributions on R. For each u ∈ S per (R, C), we denote by ( u k ) k∈Z the Fourier coefficients of u with the following normalization convention:

∀ u ∈ L 2 per (R, C), ∀ k ∈ Z, u k := (e k , u) L 2 per = 1 √ 2π ˆ2π 0 u(x)e -ikx dx,
where e k (x) := 1 √ 2π e ikx is the L 2 per -normalized Fourier mode with wavevector k ∈ Z. Recall that the 2π-periodic Sobolev spaces are the Hilbert spaces H s per (R, C), s ∈ R, defined by

H s per (R, C) := u ∈ L 2 per (R, C) k∈Z (1 + |k| 2 ) s | u k | 2 < ∞ , (u, v) H s per := k∈Z (1 + |k| 2 ) s u k v k .
We will also use the self-explanatory notation

C k per (R, R), C k per (R, C), L p per (R, R), L p per (R, C
) for k ∈ N ∪ {∞} and 1 p ∞, all these spaces being endowed with their natural norms or topologies. We now introduce, for any A > 0, the (Hardy-like) space

H A := u ∈ L 2 per (R, C) k∈Z w A (k) | u k | 2 < ∞ where w A (k) := cosh(2Ak) = 1 2 (e 2Ak + e -2Ak
), endowed with the inner product

(u, v) A := k∈Z w A (k) u k v k .
Note that H A can be canonically identified with the space of analytic functions

H A :=    u : Ω A → C analytic [-A, A] y → u(• + iy) ∈ L 2 per (R, C) continuous, ˆ2π 0 |u(x + iA)| 2 + |u(x -iA)| 2 dx < ∞    ,
where

Ω A := R + i(-A, A) ⊂ C
is the horizontal strip of width 2A of the complex plane centered on the real axis, endowed with the inner product

(u, v) H A = 1 2 (u(• + iA), v(• + iA)) L 2 per + (u(• -iA), v(• -iA)) L 2 per .
The canonical unitary mapping H A onto H A is the analytic continuation: any function u ∈ H A has a unique analytic continuation u : Ω A → C given by

∀ z = x + iy ∈ Ω A , u(z) = k∈Z u k e ikz √ 2π = k∈Z u k e -ky e k (x).
It can be easily seen that the Fourier coefficients of u(• ± iA) are the Fourier coefficients of u rescaled by a factor e ∓kA and that the function

(-A, A) y → u(• + iy) = k∈Z u k e -ky e k (•) ∈ L 2 per (R, C) has a unique continuation to [-A, A] . Therefore, u 2 H A = 1 2 u(• + iA) 2 L 2 per + u(• -iA) 2 L 2 per = 1 2 k∈Z u k e -kA 2 + k∈Z u k e +kA 2 = k∈Z w A (k) | u k | 2 = u 2 A .
We record for future use the Proposition 1. Let B > 0. Then, for all 0 < A < B, the multiplication by a function V B ∈ H B defines a bounded operator on H A .

Proof. The proof is immediate from the analyticity of

V in Ω B , which implies that sup z,|Im(z)| A |V (z)| < ∞.
3. The linear case 3.1. The linear elliptic problem. We consider in a first stage the one-dimensional linear elliptic problem

seek u ∈ H 2 per (R, C) such that -∆u + V u = f on R, (2)
where V ∈ L 2 per (R, R) and f ∈ L 2 per (R, C) are given 2π-periodic functions. We assume in this section that V 1. It is then standard that (2) has a unique solution u satisfying the a priori bounds

u L 2 per f L 2 per α and u H 1 per f H -1 per , (3)
where

α := λ 1 (-∆ + V ) 1 is the smallest eigenvalue of the self-adjoint operator H = -∆ + V on L 2 per (R, C). By bootstrap arguments, u ∈ H s+2
per (R, C) whenever V and f are in H s per , for any s 0. The following result deals with the case of real-analytic potentials V and right-hand sides f . Theorem 1. Let B > 0 and V ∈ H B be real-valued and such that V 1 on R. Then, for all 0 < A < B and f ∈ H A , the unique solution u of (2) is in H A . Moreover, we have the following estimate

∃ C > 0 independent of f such that u A C f A . ( 4
)
As a consequence, if V and f are entire, then so is u.

Proof. For N > 0, we consider the decomposition

L 2 per (R, C) = X N ⊕ X ⊥ N where X N := Span(e k , |k| N ) = u ∈ L 2 per (R, C) | u k = 0, ∀ |k| > N . (5)
Let Π N be the orthogonal projector on X N and Π ⊥ N := 1 -Π N the orthogonal projector on X ⊥ N . Note that the restriction of Π N to the Sobolev space H s per (R, C), s > 0, is also the orthogonal projector on X N for the H s per inner product, and that the same property holds for the Hilbert spaces H A .

For a fixed N , we decompose u as

u = u 1 + u 2 with u 1 ∈ X N and u 2 ∈ X ⊥ N .
As u 1 has compact Fourier support, it obviously belongs to H A and we have the estimate

u 1 A u L 2 per w A (N ) f L 2 per α w A (N ). (6)
Let us show that, for N large enough, u 2 also belongs to

H A . Projecting -∆u + V u = f onto X ⊥ N , we get T 22 u 2 + V 22 u 2 = f 2 -V 21 u 1 , (7)
where T 22 is the restriction to the invariant subspace H A ∩ X ⊥ N of the self-adjoint operator -∆ on L 2 per (R, C), and where, in view of Proposition 1,

V 22 := Π ⊥ N V Π ⊥ N ∈ L(H A ∩ X ⊥ N ), V 21 := Π ⊥ N V Π N ∈ L(H A ∩ X N , H A ∩ X ⊥ N ) and f 2 := Π ⊥ N f ∈ H A ∩ X ⊥ N . The operator T 22 on H A ∩ X ⊥ N is bounded from below by N 2 and is therefore invertible with inverse T -1 22 bounded by N -2 . As f 2 , V 21 u 1 ∈ H A ∩ X ⊥ N ,
we can therefore rewrite [START_REF] Cancès | Numerical analysis of the planewave discretization of some orbital-free and Kohn-Sham models[END_REF] as

(1 + T -1 22 V 22 )u 2 = T -1 22 (f 2 -V 21 u 1 ). (8)
Since

T -1 22 L(H A ∩X ⊥ N ) N -2 and V 22 L(H A ∩X ⊥ N )
V L(H A ) , we can choose N large enough so that the

operator (1 + T -1 22 V 22 ) is invertible on H A ∩ X ⊥ N . It then holds u 2 = (1 + T -1 22 V 22 ) -1 T -1 22 (f 2 -V 21 u 1 ) (9)
and the result follows from

V 21 L(H A ) V L(H A )
and the bound (6) on u 1 .

Remark 1 (Additional regularity). We mention here some cases where we can obtain additional regularity on the solution u to (2). First, note that, if we additionally require that

V (• ± iB) ∈ L ∞ per (R, C), then the same reasoning leads to u ∈ H B if f ∈ H B . Next, one might wonder whether there exist 0 < B < B 0 < A, f ∈ H A and V ∈ H B such that V / ∈ H C for all B 0 C A but nevertheless u ∈ H A .
The following example shows that this can happen. Consider the 2π-periodic real-valued functions V and u 0 on R defined by

∀ x ∈ R, V (x) = 1 + 1 γ + sin 2 (x)
and u 0 (x) = γ + sin 2 (x) 2γ + sin 2 (x) .

We have that (i) V ∈ H B if and only if B < B 0 := arcsinh( √ γ) and (ii) u 0 and V u 0 both belong to H A for any A < A 0 := arcsinh( √ 2γ) (in particular for any A such that B 0 A < A 0 ). Therefore, if we set f := -∆u 0 + V u 0 , then f ∈ H A for any A such that B 0 < A < A 0 and the unique solution to (2) with data V and

f is u = u 0 ∈ H A despite V / ∈ H C for all B 0 C A.

3.2.

The linear eigenvalue problem. We now focus on the linear eigenvalue problem,

-∆u + V u = λu, u L 2 per (R,C) = 1, (10)
where V ∈ H B for some B > 0. Using the same technique as for the proof of Theorem 1, we get the Theorem 2. Let B > 0, V ∈ H B be real-valued, and (u, λ) ∈ H 2 per (R, C) × R be a normalized eigenmode of H = -∆ + V , with isolated eigenvalue (i.e. a solution to (10)). Then, u is in H A for all 0 < A < B. As a consequence, if V is entire, then so is u.

Proof. The proof proceeds along the same lines as that of Theorem 1. Let (u, λ) ∈ H 2 per (R, C) × R be a solution to [START_REF] Goedecker | Separable dual-space Gaussian pseudopotentials[END_REF]. Using the same notation as in the proof of Theorem 1, we decompose u as u = u 1 + u 2 with u 1 ∈ X N and u 2 ∈ X ⊥ N , and observe that for N large enough,

u 2 = -(1 + T -1 22 (V 22 -λ)) -1 T -1 22 V 21 u 1 , with T -1 22 (V 22 -λ) L(H A ∩X ⊥ N ) V L(H A ) + |λ| N 2 .
Therefore, choosing N large enough, we have

(1 + T -1 22 (V 22 -λ)) -1 T -1 22 L(H A ∩X ⊥ N ) < 1
, from which we deduce that u 2 ∈ H A and the result follows.

3.3.

Planewave approximation of the linear Schrödinger equation. Using X N = Span(e k , |k| N ) ⊂ H 1 per (R, C) as a variational approximation space for [START_REF] Goedecker | Separable dual-space Gaussian pseudopotentials[END_REF], we obtain the finite-dimensional problem

     seek (u N , λ N ) ∈ X N × R such that u N L 2 per (R,C) = 1 and ∀ v N ∈ X N , ˆ2π 0 ∇u N • ∇v N + ˆ2π 0 V u N v N = λ N ˆ2π 0 u N v N , (11)
which is equivalent to seeking the eigenpairs of the Hermitian matrix H N with entries

[H N ] kk := |k| 2 δ kk + V k-k , k, k ∈ Z, |k| N, |k | N.
The following theorem states that if V ∈ H B for some B > 0, the planewave discretization method has an exponential convergence rate. Note that a similar result holds for the planewave approximation of the linear problem -∆u + V u = f , whenever f ∈ H A .

Theorem 3. Let B > 0, V ∈ H B be real-valued, i ∈ N * and 0 < A < B. Let λ i be the i th lowest eigenvalue of the self-adjoint operator H = -∆ + V on L 2 per (R, C) counting multiplicities, and E i = Ker(Hλ i ) the corresponding eigenspace. For N large enough, we denote by λ i,N the i th lowest eigenvalue of [START_REF] Gonze | The Abinit project: Impact, environment and recent developments[END_REF], and by u i,N an associated normalized eigenvector. Then, there is a constant C i,A ∈ R + such that, for N large enough, 2AN ) . Proof. First, note that -∆ + V has compact resolvent so that its eigenvalues λ i are isolated. Let 0 < A < B and A = A+B 2 . We have

d H 1 per (u i,N , E i ) C i,A exp (-AN ) and 0 λ i,N -λ i C i,A exp (-
∀ v ∈ H A , v -Π N v H 1 per C A,B v A e -AN with C A,B := max k∈Z (1 + |k| 2 )e 2A|k| w A (k) 1/2 < ∞.
The operator H is self-adjoint on L 2 per (R, C) with form domain H 1 per (R, C) and we infer from Theorem 2 that all the eigenfunctions of the operator H are in H A . Theorem 3 then follows from classical arguments on the variational approximations of the eigenmodes of bounded below self-adjoint operators with compact resolvent (see e.g. [1, Theorems 8.1 and 8.2]).

3.4.

Numerical results in the linear case.

3.4.1. Computational framework. In this paper, all the numerical tests are realized with the DFTK software [START_REF] Herbst | DFTK: A Julian approach for simulating electrons in solids[END_REF]. This Julia package uses a planewave basis X N , as defined in Section 3.3, through a discretization parameter E cut := N 2 /2. Then, the numerical strategy depends on the nature of the problem:

• linear eigenproblems are solved with a LOBPCG solver (see e.g. [START_REF] Nottoli | A robust, open-source implementation of the locally optimal block preconditioned conjugate gradient for large eigenvalue problems in quantum chemistry[END_REF]);

• nonlinear eigenproblems (see Section 4.2.1) and nonlinear elliptic problems (see Section 4.2.2) with source terms are solved by direct minimization of the corresponding energy, through the Julia implementation of the LBFGS minimization algorithm [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF].

3.4.2.

Results in the linear case. In this section, we provide some numerical experiments that illustrate Theorems 2 and 3. To this end, we consider the potential V defined by

∀ x ∈ [0, 2π], V (x) = 1 γ + sin 2 (x)
, with γ = 1/500. The analytic continuation of V has branching points at πZ ± iB 0 with

B 0 = arcsinh( √ γ) ≈ 0.0447,
and a direct calculation shows that V ∈ H B for any B < B 0 . From the previous results, we therefore expect the eigenfunctions of -∆ + V to belong to H B for any B < B 0 and the planewave approximation to converge with rate proportional to exp(-2B 0 N ) for the eigenvalues and proportional to exp(-B 0 N ) for the eigenfunctions in the H 1 norm.

In Figure 1, we compute numerically the ground state u 1 of -∆ + V with discretization parameter E cut = 5•10 5 . In order to be able to properly spot the convergence of the Fourier coefficients, we use quadruple precision and the tolerance of the linear solver is set to 10 -20 . As expected, the Fourier coefficients of (the numerical approximation of) u 1 decrease with rate of order w B0 (k), which is confirmed by looking at the convergence of two successive nonzero Fourier coefficients, since

lim k→∞ log cosh (2B 0 k) cosh (2B 0 (k + 1)) = -B 0 .
In Figure 2, we solved the same problem for various N 's, and computed the error on the eigenvalues and the eigenvectors in the H 1 norm with respect to a reference solution obtained with N = √ 2E cut = 1000. This time, the tolerance of the linear solver is set to 10 -13 and we use standard double precision. The convergence rates predicted by Theorem 3 are observed, as expected. We can also remark that the error on the eigenvalues seems to decrease slightly faster (with an additional 1/N 2 factor), the asymptotic rate being still given by the exponential factor exp(-2B 0 N ).

The nonlinear case

In the nonlinear case, the analyticity strip of the solution can be much smaller than the one of the data. First, we report numerical simulations on a nonlinear periodic Gross-Pitaevskii eigenvalue problem with a cubic nonlinearity and a potential admitting an entire continuation, illustrating the fact that the eigenfunction belongs to H A for some A > 0, but does not seem to be entire. Next, we study a nonlinear Schrödinger equation with entire potential and source term and cubic nonlinearity. We show that the solution is not entire and we provide an upper bound of the width of the horizontal analyticity strip of the solution.

4.1. Analyticity on a strip. We consider here the nonlinear eigenvalue problem

-∆u + V u + |u| 2 u = λu, u L 2 per (R,C) = 1, (12)
where V ∈ H B for any B > 0 is real-valued. Similarly to the linear case, we consider the planewave variational approximation: using again

X N = Span(e k , |k| N ), we obtain the finite-dimensional problem      seek (u N , λ N ) ∈ X N × R such that u N L 2 per (R,C) = 1 and ∀ v N ∈ X N , ˆ2π 0 ∇u N • ∇v N + ˆ2π 0 V u N v N + ˆ2π 0 |u N | 2 u N v N = λ N ˆ2π 0 u N v N . ( 13
)
Then, the following theorem holds, which is a simple corollary of known results on the analyticity of elliptic partial differential equations. Note that this result only yields the existence of a finite strip of analyticity in the complex plane. However, even if the data is entire, there is a priori no reason for the solution to be entire too, as shown by the counter-examples that follow.

Theorem 4. Let V real-valued be in H B for any B > 0 and (u, λ) ∈ H 2 per (R, R) × R be the ground-state of (12), uniquely defined under the assumption that u > 0. Then, there exists

A > 0 such that u ∈ H A . Moreover, if (u N , λ N ) is the variational approximation of (u, λ) in X N × R, then ∃ C A > 0, u -u N H 1 per C A exp(-AN ).
Proof. Existence and uniqueness of u is a classical result, proved for instance in [START_REF] Cancès | Numerical Analysis of Nonlinear Eigenvalue Problems[END_REF]Appendix]. In addition, as V ∈ H s per (R, R) for any s > 0, we know from [START_REF] Cancès | Numerical Analysis of Nonlinear Eigenvalue Problems[END_REF]Theorem 2] and by an immediate bootstrap argument that u also belongs to

H s per (R, R) for any s > 0. Therefore, u ∈ C ∞ per (R, R).
Knowing that u ∈ C ∞ per (R, R) yields that u is real-analytic on a neighborhood of every point x 0 of the periodicity cell [0, 2π], according to well-known results on the analyticity of solutions to nonlinear elliptic equations, which we recall here for the reader's convenience (see for instance [4, Theorem 1.1], [START_REF] Hashimoto | A Remark on the Analyticity of the Solutions for Non-Linear Elliptic Partial Differential Equations[END_REF]Theorem 1] or [START_REF] Friedman | On the Regularity of the Solutions of Non-Linear Elliptic and Parabolic Systems of Partial Differential Equations[END_REF][START_REF] Morrey | On the Analyticity of the Solutions of Analytic Non-Linear Elliptic Systems of Partial Differential Equations: Part I. Analyticity in the Interior[END_REF][START_REF] Morrey | On the Analyticity of the Solutions of Analytic Non-Linear Elliptic Systems of Partial Differential Equations: Part II. Analyticity at the Boundary[END_REF]).

Theorem 5. Let Ω be an open subset of R n , D an open subset of R n × R × R n × R n×n sym , F : D (x, v, g, M ) → F (x, v, g, M ) ∈ R
a real analytic function, and u ∈ C ∞ (Ω) a real-valued function such that for all x ∈ Ω, (x, u(x), ∇u(x), ∇ 2 u(x)) ∈ D. Assume that u solves F (x, u(x), ∇u(x), ∇ 2 u(x)) = 0 and that this equation is elliptic in the sense that

∀(x, ξ) ∈ Ω × R n , ξ 0, ξ T ∇ M F (x, u(x), ∇u(x), ∇ 2 u(x))ξ = 0.
Then u is real-analytic in Ω.

By a compactness argument, we therefore have that u is analytic on a strip of size A > 0 around the real axis, for some A > 0.

We now consider the variational approximation (u N , λ N ) of (u, λ) in X N × R. We know from [6, Theorem 1] that there exists

C > 0 such that u -u N H 1 per C u -Π N u H 1 per .
The result then follows similarly to the proof of Theorem 3.

Remark 2. In order to establish the convergence rate of the eigenvalues of the ground-state of ( 12), one can follow the proof of [START_REF] Cancès | Numerical Analysis of Nonlinear Eigenvalue Problems[END_REF]Theorem 2]. Using in particular equations ( 50) and ( 53) from this reference, with the negative Sobolev norms replaced by the dual norms of H A , one can obtain that there exists

C A > 0 such that |λ -λ N | C A exp(-2AN ).

4.2.

Counter-examples. We analyze in this section two counter-examples of Theorems 1 and 2 based on the nonlinear Gross-Pitaevskii equation with entire potentials and source terms, for which the solutions are not entire. We also present numerical tests to support our analysis.

4.2.1.

A nonlinear Gross-Pitaevskii eigenvalue problem. We study in this section the following nonlinear eigenvalue problem:

-

ε∆u ε + γ cos(x)u ε + |u ε | 2 u ε = λ ε u ε in H 1 per (R, C), u ε L 2 per = 1, u ε > 0, (14)
with entire data. In view of the results we obtained in the linear case, we could expect the solution u ε to ( 14) (which we know to exist) to be entire. However, our numerical experiments suggest on the contrary that u ε is analytic only on a band of finite size.

The potential V being real-valued, u ε is actually real-valued too and the singular limit ε → 0 yields the algebraic equation γ cos(x)u 0 (x) + u 0 (x) 3 = λ 0 u 0 (x). If we assume that u 0 does not vanish on the unit cell [0, 2π], we can divide by u 0 (x) and then integrate over [0, 2π]: using the normalization condition, we obtain λ 0 = 1/2π. This yields

u 0 (x) = 1 2π -γ cos(x), (15)
which is indeed bounded away from zero as soon as β 0 := 1/(2πγ) > 1, which we take in the following. The analytic continuation of u 0 , still denoted by u 0 , satisfies [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF] being the continuation of the square root with branch cut R -. The maximal horizontal strip of the complex plane on which u 0 is analytic is R + i(-B 0 , B 0 ), where B 0 is such that cos(±iB 0 ) = β 0 , that is to say B 0 = log β 0 + β 2 0 -1 . More precisely, the function z → u 0 (z) admits two branching points at z ± = ±iB 0 , for which u 0 (z ± ) = 0 and

γ cos(z)u 0 (z) + u 0 (z) 3 = 1 2π u 0 (z), with √ • in
du 0 dt (tz ± ) = γ sin(tz ± ) 2 1 2π -γ cos(tz ± ) → ∞ when R t → 1 -, see Figure 3.
When ε > 0, we can approximate numerically the solution of ( 14) using the planewave approximation presented before, with discretization parameter E cut = 10 6 and solver tolerance 10 -13 . The Fourier coefficients of the numerical solution for various ε are displayed in Figure 4. Although it is not possible to rule out a transition to a different asymptotic regime in the limit of extremely small values of ε or very large values of k because of finite numerical accuracy, our numerical results are compatible with the decay of the Fourier coefficients of u ε as 1/(|k| 3/2 cosh(2B 0 k)). This leads to the conclusion that the solution u ε is analytic on a band strip of size close to B 0 , which is in contradiction with the results from the case. In the next section, we study a similar nonlinear elliptic equation where we reach similar conclusions and where, in addition, we can provide an upper bound of the width of the horizontal strip on which the solutions are analytic.

4.2.2.

A nonlinear elliptic equation. Still in the perspective of studying nonlinear elliptic problems with analytic data, we now consider the nonlinear periodic elliptic equation with cubic nonlinearity

-ε∆u ε + u ε + |u ε | 2 u ε = f in H 1 per (R, C), (16)
where ε > 0, and f : R → R is a real-analytic 2π-periodic function admitting an entire continuation, still denoted by f , to the complex plane. We will show that, in this particular case, the same kind of results we showed for the linear case are not true any more and we provide an estimation of the width of the horizontal analyticity strip of the solution, which is finite even though the source term f is entire.

The right-hand-side f being real-valued, u ε is also real-valued and the singular limit ε = 0 gives rise to the algebraic equation u 0 (x) + u 0 (x) 3 = f (x), which has a unique real solution for each x ∈ R. The latter can be computed by Cardano's formula: the discriminant of the cubic equation is

R(x) = -4 + 27f 2 (x) < 0, so that u 0 (x) = 3 1 2 f (x) + -R(x) 27 + 3 1 2 f (x) - -R(x) 27 , (17)
the other two roots being complex conjugates with nonzero imaginary parts. The function u 0 can be analytically continued from the real axis upwards in the complex plane as long as R(z) does not touch zero; it has branching points (with an exploding first derivative) at the points where R(z) = 0. In the rest of this section, we will use the function f (x) = µ sin(x), for a fixed µ > 0. The analytic continuation of u 0 has branching points at 2πZ ± iB 0 , with

B 0 = arcsinh 4 27µ 2 > 0.
In particular, although the source term f has an entire continuation, the solution u 0 given by ( 17) does not.

When ε > 0, we can approximate numerically the solution to [START_REF] Morrey | On the Analyticity of the Solutions of Analytic Non-Linear Elliptic Systems of Partial Differential Equations: Part I. Analyticity in the Interior[END_REF] with the planewave approximation introduced before, with discretization parameter E cut = 10 6 and solver tolerance 10 -13 . The Fourier coefficients of the numerical solution are then displayed in Figure 6 for various ε: the plots suggest that increasing ε increases the width B ε of the horizontal analyticity strip of u ε (z), but does not make it entire. we see that ϕ ε satisfies the following second-order ODE:

εϕ ε (y) + ϕ ε (y) + ϕ 3 ε (y) = iµ sinh(y), ϕ ε (0) = u ε (0) = 0, ϕ ε (0) = iu ε (0),
with u ε (0) ∈ R. Decomposing ϕ ε in its real part θ ε and imaginary part ψ ε , we see that θ ε satisfies the equation

εθ ε + θ ε + θ 3 ε -3θ ε ψ 2 ε = 0, θ ε (0) = 0, θ ε (0) = 0,
and therefore vanishes. As a consequence, u ε is purely imaginary along the imaginary axis (as could have been anticipated by the symmetries u ε (-z) =u(z), u ε (z) = u ε (z)). We are then left with studying the imaginary part ψ ε , which satisfies the ODE

εψ ε + ψ ε -ψ 3 ε = µ sinh, ψ ε (0) = 0, ψ ε (0) = u ε (0). (18)
If we can prove that ψ ε becomes non-analytic at a finite 0 < Y ε < ∞, this will imply that the width B ε of the horizontal analyticity strip of u ε satisfies B ε Y ε < ∞ and is therefore finite even though the data in ( 16) is entire. To this end, we prove in the Appendix B the following lemma, which gives a sufficient condition for the non-analyticity of ψ ε in finite time.

Lemma 1. Let Y ε ∈ R + ∪ {∞} be the maximum time of definition of ψ ε . If, for some η > 0 small enough, there exists y η B 0 such that ψ ε (y

η ) = 1 + η > 1 and ψ ε (y η ) > 0, then Y ε < ∞ and ψ ε (y) → ∞ when y → Y ε .
To prove that the sufficient condition from Lemma 1 is satisfied, we introduce the set

X µ = (y, υ) ∈ R 2 , µ sinh(y) -υ + υ 3 0 .
This set is such that, if (y, ψ ε (y)) lies strictly inside X µ , then ψ ε is locally strictly convex. For y < B 0 , ψ ε might oscillate on both sides of the boundary of X µ [START_REF] Vrabel | Duffing-Type Oscillator with a Bounded from above Potential in the Presence of Saddle-Center Bifurcation and Singular Perturbation: Frequency Control[END_REF]. We investigated numerically the behavior of this function for the set of parameters used in Figure 7 (ε = 0.1, µ = 0.5), and observed that 0 < ψ ε (B 0 ) < 1 √ 3

and ψ ε (B 0 ) > 0. This numerical observation can be trusted as the ODE satisfied by ψ ε on the interval [0, B 0 ] for ε = 0.1 and µ = 0.5 is not stiff. It is therefore easy to solve it numerically with high accuracy with a posteriori error estimates guaranteeing that ψ ε (B 0 ) is indeed strictly between 0 and 1 √ 3 , and ψ ε (B 0 ) is positive. Therefore, ψ ε is locally strictly convex at y = B 0 . Given the shape of X µ (see Figure 7 in appendix), ψ ε lies, for any B 0 < y < Y ε , above its tangent at y = B 0 , whose slope is ψ ε (B 0 ) > 0. Strict convexity allows one to conclude that for any η > 0, there exists y η > B 0 such that

ψ ε (y η ) = 1 + η > 1, ψ ε (y η ) > ψ ε (B 0 ) > 0.
Any η > 0 is therefore suitable to apply Lemma 1 and conclude that ψ ε blows up in finite time Y ε . We can deduce from these investigations that u ε is analytic only on a horizontal strip of finite width of the complex plane although the source term f is an entire function: our results in the linear case are therefore no longer valid in general in the nonlinear case. Then, for each 0 < A < B and i ∈ N * , there exists a constant C i,A ∈ R + such that

0 max k∈Ω * (λ i,k,N -λ i,k ) C i,A e -2AN , (19)
where Ω * is the first Brillouin zone (i.e. the Voronoi cell of the lattice L of R d containing the origin). and the maximal solution Φ ε to

Φ ε (y) = G ε (Φ ε (y)), Φ ε (y η ) = 1 + η ψ ε (y η ) . (21)
As sinh(y) 0 for all y 0, we have the following differential inequality:

Ψ ε (y) G ε (Ψ ε (y)), Ψ ε (y η ) = Φ ε (y η ).
Note that, by convexity, ψ ε (y) ψ ε (y η ) > 1 for y y η : G ε is indeed quasimonotone on the domain of interest.

We are now ready to compute an upper bound of Y ε with the use of Theorem 7, which yields ∀ y y η , Ψ ε (y) Φ ε (y).

This leads us to the study of the ODE

εφ ε = -φ ε + φ 3 ε φ ε (y η ) = 1 + η, φ ε (y η ) = ψ ε (y η ) > 0.
We where C(η) is computed from the initial conditions at y = y η and B 0 < y 0 < y η is such that ψ ε (y 0 ) = 1. Note that η → -1 4 (1 + η) 4 + 1 2 (1 + η) 2 + ε 2 (ψ ε (y 0 )) 2 is decreasing on R + * and takes the value 1 4 + ε 2 (ψ ε (y 0 )) 2 > 1 4 at η = 0. Thus, for η > 0 small enough, C(η) 1 4 and we have

ε 2 (φ ε ) 2 1 4 φ 4 ε - 1 2 φ 2 ε + 1 4 = 1 4 φ 2 ε -1 2 , hence φ ε 1 √ 2ε φ 2 ε -1 on [y η , Y ε ).
Finally, we consider the ODE

ξ ε,η = 1 √ 2ε (ξ 2 ε,η -1), ξ ε,η (y η ) = 1 + η,
whose solution is

ξ ε,η (y) = 1 + 2 η + exp y-yη √ ε/2 1 + 2 η -exp y-yη √ ε/2
, which is defined only up to Y ε,η := ε 2 log 1 + 2 η +y η . Applying again Theorem 7, we have that φ ε (y) ξ ε,η (y) for any y y η such that both functions are still finite. Putting everything together, we obtain that ψ ε is only defined up to some Y ε with B 0 < Y ε Y ε,η < ∞ and that ∀ y ∈ [y η , Y ε ), ψ ε (y) ξ ε,η (y).

(

) 22
These results are illustrated on Figure 7, where we plotted the lower bound ξ ε,η for ε = 0.1 and η = 0.5 as well as a numerical approximation of ψ ε .

Figure 1 -

 1 Figure 1 -[Linear eigenvalue problem] (Left) Fourier coefficients of (the numerical approximation of) u1: they seem to decrease like 1/(|k| α √ cosh 2B0k) with α = 4, the prefactor 2B0 in the cosh being confirmed by the plot on the right. (Right) Logarithm of the ratio of two successive nonzero Fourier coefficients.

Figure 2 -

 2 Figure 2 -[Linear eigenvalue problem] (Left) Discretization error for the first eigenvalue and the associated eigenvector in the H 1 norm. Both errors seem to converge with rates given by Theorem 3, the prefactors in the cosh being confirmed by the plot on the right. (Right) Logarithm of the ratio between discretization errors for various N 's.

Figure 3 -

 3 Figure 3 -[Nonlinear eigenvalue problem] Analytic continuation of u0 for β0 = 1.00001, for which B0 ≈ 0.0045. (Left) Imaginary and real parts of y → u0(iy). A discontinuity of the derivative appears at the expected positions. (Right) Phase of z → u0(z) where z = x + iy. Branching points, in red, appear at the expected positions.

Figure 4 -

 4 Figure 4 -[Nonlinear eigenvalue problem] (Left) Fourier coefficients of (the numerical approximation of) uε: they seem to decrease like 1/(|k| α cosh(2B0k)) with α = 3/2, the prefactor 2B0 in the cosh being confirmed by the plot on the right. (Right) Logarithm of the ratio of two successive nonzero Fourier coefficients.

Figure 5 -

 5 Figure 5 -[Nonlinear elliptic equation] Analytic continuation of u0 for µ = 10, for which B0 ≈ 0.0385. (Left) Imaginary part of y → u0(iy). Discontinuities also appear at the expected positions ±B0. (Right) Phase of z → u0(z) where z = x + iy. Branching points, in red, appear at the expected positions ±iB0.

Figure 6 -

 6 Figure 6 -[Nonlinear elliptic equation] (Left) Fourier coefficients of (the numerical approximation of) uε. (Right) Logarithm of the ratio of two successive nonzero Fourier coefficients.

Theorem 6 .

 6 Let B > 0 and V ∈ H B,L . For each k ∈ R d , the eigenfunctions of the Bloch fibersH k = (-i∇ + k) 2 + V of the periodic Schrödinger operator H = -∆ + V are in H A,L for any 0 < A < B. Let λ 1,k λ 2,k• • • be the eigenvalues of H k counted with multiplicities and ranked in non-decreasing order, and λ 1,k,N λ 2,k,N • • • λ d L,N ,k,N the eigenvalues of the variational approximation of H k in the d L,Ndimensional space X L,k,N := Span(e G , G ∈ L * , |G + k| N ).

Proof.

  It suffices to replace in the proofs of Theorem 2 and Theorem 3 X N with X L,k,N and T 22 with the restriction T 22,L,k of the operator (-i∇ + k) 2 to the invariant space X ⊥ L,k,N . The latter is invertible and such that T -1 22,L,k L(X ⊥ L,k,N )

Acknowledgements

The authors would like to thank Geneviève Dusson and Michael F. Herbst for fruitful discussions. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 810367). The authors would also like to thank the anonymous reviewers for their comments and suggestions.

Data availability

All the scripts used to generate the plots of Figures 1, 2, 3, 4, 5 and 6 are available online at https://github.com/gkemlin/analytic_potentials.

Conflict of interest

All the authors declare that they have no conflicts of interest.

Appendix A. Extension to the multidimensional case with application to Kohn-Sham models.

The goal of this section is to extend the previous results to the multidimensional case and apply them to the linear version of the Kohn-Sham equations [START_REF] Babuška | Eigenvalue problems. In Handbook of Numerical Analysis[END_REF]. To this end, consider a Bravais lattice L = Za 1 + • • • + Za d where a 1 , . . . , a d are linearly independent vectors of R d (d = 3 for KS-DFT). Up to an affine change of variables (which preserves analyticity), we can take without loss of generality L = 2πZ d . We denote by Ω = [0, 2π) d a unit cell, by L * = Z d the reciprocal lattice, by e G (x) = |Ω| -1/2 e iG•x the Fourier mode with wavevector G ∈ L * , and by

the L-periodic Sobolev spaces endowed with their usual inner products. All the arguments in Sections 3.1-3.3 can be extended to the multidimensional case by introducing the Hilbert spaces

Each u ∈ H A,L can be extended to an analytic function u(z 1 , . . . , z d ) of d complex variables defined on a neighborhood on R d , and it holds

with e 1 , . . . , e d the canonical basis vectors. The approximation space X N,L is then defined as

and the inverse T -1 22,L of the restriction T 22,L of the operator -∆ on L 2 per,L to the invariant subspace

The proofs of Theorem 1, Theorem 2 and Theorem 3 can thus be straightforwardly adapted to the multidimensional case.

Lastly, if V ∈ H B,L for some B > 0, the Schrödinger operator H = -∆ + V considered this time as a Schrödinger operator on L 2 (R d , C) with an L-periodic potential, can be decomposed by the Bloch transform [START_REF] Reed | Analysis of Operators[END_REF]Section XIII.16] and its Bloch fibers are the self-adjoint operators on L 2 per,L with domain H 2 per,L and form domain H 1 per,L defined as

The following result is concerned with the Bloch eigenmodes of the H k 's.

N -2 . The result then follows similarly.

Remark 3. As a conclusion, let us mention that, under appropriate assumptions, we can extend these results to the KS-DFT equations ( 1) with GTH pseudopotentials in the case where an analytic parametrization of the exchange-correlation functional is used. This is relevant for instance when studying condensed-phase systems, for which the periodic setting is well suited, and where the commonly used approximations of V xc,ρ are analytic on the positive density values that are of interest in such cases [START_REF] Perdew | Accurate and simple analytic representation of the electron-gas correlation energy[END_REF]. Then we can rewrite (1) as a system of elliptic PDEs:

We can then use known results for elliptic systems of PDEs [START_REF] Morrey | On the Analyticity of the Solutions of Analytic Non-Linear Elliptic Systems of Partial Differential Equations: Part I. Analyticity in the Interior[END_REF][START_REF] Morrey | On the Analyticity of the Solutions of Analytic Non-Linear Elliptic Systems of Partial Differential Equations: Part II. Analyticity at the Boundary[END_REF] to conclude, in a way similar to what we proved in Theorem 4, that the orbitals ϕ i belong to H A,L for some A > 0. In particular, this leads to the exponential convergence of planewave approximations, justifying the use of GTH pseudopotentials.

Appendix B. Proof of Lemma 1

We start by rewriting [START_REF] Nottoli | A robust, open-source implementation of the locally optimal block preconditioned conjugate gradient for large eigenvalue problems in quantum chemistry[END_REF] as a first-order ODE on Ψ ε (y) := ψ ε (y) ψ ε (y) ∈ R 2 , starting at y η B 0 :

We then use the following simplified version of more general comparison results on systems of differential inequalities [START_REF] Walter | Ordinary Differential Equations[END_REF][START_REF] Wazewski | Systèmes des équations et des inégualités différentielles ordinaires aux deuxièmes membres monotones et leurs applications[END_REF]. In the sequel, the inequality a b for two vectors a, b ∈ R d means that a i b i for all 1 i d.

and the differential inequality Ψ (y) G(Ψ(y)), Ψ(y 0 ) = Φ(y 0 ).

Then we have Ψ(y) Φ(y).

on [y 0 , y M ).

To apply this result to [START_REF] Petrovskii | Sur l'analyticité des solutions des systèmes d'équations différentielles[END_REF], we introduce the function


Christophe Piscart Ae
email: christophe.piscart@univ-lyon1.fr

Benjamin Bergerot

A E Pascal

Laffaille Ae

Pierre Marmonier

C Piscart

P Laffaille

Á P Marmonier

Are amphipod invaders a threat to regional biodiversity?

Keywords: Species distributions, Crustacean invasion, Habitat, Native and invasive species, Refuge, Colonization process

The impact of invasions on local biodiversity is well established, but their impact on regional biodiversity has so far been only sketchily documented. To address this question, we studied the impact at various observation scales (ranging from the microhabitat to the whole catchment) of successive arrivals of non-native amphipods on the amphipod assemblage of the Loire River basin in France. Amphipod assemblages were studied at 225 sites covering the whole Loire catchment. Non-native species were dominant at all sites in the main channel of the Loire River, but native species were still present at most of the sites. We found that the invaders have failed to colonize most of tributaries of the Loire River. At the regional scale, we found that since the invaders first arrived 25 years ago, the global amphipod diversity has increased by 33% (from 8 to 12 species) due to the arrival of non-native species. We discuss the possibility that the lack of any loss of biodiversity may be directly linked to the presence of refuges at the microhabitat scale in the Loire channel and in the tributaries, which invasive species have been unable to colonize. The restoration of river quality could increase the number of refuges for native species, thus reducing the impact of invaders.

Introduction

The introduction of non-native species is one of the most important factors endangering native biodiversity in aquatic ecosystems [START_REF] Rahel | Assessing the effects of climate change on aquatic invasive species[END_REF], and is now considered to be the third biggest cause of biodiversity decline in aquatic ecosystems [START_REF] Sala | Global biodiversity scenarios for the year 2100[END_REF]. The arrival of invasive species is also viewed as an unprecedented form of global change [START_REF] Ricciardi | Are modern biological invasions an unprecedented form of global change?[END_REF]. Even in cases where biodiversity does not decrease, aquatic fauna can suffer from homogenization due to the arrival of non-indigenous species, the elimination of native species, and habitat alteration, which can facilitate both these processes [START_REF] Rahel | Homogenisation of freshwater faunas[END_REF][START_REF] Piscart | Biodiversity and structure of macroinvertebrate communities along a small permanent salinity gradient (Meurthe River, France)[END_REF][START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF]. Consequently, studies of the extension of invasive species have become commonplace, especially in the freshwater ecosystem literature, and numerous studies illustrating the extent of the phenomenon have recently been published [START_REF] Richardson | Fifty years of invasion ecologythe legacy of Charles Elton[END_REF]. Most of these studies have focused either on the interaction between native and non-native species, or on the spread of invasive species throughout different geographical areas and along the main channels of rivers. As a result, although the impact of invasions on local biodiversity had been well established, their consequences for regional biodiversity were still poorly documented, and this is certainly because it has proved difficult to disentangle the respective roles of environment and interspecies interactions in the success of invasive species over native ones [START_REF] Leppa ¨koski | Invasive aquatic species of Europe: distribution, impacts and management[END_REF][START_REF] Piscart | Environmental mediation of intraguild predation between the freshwater invader Gammarus pulex and the native G. duebeni celticus[END_REF]).

Among freshwater invertebrates, amphipods provide very good models for investigating the relationship between ecosystems, human disturbances and invasive species [START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF][START_REF] Piscart | Environmental mediation of intraguild predation between the freshwater invader Gammarus pulex and the native G. duebeni celticus[END_REF]. This group of species is widely distributed worldwide throughout freshwater ecosystems, and is highly diversified. Moreover, in Western Europe, non-native amphipods have met with more success than most other taxa [START_REF] Devin | Biological and ecological characteristics of invasive species: a gammarid study[END_REF][START_REF] Grabowski | How to be an invasive gammarid (Amphipoda: Gammaroidea)-comparison of life history traits[END_REF]. The ecology of European amphipod species is well documented, and several recent studies have been devoted to the dispersal of invasive species in various European catchments [START_REF] Jazdzewski | Range extensions of some gammaridean species in European inland waters caused by human activity[END_REF][START_REF] De Vaate | Geographical patterns in range extension of Ponto-Caspian macroinvertebrates species in Europe[END_REF][START_REF] Van Der Velde | Invasions by alien species in inland freshwater bodies in western Europe: the Rhine Delta[END_REF][START_REF] Jazdzewski | Recent drastic changes in the gammarid fauna (Crustacea, Amphipoda) of the Vistula River deltaic system in Poland caused by alien invaders[END_REF][START_REF] Devin | Geographic patterns in freshwater gammarid invasions: an analysis at the pan-European scale[END_REF]. For example, the Loire River basin corresponds to the largest catchment in France, and constitutes one of the most diversified large river systems for amphipods, with 12 species currently reported [START_REF] Goedmakers | Les Gammaridea (Crustaces, Amphipodes) du Massif Central[END_REF][START_REF] Chovet | Re ´partition des Gammaridae e ´pige ´s (Crustace ´s, Amphipodes) dans la Loire et les rivie `res de la Re ´gion Centre (France)[END_REF][START_REF] Bollache | Rapid range extension of the Ponto-Caspian amphipod Dikerogammarus villosus in France: potential consequences[END_REF][START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF]. Because of its central position in France, the Loire River has suffered several successive arrivals that began with Gammarus roeselii (Gervais 1835), a species that is now considered to have been naturalized in France since it arrived during the 19th century [START_REF] Jazdzewski | Bioge ´ographie de Gammarus roeselii Gervais en Europe, en particulier en France et en Pologne[END_REF], and continued with the American species Crangonyx pseudogracilis Bousfield, 1958, which arrived in 1980s (Le ´cureuil and Chovet 2003), and which has locally colonized several canals connected to the Loire River. These species were followed in 1999 by the Ponto-caspian Chelicorophium curvispinum (G.O. Sars 1895) (Le ´cureuil and [START_REF] Le ´cureuil | Connexion entre les bassins hydrographiques europe ´ens et disse ´mination des espe `ces aquatiques: le cas de la Loire moyenne[END_REF], and in 2001 by the killer shrimp Dikerogammarus villosus (Sowinskyi 1894), which has colonized most of the main channel of the Loire River (Le ´cureuil [START_REF] Le ´cureuil | Connexion entre les bassins hydrographiques europe ´ens et disse ´mination des espe `ces aquatiques: le cas de la Loire moyenne[END_REF][START_REF] Bollache | Rapid range extension of the Ponto-Caspian amphipod Dikerogammarus villosus in France: potential consequences[END_REF]. Finally, another American invasive species, Gammarus tigrinus Sexton, 1939, has recently arrived in the Loire River estuary, where it was first recorded in 2005 [START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF], and is currently expanding into the Lower Loire and adjacent canals [START_REF] Piscart | Range expansion of the North American alien amphipod Gammarus tigrinus Sexton, 1939 (Crustacea: Gammaridae) in Brittany, France[END_REF].

To understand the local and regional effects of the arrival of non-native species on the biodiversity of a large watershed, the following two questions must be addressed: (1) What is the threat to native biodiversity along the longitudinal (upstream-downstream) and lateral (from the main channel to tributaries) gradients? Or, in other terms, can native species find refuges in the Loire River or in its tributaries? (2) Does the arrival of non-native species modify the biodiversity at the regional scale? To answer these questions, we studied the amphipod assemblage of the Loire River, and the consequences of the successive arrivals of non-native amphipods at various scales of observation: from microhabitats in the main channel to the regional scale.

Materials and methods

The Loire River is the longest river in France, with a total length of 1,020 km from multiple sources in the Mont Jerbier de Jonc area (44°50 0 38 00 N, 4°13 0 12 00 E) to the estuary at St Nazaire, on the Atlantic coast (47°16 0 50 00 N, 2°12 0 31 00 W). The Loire River catchment is localized in central France. It is composed of more than 20,000 tributaries, and spreads over one fifth of the French territory (around 117,000 km 2 ). Its maximum altitude is 1,500 m, and it has a very patchy geology and landscape, consisting of woodlands, pastures, agricultural plots, and large cities downstream. The main axis of the catchment (i.e. the Loire River and its main tributary, the Allier River) is less regulated and has better longitudinal connectivity than the other tributaries (including large rivers), some of which are now substantially impounded (Lasne et al. 2007a). The Loire basin as a whole is hence usually considered to have suffered less impact than other large river systems in Europe [START_REF] Pont | Modelling habitat requirement of European fishes: do species have similar responses to local and regional environmental constraints?[END_REF].

We studied the micro-distributions of native and non-native species in three tributaries (the Beuvron River, the Cisse River and the Thouet River-Table 1). For each tributary, one site was located at the confluence with the Loire River, one site between 100 and 150 m upstream from the confluence, and one site in the Loire River channel opposite the confluence. In these nine sites, amphipods were collected separately from four microhabitats (litter, macrophyte, roots, and pebbles). Triplicate samples of 0.25 m 2 were collected at each site using a quantitative Surber net sampler (500 lm mesh size). Differences in species density with respect to the microhabitat were tested using variance nested analysis with 'microhabitat' as the fixed-effect factor, and a two-level nested analysis (microhabitat nested in position according to the confluence nested in tributaries). Tukey's HSD tests were used for multiple comparisons within the sites.

The large scale distribution of amphipods was examined at 225 sites covering the whole Loire catchment (Fig. 3). Amphipods were sampled in 2006 and 2007 using a hand net (500-lm mesh size). We used presence-absence data for each site. A total of 180 sites that harbored amphipods were used in this analysis. Environmental parameters, such as land-use type (urban, cultivated land, meadow and forest), habitat (riffles, pools, fast running channels, and slow running channels), substratum type (pebble, gravel, leaf litter, vegetation and root), mean flow velocity (m s -1 ), and river channel width (m) were measured at each sampling site. Amphipod assemblage classification was performed using self-organizing maps (SOM) according the Lasne et al. (2007a) protocol. This unsupervised artificial neural network method is used to analyze complex data sets with non-linear relationships [START_REF] Kohonen | Self-organizing map[END_REF][START_REF] Lek | Modelling community structure in freshwater ecosystems[END_REF]. In this study, the SOM procedure was used to organize the 180 sample sites into 64 neural network cells. Samples with similar species composition were classified as belonging to the same cell. By using weight vectors of trained SOMs, clustering techniques (Ward's method) were used to subdivide the SOM cells into three clusters, i.e., subgroups of amphipod assemblages. The probability that each species would be present and the diversity (number of species per sampling site) in the different clusters identified were compared using Kruskall-Wallis tests and Dunn's post test. In order to characterize clusters by testing whether they could be characterized by one or more indicator species, we used the indicator value (IndVal) according to the method developed by Dufre ˆne and Legendre (1997). The IndVal, expressed as a percentage, is based on both the fidelity and the specificity of the species for each cluster. Fidelity is highest when the species is present in all the sites in a cluster. Specificity is highest when all the individuals of a species are found in the same cluster. A Stepwise discriminant analysis was used to determine which environmental variables discriminated between the three clusters obtained with the SOM procedures.

Results

The same microdistribution pattern of the native and non-native species was observed in the three sites in relation to the distance from confluence (upstream in the tributary, confluence between the tributary and the Loire River, and downstream from the confluence in the Loire River; Fig. 1). The invader D. villosus was never found in the upstream zones in the tributary, whereas its densities had increased markedly further downstream, and it was dominant in the Loire River (F 3,8 = 5.62, P \ 0.0001). In the main channel of the Loire River, D. villosus was widely distributed in all the microhabitats without showing any clear substratum preference (F 3,8 = 2.1, P = 0.107). At the confluence, this species was mainly found in pebbles (Fig. 1; P = 0.005).

Gammarus pulex (Linnaeus 1758) was the most abundant native species found in the tributaries. This species was not randomly distributed among the microhabitats but displayed a significant substratum effect (F 9,90 = 7.57, P \ 0.0001). The microdistribution of G. pulex was considerably modified when D. villosus was present. Upstream from the confluence, the abundance of G. pulex had increased considerably from pebble to vegetation roots (Fig. 1; P \ 0.0015), whereas it had disappeared from the pebble habitat in the area of the confluence, and was restricted to the leaf litter and roots in the Loire River. We also observed that the microhabitat had a major impact on the distribution of both Echinogammarus berilloni Catta, 1878 and Echinogammarus spinulicornis Pinkster and Stock 1971 (F 9,90 = 4.33, P = 0.0001 and F 9,57 = 15.51, respectively, P \ 0.0001). Upstream of the confluence, E. berilloni was found in all microhabitats, but at a lower density in the leaf litter (P \ 0.013). In contrast, E. spinulicornis was more abundant in the leaf litter (P B 0.0001), but absent from the vegetation in the tributaries. The microdistributions of the Echinogammarus species were also considerably disrupted when D. villosus was present. The densities of E. berilloni were significantly lower at the confluence with the Beuvron River (Fig. 1a; P = 0.0001), and at that with the Thouet River (Fig. 1c; P = 0.011). In the Loire River where D. villosus was very abundant, the Echinogammarus species had disappeared near the Rivers Beuvron and Cisse, and only a few individuals were found in the leaf litter and in the roots near the Thouet River. E. spinulicornis seemed to be severely impacted by D. villosus, because had it disappeared from the confluence sites, and only few individuals remained in the pebbles and roots in the Loire River close to the Beuvron River. At the regional scale, 11 aquatic species and one terrestrial species (Orchestia cavimana Heller 1865) were found in the study area (Table 2). Among them only 7 are considered to be native (Echinogammarus berilloni, E. spinulicornis, Gammarus fossarum (Koch 1835), G. lacustris (G.O. Sars 1863), G. orinos Pinkster and Scholl 1984, G. pulex, and Orchestia cavimana), four non-native and invasive (C. pseudogracilis, C. curvispinum, D. villosus, and G. tigrinus), and one a naturalized non-native species (G. roeselii). We also observed a lack of amphipods in the upstream part of the Loire catchment basin, which is mainly composed of temporary or poorly mineralized streams (Fig. 2).

Based on the similarity of the cells in the SOM procedure, the clustering procedure identified three main clusters that were composed of 67, 92 and 21 sites, respectively. Most species varied significantly between the different clusters, except for four very scarce species G. orinos, O. cavimana, G. lacustris and C. pseudogracilis (Table 2). G. pulex was very common (mean probability of occurrence &75%), but was not present in any of the cluster-3 sites. Only D. villosus was present in all the clusters. The diversity was higher in the sites sampled in the first cluster than in those of the other two. E. berilloni was present in all the cluster-1 sites, and can be considered to be the most indicative species of this cluster (Table 2), with E. spinulicornis as the second most characteristic species of cluster 1 (Table 2). G. pulex was present in most of the cluster-2 sites, and can be considered to be indicative of the second cluster together, with G. fossarum and secondarily G. roeselii (Table 2). Finally, all the invasive species (D. villosus, C. curvispinum, G. tigrinus, and C. pseudogracilis) were found in the pooled cluster-3 sites (Table 2) and, except for C. pseudogracilis, these were the species indicative of this cluster.

The results of the stepwise discriminant analysis including all environmental variables highlighted the fact that all three clusters were significantly distinguished by river width (F 2, 172 = 32.29, P \ 0.0001) and habitat (i.e., riffles, pools, fast running channels, and slow running channels; F 2, 172 = 5.0, P = 0.0077). Sites in cluster 1 were characterized by being in small and medium size lowland rivers (mean width = 11.6 ± 11.5 m), which were all located in the north-western part of the Loire catchment and dominated by chalky geology (Fig. 2). This cluster was mainly constituted by fast running habitats (riffles, fast running channels). Cluster 2 included most of the sites also characterized by small and medium size rivers (mean width = 16.0 ± 27.0 m), but was widely distributed throughout the Loire Catchment (Fig. 2), and in most habitats (riffles, fast and slow running channels or pools). The regionalscale distribution of non-native species, corresponding mainly to the cluster 3, was characterized by large river size (mean width = 89.0 ± 103 m) with pools and slow running channels (Figs. 2,3), and was restricted to the main channel of the Loire River except for that of the naturalized G. roeselii, which has colonized some tributaries in the eastern part of the Loire catchment (Fig. 3).

An upstream-downstream pattern of alien species was found along the Loire channel, with only native species in the upstream sector, a combination of D. villosus and C. curvispinum in the first half of the river, D. villosus alone in the second half of the river and, finally, G. tigrinus in the more salty downstream area (Fig. 3).

Table 2 List, occurrence (number of sites where species were caught), origin (N native, NN non-native and NNN naturalized non-native) and mean (±SD) occurrence of each of the species found in the Loire river catchment, and the species diversity (±SD) in all the sites where amphipods were present ( 180 Although interactions between native and invasive species tend to result in reductions in the population densities of native species, some of them are able to persist in the main channel. Moreover, our findings suggest that there are refuges in the tributaries where native species have been able to resist the invaders for at least 7 years, since [START_REF] Bollache | Rapid range extension of the Ponto-Caspian amphipod Dikerogammarus villosus in France: potential consequences[END_REF] observed the arrival of D. villosus in the Loire River near the Beuvron River between 2001 and 2002. In this sector, we observed the native species E. spinulicornis and G. pulex in 2007. However, the microdistribution of native species has been severely disrupted, and they are now mainly restricted to organic substrates (i.e. roots and leaf litter). The mechanism underlying the change in the microdistribution of native species remains unclear, but it is clearly a combination of two different processes [START_REF] Dick | Post-invasion amphipod communities of Lough Neagh, Northern Ireland: influences of habitat selection and mutual predation[END_REF][START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF][START_REF] Van Riel | Interference competition among native and invader amphipods[END_REF]): a difference in microhabitat preferences that limits the overlap and interaction between native and invasive species, and competitive interaction in microhabitats that harbor both native and invasive species). This phenomenon has already been observed in Brittany and in Northern Ireland, where the microdistribution of the endemic Gammarus duebeni celticus has been disrupted by the arrival of G. pulex, which has replaced the endemic species in its optimum microhabitats [START_REF] Dick | Post-invasion amphipod communities of Lough Neagh, Northern Ireland: influences of habitat selection and mutual predation[END_REF][START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF]. In this context, the endemic G. duebeni celticus was only able to persist in microhabitats ignored by G. pulex [START_REF] Dick | Post-invasion amphipod communities of Lough Neagh, Northern Ireland: influences of habitat selection and mutual predation[END_REF]. Our results suggest that a similar process is occurring in the Loire River, where the native species E. spinulicornis and G. pulex have both resisted D. villosus in leaf litter, which is not colonized by the predatory D. villosus, and in roots, which have mainly been colonized by juveniles [START_REF] Devin | Ecological traits of the amphipod invader Dikerogammarus villosus on a mesohabitat scale[END_REF]. D. villosus juveniles probably compete with the native species less effectively than adults. The larger D. villosus, which constitute the main threat to native species [START_REF] Dick | Post-invasion amphipod communities of Lough Neagh, Northern Ireland: influences of habitat selection and mutual predation[END_REF], are indeed known to have a clear preference for hard substrata [START_REF] Devin | Ecological traits of the amphipod invader Dikerogammarus villosus on a mesohabitat scale[END_REF][START_REF] Platvoe ¨t | Dikerogammarus villosus (Sowinsky, 1894), an amphipod with bite. The relationship between morphology, behaviour, micro-distribution and impact of this invading crustacean[END_REF]MacNeil et al. 2008b), which reduces their encounters with predators, and the risk of being detected and caught by fish or invertebrate predators [START_REF] Crowder | Habitat structural complexity and the interaction between bluegills and their prey[END_REF][START_REF] Power | Habitat heterogeneity and the functional significance of fish in river food webs[END_REF].

Potential changes at the regional scale We did not find any overall pattern of extension of nonnative species in the Loire River catchment. The distribution of the non-native G. roeselii had been considerably reduced since 1994 [START_REF] Chovet | Re ´partition des Gammaridae e ´pige ´s (Crustace ´s, Amphipodes) dans la Loire et les rivie `res de la Re ´gion Centre (France)[END_REF], especially along the main channel of the Loire River, and was now restricted to a few tributaries in the middle section of the Loire catchment. The decline of G. roeselii in areas where it had been present since the 1980s [START_REF] Jazdzewski | Bioge ´ographie de Gammarus roeselii Gervais en Europe, en particulier en France et en Pologne[END_REF] [START_REF] Bollache | Rapid range extension of the Ponto-Caspian amphipod Dikerogammarus villosus in France: potential consequences[END_REF][START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF]. This was particularly true of D. villosus, which has colonized most of the Loire channel. However, the Ponto-Caspian D. villosus had failed to colonize the upstream part of the river characterized by poorly mineralized waters [START_REF] Wijnhoven | Exotic and indigenous freshwater gammarid species: physiological tolerance to water temperature in relation to ionic content of the water[END_REF]. Upstream extension of D. villosus in the future is therefore improbable due to the physico-chemical conditions (e.g. poorly mineralized waters, fast running flows), which in the Loire River are tolerated only by the native G. fossarum. Regarding the downstream extension of D. villosus, our findings suggest that this species is still extending its range beyond the distribution observed by [START_REF] Bollache | Rapid range extension of the Ponto-Caspian amphipod Dikerogammarus villosus in France: potential consequences[END_REF] in the downstream sector and in two large tributaries (the Maine and the Mayenne Rivers), both of which are located close to the downstream limit and in the northern part of the Loire catchment. It is difficult to predict the future advance of this species. However, we can predict that the extension of D. villosus has probably not yet come to an end. This species may well progress further downstream, and could eliminate G. tigrinus from the freshwater section of the Loire River, as has already been observed in several European waterways [START_REF] Dick | Invading predatory crustacean Dikerogammarus villosus eliminates both native and exotic species[END_REF][START_REF] De Vaate | Geographical patterns in range extension of Ponto-Caspian macroinvertebrates species in Europe[END_REF][START_REF] Devin | Ecological traits of the amphipod invader Dikerogammarus villosus on a mesohabitat scale[END_REF]. The continued longitudinal extension of D. villosus in the Loire River looks certain, but we did not observe any lateral extension in the tributaries, except in the downstream part of the Maine and Mayenne Rivers. This species appeared to be unable to colonize tributaries, and it was restricted to just a few meters upstream of the confluence with the Loire River (Fig. 1). The reasons for this limited lateral extension of D. villosus in tributaries are probably related to the local physico-chemical parameters, but remain unclear. This phenomenon has commonly been observed [START_REF] Jazdzewski | Four Ponto-Caspian and one American gammarid species (Crustacea, Amphipoda) recently invading Polish water[END_REF][START_REF] Fu ¨reder | (ed) Biological invaders in inland waters: profiles, distribution, and threats[END_REF][START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF]) along invasion corridors, but is surprisingly rarely discussed in the literature [START_REF] Pinkster | The introduction of the alien amphipod Gammarus tigrinus Sexton, 1939, in the Netherlands and its competition with indigenous species[END_REF][START_REF] Dick | Post-invasion amphipod communities of Lough Neagh, Northern Ireland: influences of habitat selection and mutual predation[END_REF][START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF]. Recently, [START_REF] Piscart | Environmental mediation of intraguild predation between the freshwater invader Gammarus pulex and the native G. duebeni celticus[END_REF] proposed an interesting hypothesis combining the effects of environmental conditions and intraguild predation of native species by non-native species. These authors have shown that G. pulex, an invader in Ireland, has failed to eliminate the native G. duebeni celticus as a result of intraguild predation in good quality water. This experimental finding confirmed the observation made by MacNeil et al. (2008a) in the Isle of Man that successive field introductions of G. pulex into sites previously dominated by Gammarus duebeni celticus failed in good quality rivers, whereas increased river pollution could favor the replacement of native by introduced species.

More generally, these results provided in situ evidence that environmental parameters (current velocity, microhabitats or chemical parameters) can influence the competitive equilibrium between nonnative invasive and native amphipods. This process can occur either to the detriment of the native species (which is usually what happens) or to that of nonnative species by enhancing the ability of the native species to compete with the invader.

Finally, if we look at the consequences of successive arrivals of non-native species on the biodiversity at the scale of the Loire catchment, our results do not shown that non-native species are having any strongly detrimental impact on amphipod diversity in the Loire River catchment. We only observed a major impact of invaders (i.e., D. villosus and G. tigrinus) in the main channel of the Loire River, which may have important consequences for amphipod diversity [START_REF] Dick | Post-invasion amphipod communities of Lough Neagh, Northern Ireland: influences of habitat selection and mutual predation[END_REF][START_REF] Dick | Invading predatory crustacean Dikerogammarus villosus eliminates both native and exotic species[END_REF], for the macroinvertebrate community [START_REF] Dick | Predatory impact of the freshwater invader Dikerogammarus villosus (Crustacea: Amphipoda)[END_REF][START_REF] Platvoe ¨t | Dikerogammarus villosus (Sowinsky, 1894), an amphipod with bite. The relationship between morphology, behaviour, micro-distribution and impact of this invading crustacean[END_REF], or for the fish community [START_REF] Pinkster | The introduction of the alien amphipod Gammarus tigrinus Sexton, 1939, in the Netherlands and its competition with indigenous species[END_REF][START_REF] Casellato | The predatory impact of Dikerogammarus villosus on fish[END_REF]; Platvoe ¨t 2007). In fact, at the regional scale, since the first arrival of non-native species 25 years ago, there has been an increase of 33% in amphipod diversity (from 8 to 12 species). This increase has been directly linked to the presence of refuges in the Loire channel at the microhabitat scale, and in the tributaries at the catchment scale. For example, the main changes in amphipod assemblages were caused by D. villosus but, even in the main channel of the river, our results have shown that the endemic E. spinulicornis and the native species E. berilloni and G. pulex have been able to maintain populations at least over the time scale of this study (7 years). This phenomenon has also been observed in several other European Rivers [START_REF] Dick | Post-invasion amphipod communities of Lough Neagh, Northern Ireland: influences of habitat selection and mutual predation[END_REF][START_REF] Jazdzewski | Four Ponto-Caspian and one American gammarid species (Crustacea, Amphipoda) recently invading Polish water[END_REF][START_REF] Devin | Ecological traits of the amphipod invader Dikerogammarus villosus on a mesohabitat scale[END_REF][START_REF] Piscart | Distribution and microhabitats of native and non-native gammarids (Amphipoda, Crustacea) in Brittany, with particular reference to the endangered endemic sub-species Gammarus duebeni celticus[END_REF]. The Loire River is considered to be less impacted than other large river systems in Europe [START_REF] Pont | Modelling habitat requirement of European fishes: do species have similar responses to local and regional environmental constraints?[END_REF], and the role of refuges has probably been amplified by the habitat quality and diversity, especially in the middle and downstream sections (Lasne et al. 2007b). The presence of refuges at both the microhabitat scale and the regional scale justifies optimism with regard to the conservation of native amphipod species. The restoration of altered ecosystems should reduce the expansion of invaders, and may boost re-colonization by native species from refuges. However, a lack of knowledge about the long-term ecological and evolutionary feedback between native and invasive species makes it difficult to predict the long-term impact of this restoration. Future studies therefore appear to be essential for a better prediction of the long-term consequences of biological invasions for the biodiversity of rivers.

Fig. 1

 1 Fig. 1 Mean cumulative amphipod densities (individuals per 0.25 m 2 ) at the microhabitat scale in a the Beuvron River, b the Cisse River and c the Thouet River. The first column represents sites in the tributary upstream from the confluence with the

a

  Highest values (a = 0.05) according to Dunn's post test when Kruskall-Wallis tests are significant (*** \0.001, ** \0.01, * \0.05, ns not significant) Discussion Threats to biodiversity in the Loire River and in its tributaries We have observed changes in the amphipod assemblages of the Loire River over the last 30 years as a result of the arrival of non-indigenous species (C. pseudogracilis, G. roeselii) and invasive species (D. villosus, C. curvispinum, and G. tigrinus), as well as of shifts in the distribution and density of some native (G. pulex, E. berilloni) and endemic (E. spinulicornis) species. The distributions found for regionally rare native species (G. lacustris, G. orinos, and the terrestrial landhopper O. cavimana) did not allow us to draw any conclusions about changes in their distribution.

Fig. 2

 2 Fig. 2 Map showing the locations of the 225 sites (with or without amphipods) sampled in 2006 and 2007 (a), sites constituting cluster 1 (b), cluster 2 (c) and cluster 3 (d) of the SOM procedure

  may confidently be attributed to interactions with D. villosus (e.g. intraguild predation, competition), which has colonized the middle section of the Loire River since 2001 (Le ´cureuil and Chovet 2003). The invasive species C. pseudogracilis and C. curvispinum were restricted to few backwater pools and canals (Le ´cureuil and Chovet 2003; Bollache et al. 2004), and had failed to extend their distribution in the Loire catchment.However, in this study we have confirmed an extension of the range of the invasive species D. villosus and G. tigrinus along the Loire River

Fig. 3

 3 Fig. 3 Map showing the distribution of non-native species in the Loire catchment

Acknowledgments We would like to thank Vale ´rie Briand for her essential help with literature resources and the two anonymous referees for their valuable comments. This study was supported by the 2005 PAPIER project of the ECosphe ´re COntinentale program, and the InBioProcess project (ANR-06-BDIV-007-InBioProcess (2007-2010) of the 2006 Biodiversity program of the National Research Agency (Agence Nationale de la Recherche, ANR).