David Gontier

AND Mathieu Lewin

Faizan Q Nazar

THE NONLINEAR SCHR ÖDINGER EQUATION FOR ORTHONORMAL FUNCTIONS: I. EXISTENCE OF GROUND STATES

We study the nonlinear Schrödinger equation for systems of N orthonormal functions. We prove the existence of ground states for all N when the exponent p of the non linearity is not too large, and for an infinite sequence N j tending to infinity in the whole range of possible p's, in dimensions d ≥ 1. This allows us to prove that translational symmetry is broken for a quantum crystal in the Kohn-Sham model with a large Dirac exchange constant.

Introduction and statement of the main results

1.1. Ground states for the nonlinear Schrödinger equation. The nonlinear Schrödinger equation (NLS) is one of the most famous nonlinear partial differential equation and it naturally occurs in a variety of physical situations [START_REF] Malomed | Encyclopedia of Nonlinear Science[END_REF], including Bose-Einstein condensation [START_REF] Pitaevskii | Bose-Einstein condensation[END_REF][START_REF] Lieb | The mathematics of the Bose gas and its condensation[END_REF], nonlinear optics [START_REF] Zakharov | Exact Theory of Two-dimensional Self-focusing and One-dimensional Self-modulation of Waves in Nonlinear Media[END_REF][START_REF] Manakov | On the theory of two-dimensional stationary self-focusing of electromagnetic waves[END_REF], water waves [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF], Langmuir waves in plasmas [START_REF] Taniuti | Perturbation method for a nonlinear wave modulation. I[END_REF][START_REF] Fried | On the Nonlinear Schrödinger Equation for Langmuir Waves[END_REF] and many others. We quickly recall here some of its mathematical properties before turning to its generalisation to systems of orthonormal functions, which is appropriate for quantum mechanical systems.

An NLS ground state is by definition a normalised positive solution Q ∈ L 2 (R d , R + ) to the stationary focusing NLS equation in R d :

-∆ -Q 2p-2 Q = µ Q, ˆRd Q(x) 2 dx = 1, µ < 0. (1)
In this formulation, the Lagrange multiplier µ is unknown and must be adjusted to fulfill the normalisation constraint in L 2 (R d ). For

1 < p <    ∞ in dimensions d = 1, 2, d d -2 in dimensions d ≥ 3. ( 2
)
it is known [START_REF] Coffman | Uniqueness of the ground state solution for ∆uu + u 3 = 0 and a variational characterization of other solutions[END_REF][START_REF] Kwong | Uniqueness of positive solutions of ∆uu + u p = 0 in R n[END_REF][START_REF] Mcleod | Uniqueness of positive radial solutions of ∆u + f (u) = 0 in R n . II[END_REF][START_REF] Tao | CBMS Regional Conference Series in Mathematics[END_REF][START_REF] Frank | Ground states of semi-linear PDE[END_REF] that (1) admits a unique solution (Q, µ), up to space translations for the function Q. In fact, Q is a certain dilation of the unique solution to the equation with µ = -1 and without the mass constraint. If d ≥ 3 and p ≥ d d-2 , no solution with µ < 0 can exist [START_REF] Pohozaev | On the eigenfunctions of the equation ∆u + λf (u) = 0[END_REF][START_REF] Berestycki | Nonlinear scalar field equations. I. Existence of a ground state[END_REF]. At the critical exponent p = d d-2 there is also a unique solution of the equation but it has µ = 0 (which results in an additional invariance under dilations) and it belongs to L 2 (R d ) only in dimensions d ≥ 5. When p is as in (2), the unique solution is non-degenerate [START_REF]Modulational stability of ground states of nonlinear Schrödinger equations[END_REF], which plays an important role for the behaviour of the associated time-dependent Schrödinger equation

i ∂ ∂t u = -∆ -|u| 2p-2 u,
of which Q is a stationary state. Since Q > 0 with ´Rd Q(x) 2 dx = 1, the Lagrange multiplier µ must necessarily be the lowest eigenvalue of the operator -∆ -Q 2p-2 , see [LL01,Cor. 11.9] and [START_REF]Methods of Modern Mathematical Physics. IV. Analysis of operators[END_REF]Sec. 12]. By scaling we find that the (unique) solution Q λ satisfying ´Rd Q λ (x) 2 dx = λ solves the same equation (1) with µ replaced by

µ(λ) = µ λ 2 d p-1 1+ 2 d -p .
(3) Under the additional condition

1 < p < 1 + 2 d ,
it is useful to introduce the associated NLS functional

E(u) := ˆRd |∇u(x)| 2 dx - 1 p ˆRd |u(x)| 2p dx. ( 4
)
Then Q λ is the unique solution (modulo phases and space translations) to the minimisation problem

I(d, p, λ) := min E(u), u ∈ H 1 (R d ), ˆRd |u| 2 = λ . (5)
This explains the denomination ground state. When the values of d and p are clear from the context, we will omit them in our notation and write I(λ) := I(d, p, λ).

When 1 < p < 1 + 2 d , the Gagliardo-Nirenberg inequality implies that I(λ) is finite. For larger p's one has to optimise a different functional [START_REF] Weinstein | Nonlinear Schrödinger equations and sharp interpolation estimates[END_REF].

By scaling one finds that

I(λ) = I(1) λ 1+ 2 d p-1 1+ 2 d -p . ( 6
)
Since the exponent is greater than 1 with I(1) < 0, this implies immediately that λ → I(λ) is strictly concave over R + and that I(λ) < I(λλ ′ ) + I(λ ′ ), ∀0 < λ ′ < λ.

These so-called binding inequalities guarantee the existence of a minimiser and the compactness of all the minimising sequences up to translations, by the concentrationcompactness method [START_REF] Lions | Principe de concentration-compacité en calcul des variations[END_REF][START_REF]The concentration-compactness principle in the calculus of variations. The locally compact case, Part I[END_REF][START_REF]The concentration-compactness principle in the calculus of variations. The locally compact case, Part II[END_REF][START_REF] Lewin | Describing lack of compactness in Sobolev spaces[END_REF].

1.2. The case of orthonormal functions: main results. When studying fermions, like electrons, neutrons or protons, one is naturally led to deal with systems of orthonormal functions [START_REF] Lieb | The Stability of Matter in Quantum Mechanics[END_REF], that is, u 1 , ..., u N ∈ L 2 (R d , C) with u j , u k = δ jk .

In this paper we study the nonlinear Schrödinger equation for such orthonormal systems, which could be also called the fermionic NLS equation. It takes the form

  -∆ - N j=1 |u j | 2 p-1   u i = µ i u i , i = 1, ..., N. (8)
This is a system of N coupled partial differential equations, where the coupling involves only the density

ρ(x) = N j=1 |u j (x)| 2
of the N particles. Compared to (1) and (3), we see that ρ plays the role of Q 2 , so N = ´R ρ plays here the role of λ. In what follows, we use the notation λ = N ∈ N if it is integer valued. Equation (8) has already been studied in the mathematical [LW05, AC06, BW06, BWW07] and physical [Man74, [START_REF] Kevrekidis | Solitons in coupled nonlinear Schrödinger models: A survey of recent developments[END_REF][START_REF] Zhang | The n-component nonlinear Schrödinger equations: darkbright mixed n-and high-order solitons and breathers, and dynamics[END_REF] literature, but the constraint u j , u k = δ jk is often not assumed. Dealing with this constraint is the main goal of our work. Equation (8) has several invariances. The first is the invariance under the noncompact group of space translations and it corresponds to replacing all the functions u j (x) by u j (xτ ) for some τ ∈ R d . The second is due to rotations with u j (Rx) and R ∈ SO(d). On the other hand, the multiplication by a phase for N = 1 is here replaced by the action of the (compact) group U (N ) of space-independent N × N unitary matrices, in the manner

U ∈ U (N ) → U • (u 1 , ..., u N ) = N k=1 U kj u k N j=1 .
This action does not affect the orthonormality constraint, nor the density ρ, hence it preserves the mean-field operator -∆-ρ p-1 . However it has the effect of transforming the diagonal matrix µ = diag(µ 1 , ..., µ N ) of Lagrange multipliers into U µU * . Although it could seem more appropriate to start with a general hermitian matrix of multipliers (µ jk ) 1≤j,k≤N associated with the constraints u j , u k = δ jk , we have for convenience chosen in (8) a gauge in which this matrix is diagonal.

Next, we introduce the concept of ground states for (8). We ask that the Lagrange multipliers µ i are the N first eigenvalues of the corresponding operator -∆ρ p-1 . We always use the convention that the eigenvalues are ordered increasingly and are repeated in case of degeneracies.

Definition 1 (Ground state). A ground state is a system (u 1 , ..., u N ) ∈ H 1 (R d ), orthonormal in L 2 (R d ), which solves the equations (8) where

µ 1 < µ 2 ≤ • • • ≤ µ N ≤ 0
are the N first eigenvalues of the operator -∆ρ p-1 . This definition coincides with the (strict) positivity of Q in the case N = 1, since the first eigenvalue µ 1 is always non-degenerate with a positive eigenfunction, when it exists. Our definition for N ≥ 2 is further discussed in light of the N -particle problem in Remark 8 below.

It is interesting to determine the values of p and N (depending on the dimension d ≥ 1) for which ground states exist. In this article, we focus on the case 1 < p < 1+ 2 d , which can be recast into a minimisation problem. Some kind of ground states have recently been constructed in [START_REF] Hong | Global existence versus finite time blowup dichotomy for the system of nonlinear Schrödinger equations[END_REF] for p > 1 + 2 d but their corresponding density ρ ∈ L p (R d ) is not necessarily in L 1 (R d ), so N = ´Rd ρ may in fact be infinite. The critical case p = 1 + 2 d will be handled in the second part [START_REF] Frank | The nonlinear Schrödinger equation for orthonormal functions II. Application to Lieb-Thirring inequalities[END_REF] of this work, in dimensions d ≥ 3.

A more difficult question is that of the uniqueness of ground states, when they exist. We believe that in dimension d = 1, ground states are always unique up to translations whenever they exist. Numerical simulations in dimension d = 2 presented later in Figure 2 suggest that the system may break rotational symmetry, in which case minimisers are not unique modulo space translations.

From now on, we assume

1 < p < 1 + 2 d .
As in the N = 1 case, ground states naturally occur as minimisers of the associated nonlinear functional

E(u 1 , ..., u N ) = N i=1 ˆRd |∇u i (x)| 2 dx - 1 p ˆRd N i=1 |u i (x)| 2 p dx. ( 9
)
This paper is devoted to the study of the associated minimisation problem

J(N ) = inf E(u 1 , ..., u N ), u 1 , ..., u N ∈ H 1 (R d , C), u j , u k L 2 = δ jk . (10)
For N = 1, we recover J(1) = I(1). Unfortunately, there is no simple formula such as (6) for J(N ). This is because of the orthonormality constraint, which prevents us from multiplying u i by a positive constant, as one does for I(N ) to obtain (6).

The goal of the present article is to prove that J(N ) admits minimisers, for some values of p and N . Following [START_REF]Geometric methods for nonlinear many-body quantum systems[END_REF] and as is usual in the study of nonlinear elliptic minimisation problems, our main strategy is to prove the so-called binding inequalities. Indeed, as we recall in Theorem 14 below, if N ∈ N is such that the following binding inequalities hold:

J(N ) < J(N -K) + J(K), for all K = 1, ..., N -1, (11)
then J(N ) has a minimiser and all the minimising sequences are compact, up to translations. We prove later in Proposition 16 that minimisers of J(N ) are indeed ground states. Therefore, the problem boils down to proving the binding inequalities (11). We believe that the following holds.

Conjecture 2 (Binding). For every N ≥ 2 and every

1 < p < min 2, 1 + 2 d , (12)
the binding inequalities (11) hold. In particular, J(N ) admits a minimiser, which is a ground state for (8).

Let us emphasise the new condition p < 2. The critical exponent p = 2 appears naturally in our proof, when we evaluate the interaction between two ground states placed far away. But we will also explain in Theorem 6 below that binding does not occur at p = 2 in dimension d = 1, so that the condition (12) is in fact optimal.

We were not able to prove Conjecture 2 in the whole range of parameters. We prove in this paper two weaker results. The first one is that the conjecture holds for p close enough to 1.

Theorem 3 (Binding and existence of ground states for small p). For all d ≥ 1, there exists

1 < p c (d) ≤ 1 + 2 d such that, for all 1 < p < p c (d),
the binding inequalities (11) hold for all N ≥ 2. In particular, for all 1 < p < p c (d) and all N ≥ 2 there exists a minimiser (u 1 , ..., u N ) for J(N ), which is a ground state. It solves the nonlinear system (8) where the corresponding multipliers satisfy

µ 1 < µ 2 ≤ µ 3 ≤ • • • ≤ µ N < 0
and are the N first eigenvalues of the Schrödinger operator -∆ -

( N i=1 |u i | 2 ) p-1 , counted with multiplicity.
In the proof in Section 2.5 we give an explicit lower bound of the critical exponent p c (d). This lower bound depends solely on the NLS solution Q of Equation (1). It can be numerically computed with very high accuracy using Runge-Kutta numerical methods, since Q is radial hence solves a second order Ordinary Differential Equation. We find

p c (d) >      1.614 for d = 1, 1.530 for d = 2, 1.444 for d = 3. (13)
In particular, we are able to cover the important case p = 4/3 in dimension d = 3, which is the object of Section 3. These bounds show that the critical p c (d) is not so close to 1 and let us think that p c (d) = min(2, 1 + 2 d ) should hold. See Remark 23 below for more comments on p c (d).

Our second result can cover the whole range 1 < p < min(2, 1 + 2 d ), but is valid only for an infinite sequence N j → ∞, including the two-particle case N = 2.

Theorem 4 (Binding and existence of ground states for all p). Let d ≥ 1 and

1 < p < min 2 , 1 + 2 d .
There exists an infinite increasing sequence of integers

N 1 = 1, N 2 = 2 < N 3 < • • • < N j < • • •
for which the binding inequalities (11) hold. For any such N = N j , J(N ) has a minimiser (u 1 , ..., u N ) which satisfies the same properties as in Theorem 3.

In Section 3 we use Theorem 4 to prove translational symmetry breaking of the Kohn-Sham model for a crystal with a large Dirac exchange coefficient. This result is in the same spirit as the recent work [START_REF] Ricaud | Symmetry breaking in the periodic Thomas-Fermi-Dirac-von Weizsäcker model[END_REF] by Ricaud on the Thomas-Fermivon Weizsäcker-Dirac model, and it was indeed our first motivation for studying the fermionic NLS equation ( 8). This problem naturally brings the case d = 3 and p = 4/3 which is covered by Theorem 4.

When p = 1 + 2/d the system has an additional invariance and it is not appropriate to fix the constant in front of the nonlinear term ´Rd ( N i=1 |u i | 2 ) p in (9) to be 1/p. If we study the minimisation problem similar to (10) with a constant α > 0 in front of the nonlinear term, we obtain that there exists a critical α c (N ) > 0 such that J(N ) = 0 for α ≤ α c (N ) and J(N ) = -∞ for α > α c (N ). There are no minimisers for 0 ≤ α < α c (N ). In the second part [START_REF] Frank | The nonlinear Schrödinger equation for orthonormal functions II. Application to Lieb-Thirring inequalities[END_REF] of this work, R.L. Frank and the first two authors prove a result similar to Theorem 4 for α = α c (N j ) with

N j → ∞, in dimensions d ≥ 3.
As we already mentioned, the threshold p = 2 appears naturally in our proof (Proposition 20). So our results do not cover the case 2 ≤ p < 3 in dimension d = 1. Actually, binding probably never holds for p ∈ [2, 3) in dimension d = 1.

Conjecture 5 (Absence of binding in 1D for p ∈ [2, 3)). In dimension d = 1, for all 2 ≤ p < 3 and all N ≥ 2, J(N ) does not have minimiser, and we have

J(N ) = N J(1) = N I(1).
It is explained in [FGL20a, Remark 14] that Conjecture 5 follows from the onedimensional Lieb-Thirring conjecture in [START_REF]Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities[END_REF]. Without entering into the details and using the notation of [START_REF]Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities[END_REF][START_REF]The Lieb-Thirring inequality: Recent results and open problems[END_REF], this conjecture states that the best Lieb-Thirring constant L γ,1 in dimension d = 1 coincides with the one-bound-state (Gagliardo-Nirenberg) constant L

(1) γ,1 for all 1/2 ≤ γ ≤ 3/2. In fact, Lieb and Thirring proved in [START_REF]Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities[END_REF] that L 3/2,1 = L (1) 3/2,1 for γ = 3/2 and this implies the following result, proved in [START_REF] Frank | The nonlinear Schrödinger equation for orthonormal functions II. Application to Lieb-Thirring inequalities[END_REF].

Theorem 6 (Non-existence for d = 1, p = 2 [START_REF] Frank | The nonlinear Schrödinger equation for orthonormal functions II. Application to Lieb-Thirring inequalities[END_REF]). Let d = 1 and p = 2. Then we have J(N ) = N J(1) = N I(1) for all N ∈ N. In addition, J(N ) admits no minimiser for N ≥ 2. 1.3. Comments and ideas of proof. This section contains additional comments about Theorems 3 and 4 as well as some proof ideas.

First, we place our result in a more general context. Several tools of nonlinear analysis have been generalised to systems of orthonormal functions, which can also be seen as random fields [de 15, Cd18]. The most celebrated example is the Lieb-Thirring inequality [LT75, LT76, LS10, FHJN19] which states that

N n=1 ˆRd |∇u n (x)| 2 dx ≥ c LT (d) ˆRd N n=1 |u n (x)| 2 1+ 2 d dx ( 14
)
where the positive constant c LT (d) > 0 is independent of N . This important inequality replaces the Gagliardo-Nirenberg inequality for large orthonormal systems.

It will play a role in our analysis of the large-N behaviour of J(N ) later in Section 1.4. Other Gagliardo-Nirenberg-type inequalities were considered in [START_REF]An L p bound for the Riesz and Bessel potentials of orthonormal functions[END_REF]. The fermionic NLS time-dependent equation with N = +∞ has been studied in [START_REF]On the scattering problem for infinitely many fermions in dimensions d ≥ 3 at positive temperature[END_REF]. After all these works on systems of orthonormal functions, investigating ground states of the fermionic NLS equation (8) seems a natural next step. There are many open questions and we hope that our paper will stimulate more work on the problem.

Next, we briefly explain our strategy of proof for Theorems 3 and 4 and discuss the main message stemming from this analysis. Assuming that J(N ) and J(M ) have minimisers, the proof of the binding inequality J(N +M ) < J(N )+J(M ) goes by evaluating the nonlinear interaction of these two minimisers placed far away, as is classical in such variational problems. The main difficulty here is that the functions u i all decay exponentially fast at infinity so that this interaction is exponentially small. Our main result follows from a careful evaluation of this term. In informal words, our main message is thus that quantum tunnelling can induce binding (hence existence of ground states) in a nonlinear model, provided that the nonlinearity is 'sufficiently strong'. Here 'strong' is precisely the condition that p < 2. Our simulations presented below in Section 1.4 and in [START_REF]The periodic Lieb-Thirring inequality[END_REF] suggest that this is in fact a real physical effect in dimensions d = 1, 2, not just a mathematical argument. Namely, we have numerically found ground states for J(N ) which are very close to being a combination of N copies of the NLS solution Q (see Figure 2), with an energy very close to N J(1).

To our knowledge, this is the first result of this type for a nonlinear translationinvariant problem. A similar evaluation of exponentially small nonlinear interactions can be for instance found in [Alb88, CL92, CL93a, CL93b, BL90, BL97, Dau94, OR20] but for completely different models which are not translation-invariant and have no orthonormal constraint. In [START_REF] Frank | The nonlinear Schrödinger equation for orthonormal functions II. Application to Lieb-Thirring inequalities[END_REF] our strategy is used to disprove part of the Lieb-Thirring conjecture [START_REF]Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities[END_REF].

Let us now explain where the condition p < 2 arises in the proof, which can be found in Section 2.4. We assume that there are two orthonormal systems (u n ) N n=1 and (v m ) M m=1 which are minimisers for J(N ) and J(M ), respectively. We denote the two densities by ρ = M in absolute value. Next we add to the u n 's a translation of the v m 's by a large amount R in the direction of e 1 = (1, 0, ..., 0). The resulting family of N + M functions is not orthonormal but it is exponentially close to being so. Orthonormalising it generates an error of the order O(e 2 R ) in the energy, where

e R := max n,m ˆRd |u n (x)| |v m (x -Re 1 )|dx
is the largest possible overlap between these functions. In the limit R → ∞, e R essentially behaves like exp(-min

(|µ N | 1 2 , |µ ′ M | 1
2 )R) due to the exponential decay of the eigenfunctions. On the other hand, the nonlinear term is always attractive and a careful evaluation in the intermediate region gives that it is at least of the order

exp - |µ N | |µ ′ M | |µ N | + |µ ′ M | pR .
Therefore, the nonlinearity wins over the orthonormalisation error under the assumption that

1 < p < 1 + min(|µ N |, |µ ′ M |) max(|µ N |, |µ ′ M |) (15)
where the right side is always less or equal than 2. This is how the condition p < 2 occurs. The difficulty with (15) is that the Lagrange multipliers µ N and µ ′ M are unknown. To prove Theorem 3, we derive universal lower and upper bounds (independent of N ) on the last eigenvalue µ N . This gives a critical exponent p c (d) below which binding holds. To prove Theorem 4, we observe that if N = M , then we can choose the same minimiser for J(N ) and J(M ) so that µ N = µ M . In this simpler case the interaction is attractive whenever 1 < p < min(2, 1 + 2/d). So if J(N ) has a minimiser, then

J(2N ) < 2J(N ). (16)
Since there are ground states for N = 1, we deduce that there are ground states for N = 2. By a simple pigeon-hole principle, we are then able to deduce that binding holds for an infinite sequence N j → ∞, using only (16).

We end this section with two additional comments about Theorem 3 and 4. First, we link our NLS equation for fermionic systems to a kind of Gagliardo-Niremberg-Sobolev inequality for orthonormal functions.

Remark 7 (Gagliardo-Niremberg-Sobolev for orthonormal systems). If we rescale all the u n in the manner α d/2 u n (αx) and optimise over α we obtain the inequality

N 2 d(p-1) -1 N n=1 ˆRd |∇u n (x)| 2 dx ≥ c(d, p, N ) ˆRd N n=1 |u n (x)| 2 p dx 2 d(p-1) (17) with the best constant c(d, p, N ) = N -J(N ) 1+ 2 d -p p-1 d 2p 2 d(p-1) (p -1) 1 + 2 d -p 1+ 2 d -p p-1 . ( 18
)
The constant c(d, p, N ) has a finite limit when N → ∞, as we will prove in the next section. Our theorems give the existence of optimisers for this inequality (either for small p or for a subsequence N j → ∞). The (non-sharp) inequality (17) easily follows from the Lieb-Thirring inequality, using ´Rd N n=1 |u n | 2 = N together with Hölder's inequality. In the critical case p = 1 + 2 d , the two terms in the energy scale similarly, so one cannot deduce an equality similar to (17). See [START_REF] Frank | The nonlinear Schrödinger equation for orthonormal functions II. Application to Lieb-Thirring inequalities[END_REF] for more about (17) and its link with Lieb-Thirring inequalities.

Our next remark is about how to interpret our result within the framework of N -particle anti-symmetric wave functions. The idea is that without a two-body interaction term and with a concave nonlinearity depending only on the density, a minimisation problem set on anti-symmetric wave functions can be restricted to Slater determinants.

Remark 8 (Interpretation in terms of N -particles). For a wave function Ψ ∈

L 2 ((R d ) N , C) with Ψ L 2 ((R d ) N ,C) = 1, consider the energy functional E QM (Ψ) := ˆ(R d ) N |∇Ψ(x 1 , ..., x N )| 2 dx 1 • • • dx N - 1 p ˆRd ρ Ψ (x) p dx ( 19
)
where the density ρ Ψ is defined by

ρ Ψ (x) = ˆ(R d ) N -1 |Ψ(x, x 2 , ..., x N )| 2 dx 2 • • • dx N + • • • • • • + ˆ(R d ) N -1 |Ψ(x 1 , x 2 , ..., x)| 2 dx 1 • • • dx N -1 .
Minimisation problems involving functionals of the type (19) (posed on the Nparticle space with a nonlinear term depending on ρ Ψ ) have been studied in [START_REF]Geometric methods for nonlinear many-body quantum systems[END_REF].

Here, without further constraints on Ψ, the minimum of E QM on the unit sphere is attained for a symmetric (that is, bosonic) wave function, which forms a Bose-Einstein condensate on the NLS ground state Q N / √ N defined in (5):

Ψ(x 1 , • • • , x N ) = e iθ N N 2 Q N (x 1 ) • • • Q N (x n ).
This is a simple consequence of the Hoffmann-Ostenhof inequality [HH77]

ˆ(R d ) N |∇Ψ(x 1 , ..., x N )| 2 dx 1 • • • dx N ≥ ˆRd |∇ √ ρ Ψ (x)| 2 dx which implies that E QM (Ψ) ≥ E( √ ρ Ψ )
where we recall that E is the NLS energy (4). In other words, the unconstrained N -particle problem is the same as the NLS problem for one function (5). Note that when N → ∞ the system collapses since Q N is a rescaling of Q by the factor N -(p-1)/(d+2-dp) . The situation is different if we restrict the minimisation to anti-symmetric (that is, fermionic) wave functions. From the arguments in [START_REF]Geometric methods for nonlinear many-body quantum systems[END_REF] and in Lemma 11 below, it follows that minimisers are Slater determinants (also called Hartree-Fock states), that is, of the form

Ψ(x 1 , ..., x N ) = 1 √ N ! det(u j (x k )) 1≤j,k≤N
where u 1 , ..., u N form an orthonormal system in L 2 (R d ). Slater determinants are the least correlated wave functions compatible with the anti-symmetric constraint. These wave functions satisfy that E QM (Ψ) = E(u 1 , ..., u N ), our NLS functional in (9). The N -particle interpretation of J(N ) in ( 10) is therefore that it corresponds to minimising E QM over anti-symmetric wave functions. In this light our Definition 1 of a 'ground state' is justified since the corresponding N -particle wave function Ψ is indeed a minimiser of E QM , in the anti-symmetric subspace. On the contrary to the bosonic case, the fermionic model is H-stable [Rue99] since J(N ) behaves linearly in N in the limit N → ∞, as discussed in the next section.

1.4. The large-N limit. Our main results, Theorems 3 and 4, imply that J(N ) admits a minimiser, for all N or for a subsequence. It then seems natural to ask what is happening in the limit N → ∞. This section contains results and comments in this direction. Several possible scenarios come to mind. A very natural possibility is that a sequence of minimisers ρ N = N n=1 |u n | 2 would converge (e.g. in L ∞ loc (R d ) after an appropriate translation) to some limit ρ ∞ which is extended over the whole space. This density ρ ∞ could for instance be constant (semi-classical or fluid phase), or a non-trivial periodic function (crystallisation, or solid phase [START_REF] Blanc | The crystallization conjecture: A review[END_REF]). We think that these are the only two possibilities and we explain in this section which one we expect depending on the values of p and d. We refer to [START_REF] Frank | The nonlinear Schrödinger equation for orthonormal functions II. Application to Lieb-Thirring inequalities[END_REF][START_REF]The periodic Lieb-Thirring inequality[END_REF] for a similar discussion in the context of Lieb-Thirring inequality.

1.4.1. The NLS Thomas-Fermi problem. We first introduce the NLS Thomas-Fermi problem [START_REF]Density functionals for Coulomb systems[END_REF], which will give the value of ρ ∞ in case it is constant over the whole space. For d ≥ 1 and 1 < p < 1 + 2 d we introduce the minimisation problem min

ρ≥0 ´Rd ρ=N ˆRd Cρ(x) 1+ 2 d - 1 p ρ(x) p dx ( 20
)
where C > 0 is a constant to be chosen later. This problem can be solved explicitly.

The next statement will be used several times in the paper and it states that the optimisers are for all N exactly equal to some constant ρ * depending only on C, d and p, on a set Ω of measure |Ω| = N/ρ * .

Lemma 9 (Thomas-Fermi has constant optimisers).

Let d ≥ 1, 1 < p < 1 + 2 d and C > 0. We have min ρ≥0 ´Rd ρ=N ˆRd Cρ(x) 1+ 2 d - 1 p ρ(x) p dx = -N 1 + 2 d -p d 2p d(p -1) 2pC p-1 1+ 2 d -p
with equality if and only if ρ(x) = ρ * 1 Ω (x) for some Borel set Ω of measure |Ω| = N/ρ * , where

ρ * = argmin ρ>0 Cρ 2 d - 1 p ρ p-1 = d(p -1) 2pC 1 1+ 2 d -p . ( 21
)
Proof. For α * := Cρ

2 d * -1 p ρ p-1 * , the map ρ ∈ R + → Cρ 1+ 2 d - 1 p ρ p -α * ρ
is non-negative over R + and admits exactly the two zeros 0 and ρ * . Thus we have

ˆRd Cρ(x) 1+ 2 d - 1 p ρ(x) p dx -α * N = ˆRd Cρ(x) 1+ 2 d - 1 p ρ(x) p -α * ρ(x) dx ≥ 0
with equality if and only if ρ takes only the two values 0 and ρ * , hence is of the form ρ(x) = ρ * 1 Ω (x), where Ω is a Borel set with |Ω| = N/ρ * . The minimum in the statement is thus equal to α * N and its value follows after a computation.

By Lemma 9 we see that minimisers of the Thomas-Fermi problem (20) are not unique at all, even up to translations, since any Borel set Ω is allowed as soon as |Ω| = N/ρ * . In addition, even if minimisers exist, the binding inequalities fail in NLS Thomas-Fermi theory and the energy is perfectly additive. 1.4.2. Limit N → ∞ and link with the Lieb-Thirring conjecture. Next, we prove that J(N )/N admits a limit when N → ∞ and derive upper and lower bounds on its value, in terms of the Thomas-Fermi problem (20) for two possible C's. Our main result in this section is the following Theorem 10 (Large-N limit). Let d ≥ 1 and 1 < p < 1 + 2/d. The limit e(d, p) := lim where

N →∞ J(N ) N = inf N ≥1 J(N ) N (
e LT (d, p) := -1 + 2 d -p d 2p d(p -1) 2p c LT (d) p-1 1+ 2 d -p (24)
and

e sc (d, p) := -1 + 2 d -p d 2p d(p -1) 2p c sc (d) p-1 1+ 2 d -p (25)
are the Thomas-Fermi energies obtained for C respectively equal to the Lieb-Thirring constant c LT (d) in ( 14) and to the semi-classical constant

c sc (d) = 4π 2 d (d + 2) d |S d-1 | 2 d . ( 26
)
This theorem is a corollary of our other results in this paper. Its detailed proof is provided later in Section 2.7. The lower bound in ( 23) is an immediate consequence of the Lieb-Thirring inequality (14). The existence of the limit (22) follows from the fact that J(N ) is subadditive. To understand the meaning of the semi-classical number e sc (d, p), take for instance the N first eigenfunctions u i of the Dirichlet Laplacian in a large domain Ω of volume |Ω| = N/ρ ∞ , for some constant ρ ∞ to be determined. Plugging them in the energy E, we obtain an upper bound on J(N ). By semi-classical analysis, the corresponding density N i=1 |u i | 2 is almost equal to the constant ρ ∞ inside the domain Ω. On the other hand the total kinetic energy of the u i is approximately given by Weyl's formula c sc (d)ρ

1+ 2 d
∞ |Ω| with the semi-classical constant in (26). Thus, to leading order the energy behaves as

c sc (d)ρ 1+ 2 d ∞ - ρ p ∞ p |Ω| = c sc (d)ρ 2 d ∞ - ρ p-1 ∞ p N.
From the proof of Lemma 9, this is minimized for

ρ ∞ = ρ * in (21) with C = c sc (d).
In other words, whenever ρ N converges to a constant, we expect it to be this ρ * and the limit e(d, p) to be equal to e sc (d, p). Note that the behaviour in N found in ( 22) is very different from the 'bosonic' case recalled above in (6) and Remark 8. In the later case, I(N ) behaves superlinearly in N and the exponent depends on p and d. On the contrary, for systems of orthonormal functions (fermions), J(N ) behaves linearly for all admissible p and d.

The Lieb-Thirring conjecture [START_REF]Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities[END_REF][START_REF]The Lieb-Thirring inequality: Recent results and open problems[END_REF] states that in dimensions d ≥ 3, one has c LT (d) = c sc (d). Should this conjecture be true, the upper and lower bounds would coincide in (23) and we would thus deduce that

e(d, p) = e sc (d, p) for all 1 < p < 1 + 2 d and d ≥ 3. ( 27
)
We therefore expect that the density ρ N of a ground state for J(N ) should converge to the constant ρ * given by Lemma 9 with C = c sc (d), for all [START_REF]The Lieb-Thirring inequality: Recent results and open problems[END_REF] and the two bounds in (23) do not coincide. We expect the limiting density to be a non-trivial periodic function for all 1 < p < 2 in dimensions d = 1, 2. By Theorem 6, in dimension d = 1 we think that the period will depend on p and increase when p → 2 -, whereas the density converges to Q 2 in each unit cell. In dimension d = 1, the periodicity of minimisers is confirmed by a numerical simulation 1 reported on in Figure 1. Showing such a fact is an interesting open problem [START_REF] Blanc | The crystallization conjecture: A review[END_REF]. In dimension d = 2, we could not run the computations for a too large value of N but the numerical results for N ≤ 7 presented in Figure 2 seem to suggest that the particles crystallise on a triangular lattice, as is often the case in two dimensions [START_REF] Blanc | The crystallization conjecture: A review[END_REF].

1 < p < 1 + 2/d in all dimensions d ≥ 3. On the contrary, in dimensions d = 1, 2 it is known that c LT (d) < c sc (d) [LT76,
n i P j u i 2 L 2 = i≥1 n i ˆRd |∇u i (x)| 2 dx,
and we assume it to be finite. The corresponding density is defined by

ρ γ (x) := i≥1 n i |u i (x)| 2 .
This is a non-negative integrable function with ´Rd ρ γ (x) dx = Tr(γ) = λ. The Lieb-Thirring inequality for operators [LT75, LT76, LS10, Fra20] states that

ρ γ ∈ L 1+2/d (R d ), with γ 2 d Tr(-∆γ) ≥ c LT (d) ˆRd ρ γ (x) 1+ 2 d ( 28
)
1 The code is available upon request to the authors. where γ = max(n i ) is the operator norm of γ on L 2 (R d ). When γ is an orthogonal projection,

γ = N i=1 |u i u i |, (29)
this reduces to the inequality mentioned in (14). The optimal constant c LT (d) is in fact the same as in ( 14). On the other hand, the Hoffmann-Ostenhof [START_REF] Hoffmann-Ostenhof | Schrödinger inequalities and asymptotic behavior of the electron density of atoms and molecules[END_REF] inequality states that

Tr(-∆γ) ≥ ˆRd |∇ √ ρ γ (x)| 2 dx, (30)
which implies, by the Sobolev inequality, that

ρ γ ∈      (L 1 ∩ L ∞ )(R) when d = 1, L q (R 2 ) for all 1 ≤ q < ∞ when d = 2, (L 1 ∩ L d d-2 )(R d ) when d ≥ 3.
In what follows, we assume that γ ≤ 1, and we introduce the NLS energy of any such operator γ by

E(γ) := Tr(-∆γ) - 1 p ˆRd ρ γ (x) p dx. ( 31
)
We use the same notation E as we did for the case of one or N functions, since we think that there cannot be any confusion. When γ = |u u| is a rank-one operator, we recover the usual NLS energy of one function. When γ is an orthonormal projection as in ( 29) we obtain the energy E(u 1 , ..., u N ) introduced in (9). We introduce the minimisation problem

J(λ) := inf {E(γ) : 0 ≤ γ = γ * ≤ 1, Tr(γ) = λ, Tr(-∆γ) < ∞} . ( 32
)
The following well-known result in the spirit of [START_REF] Lieb | Variational principle for many-fermion systems[END_REF] states that J(N ) coincides with our previously defined problem for orthonormal functions. This is based on the important fact that the energy is concave in γ. In fact, using that the map γ → Tr(-∆γ) is linear we find

E(tγ 0 + (1 -t)γ 1 ) =tE(γ 0 ) + (1 -t)E(γ 1 ) - 1 p ˆRd (tρ 0 + (1 -t)ρ 1 ) p -tρ p 0 -(1 -t)ρ p 1 ,
where the last integrand is negative by convexity of x → x p . Hence, the map E always attains its minimum on the extreme points of a convex set.

Lemma 11. Let d ≥ 1 and 1 < p < 1 + 2/d. The infimum J(λ) defined in (32) is finite for every λ > 0. Let N be the smallest integer such that N ≥ λ. Then we have

J(λ) = inf E(γ) : γ = N -1 i=1 |u i u i | + (λ -N + 1)|u N u N |, (u 1 , ..., u N ) ∈ H 1 (R d ) N is orthonormal in L 2 (R d ) . ( 33
)
In particular, (i) for every integer N ∈ N, J(N ) coincides with (10), (ii) for every 0 ≤ λ ≤ 1, we have J(λ) = I(λ), given by the NLS formula (6).

Note that Equation (33) is in fact true for all p > 1, even in the case J(λ) = -∞.

Proof. From the Lieb-Thirring inequality (28) and Lemma 9 we have

J(λ) ≥ inf ρ≥0 ´Rd ρ=λ ˆRd c LT (d)ρ(x) 1+ 2 d - 1 p ρ(x) p dx = e LT (d, p) λ, (34)
where e LT (d, p) is given by (24). In particular, J(λ) is finite for all λ ≥ 0.

Let us prove that we can restrict the minimisation problem to finite rank operators. Let γ = i≥1 n i |u i u i | be any admissible operator of infinite rank. Upon relabelling of the indices, we may assume 0 < n 1 < 1. Set 0 < ε < 1n 1 , and K 0 ∈ N so that i≥K0+1 n i < ε. Then, for all K ≥ K 0 , the operator

γ K :=   n 1 + i≥K+1 n i   |u 1 u 1 | + K i=2 n i |u i u i |
is also admissible and converges to γ when K → ∞, in the trace norm. We even have

Tr(-∆)γ K =   n 1 + i≥K+1 n i   ˆRd |∇u 1 | 2 + K i=2 n i ˆRd |∇u i | 2 -→ K→∞ Tr(-∆γ).
The trace-class convergence implies

ρ K → ρ in L 1 (R d ) and since ρ K is bounded in L 1+2/d (R d ) by the Lieb-Thirring inequality (28), we have ρ K → ρ in L p (R d ).
Hence E(γ K ) → E(γ) and the infimum can be restricted to finite-rank operators, as claimed.

Let now γ = i≥1 n i |u i u i | be any finite-rank admissible operator with 0 ≤ n i ≤ 1 and i≥1 n i = λ, and assume that there are two n i , say n 1 and n 2 , which belong to the open interval (0, 1). Consider the new operator obtained by varying these two occupation numbers

γt := γ + t |u 1 u 1 | -|u 2 u 2 | ,
which is admissible as soon as max(-n 1 , n 2 -1) ≤ t ≤ min(1n 1 , n 2 ). As we have explained before the lemma, the energy of γt is concave in t, hence the minimum over t must be attained at the boundary of the interval, where one of the two occupation numbers of γt is equal to either 0 or 1. The new operator γt has an energy which is lower than or equal to that of γ and it has at least one occupation number in the open interval (0, 1) less. Arguing by induction we can therefore find an operator γ ′ of the form in (33) so that E(γ ′ ) ≤ E(γ). So the minimisation problem can be restricted to such γ's.

If λ = N is an integer, the minimisation set contains only orthogonal projections, and we recover the minimisation problem introduced in (10). If 0 < λ ≤ 1 then we can take γ = λ|u u|. Its energy equals the NLS energy of √ λu. Therefore the optimum is for u = Q λ / √ λ and the minimal energy is equal to I(λ). If λ = 0 then the set is reduced to γ = 0 and we find J(0) = 0. 2.2. Binding inequalities and existence of a minimiser. Here are some standard observations. Lemma 12 (Properties of J(λ). Let d ≥ 1 and 1 < p < 1 + 2/d. Then (i) (lower bound) We have e LT (d, p) λ ≤ J(λ) < 0 for all λ > 0, where e LT (d, p) is the constant in (24). (ii) (sub-additivity) For all 0 ≤ λ ′ ≤ λ, we have

J(λ) ≤ J(λ ′ ) + J(λ -λ ′ ).
(iii) (monotonicity and continuity) The function λ ∈ R + → J(λ) is decreasing and Lipschitz, hence continuous. (iv) (concavity) It is concave on each interval (N -1, N ) with N ∈ N * .

Proof. We have already seen the inequality of (i) in (34). We now show that J(λ) is negative for λ > 0. Let γ be any admissible operator for J(λ) and set γ a = U a γU * a where U a is the dilation unitary operator by the scaling factor a, that is, in terms of operator kernels γ a (x, y) := a d γ(ax, ay). We then have 0 ≤ γ a ≤ 1, Tr(γ a ) = λ, and

E(γ a ) = a 2 Tr(-∆γ) - a d(p-1) p ˆRd ρ p γ , (35)
which is negative for a small enough since d(p -1) < 2. We turn to the proof of (ii). Take any two operators of the special form in ( 33)

γ 1 = N -1 i=1 |u i u i | + (λ ′ -N + 1)|u N u N |,
and

γ 2 = M -1 j=1 |v j v j | + (λ -λ ′ -M + 1)|v M v M |
with respectively λ ′ and λλ ′ particles. We then place γ 2 far away, at a distance R ≫ 1 along the first axis e 1 = (1, 0, ..., 0). This can be done by first translating all the functions v j into v j (• -Re 1 ) but then we have to orthonormalise the functions

u 1 , • • • , u N , v 1 (•-Re 1 ), • • • , v M (•-Re 1 )
. Although we could have chosen the u j and v j with compact support by a density argument, we use here a different reasoning which will be useful later on. We consider the Gram matrix S R of the family

u 1 , • • • , u N , v 1 (• -Re 1 ), • • • , v M (• -Re 1 ), namely S R = I N E R (E R ) * I M , E R ij = u i , v j (• -Re 1 ) .
Since S R is a Gram matrix, it is hermitian and positive. For R large enough, it is invertible, and we can set

u (R) 1 • • • u (R) N v (R) 1 • • • v (R) M := u 1 • • • u N v 1 (• -Re 1 ) • • • v M (• -Re 1 ) (S R ) -1 2 .
This family forms an (N + M )-orthonormal system. This follows from the classical fact that if Φ := (φ 1 , • • • , φ L ) is an independent family of L vectors, then the L × L matrix S with coefficients S ij := φ i , φ j is hermitian positive definite, and the family Ψ :

= (ψ 1 , • • • , ψ L ) defined by Ψ := ΦS -1/2 , that is ψ i := L k=1 φ k (S -1/2 ) ki is orthonormal, since ψ i , ψ j = L k,l=1 (S -1/2 ) ik (S -1/2 ) lj φ k , φ l =S kl = S -1/2 SS -1/2 ij = δ ij .
We then introduce the admissible operator

γ (R) = N i=1 |u (R) i u (R) i |+(λ ′ -N )|u (R) N u (R) N |+ M j=1 |v (R) j v (R) j |+(λ-λ ′ -M )|v (R) M v (R) M |
which has the required trace Tr γ (R) = λ. Since the matrix S R tends to I N +M when R → ∞, we have

u (R) i -u i H 1 (R d ) → 0, v (R) j -v j (• -Re 1 ) H 1 (R d ) → 0.
This already proves that Tr(-∆γ (R) ) -Tr(-∆γ 1 ) -Tr(-∆γ 2 ) ----→ R→∞ 0.

Similarly, we have

ρ (R) -ρ 1 -ρ 2 (•-Re 1 ) → 0 in L 1 (R d
) and since the three functions are bounded in L 1+2/d (R d ), we deduce by Hölder's inequality that the same holds in L p (R d ). For the nonlinear term, we thus have

lim R→∞ ˆRd (ρ (R) ) p = lim R→∞ ˆRd (ρ 1 + ρ 2 (• -Re 1 )) p = ˆRd ρ p 1 + ρ p 2 ,
where the last equality can be proved using the density of

C ∞ c (R d ) in L p (R d ). This proves that lim |R|→∞ E γ (R) = E(γ 1 ) + E(γ 2 ).
Hence J(λ) ≤ E(γ 1 )+E(γ 2 ). After optimising over γ 1 and γ 2 using (33) we conclude that J(λ) ≤ J(λ ′ ) + J(λλ ′ ).

The monotonicity in (iii) follows from the fact that J(λλ ′ ) < 0 by (i), hence J(λ) < J(λ ′ ) for 0 < λ ′ < λ.

Let λ, λ ′ ≥ 0 and let us prove that |J(λ) -

J(λ ′ )| ≤ C|λ -λ ′ |, that is, J is Lipschitz. If either λ or λ ′ vanishes, this is (i).
Hence, without loss of generality we may assume that 0 < λ ′ ≤ λ. Then from the monotonicity of J we have J(λ) ≤ J(λ ′ ) < 0. To get a bound in the other direction we take any trial state γ with trace Tr(γ) = λ and let γ ′ = λ ′ λ γ. We have

J(λ ′ ) ≤ E(γ ′ ) = E λ ′ λ γ = E(γ) - λ -λ ′ λ Tr(-∆γ) + 1 p 1 - λ ′ λ p ˆRd ρ p γ .
Using (1 -(1x) p ) ≤ px we deduce that

J(λ ′ ) ≤ E(γ) - λ -λ ′ λ Tr(-∆γ) - ˆRd ρ p γ . (36)
The last term takes the same form as our original problem but with no 1/p in front of the nonlinear term. Using the same argument as in (i) based on the Lieb-Thirring inequality (28) we have inf

Tr(γ)=λ
Tr(-∆γ) -

ˆRd ρ p γ ≥ -Cλ.
We have thus shown that

J(λ ′ ) ≤ E(γ) + C(λ -λ ′ ),
for a constant C depending only on p and d. Thus we obtain

J(λ) ≤ J(λ ′ ) ≤ J(λ) + C(λ -λ ′ )
after optimising over γ.

Finally we prove concavity on each interval (N -1, N ) with N ∈ N. Let N -1 ≤ λ 1 < λ < λ 2 ≤ N , and let t ∈ (0, 1) be such that λ = tλ 1 + (1t)λ 2 . For any admissible γ of the form in (33) with Tr(γ) = λ, we can write Since E is concave, this implies E(γ) ≥ tE(γ 1 ) + (1-t)E(γ 2 ) ≥ tJ(λ 1 ) + (1-t)J(λ 2 ). Minimising over γ yields the desired concavity.

γ = tγ 1 + (1 -t)γ 2 with γ 1,2 = N -1 i=1 |u i u i | + (λ 1,2 -N + 1)|u N u N |.
Remark 13 (Concavity on R + ). The concavity over R + is not expected to hold in general. Intuitively, for λ ∈ (N -1, N ) the derivative J ′ (λ) should be equal to the last (partially) filled eigenvalue µ N (see (43) below for a one-sided estimate). However at λ = N ∈ N we expect that J ′ (N ) -= µ N whereas J ′ (N ) + = µ N +1 which respectively correspond to the last filled eigenvalue when we decrease the mass or to the next eigenvalue to be filled when we increase it. Since µ N +1 ≥ µ N , J is not expected to be concave except when µ N = µ N +1 .

A numerical computation in dimension d = 1 in Figure 4 below confirms that J is not concave over R + . Also λ → J(λ)/λ is not decreasing, except when restricted to integers.

The next theorem follows from using the concentration-compactness technique for operators [Fri03, FLSS07, Lew11, LL10] and from the concavity of J(λ) on each (N -1, N ). It also follows from a profile decomposition similar to that used in [START_REF] Hong | Global existence versus finite time blowup dichotomy for the system of nonlinear Schrödinger equations[END_REF]. We will not write the proof in this paper, since the result is in fact also contained in [Lew11, Thm. 27] (with W = 0 there) using Remark 8 about the link with the N -particle problem.

Theorem 14 (Existence under the binding condition). Let d ≥ 1 and 1 < p < 1 + 2/d. Let N ∈ N be such that the following binding inequalities hold:

∀k = 1, • • • , N -1, J(N ) < J(k) + J(N -k). ( 37
)
Then, the problem J(N ) has a minimiser.

In usual concentration compactness theory, one requires the continuous binding inequality

∀0 < λ ′ < N, J(N ) < J(λ ′ ) + J(N -λ ′ ).
These inequalities automatically follow from the integer case (37) because the function λ ′ → J(λ ′ ) + J(Nλ ′ ) is concave over each interval (k, k + 1) (Lemma 12), hence its minimum is attained either at λ ′ = k or at λ ′ = k + 1.

Remark 15 (Non-integer case). Let N ∈ N and α ∈ (0, 1) be such that

∀k = 1, • • • , N, J(N + α) < J(k) + J(N -k + α).
By the concentration-compactness method one can conclude that J(N + α) has a minimiser. Note that k = N is included in the above conditions.

Properties of minimisers.

Here we state some general properties of minimisers, assuming they exist.

Proposition 16 (Euler-Lagrange equations). Let d ≥ 1 and 1 < p < 1 + 2/d. If J(λ) admits minimisers then it possesses one which is of the form in (33), with (orthonormal) real-valued eigenfunctions:

γ = N -1 i=1 |u i u i | + (λ -N + 1)|u N u N |, ( 38
)
where N is the smallest integer such that N ≥ λ. The u i are the N first eigenfunctions of the operator -∆ρ p-1 γ , counted with multiplicity:

-∆ -ρ p-1 γ u i = µ i u i , i = 1, ..., N (39)
with

µ 1 < µ 2 ≤ • • • ≤ µ N < 0.
In particular, -∆ρ p-1 γ has at least N negative eigenvalues. The functions u i are real-analytic and tend to zero at infinity. Finally, we have the estimate

2p -d(p -1) 2 -d(p -1) J(λ) λ ≤ µ N ≤ J(1)(λ -N + 1) 2 d p-1 1+ 2 d -p < 0 ( 40
)
on the last filled eigenvalue.

Equation (38) can also be written in the compact operator form

γ = 1(H γ < µ N ) + δ, with H γ := -∆ + ρ p-1 γ ,
where δ is a self-adjoint operator on Ker(H γµ N ), which is the sum of a projector plus a rank-one operator.

In the proof of the proposition, we are going to use the following remark which follows from the concavity of E.

Lemma 17 (A general inequality). Let d ≥ 1 and 1 < p < 1 + 2/d. Let 0 ≤ γ = γ * ≤ 1 and 0 ≤ γ ′ = (γ ′ ) * ≤ 1 be two admissible operators. Then we have

E(γ ′ ) ≤ E(γ) + Tr H γ (γ ′ -γ) (41)
where

H γ := -∆ -ρ p-1 γ is self-adjoint on H 2 (R d
) and the trace in (41) is understood in the quadratic form sense.

Proof of Lemma 17. The proof of (41) follows from the equality

E(γ ′ ) = E(γ) + Tr(-∆ -ρ p-1 γ )(γ ′ -γ) - 1 p ˆRd ρ p γ ′ -ρ p γ -pρ p-1 γ (ρ γ ′ -ρ γ ) ,
and the fact that the last integrand is non-negative by the convexity of x → x p . Let us prove that H γ is self-adjoint. From (30) we have

ρ γ ∈ L p (R d ) for 1 < p < ∞ in dimensions d = 1, 2 and for 1 < p < 1 + 2/(d -2) in dimensions d ≥ 3. In particular, ρ p-1 γ ∈ L r (R d ) for all 1/(p -1) < r < ∞ in dimensions d = 1, 2 and for 1 p -1 < r < 1 + 2 d-2 p -1 , where 1 + 2 d-2 p -1 > d 2 1 + 2 d -2 > max 2, d 2
in dimensions d ≥ 3. From the Rellich-Kato and Weyl theorems [START_REF]Methods of Modern Mathematical Physics. IV. Analysis of operators[END_REF], this shows that the operator

H γ = -∆ -ρ p-1 γ is self-adjoint on H 2 (R d
) and that its essential spectrum equals [0, ∞).

With Lemma 17 at hand we can write the Proof of Proposition 16. We split the proof into several steps.

Step 1: Equation. Let γ be a minimiser for J(λ). Defining the complex conjugate γ := i n i |u i u i | which has the same density ρ γ = ρ γ , we notice that

E(γ) = E γ + γ 2 ,
so we may assume that γ = γ. The eigenfunctions of γ can then be chosen real. Let γ ′ be any other admissible operator of trace λ = Tr(γ ′ ). Using (41) and the fact that E(γ ′ ) ≥ E(γ) we deduce that

Tr H γ (γ ′ -γ) ≥ 0.
In other words, γ also solves the linear minimisation problem inf

Tr(γ ′ )=λ Tr H γ γ ′ . ( 42
)
Minimisers of the linear problem (42) exist only when H γ has at least N non-positive eigenvalues, and are all of the form

γ ′ = 1 (-∞,µ N ) H γ + δ with 0 ≤ δ = δ * ≤ 1 {µ N }
H γ and δ = 0. This is called the aufbau principle in quantum chemistry. The eigenvalues are filled starting from the bottom and only the last eigenvalue can be partially filled. Our minimiser γ must therefore be of this form. Note that H γ is real since ρ γ is a real function, therefore 1 (-∞,µ N ) (H γ ) is real as well. We conclude that δ = δ. It remains to show that δ is a projection plus a rank-one operator. In the next step we prove that µ N < 0, which already implies that δ must be finite rank.

Step 2: Estimates on µ N . We first show that µ N < 0 = min σ ess (H γ ). We consider

γ ′ = γ -t|u N u N |, which is admissible for 0 ≤ t ≤ (λ -N + 1) ≤ 1. From (41) we have 2 J(λ -t) ≤ E(γ ′ ) ≤ E(γ) -µ N t = J(λ) -µ N t.
(43) Using J(λt) ≥ J(λ) -J(t) and the explicit formula for J(t) = I(t) in (6) we obtain the inequality

µ N ≤ J(λ -N + 1) λ -N + 1 = J(1)(λ -N + 1) 2 d p-1 1+ 2 d -p < 0.
Next we derive the lower bound (40) on µ N . To this end we use the virial (also called Pohozaev) identity. Let γ a be the rescaled operator as in (35). Then the function a → E(γ a ) = a2 Tr(-∆γ) -a d(p-1) p ˆRd ρ p γ must attain its minimum at a = 1. Writing that the derivative vanishes at this point we find the virial identity

Tr(-∆γ) = d(p -1) 2p ˆRd ρ p γ .
This gives

J(λ) = E(γ) = d(p -1) -2 2p ˆRd ρ p γ and Tr(H γ γ) = d(p -1) -2p 2p ˆRd ρ p γ , so that 2p -d(p -1) 2 -d(p -1) J(λ) = Tr(H γ γ) = N -1 i=1 µ i + (λ -N + 1)µ N ≤ λµ N , (44)
where the last inequality comes from the fact that µ i ≤ µ N for all i. We obtain as claimed

µ N ≥ 2p -d(p -1) 2 -d(p -1) J(λ) λ .
Step 3: Regularity and decay. Note that the first eigenfunction u 1 of H γ is always positive and non-degenerate. Therefore ρ γ > 0. Since we now have a system of finitely many coupled Partial Differential Equations, the real-analyticity of the u i 's follows from classical results [START_REF] Morrey | On the analyticity of the solutions of analytic non-linear elliptic systems of partial differential equations. I. Analyticity in the interior[END_REF][START_REF] Kato | New idea for proof of analyticity of solutions to analytic nonlinear elliptic equations[END_REF].

Next we show that the functions u i tend to 0 at infinity. In dimension d ≤ 3, this follows from the fact that u i ∈ H 2 (R d ) as we have seen in Lemma 17. In dimensions d ≥ 4 we need to employ a simple boot-strap argument. Assuming that u i ∈ L r (R d ) for all i = 1, ..., N , we infer that ρ ∈ L r/2 (R d ) and then (-

∆+|µ i |)u i = ρ p-1 u i ∈ L r 2p-1 (R d ). When r < d(2p -1)/2, the Sobolev embedding shows that u i ∈ L f (r) (R d ) with f (r) = dr/(2dp -d -2r
). This function f has the two fixed points r = 0 and r = d(p -1) < 2, the latter being unstable. Starting from any r 0 > 2 we obtain after iterating f finitely many times an r > d(2p -1)/2. The Sobolev-Morrey embeddings now prove that the u i are continuous and tend to 0 at infinity.

Step 4. Form of δ. To prove that δ = δ is a finite rank projection plus a rank-one operator, we assume by contradiction that δ has two eigenvalues δ 1 , δ 2 ∈ (0, 1) with corresponding orthonormal real-valued eigenfunctions u i , u j and we vary the corresponding eigenvalues linearly like δ 1 + t, δ 2t, as we did in the proof of Lemma 11. The energy is concave in t, hence must be constant since γ is a minimiser. The nonlinear term is even strictly concave, unless |u i | = |u j |. Since these are real-analytic real-valued functions, it would imply u i = ±u j everywhere, a contradiction. Therefore, at most one eigenvalue of δ can be in (0, 1) and this concludes the proof of Proposition 16.

Remark 18. The upper bounds on µ N in (40) deteriorates when λ → (N -1) + . We were not able to bound J ′ (N ) + , as it probably requires the evaluation of µ N +1 (λ), which may vanish as λ → N + .

Using the Euler-Lagrange equation (39), we can prove that the functions u i (and therefore the density ρ) are exponentially decaying. Actually, we may provide lower bounds as well. This is not obvious, because the functions u i have non trivial nodal sets for i ≥ 2. Only u 1 is positive everywhere. Following [BM77, HOHOS85], we introduce, for f ∈ L q loc (R d , C), the q-spherical average

[f ] q (x) := 1 |S d-1 | ˆSd-1 f |x|ω q dσ(ω) 1 q = ˆSO(d) f Rx q dR 1 q . ( 45
)
In the second integral we use the normalised Haar measure on SO(d).

Lemma 19 (Decay of minimisers at infinity). Let d ≥ 1 and 1 < p < 1 + 2/d. Let γ be a real minimiser of J(λ) of the form (33), with density ρ = ρ γ . Then we have the bounds

1 C e -2 √ |µ N ||x| 1 + |x| d-1 ≤ [ρ] 1 (x) and ρ(x) + |∇ρ(x)| ≤ C e -2 √ |µ N ||x| 1 + |x| d-1 (46)
for some constant C > 0. Similarly, for the eigenfunctions u i of -∆ρ p-1 with eigenvalue µ i as in (39), we have

1 C i e - √ |µi||x| 1 + |x| d-1 2 ≤ [u i ] 2 (x) and |u i (x)| + |∇u i (x)| ≤ C i e - √ |µi||x| 1 + |x| d-1 2 . ( 47
)
The slowest exponential is the one corresponding to the N th eigenvalue µ N and it is the leading term in the density ρ.

Proof of Lemma 19. First we derive a non-optimal exponential bound on ρ. We have

-∆ρ = 2 N i=1 u i (-∆u i ) -|∇u i | 2 = 2 N i=1 µ i |u i | 2 + ρ p-1 |u i | 2 -|∇u i | 2 ≤ 2(µ N + ρ p-1 )ρ,
where we used the fact that µ i ≤ µ N < 0 in the last inequality. Since ρ goes to 0 at infinity and p > 1, there exists a large enough R > 0 such that ρ p-1 (x) < 1 2 |µ N | for all |x| ≥ R. In particular, we have

(-∆ + |µ N |) ρ(x) ≤ 0, ∀|x| ≥ R. Let Y m (x) := √ m |x| 1-d/2 K d-2 2
(m|x|) > 0 be the Yukawa potential, solution to

(-∆ + m 2 )Y m = 0 on R d \ {0}.
Here K α is the modified Bessel function of the second kind. From the asymptotic behavior of K α [AS64, 9.7.2] we have

Y m (x) ∼ |x|→∞ π 2 e -m|x| |x| d-1 2 . ( 48
)
On the sphere RS d-1 of radius R, the function ρ is bounded, so we have ρ

(x) ≤ CY m (x) for some C > 0 with m = |µ N |. Since (-∆ + |µ N |)(ρ -CY m ) = (-∆ + |µ N |)ρ ≤ 0 on R d \ B R ,
we deduce from the maximum principle [LL01, Chapter 9.4] that ρ(x)

≤ CY m (x) on R d \ B R .
The function ρ is bounded on B R and therefore we have proved the pointwise upper bound

∀x ∈ R d , ρ(x) ≤ C ′ e - √ |µ N ||x| 1 + |x| d-1 2 . ( 49
)
Compared with (46), we see that a factor 2 is missing in the exponential. However we have learned that the potential -ρ p-1 is exponentially decaying at infinity and for proving (47) we can now rely on existing results for the eigenfunctions u i of the linear Schrödinger operator -∆ρ p-1 . The upper bound in (47) is well known but we give a detailed proof for completeness. Let for instance

F i (x) := √ m i |x| 1-d 2 K d-2 2 +ε (m i |x|) with m i = |µ i | which solves, this time, -∆ + |µ i | + ε(d -2 + ε) |x| 2 F i = 0 on R d \ {0} (50)
see for instance [START_REF] Dereziński | On Schrödinger operators with inverse square potentials on the half-line[END_REF]. We choose ε small enough so that the inverse-square potential is attractive: ε(d -2 + ε) < 0. In other words, we take ε > 0 in dimension d = 1 and ε < 0 in dimension d ≥ 3. The function F i is then positive on R d \{0}. In dimension d = 2 we have to take ε ∈ iR and the resulting radial function F i vanishes infinitely many times close to the origin but it is positive for |x| ≥ |ε| [START_REF] Ferreira | Zeros of the modified Hankel function[END_REF]. In all cases we have again

F i (x) ∼ |x|→∞ π 2 e - √ |µi||x| |x| d-1 2
by [AS64, 9.7.2]. Using Kato's inequality [RS75, Theorem X.27] and the exponential decay of ρ p-1 we find

-∆ + ε(d -2 + ε) |x| 2 + |µ i | |u i | ≤ 0 on R d \ B R
for R large enough. The maximum principle then gives as before |u i | ≤ CF i on R d \ B R , and thus the upper bound on |u i | in (47). For the estimate on ∇u i , we can use that

u i = (-∆ + |µ i |) -1 (ρ p-1 u i ) = CY mi * (ρ p-1 u i ) with m i := |µ i |, so that |∇u i | ≤ C|∇Y mi | * (ρ p-1 |u i |).
The function |∇Y mi | is integrable in a neighborhood of the origin and it behaves like |µ i |Y mi at infinity, due to the fact that

K ′ α (r) ∼ r→∞ -K α (r), (51)
see [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]9.7.4]. Thus, using our bound on |u i | and the exponential bound (49) for ρ p-1 , we find

|∇u i | ≤ C(|∇Y mi |1 B1 + Y mi ) * (ρ p-1 |u i |) ≤ C(|∇Y mi |1 B1 ) * Y Mi + CY mi * Y Mi with M i := (p -1) |µ N | + |µ i | > m i .
(|k| 2 + m 2 i )(|k| 2 + M 2 i ) = 1 M 2 i -m 2 i 1 |k| 2 + m 2 i - 1 |k| 2 + M 2 i and hence Y mi * Y Mi = C Y mi -Y Mi M 2 i -m 2 i ≤ C M 2 i -m 2 i Y mi . (52)
The pointwise bound on |∇u i | in (47) follows from the behaviour (48) of Y mi at infinity and the fact that u i is C ∞ , hence |∇u i | is bounded on compact sets.

For the lower bound in (47) this is more complicated and we just apply [HOHOS85, Corollary 2.2]. This result exactly states that if u i ∈ H 1 (R d ) is solution to (-∆ + V )u i = µ i u i with µ i < 0 and |V | ≤ C|x| -1-ε for some ε > 0 at infinity, then the lower bound in (47) holds. Thus we have proved all the estimates in (47).

Finally, we come back to (46) and notice that 2.4. Proof of binding. We now focus on the proof of the binding inequality J(λ + λ ′ ) < J(λ) + J(λ ′ ). The usual proof is to consider minimisers for J(λ) and J(λ ′ ), and to construct from them a good candidate for J(λ + λ ′ ) by putting these two minimisers far from each other. In our case, all quantities are exponentially decaying, which makes the evaluation of the interaction quite delicate. The following is the heart of the paper.

ρ(x) = N i=1 u i (x) 2 , [ρ] 1 (x) = N i=1 [u i ] 2 (x) 2 ≥ [u N ] 2 (x) 2 . Since 0 < |µ N | ≤ |µ N -1 | ≤ • • • < |µ 1 |,
Proposition 20 (Exponentially small binding). Let d ≥ 1, 1 < p < 1 + 2/d and λ, λ ′ > 0. Assume J(λ) and J(λ ′ ) admits the respective minimisers γ and γ ′ , satisfying the properties in Proposition 16. Let µ = µ N and µ ′ = µ ′ M be the associated last filled eigenvalue of H γ and H γ ′ . Then, under the additional condition

1 < p < 1 + min{|µ|, |µ ′ |} max{|µ|, |µ ′ |} , ( 53
)
we have the binding inequality J(λ + λ ′ ) < J(λ) + J(λ ′ ).

Proof. We only write the proof in the integer case λ = N ∈ N and λ ′ = M ∈ N for clarity. The arguments are exactly the same in the non-integer case, but the notation is a bit more heavy due to the additional rank-one operator.

Let γ := N i=1 |ũ i ũi | and γ′ := M j=1 |ṽ j ṽj | be two real minimisers for J(N ) and J(M ) respectively. Recall that our problem is invariant under rotations. Thus we can introduce u i (x) = ũi (Rx) and v j (x) = ṽj (R ′ x) for some R, R ′ ∈ SO(d) and deduce that these are minimisers of J(N ) and J(M ) as well. We will need to choose R and R ′ appropriately. For clarity, we prefer to postpone this discussion to Lemma 21 below. Note that the functions u i , v j satisfy the exponential bounds (47) uniformly in R, R ′ since those are invariant under rotations.

Next, we place the second system far away. For R > 0, we set v j,R (x) := v j (x -Re 1 ) where e 1 = (1, 0, ..., 0), and we introduce the Gram matrix

S R := I N E R (E R ) * I M , with E R ij := u i , v j,R = ˆRd u i (x)v j (x -Re 1 )dx,
as we did in the proof of Lemma 12. Since the functions u i and v j are real-valued and exponentially decaying, E R is real and goes to 0 exponentially fast. So the Gram matrix S R is real symmetric positive definite for R large enough. As before, the frame

ψ 1,R • • • ψ N +M,R := u 1 • • • u N v 1,R • • • v M,R (S R ) -1 2
is orthonormal. Our trial state is the orthogonal projection onto this frame, given by

γ R = N i=1 |ψ i,R ψ i,R | + M k=1 |ψ N +k,R ψ N +k,R | = N i,j=1 (S -1 R ) ij |u i u j | + M k,ℓ=1 (S -1 R ) N +k,N +ℓ |v k,R v ℓ,R | + N i=1 M k=1 (S -1 R ) i,N +k |u i v k,R | + (S -1 R ) N +k,i |v k,R u i | .
To compute E(γ R ), we consider the Taylor expansion with respect to the largest overlap e R := max i,j ˆRd |u i (x)| |v j (x -Re 1 )|dx.

Note that E R ≤ e R . We compute all quantities to the order O(e 2 R ).

First, from (1 + E) -1 = 1 -E + O(E 2 ), we have (S R ) -1 = I N 0 0 I M - 0 E R (E R ) * 0 + O(e 2 R ).
This gives, to first order and with

γ ′ R (x, y) = γ ′ (x -Re 1 , y -Re 1 ) the translation of γ ′ , γ R = γ + γ ′ R - N i=1 M j=1 E R ij (|u i v j,R | + |v j,R u i |) + O • 1,1 (e 2 R ) (54)
where γ 1,1 = Tr| √ 1 -∆γ √ 1 -∆| is the Sobolev-type trace norm. Let us evaluate the different terms in the energy. For the kinetic energy, we obtain (recall that everything is real-valued and that N and M are finite)

Tr(-∆γ R ) -Tr(-∆γ) -Tr(-∆γ ′ ) = -2   N i=1 M j=1 E R ij ˆRd ∇u i • ∇v j,R   + O(e 2 R ).
Using the Euler-Lagrange equations (8), we see that

-∆ -ρ p-1 u i = µ i u i and -∆ -(ρ ′ R ) p-1 v j,R = µ ′ j v j,R , with ρ = ρ γ and ρ ′ R = ρ γ ′ (• -Re 1 )
. This gives for instance

ˆRd ∇u i • ∇v j,R = µ i E R ij + ˆRd ρ p-1 u i v j,R .
Since ρ is bounded, together with the definition of e R , we deduce that

Tr(-∆γ R ) -Tr(-∆γ) -Tr(-∆γ ′ ) = O(e 2 R
). We now compute the difference for the term ´ρp . We first find an expression for ρ R , the density of γ R . From (54), we get

ρ R = ρ + ρ ′ R -2 N i=1 M j=1 E R ij u i v j,R + O L p (R d ) (e 2 R ).
This gives

ˆRd ρ p R = ˆRd (ρ + ρ ′ R ) p -2p ˆRd (ρ + ρ ′ R ) p-1 N i=1 M j=1 E R ij u i v j,R + O(e 2 R ).
Again, using that ρ and ρ ′ R are bounded functions, the last integral is of order O(e 2 R ). Altogether, this proves that

J(N + M ) -J(N ) -J(M ) ≤ E(γ R ) -J(N ) -J(M ) = - 1 p ˆRd (ρ + ρ ′ R ) p -ρ p -(ρ ′ R ) p + O(e 2 R ).
The orthonormalisation procedure generates an error of the order O(e 2 R ) in the energy. The first term of the second line is the nonlinear interaction and it is always negative, from the concavity of x → -x p . The question is whether it wins over the error term O(e 2 R ). This is the topic of the next result. We recall that our functions u i and v j depend on the two rotations R and R ′ which we have introduced at the beginning of the proof and which we choose now, in terms of R.

Lemma 21. Let ε := |µ N | and ε ′ := |µ M | and assume, without loss of generality, that ε ′ ≤ ε. Then, there is C ≥ 0 so that for all R large enough

e R ≤ CR 3-d 2 e -ε ′ R , ( 55
)
uniformly in R, R ′ ∈ SO(d).
On the other hand, for every R large enough there exists R, R ′ ∈ SO(d) such that

I R := ˆRd (ρ + ρ ′ R ) p -ρ p -(ρ ′ R ) p ≥ CR -p(d-1) e -2p εε ′ ε+ε ′ R . ( 56
)
In particular, if

1 < p < 1 + ε ′ ε , then e 2 R = o(I R ).
The polynomial factors R

3-d 2

and R -p(d-1) in ( 55) and (56) are not necessarily optimal but will suffice for our argument. Before we provide the proof of Lemma 21, we remark that it immediately gives the strict inequality

J(N + M ) -J(N ) -J(M ) < 0 after taking R large enough, whenever 1 < p < 1 + ε ′
ε . This condition is equivalent to the one in (53). It thus only remains to provide the proof of Lemma 21.

ˆRd |u i ||v j,R | ≤ CY mi * Y m ′ j (Re 1 ) ≤ CY ε * Y ε ′ (Re 1 ). When ε ′ < ε we have Y ε * Y ε ′ ≤ C(ε -ε ′ ) -1 Y ε ′ by (52). When ε ′ = ε we can use that Y ε * Y ε = Cε -1 ∂ m Y m |
m=ε which behaves at infinity like rY ε (r), by (51). This gives the bound (55).

Let us now bound I R from below. Recall that ρ(x) = ρ(Rx) and ρ ′ (x) = ρ′ (R ′ x) From Lemma 19, the spherical average of ρ is bounded below in terms of (1 + |x|) 1-d e -2ε|x| . This means that for any large enough r, there exists at least one x r on the sphere of radius r such that ρ(x r ) ≥ cr 1-d e -2εr . The same lemma provides the pointwise upper bound |∇ρ(x)| ≤ Cr 1-d e -2εr and this implies that ρ(x) ≥ c e -2εr 2r d-1 on the ball B(x r , η) with η = c/(2C). The same property holds for ρ′ at a point x ′ r . We choose the rotations R and R ′ to align the points x r and x ′ r ′ with the radii r + r ′ = R chosen to obtain the largest possible interaction. More precisely, we introduce the point

x * = ε ′ ε + ε ′ Re 1 where |x| 1-d e -2ε|x| and |x-Re 1 | 1-d e -2ε|x-Re1| are of the same order R 1-d e -2 εε ′ ε+ε ′ R . We then take R, R ′ such that x r = Rx * , x ′ r ′ = R ′ (x * -Re 1 ), r = ε ′ ε + ε ′ R, r ′ = ε ε + ε ′ R. We obtain ρ(x) ≥ c ′ R 1-d e -2 εε ′ ε+ε ′ R , ρ ′ R (x) ≥ c ′ R 1-d e -2 εε ′ ε+ε ′ R , for all x ∈ B(x * , η), with c ′ = c 2 ε ε + ε ′ 1-d .
Since the function (x, y) → (x + y) px py p is increasing in x and in y separately for p > 1, we deduce that

ρ + ρ ′ R p -ρ p -(ρ ′ R ) p (x) ≥ (c ′ ) p (2 p -1)R p(1-d) e -2p εε ′ ε+ε ′ R ,
for all x ∈ B(x * , η). From the positivity of the integrand on R d we thus obtain

I R ≥ ˆB(x * ,η) ρ + ρ ′ R p -ρ p -(ρ ′ R ) p ≥ |B(0, η)|(c ′ ) p (2 p -1)R p(1-d) e -2p εε ′ ε+ε ′ R
which is the claimed lower bound (56).

This concludes the proof of Proposition 20.

One immediate consequence of Proposition 20 is the case λ = λ ′ , where the condition (53) is always satisfied whenever p < 2. Since J(1) = I(1) always has a minimiser, we also conclude that J(2) always has one, for p < 2.

Corollary 22. Let d ≥ 1 and 1 < p < min{2, 1 + 2/d}. Then, if J(λ) has a minimiser, we have J(2λ) < 2J(λ). In particular, J(2) < 2J(1) and J(2) has a minimiser.

2.5. Proof of Theorem 3. Assume that J(N ) and J(M ) have minimisers. Let µ N and µ ′ M be the corresponding last filled eigenvalues. From (40) in the integer case λ = N , and the fact that 0 > J(N )/N ≥ e LT (d, p) by (i) in Lemma 12, we have min{|µ

N |, |µ ′ M |} max{|µ N |, |µ ′ M |} ≥ |I(d, p, 1)| |e LT (d, p)| 2 -d(p -1) 2p -d(p -1) . (57)
Let p c (d) ∈ (1, 2) be the first zero of the function

p → 1 + |I(d, p, 1)| |e LT (d, p)| 2 -d(p -1) 2p -d(p -1) -p.
Using that e LT (d, p) and I(d, p, 1) have a finite limit when p → 1 + we deduce that p c (d) > 1. Hence, if J(N ) and J(M ) have minimisers and if 1 < p < p c (d), then J(N + M ) < J(N ) + J(M ) by Proposition 20. Since we already know that J(1) = I(1) has a minimiser, we can deduce by induction on N that binding holds and that J(N ) has a minimiser for all N , when 1 < p < p c (d). This concludes the proof of Theorem 3. An exact upper bound on I(d, p, 1) can be obtained using Gaussian functions as trial states:

I(d, p, 1) ≤ - d 2 1 + 2 d -p p -1 π d(p-1) 2+d-dp p -2+d 2+d-dp . ( 59
)
This happens to be very precise in the regime of interest. In fact, using (59) and (58), we already obtain

p c (d) >      1.612 for d = 1, 1.526 for d = 2, 1.440 for d = 3. (60)
We can therefore safely claim that 4/3 < p c (3).

2.6. Proof of Theorem 4. We now prove that J(N ) has a minimiser for an infinity of integers N ∈ N. Let us call N ⊂ N the set of all the integers n which satisfies the binding inequalities J(n) < J(k) + J(nk) for all k = 1, ..., n -1. In particular, J(n) has a minimiser for all n ∈ N . For 1 < p < min(2, 1 + 2/d) we already know that 1, 2 ∈ N . Our goal is to show that N is not finite.

If N / ∈ N , then there is 1 ≤ k ≤ N -1 so that J(N ) = J(k) + J(N -k). If k / ∈ N or N -k /
∈ N , we can further decompose J(k) or J(Nk), and so on, until we obtain a decomposition of the form

J(N ) = n∈N k n J(n), N = n∈N k n n, k n ∈ N. (61)
We claim that k n ∈ {0, 1} and prove this by contradiction. Let n 0 ∈ N be so that k n0 ≥ 2 in the previous decomposition. Then, since n 0 ∈ N , we have by Proposition 20 that J(2n 0 ) < 2J(n 0 ). Together with the weak-binding inequality, this would imply

J(k n0 n 0 ) ≤ (k n0 -2)J(n 0 ) + J(2n 0 ) < k n0 J(n 0 ).
In other words, if two of more "bubbles" have the same number of particles, it is energetically favourable for these bubbles to merge. This gives

J(N ) = J    n∈N n =n 0 k n n + k n0 n 0    < n∈N n =n 0 k n J(n) + k n0 J(n 0 ) = n∈N k n J(N ),
which is in contradiction with (61). So the coefficients k n in (61) must all be equal to 0 or 1. Now N cannot be finite otherwise we would not be able to write all the N / ∈ N as in (61) with k n ∈ {0, 1}. This concludes the proof that N is infinite.

2.7. Proof of Theorem 10 on the large-N limit. Using the Lieb-Thirring inequality (14), we have already seen in (i) in Lemma 12 that

J(N ) ≥ e LT (d, p) N, ∀N ∈ N.
In particular, we see that J(N )/N is bounded from below.

In Lemma 12 we have also shown in (ii) that J(N ) ≤ J(N -K) + J(K) for every integer K = 1, ..., N -1. We deduce first that J(N )/N ≤ J(1) = I(1) < 0 which appears in the upper bound (23). The first inequality is in fact strict when p < 2, by Corollary 22. Since the function N → J(N ) is subadditive and bounded from below by a constant times N we conclude that J(N )/N converges to its infimum e(d, p) as in ( 22), by Fekete's subadditive lemma [START_REF] Fekete | Über die Verteilung der Wurzeln bei gewissen algebraischen Gleichungen mit ganzzahligen Koeffizienten[END_REF].

It remains to prove the upper bound on e(d, p) in terms of the semi-classical constant c sc (d). Instead of using the Dirichlet eigenfunctions of a large domain Ω, as we mentioned after Theorem 10, we rather localise the periodic eigenfunctions in a cube (plane waves), which gives more explicit formulas. Let C L be the cube of side length L centered at the origin and let

χ ∈ C ∞ c (R d , R + ) be such that ´Rd χ = 1. Denote u k (x) = L -d 2 1 C L * χ e -ik•x for k ∈ (2π/L)Z d . It turns out that these functions are orthonormal, since u k , u k ′ = L -d ˆRd 1 C L * χ e i(k-k ′ )•x dx = (2π) d 2 1 C L * χ(k -k ′ ) = 0 due to the fact that 1 C L (ℓ) = 0 for ℓ ∈ (2π/L)Z d \ {0}.
A computation gives that the kinetic energy of each such function is equal to

ˆRd |∇u k (x)| 2 = |k| 2 + 1 L d ˆRd |∇ 1 C L * χ(x)| 2 dx.
The second term is a O(1/L) since the function 1 C L * χ is equal to 1 inside C L , at a distance of order one to its boundary, and vanishes outside at a similar distance. We take N ≈ |C L |ρ * = L d ρ * such functions, each with a different k, with ρ * given by (21) in Lemma 9 and C = c sc (d). We find the energy

E(u k1 , ..., u k N ) = N j=1 k 2 j - N p pL dp ˆRd (1 C L * χ) p =L d +O(L d-1 ) +O N L .
The first term is minimum when we take for the k j all the points of (2π/L)Z d in a ball of fixed radius R, where R is chosen so that there are N points, that is,

R ≈ 1 (2π) dρ * |S d-1 | 1/d
. Taking the limit L → ∞ gives the upper bound in (23).

Using that J is Lipschitz by (iii) in Lemma 12, we infer

J(N + α) N + α - J(N ) N ≤ |J(N + α) -J(N )| N + α + α|J(N )| N (N + α) ≤ α N + α C + e(d, p) -→ N →∞ 0,
for any α ∈ [0, 1]. This proves that J(λ)/λ has the same limit as when it is restricted to integers:

lim λ→∞ J(λ) λ = e(d, p). ( 62
)
Using similar arguments as in the proof of Theorem 4, we can then prove that J(λ)/λ is always strictly above its limit, independently of whether it admits a minimiser or not.

Corollary 24. Let d ≥ 1 and 1 < p < min{2, 1 + 2 d }. For all λ > 0 and all m ∈ N \ {1} we have J(mλ) < mJ(λ). In particular, we have

J(λ) λ > e(d, p) (63)
for all λ > 0.

Proof. If J(λ) has a minimiser, this was already proved in Corollary 22. This covers in particular the case 0 ≤ λ ≤ 1. If J(λ) does not have a minimiser, then, according to Remark 15, there is an integer 1 ≤ k ≤ λ so that J(λ) = J(k) + J(λk). By further decomposing J(k), we can therefore write as in (61),

J(λ) = n∈N k n J(n) + J(λ -k), λ = n∈N k n n + (λ -k), k n ∈ {0, 1},
and at least one k n has value 1. So we have, as before

mJ(λ) = n∈N mk n J(n) + mJ(λ -k) > n∈N k n J(mn) + J(m(λ -k)) ≥ J(mλ).
This concludes the proof of Corollary 24.

3. Application: Symmetry breaking for a crystal in the Kohn-Sham model with large Dirac exchange

In this section we explain how the previous results can be used to prove symmetry breaking for an infinite periodic system, within a simple Kohn-Sham model with a Dirac (a.k.a. Slater) term. The results of this section are similar to a recent work by Ricaud [Ric18] on the Thomas-Fermi-von Weisäcker-Dirac model, so some technical details will be omitted for shortness. The main difference is that we deal with operators instead of functions. Our results can be generalised to other contexts, such as the symmetry breaking in the dissociation of the hydrogen molecule considered in [HHL + 19].

3.1. Notation and main results. Everywhere in this section we fix the dimension d = 3. Let R be a lattice of R 3 , with unit cell denoted by K and dual lattice by R * . We consider the infinite system obtained by placing one point nucleus of charge Z = N at each site of this lattice, 3 together with an infinite sea of quantum electrons in a periodic state. We assume that the system is locally neutral, which means that the number of electrons per unit volume must be equal to N |K| -1 . Our goal is to determine whether these electrons will have the same periodicity R as the lattice of the nuclei or whether it is more favourable energetically to place them with a 3 More generally we could place several nuclei of charges z 1 , ..., z M in each unit cell, so that the total charge is M m=1 zm = N .

different period. In the latter case we say that there is spatial symmetry breaking.

More specifically, we will study whether the (ℓR)-periodic electronic ground state is R-periodic or not, for ℓ ≥ 2. We recall that an R-periodic density matrix γ is a self-adjoint operator 0 ≤ γ = γ * ≤ 1 on L 2 (R 3 ) (we neglect the spin for simplicity) which commutes with all the translations of the lattice R:

∀R ∈ R, τ R γ = γτ R .
Here τ R is the unitary operator on L 2 (R 3 ) defined by (τ R f )(x) := f (x -R). We restrict ourselves to density matrices which have a finite trace and a finite kinetic energy per unit volume, which means that γ and √ -∆γ √ -∆ are locally traceclass. The density of γ is the unique R-periodic function

ρ γ ∈ L 1 loc (R 3 , R + ) such that Tr(χγχ) = ˆR3 χ(x) 2 ρ γ (x) dx
for every χ ∈ L ∞ (R 3 ) of compact support. Any such density matrix γ represents an infinite periodic system of electrons. The number of electrons in each unit cell is defined by

Tr R (γ) := Tr (1 K γ1 K ) = ˆK ρ γ (x) dx.
In this section we work with electronic density matrices which are (ℓR)-periodic for some ℓ ≥ 1. All the previous definitions are easily extended to the case ℓ ≥ 2. Our main goal is to determine whether an (ℓR)-periodic minimiser is necessarily R-periodic or not. In what follows, we enforce neutrality of the system. So, in any supercell of the type ℓK, we impose

Tr ℓR (γ) = ˆℓK ρ γ = ℓ 3 N.
The functional to minimise is the Kohn-Sham energy per unit cell which is defined by

E KS c,ℓ (γ) := Tr ℓR (-∆γ) -N ˆℓK G R (x)ρ γ (x) dx + 1 2 D ℓR (ρ γ , ρ γ ) - 3c 4 ˆℓK ρ γ (x) 4 3 dx (64)
for any (ℓR)-periodic density matrix γ. The first term is the kinetic energy per unit cell ℓK, interpreted in the sense of quadratic forms. The second term is the interaction between the (ℓR)-periodic electrons and the lattice R of the nuclei of charge Z = N . The function G R is the R-periodic Green's function, solution to the periodic Laplace equation (66) The existence of minimisers easily follows from the direct method of the calculus of variations, since the problem is posed on the compact set ℓK [START_REF] Catto | On the thermodynamic limit for Hartree-Fock type models[END_REF][START_REF] Cancès | A new approach to the modelling of local defects in crystals: the reduced Hartree-Fock case[END_REF]. An R-periodic state is of course (ℓR)-periodic and its Kohn-Sham energy is found to be equal to E KS c,ℓ (γ) = ℓ 3 E KS c,1 (γ). In particular we deduce that

-∆G R = 4π R∈R δ R -|K| -1 and ˆK G R = 0. ( 65
E KS (c, ℓ) ≤ ℓ 3 E KS (c, 1)
for every ℓ ∈ N and every c ≥ 0.

Definition 25 (Symmetry breaking). We say that there is spatial symmetry breaking for the Dirac-Kohn-Sham model with parameter c ≥ 0 if there exists ℓ ∈ N \ {1} such that

E KS (c, ℓ) ℓ 3 < E KS (c, 1).
The definition means that an ℓ-periodic minimiser has a lower energy per unit volume than the 1-periodic state. It does not mean that the electrons will necessarily be in this ℓ-periodic state. But at least we can deduce that they will not be 1-periodic.

The case c = 0 is studied at length in [START_REF] Catto | On the thermodynamic limit for Hartree-Fock type models[END_REF] and in [CDL08a, App. A]. In this situation the energy γ → E KS 0,ℓ (γ) is convex and the problem E KS (0, ℓ) admits a unique minimiser γ 0,ℓ for every ℓ ≥ 1. This state solves the nonlinear equation

γ 0,ℓ = 1 -∆ -N G R + ρ γ 0,ℓ * G ℓR ≤ ε 0,ℓ (67)
where ε 0,ℓ is a Lagrange multiplier chosen to enforce the constraint that γ 0,ℓ has N ℓ 3 electrons per unit cell. In addition, it is unique in the sense that any (ℓR)periodic solution to equation (67) for some ε 0,ℓ with the right number of electrons N ℓ 3 must be equal to γ ℓ . Since the R-periodic state γ 0,1 with ℓ = 1 is a solution for all ℓ ≥ 2, it follows that γ 0,ℓ = γ 0,1 and ε 0,ℓ = ε 0,1 , for all ℓ ∈ N, when c = 0 and therefore that E KS (0, ℓ) = ℓ 3 E KS (0, 1), for all ℓ ∈ N, when c = 0.

No symmetry breaking occurs for c = 0. For latter purposes, we mention that the system is called an insulator when ε 0,1 can be chosen in a spectral gap of the operator -∆ -N G R + ρ γ0,1 * G R in (67) and that it is a metal otherwise. Which of the two cases occurs depends on the shape of the lattice R and on the number of particles N per unit cell. The Dirac term is not convex when c > 0. It is natural to expect that symmetry will not be broken for c small enough whereas it could be broken for large c. This is confirmed by the following result.

Theorem 26 (Occurrence of symmetry breaking). Let R be a lattice in R 3 and N ∈ N. There is a critical c * = c * (R, N ) ∈ [0, ∞) such that, for all c > c * , the system breaks spatial symmetry. In addition, if the system is insulating at c = 0, then c * > 0: there exists c * 0 > 0 such that E KS (ℓ, c) = ℓ 3 E KS (1, c) for all ℓ ∈ N and all c < c * 0 .

The spirit of the result is exactly the same as [START_REF] Ricaud | Symmetry breaking in the periodic Thomas-Fermi-Dirac-von Weizsäcker model[END_REF] in the Thomas-Fermivon Weisäcker-Dirac case. For c very large the kinetic energy and the Dirac term dominate, the other terms being of lower order. The very large constant c has the effect of concentrating the electrons at the scale 1/c. After rescaling length by a factor 1/c about a blow-up point, in the limit the problem converges to the fermionic NLS problem in the whole space with p = 4/3. This is the content of the following result.

Proposition 27 (Convergence to the NLS problem in R 3 ). Let R be a lattice in R 3 and N ∈ N. For all ℓ ∈ N, we have

lim c→∞ E KS (c, ℓ) c 2 = J(ℓ 3 N )
where J(ℓ 3 N ) is defined as in (10) with d = 3 and p = 4/3.

In Theorem 4 and in Corollary 24 we have proved that J(ℓ 3 N ) < ℓ 3 J(N ) for all ℓ ∈ N \ {1}. This shows that for c large enough E KS (c, ℓ) < ℓ 3 E KS (c, 1), hence that there is symmetry breaking. The intuitive picture is that it is more favourable to concentrate ℓ 3 N particles at one point rather than having ℓ 3 bumps of N concentrated electrons, as is the case for the R-periodic minimiser placed in the (ℓR)-periodic energy. This is how we can prove the first part of Theorem 26.

Remark 28. The previous result does not use that p = 4/3 < p c (3). However, since this inequality has been numerically found to hold (see Remark 23), minimisers for J(N ) always exist, and one can say more. Following the approach of [START_REF] Ricaud | Symmetry breaking in the periodic Thomas-Fermi-Dirac-von Weizsäcker model[END_REF], it is possible to prove that minimisers γ c for E KS (c, 1) satisfy U c γ c U c -1 ⇀ γ weakly- * locally in the trace class, where γ minimises J(N ) and U c is the dilation operator defined by (U c f )(x) = c 3/2 f (cx). In other words, the electrons concentrate at the origin where the nucleus is placed, in the unit cell K. For E KS (c, ℓ) the result is similar but the ℓ 3 N electrons concentrate at one of the ℓ 3 nuclei of the larger unit cell ℓK. Finally, we have the expansion

E KS (c, ℓ) = J(ℓ 3 N )c 2 + c min γ min. for J(ℓ 3 N ) - ˆR3 ρ γ (x) |x| dx + 1 2 ¨R3 ×R 3 ρ γ (x)ρ γ (y) |x -y| dx dy + o(c).
The literature contains several results in the same spirit as Theorem 26 and Proposition 27. The closest to our work is [HHL + 19] which studies the case of the Kohn-Sham hydrogen molecule (two electrons in the field of two nuclei separated by a distance R). This corresponds to N = 2 but since the spin is taken into account there is no orthogonality constraint between u 1 and u 2 . In this model, spin symmetry breaking arises in the limit c → ∞ because each electron has to concentrate about one of the two nuclei, where it asymptotically solves the NLS problem I(1). This is therefore a completely different phenomenon from this present work, where the two particles concentrate at the same point. Other works in the same spirit include for instance [AFG + 02, GS14, GZZ16] for the Hartree model in multiple well potentials, which also has no orthogonality constraint.

In the next section we outline the proof of Proposition 27 whereas in Section 3.3 we quickly discuss the absence of symmetry breaking for c small enough, under the additional assumption that the system is an insulator at c = 0 (second part of Theorem 26).

3.2. Sketch of the proof of Proposition 27 and of the first part of Theorem 26. The symmetry breaking stated in the first part of Theorem 26 follows immediately from Proposition 27, the proof of which we outline in this section.

We set for simplicity ℓ = 1 (the proof is similar in the general case). For γ an R-periodic density matrix, we call γ c = U c γU c -1 the rescaled operator whose kernel is γ c (x, y) := c -3 γ(x/c, y/c).

Using that G cR (x) = c -1 G R (c -1 x), we obtain the following scaling relations:

Tr cR (γ c ) = Tr R (γ), ˆcK G cR ρ c = 1 c ˆK G R ρ, D cR (ρ c , ρ c ) = 1 c D R (ρ, ρ), Tr cR (-∆γ c ) = 1 c 2 Tr R (-∆γ), ˆcK ρ 4/3 c = 1 c 2 c ˆK ρ4/3
. We deduce that the energy of γ can be re-expressed as

E KS c,1 (γ) := c 2 E cR (γ c ) + cF cR (γ c ), (69)
with

E cR (γ c ) := Tr cR (-∆γ c ) - 3 4 ˆcK ρ 4/3 c , F cR (γ c ) := -N ˆcK G cR ρ c + D cK (ρ c , ρ c ).
The energy E cR is similar to the NLS energy E in (31) except that the problem is restricted to the flat torus of size c, instead of being posed over the whole of R 3 .

Step 1. Let us first prove that

lim sup c→∞ E KS (c, 1) c 2 ≤ J(N ) (70)
Let γ be a smooth rank-N projector of compact support such that 4 E(γ) ≤ J(N )+ε. This state can be used as a trial state in the rescaled box cK, as soon as its support is strictly included in cK. This amounts to (cR)-periodising γ in the manner R∈R τ * cR γτ cR . Then E cR (γ) = E(γ) whereas

lim c→∞ F cR (γ) = -N ˆR3 ρ γ (x) |x| dx + 1 2 ¨R3 ×R 3 ρ γ (x)ρ γ (y) |x -y| dx dy. Hence E KS (c, 1) ≤ c 2 (J(N ) + ε) + O(c).
The claimed bound (70) follows after taking c → ∞ and then ε → 0.

Step 2. To prove the other inequality

lim inf c→∞ E KS (c, 1) c 2 ≥ J(N ) (71)
we consider a minimiser γc for E KS (c, 1) and call γ c the rescaled operator as in (68).

From the previous step and the positivity of the Hartree term, we have for c large enough

Tr R (-∆ γc ) -N ˆK G R ρc - 3c 4 ˆK ρc 4/3 ≤ c 2 J(N ) 2 . ( 72
)
Using the Gagliardo-Nirenberg and Hoffmann-Ostenhof periodic inequalities, we have

ˆK ρc 4 3 ≤ C 1 N 5 6 ˆK ρ c + ˆK |∇ √ ρ c | 2 1 2 ≤ C 2 N 5 6 N + Tr R (-∆ γc ) 1 2 .
4 To obtain such a γ one can start with a trial state γ = N i=1 |u i u i | and then truncate and regularise the u i 's. The new functions can be orthonormalised using the same procedure as in Lemma 12.

Similarly, to control the potential energy, we use that G R ≤ |x| -1 + C, and obtain by Hardy's inequality

ˆK G R ρc ≤ C 2 √ N N + Tr R (-∆ γc ) 1 2 .
Inserting in (72) this gives Tr(-∆ γc ) = O(c 2 ) and hence after scaling we obtain

Tr(-∆γ c ) = O(1), ˆcK G cR ρ c = O(1).
This gives

E KS (c, 1) ≥ c 2 E cR (γ c ) + O(c). The last step is to show that lim inf c→∞ E cR (γ c ) ≥ J(N ). ( 73
)
To prove (73) we decompose γ c into bubbles. We use the operator version of the bubble decomposition, which has implicitly appeared several times in the literature and can be read with full details in the recent work [HKY19, Theorem 3.1]. The present setting is slightly different from [START_REF] Hong | Global existence versus finite time blowup dichotomy for the system of nonlinear Schrödinger equations[END_REF] due to the periodic boundary condition but the proof is similar, see, e.g., [START_REF] Ricaud | Symmetry breaking in the periodic Thomas-Fermi-Dirac-von Weizsäcker model[END_REF] in the case of functions. For operators the result is that there exists a sequence of density matrices {γ

(1) , γ (2) , • • • } over R 3 with Tr(-∆γ (i) ) < ∞ such that N ≥ i Tr(γ (i) ), lim inf c→∞ Tr cR (-∆γ c ) ≥ i Tr(-∆γ (i) ) and lim c→∞ ˆcK ρ 4 3 γc = i ˆR3 ρ 4 3 γ (i) .
The Dirac term decomposes exactly since 4/3 is a sub-critical power, whereas for the mass and the kinetic energy one only obtains lower bounds. The missing mass and kinetic energy are contained in the vanishing part of γ c , to employ the vocabulary of the concentration-compactness method. Each γ (i) is constructed as the strong local limit of χ i,c (• + x i,c )γ c χ i,c (• + x i,c ) for some translation x i,c and some localisation function χ i,c , with |x i,cx j,c | → ∞ when i = j, up to subsequences. Using the subadditivity of J proved in Lemma 12, we deduce that

lim inf c→∞ E cR (γ c ) ≥ i E(γ (i) ) ≥ i J Tr(γ (i) ) ≥ J i Tr(γ (i) ) ≥ J(N ). ( 74
)
This concludes our sketch of the proof of Proposition 27.

3.3. Proof of the second point of Theorem 26: stability for small c. When c = 0, we have recalled from [CDL08a, App. A] that the minimisation problem E KS (0, 1) admits a unique minimiser, which we denote here by γ 0 (it was called γ 0,1 above). It solves the nonlinear operator equation

γ 0 = 1(H 0 ≤ ε 0 ) (75)
where ε 0 is a Lagrange multiplier and

H 0 := -∆ -N G R + ρ 0 * G R .
The assumption that the system is an insulator means that ε 0 belongs to a spectral gap of the operator H 0 . For simplicity, we denote by a := max σ(H 0 ) ∩ (-∞, ε 0 ), b := min σ(H 0 ) ∩ (ε 0 , ∞)

and, without loss of generality, we can choose

ε 0 = 1 2 (a + b).
The length of the gap is g := ba > 0. Let us prove that ρ 0 > 0. After a Bloch-Floquet transform [START_REF]Methods of Modern Mathematical Physics. IV. Analysis of operators[END_REF], the equation (75) means that the kernel of γ 0 is given by γ 0 (x, y) = n≥0 ˆB 1(λ n (ξ) ≤ ε 0 ) u n (ξ, x)u n (ξ, y) dξ with the density

ρ 0 (x) = n≥0 ˆB 1(λ n (ξ) ≤ ε 0 ) |u n (ξ, x)| 2 dξ.
Here B is the Brillouin zone (the unit cell of the dual lattice R * ) and (u n , λ n ) are the Bloch eigenfunctions and (ordered) eigenvalues, which solve

| -i∇ x + ξ| 2 -N G R + ρ 0 * G R u n (ξ, •) = λ n (ξ) u n (ξ, •)
with periodic boundary conditions on ∂K. By Perron-Frobenius we have u 0 > 0 and λ 0 (0) < λ 1 (0). By perturbation theory we then deduce that λ 0 (ξ) is nondegenerate with a positive eigenfunction u 0 (ξ, •), for ξ small enough. Then ρ 0 > 0 and in the following we denote by α := min K ρ 0 > 0 the minimal value of the periodic density.

The following shows that the gap does not close and the density stays strictly positive for c small enough.

Lemma 29 (Stability of the gap).

There is c 1 > 0 such that, for all 0 ≤ c < c 1 , any minimiser γ c for E KS (c, 1) satisfies

min K ρ c > α 2
, and dist (σ(H c ), ε 0 ) > g 4 ,

where we set ρ c := ρ γc , and

H c := -∆ -N G R + ρ c * G R -cρ 1/3
c . Finally, we have γ c = 1(H c < ε 0 ) and there is C > 0 independent of c < c 1 so that the following operator inequality holds:

C -1 (1 -∆) ≤ |H c -ε F | ≤ C(1 -∆). ( 76
)
Proof. Let c n → 0 + . The energy c → E KS (c, 1) is continuous at c = 0 and any minimiser γ cn is a minimising sequence for E KS (0, 1). Hence it must converge to the unique minimiser γ 0 weakly and lim cn→0 + Tr R (-∆γ cn ) = Tr R (-∆γ 0 ). This implies that ρ cn → ρ 0 strongly in L 1 ∩ L 3 (K). Then we write the associated mean-field operator in the form

H cn = -∆ -N G R + ρ cn * G R -c n ρ 1/3 cn = H 0 + (ρ cn -ρ 0 ) * G R -c n ρ 1/3
cn and estimate the operator norms of the last two terms by

(ρ cn -ρ 0 ) * G R (1 -∆) -1 ≤ ||ρ cn -ρ 0 || L 1 (K) G R (1 -∆) -1 → 0 and ρ 1/3 cn (1 -∆) -1 ≤ C ||ρ cn || 1/3 L 1 (K)
. With similar estimates we know that (H 0 + C)(1 -∆) -1 and (H 0 + C) -1 (1 -∆) are bounded for C large enough, locally uniformly in c, see [START_REF] Cancès | A new approach to the modelling of local defects in crystals: the reduced Hartree-Fock case[END_REF]Lem. 1]. By the Rellich-Kato theorem, this proves that the spectrum of H cn converges to that of H 0 . In particular, H cn has a gap around ε 0 , independent of c n for c n small enough. To conclude we have therefore shown that there exists c 1 > 0 so that any minimiser γ c for E KS (c, 1) has a mean-field operator H c with the gap g/4 around ε 0 . This implies (76) by [START_REF] Cancès | A new approach to the modelling of local defects in crystals: the reduced Hartree-Fock case[END_REF]Lem. 3].

Let then γ c be any such minimiser for c < c 1 . Since the family H(t) = H 0 + t(ρ cρ 0 ) * G Rtcρ 1/3 c has a gap for all t ∈ [0, 1] and the rank of a continuous family of orthogonal projectors is always constant, we obtain

Tr R 1(H c ≤ ε 0 ) = N.
By [CDL08a, App. A] we know that γ c = 1(H c ≤ ε c ) where ε c is the unique Lagrange multiplier chosen such that Tr(γ c ) = N , and we conclude that ε c = ε 0 is independent of c. In particular γ c = 1(H c ≤ ε 0 ).

Finally, we have

Tr(C + H c )γ c (C + H c ) ≤ (C + ε 0 ) 2 N.
Since (H 0 + C)(1 -∆) -1 and (H 0 + C) -1 (1 -∆) are bounded, this shows that Tr(1 -∆)γ c (1 -∆) ≤ C uniformly in c < c 1 . This implies that ρ c is bounded in W 2,1 (K) and therefore we have ρ c → ρ 0 in L ∞ (K) when c → 0 + . In particular ρ c ≥ α/2 > 0 for c small enough.

Next we use the properties of minimisers for c < c 1 in Lemma 29 to show that there is indeed only one, for every ℓ ≥ 1.

Lemma 30. There is c 2 > 0 so that, for all 0 ≤ c < c 2 , E KS (c, 1) has a unique minimiser γ c , satisfying the properties of Lemma 29. This minimiser is also the unique minimiser for E KS (c, ℓ) for all ℓ ≥ 1, hence there is no symmetry breaking for c < c 2 .

Proof. We use the framework developed in [BBHS99, [START_REF] Hainzl | Existence of a stable polarized vacuum in the Bogoliubov-Dirac-Fock approximation[END_REF][START_REF] Cancès | A new approach to the modelling of local defects in crystals: the reduced Hartree-Fock case[END_REF][START_REF]Non-perturbative embedding of local defects in crystalline materials[END_REF][START_REF] Frank | A positive density analogue of the Lieb-Thirring inequality[END_REF]. Let γ c be any minimiser for E KS (c, 1) with c < c 1 and let γ be any other (ℓR)-periodic density matrix. Using Tr ℓR γ = Tr ℓR γ c = ℓ 3 N , we can rewrite and estimate the difference of the two energies as

E KS c,ℓ (γ) -E KS c,ℓ (γ c ) = Tr ℓR (H c -ε 0 )Q + 1 2 D ℓR (ρ Q , ρ Q ) - 3 4 c ˆℓK (ρ c + ρ Q ) 4 3 -ρ 4 3 c - 4 3 ρ 1 3 c ρ Q ≥ Tr ℓR (H c -ε 0 )Q + 1 2 D ℓR (ρ Q , ρ Q ) -cK ˆℓK min ρ 2 Q , ρ 4 3 Q ( 77
)
where Q := γγ c . In the second line we have used that

(1 + t) 4/3 -1 -4 3 t ≤ C min(t 4/3 , t 2 ) for all t ≥ -1 and that ρ c ≥ α/2. Our goal is to show that (77) is non-negative and vanishes only at γ = γ c . We claim that ˆℓK min ρ 2 Q , ρ

4 3 Q ≤ CTr ℓR (H c -ε 0 )Q. ( 78
)
The result then follows under the assumption that c < c 2 := min(c 1 , (2CK) -1 ).

To prove (78) we introduce

Q --:= γ c Qγ c , Q -+ := γ c Q(1 -γ c ), Q +-:= (1 -γ c )Qγ c , Q ++ := (1 -γ c )Q(1 -γ c )
and note that

Tr ℓR (H c -ε 0 )Q = Tr ℓR |H c -ε 0 |(Q ++ -Q --) ≥ CTr ℓR (1 -∆)(Q ++ -Q --) ≥ CTr ℓR (1 -∆)Q 2 .
We have used Bach's inequality Q 2 ≤ Q ++ -Q --from [BBHS99, Eq. ( 18)-( 19)].

For q = Q ++ , Q --we use the Lieb-Thirring inequality which implies

Tr ℓR (1 -∆)q ≥ C ˆℓK ρ q + ρ 5 3 q ≥ 2C ˆℓK ρ 4 3 q and provides the desired bound on the two densities ρ Q ++ and ρ Q --. For Q +-and Q -+ the argument is slightly more involved. Following [CDL08a, Prop. 1] we claim that

ˆℓK ρ 2 Q +-+ ρ 2 Q -+ ≤ C Tr ℓR Q 2 (79)
where the constant C is independent of ℓ. The argument goes by duality in the form ˆℓK

ρ Q +-V = Tr ℓR γ c V (1 -γ c )Q ≤ ||Q|| S 2 (L 2 (ℓK)) ||γ c V || S 2 (L 2 (ℓK)) ≤ ||Q|| S 2 (L 2 (ℓK)) ||γ c (1 -∆)|| (1 -∆) -1 V S 2 (L 2 (ℓK)) ,
with S 2 (H) the Hilbert-Schmidt norm on a Hilbert space H. We have

||γ(1 -∆)|| ≤ ||γ(H c -i)|| (H c + i) -1 (1 -∆) ≤ K,
for a constant K independent of c. We obtain an upper bound involving

(1 -∆) -1 V 2 S 2 (L 2 (ℓK)) . We compute this Hilbert-Schmidt norm in the Fourier basis e k (x) := ℓ -3/2 e ik•x , which gives

(1 -∆) -1 V 2 S 2 (L 2 (ℓK)) = k1,k2∈R * /ℓ | e k1 , (1 -∆) -1 V e k2 | 2 = k1,k2∈R * /ℓ 1 (1 + |k 1 | 2 ) 2 | e k1 , V e k2 | 2 = 1 ℓ 3 k∈R * /ℓ 1 (1 + |k| 2 ) 2 ˆℓK V 2 ≤ C ˆℓK V 2 .
This concludes our sketch of the proof of (79), hence of (78) and of Theorem 26.

N

  n=1 |u n | 2 and ρ ′ = M m=1 |v m | 2 , the decay of which is dictated by the last eigenfunctions u N and v M which have the smallest eigenvalues µ N and µ ′

  22) exists and satisfies e LT (d, p) ≤ e(d, p) ≤ min e sc (d, p), I(1) < 0 (23)

Figure 1 .

 1 Figure 1. Numerical computation of the density ρ N = N i=1 |u i | 2 of the minimiser for J(N ) in dimension d = 1 with N = 15 and p = 1.3. The system exhibits a crystallised phase with 15 local maxima.

Figure 2 .

 2 Figure 2. Numerical computation of the density ρ N = N i=1 |u i (x)| 2 of the minimiser of J(N ) for N varying between 1 and 7 in dimension d = 2 with p = 1.5.

Figure 3 .

 3 Figure 3. Plot of J(N )/N as a function of N in dimension d = 2for p = 1.5. It is unclear whether J(7)/7 is smaller than J(6)/6, due to numerical errors.

Figure 4 .

 4 Figure 4. Numerical computation of the functions λ → J(λ) (left) and λ → J(λ)/λ (right) in dimension d = 1 with p = 1.3. The function has the behaviour described in Remark 13.

  the leading term in the first sum is the exponential involving |µ N | and (46) follows.

  Proof of Lemma 21. Let us first bound e R . By Lemma 19 we have |u i | ≤ CY mi and |v j | ≤ CY m ′ j with m i = |µ i | and m ′ j = |µ ′ j |. Thus we find

Remark 23 (

 23 Numerical evaluation of p c (d)). Any lower bound on |I(d, p, 1)| and on the Lieb-Thirring constant c LT (d) appearing in the definition (24) of e LT (d, p) yields a lower bound on p c (d). The NLS energy I(d, p, 1) can easily be computed numerically to a high precision since this amounts to solving an ordinary differential equation for the radial function Q. Using a Runge-Kunta solver for Q and the recent lower bound c LT (d) ≥ (0.471851) 1 d c sc (d) (58) from [FHJN19, Prop. 10] (see also [Fra20, Thm. 5]), we obtained the lower bounds mentioned in (13) for p c (d) in dimensions d = 1, 2, 3. Should c LT (3) be equal to its conjectured value c sc (3), we would obtain the slightly better bound p c (3) ≥ 1.494.

)

  In other words, N G R is the Coulomb potential of the infinite lattice of nuclei, screened by a uniform background. The third term in (64) is the Coulomb interaction between the electrons in the Hartree approximation and it readsD ℓR (f, g) := ˆℓK ˆℓK G ℓR (xy)f (x)g(y)dx dywhere G ℓR is defined similarly as in (65) with R replaced by ℓR and K replaced by ℓK. Finally, the last term of (64) is the Dirac or Slater term. This term is an approximation of the exchange-correlation energy of γ, in terms of the density ρ γ only. The parameter c ≥ 0 usually has a fixed value given by physical considerations (for the exchange part of the energy the constant is c ≃ 1.24 without spin [LS10, Sec. 6.2]). Here, we change the value of c, and compare the resulting energies for different values of the periodicity ℓ of the electrons.The minimisation problem to be considered reads E KS (c, ℓ) := min E KS c,ℓ (γ), γ (ℓR)-periodic density matrix, ˆℓK ρ γ = ℓ 3 N .

  More recently, the Strichartz inequality has been extended to orthonormal systems in [FLLS14, FS17, BHL + 17] and it has played a central role for the existence and the long time behaviour of infinite systems [LS15, LS14, CHP17a, CHP17b, Cd18].

  The first term on the right behaves like Y Mi = o(Y mi ) at infinity, since |∇Y mi |1 B1 is integrable and has compact support. For the second term we remark that the Fourier transform of Y mi * Y Mi is proportional to 1

This inequality also implies J ′ (λ) -= lim t→0 + J(λ-t)-J(λ) -t ≥ µ N , see Remark 13.

Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org)

Acknowledgement. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement MDFT 725528 of M.L.).