Mathieu Lewin
email: mathieu.lewin@math.cnrs.fr

Séverine Paul
email: severine.paul@u-cergy.fr

A Numerical Perspective on Hartree-Fock-Bogoliubov Theory

The method of choice for describing attractive quantum systems is Hartree-Fock-Bogoliubov (HFB) theory. This is a nonlinear model which allows for the description of pairing effects, the main explanation for the superconductivity of certain materials at very low temperature.

This paper is the first study of Hartree-Fock-Bogoliubov theory from the point of view of numerical analysis. We start by discussing its proper discretization and then analyze the convergence of the simple fixed point (Roothaan) algorithm. Following works by Cancès, Le Bris and Levitt for electrons in atoms and molecules, we show that this algorithm either converges to a solution of the equation, or oscillates between two states, none of them being a solution to the HFB equations. We also adapt the Optimal Damping Algorithm of Cancès and Le Bris to the HFB setting and we analyze it.

The last part of the paper is devoted to numerical experiments. We consider a purely gravitational system and numerically discover that pairing always occurs. We then examine a simplified model for nucleons, with an effective interaction similar to what is often used in nuclear physics. In both cases we discuss the importance of using a damping algorithm.

Introduction

Hartree-Fock-Bogoliubov (HFB) theory is the method of choice for describing some special features of attractive fermionic quantum systems [START_REF] Ring | The nuclear many-body problem[END_REF]. It is a generalization of the famous Hartree-Fock (HF) method [START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF] used in quantum chemistry. It also generalizes the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity [START_REF] Bardeen | Theory of superconductivity[END_REF] which was invented in 1957 to explain the complete loss of resistivity of certain materials at very low temperature. In 1958 Bogoliubov realized in [START_REF] Bogoliubov | On a New Method in the Theory of Superconductivity[END_REF] that the BCS theory was actually very similar to his previous works [START_REF] Bogoliubov | About the theory of superfluidity[END_REF][START_REF]On the theory of superfluidity[END_REF][START_REF]Energy levels of the imprefect Bose gas[END_REF] on the superfluidity of certain bosonic systems. He adapted it to fermions, leading to a model that is now called Hartree-Fock-Bogoliubov, and which is the main interest of this article.

In Hartree-Fock-Bogoliubov theory the state of the system is completely determined by two operators [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF]. The one-particle density matrix γ is the same as in Hartree-Fock theory [START_REF] Lieb | The Stability of Matter in Quantum Mechanics[END_REF], whereas the pairing density matrix α describes the Cooper pairing effect. This effect can only hold in attractive systems. When the interaction potential is positive, the energy is decreased by replacing α by 0.

One of the most interesting questions in HFB theory is precisely the existence of pairing, that is, the non-vanishing of the matrix α for a minimizer. This question has been settled in the simpler translation-invariant BCS theory [START_REF] Billard | An existence proof for the gap equation in the superconductivity theory[END_REF][START_REF] Vansevenant | The gap equation in superconductivity theory[END_REF][START_REF] Yang | On the Bardeen-Cooper-Schrieffer integral equation in the theory of superconductivity[END_REF][START_REF] Mcleod | The uniqueness and approximation of a positive solution of the Bardeen-Cooper-Schrieffer gap equation[END_REF][START_REF] Frank | The critical temperature for the BCS equation at weak coupling[END_REF][START_REF] Hainzl | The BCS functional for general pair interactions[END_REF][START_REF]The BCS critical temperature for potentials with negative scattering length[END_REF] and in some cases for the translation-invariant Hubbard model [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF]. But it remains completely open for general attractive systems with a few particles, like those encountered in nuclear physics. In [START_REF] Lenzmann | Minimizers for the Hartree-Fock-Bogoliubov theory of neutron stars and white dwarfs[END_REF], the first author of this paper has shown with Lenzmann the existence of HFB minimizers for a purely Newtonian system of N fermions, but it is not yet known if α = 0. One of the purpose of this work is to answer this question numerically.

The HFB energy is a nonlinear function of γ and α. Because of nonlinearity, minimizing this functional on a computer is not an easy task. The simpler Hartree-Fock model in which α = 0 is now well understood from the point of view of numerical analysis [START_REF] Bris | Computational chemistry from the perspective of numerical analysis[END_REF][START_REF] Cancès | Computational quantum chemistry: a primer[END_REF], even though most authors have concentrated their attention to the special case of electrons in an atom or a molecule. Cancès and Le Bris have studied in [START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF][START_REF] Cancès | Can we outperform the DIIS approach for electronic structure calculations?[END_REF] the simpler fixed point algorithm called the Roothaan algorithm [START_REF] Roothaan | New developments in molecular orbital theory[END_REF] and they have shown that this algorithm either converges or oscillates between two points, none of them being the solution of the HF equation. This result was recently improved by Levitt [START_REF] Levitt | Convergence of gradient-based algorithms for the Hartree-Fock equations[END_REF]. Cancès and Le Bris have also proposed a new algorithm called Optimal Damping, which is now used in several chemistry programs. It is based on the fact that one can freely minimize the energy over mixed HF states instead of pure HF states, by Lieb's variational principle [START_REF] Lieb | Variational principle for many-fermion systems[END_REF].

Because the HFB model is an extension of HF theory, it is natural to believe that these ideas can be applied to the HFB case as well. In particular, under appropriate assumptions on the interaction potential, Bach, Fröhlich and Jonsson have recently shown in [START_REF] Bach | Bogolubov-Hartree-Fock mean field theory for neutron stars and other systems with attractive interactions[END_REF] an HFB equivalent of Lieb's variational principle. Up to some difficulties that will be explained later, we will show in this paper that the previously mentioned results can indeed be transposed to the Hartree-Fock-Bogoliubov model.

The paper is organized as follows. In the next section we quickly recall the basic formulation of Hartree-Fock-Bogoliubov theory. Then, in Section 3, we derive the discretized HFB equations and we prove that, in the limit of a large Galerkin basis set, the discretized solution converges to the true solution. We also discuss at length the possible symmetries of the system and we formulate the theory when these symmetries are taken into account.

In Section 4 we study the HFB Roothaan algorithm and we prove that it either converges to a solution of the HFB equation, or oscillate between two points, none of them being a solution of the equations. We then introduce an equivalent of the Optimal Damping Algorithm of Cancès and Le Bris, which is based on an optimization in the set of mixed HFB states.

Section 5 is devoted to the presentation of some numerical results. We first consider the purely gravitational model studied by Lenzmann and Lewin [START_REF] Lenzmann | Minimizers for the Hartree-Fock-Bogoliubov theory of neutron stars and white dwarfs[END_REF] and we numerically discover that there is always pairing. Then, we introduce a simplified model for nucleons, for which we as well present some preliminary numerical results. We particularly discuss the importance of using a damped algorithm instead of a simple fixed point method, a fact which has already been noticed in nuclear physics [START_REF] Dechargé | Hartree-Fock-Bogolyubov calculations with the D1 effective interaction on spherical nuclei[END_REF]. Our approach could help in improving the existing numerical techniques.

Acknowledgement. The authors would like to thank Laurent Bruneau and Julien Sabin for useful discussions. They acknowledge financial support from the French Ministry of Research (ANR-10-BLAN-0101) and from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement MNIQS 258023).

2.

A quick review of Hartree-Fock-Bogoliubov theory 2.1. Hartree-Fock-Bogoliubov states and their energy We consider a system composed of N identical fermions, described by the many-body Hamiltonian

H(N ) = N j=1 T j + 1 k<ℓ N W kℓ , (2.1)
acting on the fermionic N -body space H N = N 1 H, where H is the space for one particle. Here, T : H → H is a one-body operator and W : H 2 → H 2 accounts for the interactions between the particles. We use the notation T j for the operator T which acts on the jth component of the tensor product

H N = N 1 H, that is T j = 1 ⊗ • • • ⊗ T ⊗ • • • ⊗ 1,
and a similar convention for W kℓ .

Most of what follows is valid in an abstract setting. However, for the sake of simplicity, in the whole paper we will restrict ourselves to the special case of nonrelativistic fermions with q internal degrees of freedom, moving in R 3 (q = 2 for spin-1/2 particles like electrons). We also assume that no external force is applied, and that their interaction is translationinvariant. Then, in units where m = 1/2 and = 1, we have

H = L 2 (R 3 , C q ), T = -∆, W kℓ = W (x k -x ℓ ).
The N -body space H N = N 1 L 2 (R 3 , C q ) consists of wave functions Ψ(x 1 , σ 1 , ..., x N , σ N ) which are antisymmetric with respect to exchanges of the variables (x i , σ i ). In principle W (x kx ℓ ) is also a function of the two internal variables σ k , σ ℓ ∈ {1, ..., q} of the particles k and ℓ. Again for simplicity, we will assume that W only depends on the space variable x kx ℓ . Finally, we make the assumption that W is smooth and decays fast enough at infinity to ensure that H is bounded from below. To make this more explicit, we assume in the whole paper that W = W 1 + W 2 ∈ L p (R 3 ) + L q (R 3 ) for some 2 p q < ∞.

(

Sometimes we will make more precise assumptions on W . We are interested in the case where W is attractive (W 0), or at least partially attractive (W 0 on a set of measure non zero). By translation invariance, the Hamiltonian H(N ) has no ground state (that is, the bottom of its spectrum cannot be an eigenvalue). But it may have one once the center of mass is removed, if W is sufficiently negative.

In a nonlinear model approximating the many-body problem above, there could be a ground state, even if the system is translation-invariant. Of course, translation invariance is not lost and there are then infinitely many ground states, obtained by translating the system arbitrarily. In Hartree-Fock theory [START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF][START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF], such breaking of symmetry is known to occur for instance when W (x) = -1/|x| is a purely gravitational interaction and T = -∆ (nonrelativistic), or T = √ 1 -∆ -1 (pseudo-relativistic), see [START_REF] Lenzmann | Minimizers for the Hartree-Fock-Bogoliubov theory of neutron stars and white dwarfs[END_REF][START_REF] Lewin | Geometric methods for nonlinear many-body quantum systems[END_REF] and Theorem 2.3 below.

For attractive systems, it is often convenient to allow for another symmetry breaking, namely that of particle number. This means that the fixed particle number N is replaced by an operator N whose eigenvalues are 0, 1, 2, .... Only the average particle number is well defined for a quantum state. The classical way to define N is to introduce the fermionic Fock space

F = C ⊕ n 1 H n ,
which gathers all the possible n-particle subspaces in a direct sum. A (pure) quantum state in

F is a vector Ψ = ψ 0 ⊕ ψ 1 ⊕ • • • which is normalized in the sense that ||Ψ|| 2 F = |ψ 0 | 2 + n 1 ||ψ n || 2 H n = 1.
The average particle number is the diagonal operator

N := 0 ⊕ n 1 n,
such that the average number of particles in a state Ψ is given by the formula

Ψ, N Ψ = n 1 n | |ψ n || 2 H n .
Instead of imposing that Ψ ∈ H N which is equivalent to Ψ being an eigenvector of N , N Ψ = N Ψ, we will only fix the average particle number of Ψ:

Ψ, N Ψ = N.
Allowing to have ψ n = 0 for n = N is useful to describe some physical properties of attractive systems. In most practical cases it is expected that the variance n 0 (n -N ) 2 ||ψ n || 2 H n will be quite small, i.e. that Ψ will live in a neighborhood of H N .

Similarly to the particle number operator N , the many-body Hamiltonian H(N ) is now replaced by a many-body Hamiltonian H on Fock space

H := 0 ⊕ n 1 H(n) (2.3)
which is nothing else but the diagonal operator which coincides with H(n) on each n-particle subspace. We will not discuss here the problem of defining H as a self-adjoint operator on F .

The Hartree-Fock-Bogoliubov (HFB) model generalizes the well-known Hartree-Fock (HF) method and it allows for breaking of particle number in a very simple fashion. The method consists in restricting the many-body Hamiltonian H on F to a special class of states called Hartree-Fock-Bogoliubov states (or quasi-free states), which are completely characterized by their one-particle density matices [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF].

Let us recall that a state in Fock space has two one-particle density matrices, instead of one in usual HF theory. These are two operators γ : H → H and α : H → H, which are defined by means of creation and annihilation operators by the relations [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF]

Ψ, a † (f )a(g)Ψ F = g, γf H , Ψ, a(f )a(g)Ψ F = g, αf H .
When Ψ lives in a particular N -particle subspace H N , then a(f )a(g)Ψ ∈ H N -2 hence Ψ, a(f )a(g)Ψ F = 0 for all f, g ∈ H and the matrix α vanishes. However for a general state Ψ ∈ F , one can have α = 0. The two operators γ and α satisfy several constraints. First, we have γ * = γ, 0 γ 1 (in the sense of operators) and Tr γ = Ψ, N Ψ = N , for the one-particle matrix γ. On the other hand, the so-called pairing matrix α satisfies α T = -α. Its kernel α(x, σ, x ′ , σ ′ ) = -α(x ′ , σ ′ , x, σ) can thus be seen as an antisymmetric two-body wavefunction in H 2 . It is interpreted as describing pairs of virtual particles, called Cooper pairs.

It is well known that a quantum state Ψ ∈ H N such that (γ Ψ ) 2 = γ Ψ is necessarily a Slater determinant (that is, a Hartree-Fock state). The same is true for states in Fock space. Consider a pair (γ, α) which is such that

Γ 2 = Γ, with Γ := γ α α * 1 -γ (2.4)
on H ⊕ H. Hence we have for instance αα * = γγ 2 . Then there exists a unique state Ψ in F which has γ and α as density matrices. This state has the property that any observable can be computed using only γ and α, by Wick's formula (see Thm 2.3 in [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF]). The quantum states obtained by considering projections Γ are called Hartree-Fock-Bogoliubov states and they generalize usual Hartree-Fock states. When α ≡ 0, then γ = N j=1 |ϕ j ϕ j | is a rank-N projection and the corresponding state is the usual Slater determinant

Ψ = ϕ 1 ∧ • • • ∧ ϕ N = 1 √ N ! det ϕ i (x j , σ j ) .
When α = 0, Ψ can be obtained by applying a Bogoliubov rotation to the vacuum but we will not explain this further. For the present work, we will only need the formula of the total energy, in terms of γ and α:

Ψ, HΨ F = n 0 ψ n , H(n)ψ n H n = Tr(-∆)γ + 1 2 R 3 R 3 W (x -y) ρ γ (x)ρ γ (y) -|γ(x, y)| 2 + |α(x, y)| 2 dx dy := E(γ, α) (2.5)
where ρ γ (x) = Tr C q (γ(x, x)) is the density of particles in the system. The terms in the double integral are respectively called the direct, exchange and pairing terms. Taking α ≡ 0 one recovers the usual Hartree-Fock energy which has been studied by many authors [START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF][START_REF] Lions | Solutions of Hartree-Fock equations for Coulomb systems[END_REF][START_REF] Bach | Error bound for the Hartree-Fock energy of atoms and molecules[END_REF][START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF]. Our main goal in this paper is to investigate the minimization of the more complicated nonlinear functional E(γ, α), when γ and α are submitted to the above constraints, and its numerical implementation. We will show below that the energy E is well defined in an appropriate function space, under our assumption (2.2) on W . Note that the variance of the particle number for a HFB state Ψ in Fock space can be expressed only in terms of α by

Ψ, N -Ψ, N Ψ F 2 Ψ F = n 0 (n -N ) 2 ||ψ n || 2 H n = 2 Tr H (αα * ),
see Lemma 2.7 in [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF]. The spreading of the HFB state over the different spaces H n is therefore determined by the Hilbert-Schmidt norm of the pairing matrix α. We recover the fact that an HFB state has a given particle number if and only if its pairing matrix α vanishes.

Pure vs mixed states

In our previous description, we have only considered pure states in Fock space, that is states given by a normalized vector Ψ ∈ F . For practical purposes, it is very convenient to extend the model to mixed states, which are nothing else but convex combinations of pure states, given by a (many-body) density matrix in

F D = j λ j |Ψ j Ψ j | with λ j 0, j λ j = 1, Ψ i , Ψ j F = δ ij .
The average particle number and energy are then given by the formulas

Tr F N D) = j λ j Ψ j , N Ψ j , Tr F HD) = j λ j Ψ j , HΨ j .
Resorting to mixed states is mandatory at positive temperature, when the equilibrium state of the system will actually always be a mixed state. But it is also very useful at zero temperature, even if the true ground state is a pure state. We will recall later the practical advantages of using mixed HFB states.

Similarly to what we have explained in the previous section, there is a class of mixed Hartree-Fock-Bogoliubov states, which are completely characterized by their one-particle density matrices γ and α. The latter now satisfy the constraint

0 0 0 0 Γ := γ α α * 1 -γ 1 0 0 1 (2.6)
on H ⊕ H, which is nothing else but the relaxation of the constraint Γ 2 = Γ of pure states.

The set of one-particle density matrices of mixed HFB states is therefore a convex set, whose extremal points are the density matrices of pure HFB states. One should remember that the set of mixed HFB states is not the convex hull of HFB pure states, however. The relation between the density matrices (γ, α) and the corresponding HFB states in F is highly nonlinear.

The energy of a mixed HFB state described by the density matrices (γ, α) is given by the same formula (2.5) as for pure states. Hence, minimizing this energy under the relaxed constraint (2.6) is equivalent to minimizing the full quantum energy over all mixed HFB states. The natural question arises whether a minimizer, when it exists, is automatically a pure state. The answer to this question is positive in many situations, as we will see below.

Before turning to the comparison between the minimization among pure and mixed states, we first introduce the variational sets on which the energy is well defined. The sets of all pure and mixed HFB states with finite kinetic energy are respectively given by

P := (γ, α) ∈ S 1 (H) × S 2 (H) : α T = -α, Γ = Γ * = Γ 2 , Tr(-∆)γ < ∞ (2.7)
and

K := (γ, α) ∈ S 1 (H) × S 2 (H) : α T = -α, 0 Γ = Γ * 1 H⊕H , Tr(-∆)γ < ∞ , (2.8)
(the matrix Γ is the one appearing in (2.4)). Here S 1 (H) and S 2 (H) denote the spaces of trace-class and Hilbert-Schmidt operators [START_REF] Simon | Geometric methods in multiparticle quantum systems[END_REF]. The expression Tr(-∆)γ is to be understood in the sense of quadratic forms, that is

Tr(-∆)γ = 3 k=1 Tr(p k γp k ) ∈ [0, +∞], with p k = -i∂ x k .
In practice we want to fix the total average number of particles. For this reason we also define the constrained sets

P(N ) := {(γ, α) ∈ P : Tr γ = N } (2.9)
and

K(N ) := {(γ, α) ∈ K : Tr γ = N } , (2.10)
of pure and mixed states with average particle number N . In practice N is an integer but it is convenient to allow any non-negative real number.

The following lemma says that the energy is a well-defined functional on the largest of the above sets K, and that it is bounded from below on K(N ) for any N 0.

Lemma 2.1 (The HFB energy is bounded-below on K(N )).

When W = W 1 + W 2 ∈ L p (R 3 ) + L q (R 3 ) with 2 p q < ∞, then E(γ, α
) is well defined for any (γ, α) ∈ K. It also satisfies a bound of the form

∀(γ, α) ∈ K, E(γ, α) 1 2 Tr(-∆)γ -C(N ) (2.11)
for some constant C(N ) depending only on N = Tr(γ).

Proof.

The assumption that W = W 1 + W 2 ∈ L p (R 3 )+ L q (R 3 ) with 2 p q < ∞ implies that W is relatively form-bounded with respect to the Laplacian, with relative bound as small as we want [START_REF] Davies | Spectral theory and differential operators[END_REF]. This means |W | ǫ(-∆) + C ǫ in the sense of quadratic forms, for all ǫ > 0 and for some constant C ǫ . This can now be used to verify that the energy is well defined under the assumption that Tr(-∆)γ < ∞. First, we have for the direct term

R 3 R 3 |W (x -y)| ρ γ (x)ρ γ (y) dx dy ǫN R 3 |∇ √ ρ γ | 2 + C ǫ N 2 ǫN Tr(-∆)γ + C ǫ N 2 ,
where in the last line we have used the Hoffmann-Ostenhof inequality [START_REF] Hoffmann-Ostenhof | Schrödinger inequalities and asymptotic behavior of the electron density of atoms and molecules[END_REF],

R 3 |∇ √ ρ γ | 2 Tr(-∆)γ.
(2.12)

The exchange term is bounded similarly by applying the inequality |W | ǫ(-∆) + C ǫ in x with y fixed:

R 3 R 3 |W (x -y)| |γ(x, y)| 2 dx dy ǫ R 3 R 3 |∇ x γ(x, y)| 2 dx dy + C ǫ R 3 R 3 |γ(x, y)| 2 dx dy = ǫ Tr(-∆)γ 2 + C ǫ Tr γ 2 ǫ Tr(-∆)γ + C ǫ N,
since γ 2 γ. Similarly we have, since αα * γγ 2 γ,

R 3 R 3 |W (x -y)| |α(x, y)| 2 dx dy Tr ǫ(-∆) + C ǫ αα * ǫ Tr(-∆)γ + C ǫ N.
All this shows that all the terms in the energy are well defined when (γ, α) ∈ K(N ). Also, we have

E(γ, α) 1 -ǫ -ǫN/2 Tr(-∆)γ -C ǫ N -C ǫ N 2 /2. (2.13)
Taking ǫ = 1/(2 + N ) finishes the proof.

Lemma 2.1 allows us to define the minimization problems for pure and mixed states as follows:

I(N ) := inf (γ,α)∈K(N ) E(γ, α), (2.14)
J(N ) := inf (γ,α)∈P(N ) E(γ, α). (2.15)
Of course we have J(N ) I(N ) since P(N ) ⊂ K(N ). In many cases, we have that I(N ) = J(N ) and that any minimizer, when it exists, is automatically a pure HFB state. We give two results in the literature going in this direction. The first deals with purely repulsive interactions and it is Lieb's famous variational principle [START_REF] Lieb | Variational principle for many-fermion systems[END_REF] (see also Thm. 2.11 in [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF]).

Theorem 2.1 (Lieb's HF Variational Principle [START_REF] Lieb | Variational principle for many-fermion systems[END_REF]). Assume that W 0 and let N be an integer. Then for any (γ, α) ∈ K(N ), there exists (γ ′ , 0) ∈ P(N ) such that

E(γ, α) E(γ ′ , 0).
In particular, we have I(N ) = J(N ).

If W > 0 a.e., then any minimizer for I(N ), when it exists, is necessarily of the form (γ ′ , 0) with (γ ′ ) 2 = γ ′ .

We see that for repulsive interactions, W > 0, there is never pairing (α ≡ 0) and the ground state is always a pure HF state, that is, a Slater determinant. The fact that, in HF theory, one can minimize over mixed states and get the same ground state energy is very important from a numerical point of view. This was used by Cancès and Le Bris [START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF][START_REF] Cancès | Can we outperform the DIIS approach for electronic structure calculations?[END_REF] to derive well-behaved numerical strategies, to which we will come back later in Section 4.2.

For purely attractive interactions, the following recent result of Bach, Fröhlich and Jonsson [START_REF] Bach | Bogolubov-Hartree-Fock mean field theory for neutron stars and other systems with attractive interactions[END_REF] is relevant. Theorem 2.2 (HFB Constrained Variational Principle [START_REF] Bach | Bogolubov-Hartree-Fock mean field theory for neutron stars and other systems with attractive interactions[END_REF]). Assume that the number of spin states is q = 2 (spin-1/2 fermions), and that W can be decomposed in the form

W (x -y) = - Ω dµ(ω)1 Aω (x) 1 Aω (y) (2.16)
on a given measure space (Ω, µ), with A ω a measurable family of bounded domains in R 3 .

Let N be any positive real number. Then for any (γ, α) ∈ K(N ), we have

E(γ, α) E(γ ′ , α ′ ), (2.17)
with

γ ′ = g ⊗ 1 0 0 1 , α ′ = g(1 -g) ⊗ 0 1 -1 0 (2.18)
(the matrices act on the spin variables), and

g = g T = g = γ ↑↑ + γ ↓↓ + γ ↑↑ + γ ↓↓ 4 .
This HFB state is pure: (γ ′ , α ′ ) ∈ P(N ). In particular, we have I(N ) = J(N ). Furthermore, if W < 0 a.e., then any minimizer for I(N ), when it exists, is necessarily of the previous form.

Remark 2.1. Note that N does not have to be an even integer in this result. Since Tr L 2 (R 3 ) (g) = N/2, the operator g must have one eigenvalue different from 1 when N is odd, and it follows that α = 0 in this special case.

In Theorem 2.2, we have decomposed the operator γ acting on L 2 (R 3 × {↑, ↓}, C) according to the spin variables as follows:

γ = γ ↑↑ γ ↓↑ γ ↑↓ γ ↓↓ .
Theorem 2.2 says that when W satisfies (2.16), one can minimize over states which are pure, real, and have a simple spin symmetry. The antisymmetry of α is only contained in the spin variables, hence the Cooper pairs are automatically in a singlet state. Of course, one can express the total energy only in terms of the real operator g, as follows

E(γ ′ , α ′ ) = 2 Tr L 2 (R 3 ) (-∆)g + R 3 R 3 W (x -y) 2ρ g (x)ρ g (y) -|g(x, y)| 2 + | g(1 -g)(x, y)| 2 .
In practice it will be more convenient to keep a pairing term a(x, y) not a priori related to g and to optimize over both g and a, that is, to consider mixed states. When W satisfies the assumptions of the theorem, any ground state will automatically lead to a = ± g(1g).

Let us conclude our comments on Theorem 2.2, by noticing that several simple attractive potentials W can be written in the form (2.16). For instance the Fefferman-de la Llave formula [START_REF] Fefferman | Relativistic stability of matter. I[END_REF]

1 |x -y| = 1 π ∞ 0 dr r 5 R 3 dz 1 B(z,r) (x) 1 B(z,r) (y)
shows that a simple Newtonian interaction W (xy) = -|x -y| -1 is covered (here 1 B(z,r) is the characteristic function of the ball centered at z, of radius r). Hainzl and Seiringer showed in [START_REF] Hainzl | General decomposition of radial functions on R n and applications to N -body quantum systems[END_REF] that any smooth enough radial function W can be written in the form

W (x -y) = ∞ 0 dr g(r) R 3 dz 1 B(z,r) (x) 1 B(z,r) (y)
for some explicit function g, whose sign can easily be studied.

Existence results and properties of minimizers

Before turning to the discretization and the numerical study of the HFB minimization problem, we make some comments on the existence of minimizers. In addition to the infinite dimension and the nonlinearity of the model, an important difficulty is the invariance under translations. For instance, there are always minimizing sequences which do not converge strongly (assuming there exists a minimizer, one can simply translate it far away). Consider the electrons in an atom or in a molecule, with fixed classical nuclei (Born-Oppenheimer approximation). From the point of view of the electrons, the problem is no more translation-invariant, once the nuclei have been given a fixed position. Since the Coulomb interaction between the electrons is repulsive, W (xy) = |x -y| -1 , Theorem 2.1 tells us that there is never pairing, α ≡ 0. In this case there are several existence results, starting with the fundamental works of Lieb and Simon [START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF] and continuing with works by Lions [START_REF] Lions | Solutions of Hartree-Fock equations for Coulomb systems[END_REF], Bach [START_REF] Bach | Error bound for the Hartree-Fock energy of atoms and molecules[END_REF], Solovej [START_REF] Solovej | Proof of the ionization conjecture in a reduced Hartree-Fock model[END_REF][START_REF]The ionization conjecture in Hartree-Fock theory[END_REF].

For interactions W which have no particular sign, the pure HF problem was studied by Friesecke in [START_REF] Friesecke | The multiconfiguration equations for atoms and molecules: charge quantization and existence of solutions[END_REF] and by the first author of this paper in [START_REF] Lewin | Geometric methods for nonlinear many-body quantum systems[END_REF]. There is an HVZ-type theorem for HF wavefunctions which states that some binding conditions imply the existence of minimizers (Theorem 22 in [START_REF] Lewin | Geometric methods for nonlinear many-body quantum systems[END_REF]).

After the fundamental paper of Bach, Lieb and Solovej [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF] with its study of the Hubbard model, to our knowledge the existence of ground states for the HFB model with pairing was only studied recently by Lenzmann and the first author of this paper in [START_REF] Lenzmann | Minimizers for the Hartree-Fock-Bogoliubov theory of neutron stars and white dwarfs[END_REF]. Some caricatures of HFB in nuclear physics had been previously considered by Gogny and Lions in [START_REF] Gogny | Hartree-Fock theory in nuclear physics[END_REF].

We give here an existence result for the variational problem I(N ). In some cases we have I(N ) = J(N ) (see the previous section) and for this reason, we only concentrate on I(N ). The next theorem can be proved by following the method of [START_REF] Lenzmann | Minimizers for the Hartree-Fock-Bogoliubov theory of neutron stars and white dwarfs[END_REF], which dealt with the more complicated case of a pseudo-relativistic kinetic operator √ 1 -∆ -1.

Theorem 2.3 (Existence of minimizers and compactness of minimizing sequences). We assume as before that W = W 1 +W 2 ∈ L p (R 3 )+L q (R 3 ) with 2 p q < ∞.

Let λ > 0. Then the following assertions are equivalent:

(1) All the minimizing sequences (γ n , α n ) ⊂ K(λ) for I(λ) are precompact up to translations, that is there exists a sequence (x k ) ⊂ R 3 and (γ, α) ∈ K(λ) such that, for a subsequence,

lim k→∞ (1 -∆) 1/2 (τ x k γ n k τ -x k -γ)(1 -∆) 1/2 S1 = lim k→∞ (1 -∆) 1/2 (τ x k α n k τ -x k -α) S2 = 0.
In particular (γ, α) is a minimizer for I(λ).

(2) The binding inequalities

I(λ) < I(λ -µ) + I(µ) for all 0 < µ < λ (2.19) are satisfied.
Furthermore, if W is Newtonian at infinity, that is

W (x) - a |x| for a > 0 and |x| R, (2.20)
then the previous two equivalent conditions are verified.

The assumption that the interaction is Newtonian at infinity is a big simplification, as it means that two subsystems receiding from each other always attract at large distances. One can expect that minimizers exist even if the potential is not attractive at infinity, as soon as it has a sufficiently large negative component. A typical effective potential W (x) used in nuclear physics is nonnegative for small |x|, and has a negative well at intermediate distances [START_REF] Ring | The nuclear many-body problem[END_REF]. At infinity it typically decays like +κ|x| -1 for two protons, and exponentially fast when one of the two particles is a neutron. Even in HF theory, we are not aware of any existence result dealing with such potentials, however. We will numerically investigate a model of this form in Section 5.3.

The form of the nonlinear equation solved by minimizers is well-known in the physics literature, and it was re-explained in [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF]. The following result summarizes some known properties.

Theorem 2.4 (HFB equation and properties of minimizers [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF][START_REF] Lenzmann | Minimizers for the Hartree-Fock-Bogoliubov theory of neutron stars and white dwarfs[END_REF]). A HFB minimizer on K(N ) solves the nonlinear equation

Γ = 1 (-∞,0) F Γ -µN + δ (2.21)
where 0 δ 1 {0} (F Γ -µN ) has the same form as Γ, and where

N := 1 0 0 -1 , F Γ = h γ π π * -h γ (2.22) with h γ = -∆ + ρ γ * W -W (x -y)γ(x, y) and π(x, y) = α(x, y)W (x -y).
If W (xy) = -κ|x -y| -1 (Newtonian interaction) and N is an integer, then all the minimizers are of the special form (2.18). In this case, we have either α ≡ 0 and γ is a projector of rank N , or α = 0 and γ has an infinite rank.

The nonlinear equation (2.21) is in principle similar to the usual equation obtained in HF theory,

γ = 1 (-∞,µ) (h γ ) + δ (2.23)
Indeed, (2.21) reduces to (2.23) when α ≡ 0. Let us however emphasize that the mean-field operator F Γ has a spectrum which is symmetric with respect to 0. Hence F Γ is usually not even semi-bounded, on the contrary to h γ which is always bounded from below. Furthermore, the operator N does not commute with F Γ (except when α ≡ 0) and the equation cannot be written in a simple form as in HF theory. This will cause several difficulties to which we will come back at length later. The fact that γ has an infinite rank when there is pairing, α = 0, is a dramatic change of behavior compared to the simple HF case. However, no information on the decay of the eigenvalues of γ seems to be known.

An important open question is to show that minimizers actually exhibit non-vanishing pairing α = 0, at least for a sufficiently strong attractive potential W . On heuristic grounds, one expects such a phenomenon of "Cooper pair formation" to be energetically favorable due to the attractive interaction among particles. However, it seems to be a formidable task to find mathematical proof for this claim. The existence of pairing is known in some particular situations (when N is odd and W is Newtonian, see Remark 2.1, for the Hubbard model [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF], or in translation-invariant BCS theory [START_REF] Billard | An existence proof for the gap equation in the superconductivity theory[END_REF][START_REF] Vansevenant | The gap equation in superconductivity theory[END_REF][START_REF] Yang | On the Bardeen-Cooper-Schrieffer integral equation in the theory of superconductivity[END_REF][START_REF] Mcleod | The uniqueness and approximation of a positive solution of the Bardeen-Cooper-Schrieffer gap equation[END_REF][START_REF] Frank | The critical temperature for the BCS equation at weak coupling[END_REF][START_REF] Hainzl | The BCS functional for general pair interactions[END_REF][START_REF]The BCS critical temperature for potentials with negative scattering length[END_REF]), but for the model presented here, we are not aware of any result of this sort. One of the purposes of this paper is to investigate this question numerically.

Discretized Hartree-Fock-Bogoliubov theory

In this section, we write and study the Hartree-Fock-Bogoliubov model in a discrete basis.

Convergence analysis

Here we show that the ground state HFB energy in a finite basis converges to the true HFB ground state energy when the size of the basis grows. We consider a sequence of finitedimensional spaces

V h ⊂ H 1 (R 3 , C q ) for h → 0. We assume that any function f ∈ H 1 (R 3 , C q ) can be approximated by functions from V h : ∀f ∈ H 1 (R 3 , C q ), ∃f h ∈ V h such that ||f -f h || H 1 -→ h→0 0. (3.1)
We typically think of a sequence V h given by the Finite Elements Method. Let π h denote the orthogonal projection on V h in L 2 (R 3 , C q ). We define the set of density matrices living on V h (with average particle number N ) as follows

K h (N ) = (γ, α) ∈ K(N ) : π h γπ h = γ, π h α π h = α . (3.2)
The corresponding minimization problem is now

I h (N ) = inf (γ,α)∈K h (N ) E(γ, α). (3.3) Since K h (N ) ⊂ K(N ) by definition, it is obvious that I h (N ) I(N ).
The following result is a consequence of Theorem 2.3.

Theorem 3.1 (Convergence of the approximate HFB problem).

When W = W 1 + W 2 ∈ L p (R 3 ) + L q (R 3
) with 2 p q < ∞ and under Assumption (3.1) on the sequence (V h ), we have

lim h→0 I h (N ) = I(N ). (3.4)
If the binding inequality (2.19) is satisfied, then any sequence of minimizers (γ h , α h ) ∈ K h (N ) for I h (N ) converges, up to a subsequence and up to a translation, to a minimizer

(γ, α) ∈ K(N ) of I(N ), in the sense that lim h k →0 (1 -∆) 1/2 (τ x k γ h k τ -x k -γ)(1 -∆) 1/2 S1 = lim h k →0 (1 -∆) 1/2 (τ x k α h k τ -x k -α) S2 = 0. (3.5)
Proof. We only have to show that I h (N ) → I(N ) as h → 0. Then, any sequence of exact minimizers (γ h , α h ) for I h (N ) is also a minimizing sequence for I(N ). Applying Theorem 2.3 concludes the proof. We know that finite-rank operators are dense in K(N ). Let (γ, α) ∈ K(N ) be any such finite rank operator. Let (f i ) K i=1 be an orthonormal basis of the range of γ. By Löwdin's theorem (Lemma 13 in [START_REF] Lewin | Geometric methods for nonlinear many-body quantum systems[END_REF]), we know that the two-body wavefunction α can be expanded in the same basis (f 1 , ..., f K ). Now, for every i = 1, .., K, we apply (3.1) and take a sequence

f h i ∈ V h be such that f h i → f i in H 1 (R 3 ). The system (f h i ) k i=1 is not necessarily orthonormal but we have f h i , f h j → f i , f j = δ ij .
Applying the Gram-Schmidt procedure, we can therefore replace the (f h i ) K i=1 by an orthonormal set (g h i ) K i=1 ⊂ V h having the same properties. An equivalent procedure is to take

g h i = K j=1 (S -1/2 h
) ji f h j where S h is the Gram matrix ( f h i , f h j ). Let now U h be any unitary operator on L 2 (R 3 ) which is such that U h f i = g h i for all i = 1, ..., K. We then take γ h := U h γU * h and α h := U h αU T h . In words, we just replace the f i by g h i in the decomposition of γ and α. To see that (γ h , α h ) ∈ K(N ), we just notice that

γ h α h α * h 1 -γ h = U h 0 0 U h γ α α * 1 -γ U * h 0 0 U h * .
Note also that Tr(γ h ) = Tr(γ) = N since U h is unitary. Now γ h and α h belong to K h (N ) by definition, hence we have that E(γ h , α h ) I h (N ). On the other hand, by the convergence of

f h i (hence of g h i ) towards f i in H 1 (R 3 ), we easily see that lim h→0 E(γ h , α h ) = E(γ, α)
by continuity of E. Therefore we have proved that lim sup

h→0 I h (N ) E(γ, α).
This is valid for all finite-rank (γ, α) ∈ K(N ), hence we deduce that lim sup

h→0 I h (N ) inf (γ,α)∈K(N ) E(γ, α) = I(N ).
On the other hand the inequality I h (N ) I(N ) is always satisfied, hence we have proved the claimed convergence I h (N ) → I(N ).

Discretization

In this section we compute the HFB energy E(γ, α) of a discretized state (γ, α) ∈ K h (N ) and we write the corresponding self-consistent equation. We fix once and for all the approximation space V h and we consider a basis set (χ i ) N b i=1 of V h , which is not necessarily orthonormal. We will assume that V h is stable under complex conjugation, which means that f ∈ V h ⇒ f ∈ V h (this amounts to replacing V h by Span(V h , V h )). Then we can choose a basis (χ i ) N b i=1 which is real, that is χ i = χ i for all i = 1, ..., N b . This will dramatically simplify the calculation below.

Since π h γπ h = γ and π h απ h = α, we can write the kernels of γ and α as follows

γ(x, y) σ,σ ′ = N b i,j=1 G ij χ i (x) σ χ j (y) σ ′ , α(x, y) σ,σ ′ = N b i,j=1 A ij χ i (x) σ χ j (y) σ ′ . (3.6)
The complex conjugation on χ j is superfluous but we keep it to emphasize the difference between γ and α. The matrices G and A (defined by the previous relation) satisfy the constraints G * = G and A T = -A. Note that since A is antisymmetric, we can also write

α(x, y) σ,σ ′ = 1 i<j N b A ij χ i (x) σ χ j (y) σ ′ -χ j (x) σ χ i (y) σ ′ = √ 2 1 i<j N b A ij χ i ∧ χ j (x, σ, y, σ ′ )
where

(χ i ∧ χ j )(x, σ; y, σ ′ ) := χ i (x) σ χ j (y) σ ′ -χ i (y) σ ′ χ j (x) σ / √
2 is a two-body Slater determinant. Let us also remark that (3.6) can be written in the operator form

γ = N b i,j=1 G ij |χ i χ j |, α = N b i,j=1 A ij |χ i χ j |.
Again the complex conjugation on χ j is superfluous.

Let us define the so-called overlap matrix S = S * = S by

S ij = χ i , χ j H = q σ=1 R 3 χ i (x) σ χ j (x) σ dx = q σ=1 R 3 χ i (x) σ χ j (x) σ dx. (3.7)
It is tedious but straightforward to verify that the constraint Another way to write the constraint is 0 S 1/2 ΥS 1/2 1. Let us notice that we have used everywhere the fact that S = S * = S = S T . The formulas are much more complicated when S is not real. The energy can be expressed in terms of the matrices G and A, as well. A calculation shows that

0 0 0 0 Γ := γ α α * 1 -γ 1
E(γ, α) = Tr(hG) + 1 2 Tr(G J(G)) - 1 2 Tr(G K(G)) + 1 2 Tr(A * X(A)). (3.11)
The trace here is the usual one for N b × N b matrices. As we think that there is no possible confusion with E(γ, α), we will also denote by E(G, A) this discretized energy functional. In formula (3.11),

h ij = χ i , (-∆)χ j = q σ=1 R 3 ∇χ i (x) σ • ∇χ j (x) σ dx, J(G) ij = N b k,ℓ=1 (ij|ℓk) G kℓ , K(G) ij = N b k,ℓ=1 (ik|ℓj) G kℓ , X(A) ij = N b k,ℓ=1 (ik|jℓ) A kℓ , (3.12) and
(ij|kℓ) :=

R 3 R 3 W (x -y) χ i (x) * χ j (x) χ k (y) * χ ℓ (y) dx dy. (3.13)
We use here the notation a * b = q σ=1 a σ b σ (but the complex conjugation is superfluous for our real basis, hence K = X). Similarly, we have Tr(γ) = Tr(SG).

(3.14)

We define the discretized number operator as

N = S 0 0 -S (3.15)
The constraint Tr(γ) = N can be written equivalently as

Tr(NΥ) = 2N -N b
We deduce from this calculation that the variational problem I h (N ) can be written in finite dimension as

I h (N ) = min E(G, A) : 0 ΥSΥ Υ, Tr(NΥ) = 2N -N b (3.16)
where we recall that Υ and S have been defined in (3.10). Here the infimum is always attained because the problem is finite dimensional. In this form, the discretized problem is very similar to the usual discretized Hartree-Fock problem [START_REF] Cancès | Computational quantum chemistry: a primer[END_REF][START_REF] Bris | Computational chemistry from the perspective of numerical analysis[END_REF], in dimension 2N b instead of N b . There is a big difference, however. In HF theory the constraint involves the matrix S instead of N. This difference will cause several difficulties. To understand the problem, let us introduce a new variable Υ ′ = S 1/2 ΥS 1/2 (which is the same as orthonormalize the basis (χ i )). Then the constraint 0 ΥSΥ Υ is transformed into 0 Υ ′ 1. However, the constraint on the number of particles becomes Tr(S -1/2 NS -1/2 Υ ′ ) = 2N -N b . In usual Hartree-Fock theory, the matrix S -1/2 NS -1/2 is replaced by the identity. The fact that this matrix then commutes with the Fock Hamiltonian (defined below) simplifies dramatically the self-consistent equations. Here we get

S -1/2 NS -1/2 = 1 0 0 -1
which commutes with the Fock Hamiltonian if and only if A ≡ 0. The self-consistent equation is obtained like in [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF]. The result is as follows.

Lemma 3.1 (Discretized HFB equation). Let Υ be a minimizer for the variational problem I h (N ). Then there exists µ ∈ R such that Υ solves the linear problem

min Tr(F Υ -µN) Υ : 0 ΥS Υ Υ (3.17)
where

F Υ := h G X(A) X(A) * -h G , h G := h + J(G) -K(G). (3.18)
The solution Υ can be written in the form

Υ = S -1/2 1 (-∞,0) S -1/2 (F Υ -µN)S -1/2 + δ S -1/2 (3.19)
where 0 δ 1 lives only in the kernel of

S -1/2 (F Υ -µN)S -1/2 .
The solution Υ of the self-consistent equation (3.19) may be equivalently written by considering the generalized eigenvalue problem

F Υ -µN f i = ǫ i S f i , f i , Sf j = δ ij . (3.20)
Then we have simply (assuming

ǫ i = 0 for all i = 1, ..., 2N b ) Υ = ǫi<0 f i f * i .
Again, this is similar to the Hartree-Fock solution [START_REF] Cancès | Computational quantum chemistry: a primer[END_REF][START_REF] Bris | Computational chemistry from the perspective of numerical analysis[END_REF] except that µ is unknown and N does not always commute with F Υ . Remark that although the basis functions χ i are real, the density matrix Υ is not necessarily real. In the next section, we will restrict ourselves to real-valued density matrices and impose some spin symmetry.

Using symmetries

The HFB energy E(Γ) has some natural symmetry invariances which we describe in detail in this section. Recall that since E is a nonlinear functional, it cannot be guaranteed that the HFB minimizers will all have the same symmetries as the HFB energy. The set of all minimizers will be invariant under the action of the symmetry group, but each minimizer alone does not have to be invariant.

We have already allowed for the breaking of particle-number symmetry and we hope to find an HFB ground state. It will then automatically break the translational invariance of the system. There are three other symmetries (spin, complex conjugation and rotations) which are of interest to us. We have the choice of imposing these symmetries by adding appropriate constraints, or not. Because this reduces the computational cost, it will be convenient to impose them.

Time-reversal symmetry

Let us now assume that q = 2, which means that our fermions are spin-1/2 particles. Since the Laplacian and the interaction function W do not act on the spin variable, the HFB energy has some spin symmetry, which can be written for q = 2 as

∀k = 1, 2, 3, E(Σ k ΓΣ * k ) = E(Γ) where Σ k := iσ k 0 0 iσ k ,
with σ 1 , σ 2 , σ 3 being the usual Pauli matrices. Note that Σ k has the form of a Bogoliubov transformation, hence Σ k ΓΣ * k is also an HFB state. The number operator is also invariant which means that

Σ k N Σ * k = N for all k = 1, 2, 3. Thus, the contraint Tr(γ) = N is conserved and we have Σ k ΓΣ * k ∈ K(N ) when Γ ∈ K(N ).
Another important symmetry is that of complex conjugation which means this time that

E(Γ) = E(Γ)
and which is based on the fact that the Laplacian and W are real operators. Again we have Tr(γ) = Tr(γ) hence K(N ) is invariant under complex conjugation. As was explained in [START_REF] Hainzl | On blowup for time-dependent generalized Hartree-Fock equations[END_REF] (see Remark 5 page 1032), the density matrices (γ, α) can be written in the special form In other words, Γ is invariant under the action of the group generated by Σ 1 , Σ 2 and C. This invariance is sometimes called the time-reversal symmetry. As remarked in [START_REF] Hainzl | On blowup for time-dependent generalized Hartree-Fock equations[END_REF], imposing Σ k ΓΣ * k = Γ for all k = 1, 2, 3 implies α ≡ 0 which is not interesting for us. When W can be written in the form (2.16), the Theorem 2.2 of Bach, Fröhlich and Jonsson [START_REF] Bach | Bogolubov-Hartree-Fock mean field theory for neutron stars and other systems with attractive interactions[END_REF], tells us that there is no breaking of the time-reversal symmetry. That is, we can always minimize over such special states. For other interactions W this is not necessarily true but it is often convenient to impose this symmetry anyhow.

γ ′ = g ⊗ 1 0 0 1 , α ′ = a ⊗ 0 1 -1 0 , with g = g T = g
Because it holds

F Σ k ΓΣ * k = Σ k F Γ Σ * k , F Γ = F Γ , it can
then be verified that minimizers under the additional symmetry constraint, satisfy the same self-consistent equation as when no constraint is imposed.

When we discretize the problem by imposing time-reversal symmetry, we use two real symmetric matrices G and A, related through the constraint that

0 ΥSΥ Υ, with Υ := G A A S -1 -G and S = S 0 0 S . (3.23)
The energy becomes

E(G, A) = 2 Tr(hG) + 2 Tr(G J(G)) -Tr(G K(G)) + Tr(A K(A)). (3.24)
and the associated particle number constraint is Tr(SG) = N/2. In practice we always assume that N is even for simplicity. The basis (χ i ) is now composed of (real-valued) functions in H 1 (R 3 , R), instead of functions in H 1 (R 3 , C q ) as before, and the formulas for S, h, J, K and X are the same as before.

Rotational symmetry

The group SO(3) of rotations in R 3 also acts on HFB states and it leaves the energy invariant when the interaction W is a radial function. In this section we always assume that the spin variable has already been removed according to the previous section and we denote by g and a the corresponding (real symmetric) density matrices. If the spin were still present, rotations would act on it as well.

To any rotation R ∈ SO(3) we can associate a unitary operator on L 2 (R 3 ), denoted also by R, defined by (Rϕ)(x) = ϕ R -1 x . Now we say that an HFB state Γ with density matrices (γ, α) is invariant under rotations when it satisfies

R ΓR * = Γ, where R = R 0 0 R . Note that R is a Bogoliubov rotation since R = R.
The density matrices of an invariant state satisfy g(Rx, Ry) = g(x, y), a(Rx, Ry) = a(x, y) for all x, y ∈ R 3 and any rotation R ∈ SO(3).

As the angular momentum L = x × (-i∇) generates the group of rotations, a (smooth enough) HFB state is invariant under rotations if and only if

L Γ = ΓL , where L = L 0 0 L .
The density matrices g and a can then be written in the special form

g(x, y) = 1 4π ℓ 0 g ℓ (|x|, |y|) (2ℓ+1) P ℓ ω x •ω y , a(x, y) = 1 4π ℓ 0 a ℓ (|x|, |y|) (2ℓ+1) P ℓ ω x •ω y (3.25
) where P ℓ is the Legendre polynomial of degree ℓ, which is such that P ℓ (1) = 1. The constraint on g and a are transfered in each angular momentum sector (labelled by ℓ 0), leading to 0 0 0 0

g ℓ a ℓ a ℓ 1 -g ℓ 1 0 0 1 on L 2 ([0, ∞), r 2 dr)⊕L 2 ([0, ∞), r 2
dr). However, there is no such constraint between different ℓ's. The total average particle number is given by

Tr(g) = ℓ 0 (2ℓ + 1) Tr(g ℓ ) = N/2.
This is the only constraint which mixes the different angular momentum density matrices. Now we can discretize the radial HFB problem. We choose a finite-dimensional subspace V rad in L 2 ([0, ∞), r 2 dr) with basis (χ 1 , ..., χ N b ), which we use to expand the density matrices g ℓ and a ℓ . Then we fix a maximal angular momentum ℓ max and we truncate the series in (3.25). This is the same as taking as discretization space

V = f (|x|) Y m ℓ (ω x ) : f ∈ V rad , 0 ℓ ℓ max , -ℓ m ℓ ⊂ L 2 (R 3 , R)
where Y m ℓ is the spherical harmonics of total angular momentum ℓ and azimuthal angular momentum m. We then assume that g and a are radial and live in this space. The matrices G ℓ and A ℓ of g ℓ and a ℓ in the basis (χ i ) are defined similarly as before by

g ℓ (r, r ′ ) = N b i,j=1 G ℓ ij χ i (r)χ j (r ′ ), a ℓ (r, r ′ ) = N b i,j=1 A ℓ ij χ i (r)χ j (r ′ ). (3.26)
The constraints on the matrices G ℓ and A ℓ are

0 Υ ℓ SΥ ℓ Υ ℓ , with Υ ℓ := G ℓ A ℓ A ℓ S -1 -G ℓ and S = S 0 0 S (3.27)
with The total HFB energy is now

S ij = ∞ 0 χ i (r) χ j (r)
E(G 0 , ..., G ℓmax , A 0 , ..., A ℓmax ) = 2 ℓmax ℓ=0 (2ℓ + 1) Tr(h ℓ G ℓ ) + ℓmax ℓ,ℓ ′ =0 (2ℓ + 1)(2ℓ ′ + 1) 2 Tr(G ℓ J(G ℓ ′ )) -Tr(G ℓ K ℓℓ ′ (G ℓ ′ )) + Tr(A ℓ K ℓℓ ′ (A ℓ ′ )) , (3.29)
where

h ℓ ij := ∞ 0 χ ′ i (r) χ ′ j (r) r 2 dr + ℓ(ℓ + 1) ∞ 0 χ i (r) χ j (r) dr, J(G ℓ ′ ) ij := N b m,n=0 (ij|nm) 0,0 G ℓ ′ mn , K ℓℓ ′ (G ℓ ′ ) ij := N b m,n=0 (im|jn) ℓ,ℓ ′ G ℓ ′ mn , (ij|mn) ℓ,ℓ ′ := ∞ 0 r 2 dr ∞ 0 s 2 ds χ i (r) χ j (r) χ m (s) χ n (s) w ℓ,ℓ ′ (r, s)
and

w ℓ,ℓ ′ (r, s) := 1 2 1 -1 W r 2 + s 2 -2rst P ℓ (t) P ℓ ′ (t) dt.
Any minimizer (G 0 , ..., G ℓmax , A 0 , ..., A ℓmax ) of E under the constraints (3.27) and (3.28) will be of the form

Υ ℓ = ǫ ℓ i <0 f ℓ i f ℓ i T , 0 ℓ ℓ max , (3.30)
where the vectors f ℓ i 's solve the generalized eigenvalue problem

F ℓ -µ(2ℓ + 1)N f ℓ i = ǫ ℓ i S f ℓ i , f i , Sf j = δ ij (3.31)
with

F ℓ := h ℓ 0 0 -h ℓ + ℓmax ℓ ′ =0 2J G ℓ ′ -K ℓℓ ′ G ℓ ′ K ℓℓ ′ A ℓ ′ K ℓℓ ′ A ℓ ′ -2J G ℓ ′ + K ℓℓ ′ G ℓ ′ .
The Euler-Lagrange multipler µ appearing in (3.31) is common to all the different angular momentum sectors and it is chosen to ensure the validity of the constraint (3.28).

Algorithmic strategies and convergence analysis

In this section we study the convergence of two algorithms which can be used in practice to solve the HFB minimization problem (2.15). In order to simplify our presentation, we restrict ourselves to the finite-dimensional case, that is, to the discretized problem (3.16).

We also assume that the discretization basis (χ j ) is orthonormal, such that S = I 2N b , the (2N b ) × (2N b ) identity matrix. Finally, we only consider states which are invariant under time-reversal symmetry like in Section 3.3.1. This means that the HFB state is described by real and symmetric matrices G and A such that

0 0 0 0 Υ := G A A 1 -G 1 0 0 1 (4.1)
The energy is given by (3.24),

E(Υ) = 2 Tr(hG) + 2 Tr(G J(G)) -Tr(G K(G)) + Tr(A K(A)). (4.2)
The extension to more general situations is straightforward.

The energy E is continuous (it is indeed real-analytic) with respect to Υ. Also the set K of density matrices Υ of the form (4.1) is compact in finite dimension. Hence minimizers always exist and, as we have seen, they solve the nonlinear equation

Υ = 1 (-∞,0) F Υ -µN + δ, (4.3)
where µ is a Lagrange multiplier chosen to ensure the constraint that Tr(G) = N/2. Of course we must have N/2 N b , the dimension of the (no-spin) discretization space V h , otherwise the minimization set is always empty.

Roothaan Algorithm

The most natural technique used in practice to solve the equation (4.3) is a simple fixed point method [START_REF] Quentin | Self-Consistent Calculations of Nuclear Properties with Phenomenological Effective Forces[END_REF][START_REF] Dechargé | Hartree-Fock-Bogolyubov calculations with the D1 effective interaction on spherical nuclei[END_REF]. This is usually refered to as the Roothaan algorithm in the chemistry literature [START_REF] Roothaan | New developments in molecular orbital theory[END_REF]. The iteration scheme is the following

Υ n+1 = 1 (-∞,0) F Υn -µ n+1 N + δ n+1 . (4.4)
At each step, one has to determine µ n+1 such as to satisfy the constraint Tr(G n+1 ) = N/2.

If the operator F Υnµ n+1 N has a trivial kernel, then δ n+1 ≡ 0. This is the usual situation encountered in practice. In the iteration (4.4), the state is assumed to be pure at each step (Υ is an orthogonal projection). Recall that by Theorem 2.2 we know that minimizers of E under a particle number constraint are always pure, under suitable assumptions on the interaction potential W . The algorithm is stopped when the commutator Υ n , F Υn and/or when the variation of the HFB state

||Υ n+1 -Υ n ||
are smaller than a prescribed ε.

Our purpose in this section is to study the behavior of the Roothaan algorithm (4.4). In the Hartree-Fock case, it was shown in a fundamental work of Cancès and Le Bris [START_REF] Cancès | Can we outperform the DIIS approach for electronic structure calculations?[END_REF][START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF], that the algorithm converges or oscillate between two points, none of them being a solution to the equation (4.3). This result was recently improved by Levitt in [START_REF] Levitt | Convergence of gradient-based algorithms for the Hartree-Fock equations[END_REF]. We will explain that the results of Cancès-Le Bris and Levitt can be generalized to the HFB model. Actually, in a discretization space of dimension N b , HFB is equivalent to a Hartree-Fock-like minimization problem in dimension 2N b , with additional constraints. The adaptation of the previously cited works in the HF case reduces to handling these contraints.

In order to avoid the convergence problems of the Roothaan algorithm, Cancès and Le Bris have proposed the Optimal Damping Algorithm (ODA). We will study the equivalent of this algorithm in HFB theory in Section 4.2.

To start with, we show that the Roothaan algorithm is well defined, in the sense that for any HFB state Υ n , there exists (Υ n+1 , µ n+1 , δ n+1 ) solving (4.4). To this end, we follow [START_REF] Cancès | Can we outperform the DIIS approach for electronic structure calculations?[END_REF][START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF] and introduce the auxillary functional Ẽ(Υ, Υ ′ ) := Tr(hG) + Tr(hG

′ ) + 2 Tr(G J(G ′ )) -Tr(G K(G ′ )) + Tr(A K(A ′ )) (4.5)
as well as the variational problem

I Υ (λ) := min Υ ′ Ẽ(Υ, Υ ′ ) : Tr(G ′ ) = λ (4.6)
which consists in minimizing over Υ ′ with Υ fixed. The matrix Υ ′ must be an admissible HFB state which, in our context, means that

0 0 0 0 Υ ′ := G ′ A ′ A ′ 1 -G ′ 1 0 0 1 , (G ′ ) T = G ′ = G ′ , (A ′ ) T = A ′ = A ′ . (4.7)
Recall that we have chosen an orthonormal basis and that the spin has been eliminated. It is clear that (4.6) admits at least one solution Υ ′ , as soon as 0 λ N b , where we recall that N b is the dimension of the discretization space V h . The following says that these solutions are exactly those solving the equation of the Roothaan method.

Υ ′ = 1 (-∞,0) F Υ -µ ′ N + δ ′ (4.8) where µ ′ ∈ [∂ -I Υ (λ), ∂ + I Υ (λ)] and 0 δ ′ 1 {0} (F Υ -µ ′ N). If 0 / ∈ σ(F Υ -µ ′ N), then δ ′ ≡ 0 and Υ ′ is unique, for any such µ ′ ∈ [∂ -I Υ (λ), ∂ + I Υ (λ)].
Proof. To see that I Υ (λ) is convex, let 0 λ 1 < λ 2 N b and let Υ ′ i be a minimizer for

I Υ (λ i ) with i = 1, 2. Then tΥ ′ 1 + (1 -t)Υ ′
2 is a test state for the problem I Υ (tλ 1 + (1t)λ 2 ). Therefore it holds

I Υ (tλ 1 + (1 -t)λ 2 ) Ẽ(Υ, tΥ ′ 1 + (1 -t)Υ ′ 2 ) = t Ẽ(Υ, Υ ′ 1 ) + (1 -t) Ẽ(Υ, Υ ′ 2 ) = tI Υ (λ 1 ) + (1 -t)I Υ (λ 2 ).
Then, by convexity we get that I Υ (λ ′ )

I Υ (λ) + µ(λ ′ -λ) for any λ ′ ∈ [0, N b ] and any µ ∈ [∂ -I Υ (λ), ∂ + I Υ (λ)]. Thus I Υ (λ) -µλ = min{I Υ (λ ′ ) -µλ ′ : 0 λ ′ N b } = min Υ ′ Ẽ(Υ, Υ ′ ) -µ Tr(G ′ ) = Tr(hG) + 1 2 min Υ ′ Tr(F Υ -µN)Υ ′ .
In the previous two mins, Υ ′ is varied over all possible HFB states, without any particle number constraint. It is well known that the minimizers of the problem on the right side are exactly the solutions of the equation (4.8).

Lemma 4.1 tells us that for any given Υ n , there always exists at least one solution (Υ n+1 , µ n+1 , δ n+1 ) of the equation (4.4). It is obtained by solving the minimization problem I Υn (N/2), and one has to take µ n+1 ∈ [∂ -I Υn (N/2), ∂ + I Υn (N/2)]. We can always take by convention

µ n+1 := ∂ -I Υn (N/2) + ∂ + I Υn (N/2) 2 .
However, Υ n+1 is not uniquely defined yet because of the possibility of having δ n+1 = 0. As we have seen it is unique when

0 / ∈ σ F Υn -µ n+1 N . (4.9)
In this section we always assume that it is possible to find (Υ n+1 , µ n+1 , δ n+1 ) exactly. In practice, we will only know (Υ n+1 , µ n+1 , δ n+1 ) approximately. Later in Section 4.3 we explain how to do this numerically. We will also see that the condition (4.9) is "very often" satisfied. This vague statement is made precise in Lemma 4.4 below. Following Cancès and Le Bris, we now introduce the concept of uniform well-posedness.

Definition 4.1 (Uniform well-posedness). We say that for a given initial HFB state Υ 0 , the sequence (Υ n ) generated by the Roothaan algorithm is uniformly well posed when there exists η > 0 such that

|F Υn -µ n+1 N| η (4.10)
for all n 0, where

µ n+1 = ∂ -I Υn (N/2) + ∂ + I Υn (N/2) /2.
Note that the condition |F Υn -µ n+1 N| η is equivalent to (-η, η)∩σ(F Υn -µ n+1 N) = ∅. Later in Section 4.3 we will make several comments concerning Assumption (4.10).

We have seen that the sequence generated by the Roothaan algorithm can be obtained by solving the minimization problem

I Υn = min Υ ′ Ẽ(Υ n , Υ ′ ).
Since Ẽ(Υ, Υ ′ ) = Ẽ(Υ ′ , Υ), we conclude that the Roothaan algorithm is the same as minimizing Ẽ with respect to its first and second variables one after another, inductively. This fact allows to prove the following result, which is the HFB equivalent of Theorem 7 in [START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF] and Theorem 5.1 in [START_REF] Levitt | Convergence of gradient-based algorithms for the Hartree-Fock equations[END_REF] in the HF case. • The sequence (Υ 2n , Υ 2n+1 ) converges towards a critical point (Υ, Υ ′ ) of Ẽ;

• If Υ = Υ ′ , then this state is a solution of the original HFB equation (4.3), but if Υ = Υ ′ , then none of these two states is a solution to (4.3).

Theorem 4.1 says that (provided it is uniformly well posed) the sequence Υ n will either converge to a solution of the self-consistent Equation (4.3), or oscillate between two points Υ and Υ ′ , none of them being a solution to the desired equation.

Proof. We split the proof into several steps.

Step 1: µ n is uniformly bounded. It will be very useful to know that the sequence µ n is uniformly bounded. The following says that, as soon as we fix Tr(G) = N/2 with 0 < N/2 < N b , the chemical potential µ cannot be too negative and too positive. Lemma 4.2 (Bounds on the multiplier µ). Let Υ ′ be any fixed HFB state and

Υ µ := 1 (-∞,0) F Υ ′ -µN = G µ A µ A µ 1 -G µ
the corresponding HFB ground state at chemical potential µ. There exists a constant C which is independent of Υ ′ and µ, such that

∀µ -C, Tr G µ C |µ| and ∀µ C, Tr G µ N b - C µ . (4.11)
The lemma says that the average number of particles in 1 (-∞,0) F Υ -µN tends to N b when µ → ∞ whereas it tends to 0 when µ → -∞, this uniformly with respect to the state Υ ′ used to build the mean-field operator F Υ ′ .

Proof. We first remark that there exists a constant C such that

||F Υ ′ || C (4.12)
for any HFB state Υ ′ . This follows from the fact that F Υ ′ is continuous with respect to Υ ′ and that the latter lives in a compact set since we always have 0 Υ 1. The chosen norm for ||F Υ ′ || does not matter since we are in finite dimension. Now, for µ large enough we can use regular perturbation theory and obtain that

1 (-∞,0) F Υ ′ -µN -1 (-∞,0) -µN = 1 (-∞,0) F Υ ′ |µ| - µ |µ| N -1 (-∞,0) - µ |µ| N C |µ| .
Note that

1 (-∞,0) - µ |µ| N =              1 0 0 0 for µ > 0, 0 0 0 1 for µ < 0.
Taking the trace against N gives the result.

From Lemma 4.2 we deduce that our sequence µ n is bounded. Indeed, since we have by construction Tr(G n+1 ) = N/2 with 0 < N/2 < N b , we must have

-max C , 2C N µ n+1 max C , C N b -N/2 (4.13)
otherwise Tr(G n+1 ) would be too small or too large.

Step 2: convergence of Ẽ(Υ n , Υ n+1 ). We now follow [START_REF] Cancès | Can we outperform the DIIS approach for electronic structure calculations?[END_REF][START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF][START_REF] Cancès | SCF algorithms for HF electronic calculations[END_REF]. At each step, we know from Lemma 4.1 that Υ n+1 is a solution of the minimization problem min Υ ′ Ẽ(Υ n , Υ ′ ). In particular, we deduce that

Ẽ(Υ n , Υ n+1 ) Ẽ(Υ n , Υ n-1 ) = Ẽ(Υ n-1 , Υ n ). (4.14)
Thus the sequence of real numbers Ẽ(Υ n , Υ n+1 ) is non-increasing. It is also bounded from below, hence it converges to a limit ℓ. We now use the uniform well-posedness to prove an inequality which is more precise than (4.14). We remark that

Ẽ(Υ n , Υ n-1 ) -Ẽ(Υ n , Υ n+1 ) = 1 2 Tr F Υn Υ n-1 -Υ n+1 = 1 2 Tr F Υn -µ n+1 N Υ n-1 -Υ n+1 1 2 Tr F Υn -µ n+1 N Υ n-1 -Υ n+1 2 η 2 Tr Υ n-1 -Υ n+1 2 = η 2 ||Υ n-1 -Υ n+1 || 2 . (4.15)
In the above calculation we have used that Tr N Υ n+1 = Tr N Υ n-1 = N -N b . We have also used that Υ n+1 is the negative spectral projector of F Υnµ n+1 N, such that we can write

F Υn -µ n+1 N = F Υn -µ n+1 N Υ ⊥ n+1 -Υ n+1 .
Finally, we have used that 0 γ 1 is equivalent to (γ -P ) 2 P ⊥ (γ -P )P ⊥ -P (γ -P )P , for any orthogonal projector P , thus

Υ ⊥ n+1 Υ n-1 -Υ n+1 Υ ⊥ n+1 -Υ n+1 Υ n-1 -Υ n+1 Υ n+1 Υ n-1 -Υ n+1 2 .
Since Ẽ(Υ n , Υ n+1 ) converges to a limit ℓ, we deduce that

n 1 ||Υ n+1 -Υ n-1 || 2 < ∞.
In particular, ||Υ n+1 -Υ n-1 || → 0 which is called numerical convergence in [START_REF] Cancès | Can we outperform the DIIS approach for electronic structure calculations?[END_REF][START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF].

Step 3: convergence of Υ 2n and Υ 2n+1 . In order to upgrade the numerical convergence to true convergence, we use a recent method of Levitt [START_REF] Levitt | Convergence of gradient-based algorithms for the Hartree-Fock equations[END_REF]. Namely we show that

Ẽ(Υ n , Υ n-1 ) -ℓ θ -Ẽ(Υ n , Υ n+1 ) -ℓ θ η ′ 2 ||Υ n-1 -Υ n+1 || (4.16)
for a well chosen 0 < θ 1/2. Summing over n and using the convergence of Ẽ(Υ n , Υ n-1 ), hence of ( Ẽ(Υ n , Υ n-1 )ℓ) θ , then gives the convergence of Υ 2n and Υ 2n+1 .

For the proof of (4.16), we argue as follow. Consider a (real, no-spin) pure HFB state Υ. It is possible to parametrize the manifold of pure HFB states around Υ by using Bogoliubov transformations as follows:

H → e H Υe -H
where H is assumed to be of the form

h p -p -h , h T = -h = -h, p T = p = p.
These constraints ensure that iH is a self-adjoint Hamiltonian such that e H = e -i(iH) is a Bogoliubov rotation. They also ensure that e H Υe -H stays real. That H → e H Υe -H is a local chart of the manifold of pure HFB states around Υ follows from the arguments in [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF] as well as simple considerations in linear algebra.

Let us now consider the energy Ẽ in a neighborhood of any fixed (Υ, Υ ′ ). The map

f : (H, H ′ ) → Ẽ e H Υe -H , e H ′ Υ ′ e -H ′ - µ ′ 2 Tr Ne H Υe -H - µ 2 Tr Ne H ′ Υ ′ e -H ′
is real analytic in a neighborhood of (0, 0) for any fixed µ, µ ′ ∈ R and any fixed pure HFB states (Υ, Υ ′ ). The Lojasiewicz inequality (Theorem 2.1 in [START_REF] Levitt | Convergence of gradient-based algorithms for the Hartree-Fock equations[END_REF]) then tells us that there exist 0 < θ 1/2 and a constant κ > 0 such that

||H|| + ||H ′ || κ implies |f (H, H ′ ) -f (0)| 1-θ κ -1 |∇ H f (H, H ′ )| + |∇ H ′ f (H, H ′ )| .
A simple computation shows that

∇ H f (H, H ′ ) = 1 2 F Υ ′ -µ ′ N , e H Υe -H , ∇ H ′ f (H, H ′ ) = 1 2 F Υ -µN , e H ′ Υ ′ e -H ′ .
If we rephrase all this in our setting, this means that for any fixed pure HFB states (Υ 1 , Υ ′ 1 ) and any µ, µ ′ ∈ R, there is a κ > 0 such that for any (Υ 2 , Υ ′

2 ) another pure HFB state which is at most at a distance κ from (Υ 1 , Υ ′ 1 ), we have

Ẽ Υ 1 , Υ ′ 1 -Ẽ Υ 2 , Υ ′ 2 + µ ′ Tr N(G 2 -G 1 ) + µ Tr N(G ′ 2 -G ′ 1 ) 1-θ κ -1 F Υ ′ 2 -µ ′ N , Υ 2 + F Υ2 -µN , Υ ′ 2 . (4.17)
The constants κ and θ depend on µ, µ ′ and of the reference point (Υ 1 , Υ ′ 1 ). But they stay positive as soon as µ, µ ′ and (Υ 1 , Υ ′ 1 ) stay in a compact set. By a simple compactness argument, we therefore deduce that there exists a neighborhood of the compact set

(Υ, Υ ′ ) : Ẽ(Υ, Υ ′ ) = ℓ, Tr(G) = Tr(G ′ ) = N/2
such that for any (Υ, Υ ′ ) in this neighborhood and µ, µ ′ in a compact set in R, we have

Ẽ Υ, Υ ′ -ℓ + µ ′ N/2 -Tr NG + µ N/2 -Tr NG ′ 1-θ κ -1 F Υ ′ -µ ′ N , Υ + F Υ -µN , Υ ′ (4.18)
for some 0 < θ 1/2 and some κ > 0. We recall that ℓ is by definition the limit of Ẽ Υ n , Υ n+1 .

Recall our inequality (4.13) which says that µ n is uniformly bounded. Also, we know that Ẽ(Υ n , Υ n+1 ) converges to ℓ so, for n large enough, (Υ n , Υ n+1 ) must be in the neighborhood of the level set ℓ. Choosing µ = µ n+1 and µ ′ = µ n and using that G n and G n+1 have the correct trace, we get the estimate

Ẽ Υ n , Υ n-1 -ℓ 1-θ κ -1 F Υn -µ n+1 N , Υ n-1 + F Υn-1 -µ n N , Υ n = κ -1 F Υn -µ n+1 N , Υ n-1 -Υ n+1 C ||Υ n-1 -Υ n+1 ||
for n large enough. Here we have used that Υ n commutes with F Υn-1µ n N and that Υ n+1 commutes with F Υnµ n+1 N by construction, and that ||F Υn || and µ n+1 are both uniformly bounded. In order to conclude, we use the concavity of x → x θ and (4.15) like in [START_REF] Levitt | Convergence of gradient-based algorithms for the Hartree-Fock equations[END_REF] to obtain

Ẽ Υ n , Υ n-1 -ℓ θ -Ẽ Υ n , Υ n+1 -ℓ θ θ Ẽ Υ n , Υ n-1 -ℓ 1-θ Ẽ Υ n , Υ n-1 -Ẽ Υ n , Υ n+1 η θ 2 Ẽ Υ n , Υ n-1 -ℓ 1-θ ||Υ n+1 -Υ n-1 || 2 ηθ/(2C) ||Υ n+1 -Υ n-1 ||
by (4.15). This concludes the proof of the inequality (4.16), hence the proof of the convergence of (Υ 2n , Υ 2n+1 ), towards some pure HFB states (Υ, Υ ′ ).

Step 4: the limit (Υ, Υ ′ ) of (Υ 2n , Υ 2n+1 ) is a critical point of Ẽ. Since we have Υ 2n → Υ and Υ 2n+1 → Υ ′ , we deduce that F Υ2n → F Υ and F Υ2n+1 → F Υ ′ , by continuity of the map Υ → F Υ . Extracting a subsequence, we can assume that µ 2n k → µ ′ and µ 2n k +1 → µ.

We have

Υ 2n k = 1 (-∞,0) (F Υ2n k -1 -µ 2n k N), Υ 2n k +1 = 1 (-∞,0) (F Υ2n k -µ 2n k +1 N)
and, by uniform well-posedness,

F Υ2n k -1 -µ 2n k N η, F Υ2n k -µ 2n k +1 N η.
Passing to the limit k → ∞ we get

Υ = 1 (-∞,0) (F Υ ′ -µ ′ N) and Υ ′ = 1 (-∞,0) (F Υ -µN).
This exactly means that (Υ, Υ ′ ) is a critical point of Ẽ on P h (N/2) × P h (N/2). Note that we have also F Υ ′µ ′ N η and F Υ -µN η. The remaining statements are verified exactly like in the HF case. This concludes the proof of Theorem 4.1.

Optimal Damping Algorithm

In the previous section we have studied the convergence properties of the Roothaan algorithm, which consists in solving the self-consistent equation by a fixed point method. We have seen that the algorithm can either converge or oscillate between two states, none of them being a solution to the problem.

Examples of such oscillations in quantum chemistry have been exhibited by Cancès and Le Bris [START_REF] Cancès | Can we outperform the DIIS approach for electronic structure calculations?[END_REF][START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF]. In this case the potential W is repulsive and there is no pairing. In order to cure this problem of oscillations, Cancès and Le Bris proposed in [START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF] a relaxed algorithm called the Optimal Damping Algorithm (ODA). This method makes use of the important fact that one can minimize over mixed states and get the same ground state as when minimizing over pure states only (Theorem 2.1).

The same oscillations can a priori happen in HFB with an attractive potential W . They are frequently seen with the Roothaan algorithm and we will give several numerical examples later in Section 5. Even when the sequence Υ n eventually converges towards a single state Υ, these oscillations can slow down the convergence considerably. This phenomenon is well known in nuclear physics. Dechargé and Gogny already advocate in [START_REF] Dechargé | Hartree-Fock-Bogolyubov calculations with the D1 effective interaction on spherical nuclei[END_REF] the use of a damping parameter between two successive iterations, in order to "slow down the convergence on the density matrix. In this way the average field varies slowly and we can insure the convergence on the pairing tensor step by step" (see [START_REF] Dechargé | Hartree-Fock-Bogolyubov calculations with the D1 effective interaction on spherical nuclei[END_REF] page 1574). Even in the modern computations, this damping parameter is fixed all along the algorithm (Nathalie Pillet, private communication).

We suggest to transpose the method of Cancès and Le Bris to the HFB setting by using an optimal damping parameter, chosen such as to minimize the energy. This means resorting to mixed states even if the final ground state is always a pure HFB state. This is theoretically justified when the assumptions of the Bach-Fröhlich-Jonsson Theorem 2.2 are fulfilled.

The ODA involves two density matrices Υ n and Υn . The HFB state Υ n is always pure but Υn can (and will usually) be a mixed HFB state. The starting point Υ 0 = Υ0 being chosen, the sequence is then constructed by induction as follows:

(1) One finds (Υ n+1 , µ n+1 ) solving

Υ n+1 = 1 (-∞,0) F Υn -µ n+1 N and Tr(G n+1 ) = N/2.
This is always possible, by Lemma 4.1 and we can take as before

µ n+1 := ∂ -I Υn (N/2) + ∂ + I Υn (N/2) 2 , in case 0 is in the spectrum of F Υn -µ n+1 N. (2) One lets Υn+1 = t n+1 Υn + (1 -t n+1 )Υ n+1
where the damping parameter t n+1 ∈ [0, 1] is chosen such as to minimize the (quadratic) function

t → E t Υn + (1 -t)Υ n+1 .
(3) The algorithm is stopped when ||[Υ n , F Υn ]|| and/or ||Υ n+1 -Υ n || are smaller than a prescribed ε. The general strategy of the ODA is displayed in Figure 1. By construction we see that E( Υn ) is a non-increasing sequence. This guarantees the convergence of the ODA. The result is the following Theorem 4.2 (Convergence of the ODA). Assume that 0 < N/2 < N b . Let Υ 0 = Υ0 be an initial HFB state such that the sequence (Υ n ) generated by the ODA is uniformly well posed, that is ∀n,

Υ 0 = Υ0 Υ 1 = 1 (-∞,0) (F Υ -µ 1 N) Υ1 mixed HFB states pure HFB states Υ 2 = 1 (-∞,0) (F Υ1 -µ 2 N)
F Υn -µ n+1 N η > 0. (4.19)

Then

• The sequence E( Υn ) decreases towards a critical value of E;

• The sequence Υ n numerically converges towards a critical point Υ of E, in the sense that Υ n+1 -Υ n → 0, Υ n+1 -Υn → 0 and that all the limit points Υ of subsequences of

(Υ n ) solve Υ = 1 (-∞,0) (F Υ -µN).
Proof. The proof is exactly the same as in the Hartree-Fock case [START_REF] Cancès | SCF algorithms for HF electronic calculations[END_REF][START_REF] Cancès | Méthodes mathématiques en chimie quantique. Une introduction[END_REF] and we only sketch it. First we have by definition E( Υn+1 ) E( Υn ), so E( Υn ) must converge to a limit ℓ. Now we have

E( Υn+1 ) = E (1 -t n+1 ) Υn + t n+1 Υ n+1 = E( Υn ) -t n+1 a n+1 + t 2 n+1 b n+1 where t n+1 = argmin t∈[0,1] -ta n+1 + t 2 b n+1
and with

a n+1 := Tr F Υn Υn -Υ n+1 = Tr |F Υn -µN| Υn -Υ n+1 2 η Υn -Υ n+1 2 , b n+1 = 2 Tr( Gn+1 -G n )J( Gn+1 -G n ) -Tr( Gn+1 -G n )K( Gn+1 -G n ) + Tr( Ãn+1 -A n )K( Ãn+1 -A n ).
In finite dimension we have

|b n+1 | C Υn+1 -Υ n 2 (C/η)a n+1
. This can be used to prove that

-t n+1 a n+1 + t 2 n+1 b n+1 -ǫ a n+1 -ǫ η Υn -Υ n+1 2
for some ǫ > 0 independent of n. This now proves that

n Υn -Υ n+1 2 < ∞,
hence that Υn -Υ n+1 → 0. In order to conclude the proof, we notice that

Υ n+1 -Υ n = Υ n+1 -Υn + (1 -t n ) Υn-1 -Υ n which finally implies n ||Υ n -Υ n+1 || 2 < ∞ and n Υn -Υn+1 2 < ∞.
Since Υ n+1 = 1 (-∞,0) (F Υnµ n+1 N) by definition, the proof that any limit Υ of a subsequence of (Υ n ) satisfies the self-consistent equation is elementary.

Handling constraints

Both the Roothaan algorithm and the ODA are based on Lemma 4.1 which says that for any given F Υ , there exist µ ′ , δ ′ and Υ ′ such that

Υ ′ = 1 (-∞,0) F Υ -µ ′ N + δ ′ , Tr NΥ ′ = N -N b . (4.20)
The purpose of this section is to explain how to solve this problem numerically. To simplify our notation, we consider in this section a generic matrix

F = h p p -h , with p = p = p T and h = h = h T (4.21)
and we study the problem consisting in finding Υ, µ and δ such that

Υ = 1 (-∞,0) F -µN + δ, Tr NΥ = N -N b . (4.22)
Assume first that p ≡ 0 (Hartree-Fock case). Then we have

F = h 0 0 -h
The matrix F -µN is a linear function of µ ∈ R, hence by [START_REF] Kato | Perturbation theory for linear operators[END_REF], we know that its eigenvalues form a set of real analytic functions. They cannot be constant because the matrix N does not vanish. The eigenvalues of F -µN all behave like ±µ for large µ, by perturbation theory. We conclude that 0 can be an eigenvalue of F -µN for a finite number of µ's, say µ 1 < • • • < µ K . On the other hand, the map µ → 1 (-∞,0) F-µN is real-analytic outside of the µ k 's (see [START_REF] Kato | Perturbation theory for linear operators[END_REF]). So ν F is itself real-analytic outside of this set and it can have at most a finite number of jumps.

Outside of the µ k 's, it is possible to compute the derivative of ν F by usual perturbation methods [START_REF] Kato | Perturbation theory for linear operators[END_REF]. The answer is (4.24) and the fact that dν F /dµ 0 proves that ν F is increasing with respect to µ, in between these points. That the jumps are all positive can be easily seen by an approximation argument using Lemma 4.4 below. We skip the details.

The shape of the function ν F is very different in the HF and HFB cases. For a Hartree-Fock state, the function ν F is piecewise constant and it has jumps at the eigenvalues

ǫ 1 < • • • < ǫ N b of h G .
The size of the jumps is equal to the multiplicity of the associated eigenvalue. An HFB state will most always have a very smooth ν F . Of course, the smaller p in the Hamiltonian F, the more ν F looks like a step function.

In Figure 2 below, we show the function ν F for different values of the pairing term. More precisely, we have randomly chosen two symmetric real matrices h and p of size N b = 5, and we display the function ν F when the pairing is replaced by tp for t = 0 (Hartree-Fock case), t = 0.1 and t = 1. Figure 3 is a plot of the eigenvalues of F -µN for t = 0.1, as functions of µ. Note that there are some crossings of eigenvalues above and below the real line (recall that the spectrum is symmetric with respect to 0). But, around 0 the crossings are avoided and there is a gap.

If we repeat the numerical experiment with several random matrices h and p, we never see any jump for ν F . The purpose of the next result is to clarify this observation. Lemma 4.4 (Generic behavior of ν F ). The Fock matrix F -µN is invertible if and only if h ± ipµ are invertible. More precisely,

min σ |F -µN| = min (h + ip -µ) -1 -1 , (h -ip -µ) -1 -1 . (4.25)
The set of real symmetric matrices h and p such that

σ F -µN ∩ {0} = ∅ for all µ ∈ R is open and dense in {(h, p) : h = h T = h, p = p T = p}. For h and p in this set, ν F is real-analytic on R.
It is obvious that there are matrices h and p for which F -µN has 0 as eigenvalue for some µ ∈ R. The simplest examples are HF Hamiltonians for which p ≡ 0 and F -µN is not invertible each time µ equals an eigenvalue of h. If p does not vanish but commutes with h, then we have |h + ip -µ| 2 = |hip -µ| 2 = (hµ) 2 + p 2 and we see that 0 is never in the spectrum of F -µN when the kernel of p does not contain the eigenvectors of h. However there are counterexamples with p invertible not commuting with h. For instance, F + N is not invertible for

h = -1 0 0 2 , p = 0 2 2 0 .
We now turn to the proof of Lemma 4.4. Proof. The operator F -µN is unitarily equivalent to

1 √ 2 i √ 2 i √ 2 1 √ 2 F -µN 1 √ 2 -i √ 2 -i √ 2 1 √ 2 = -i 0 h + ip -µ -(h -ip -µ) 0 .
From this we deduce that F -µN is invertible if and only if h + ipµ and hipµ are invertible. Then we have

F -µN -1 = i 1 √ 2 -i √ 2 -i √ 2 1 √ 2 0 -(h -ip -µ) -1 (h + ip -µ) -1 0 1 √ 2 i √ 2 i √ 2 1 √ 2
and

F -µN -1 = max (h + ip -µ) -1 , (h -ip -µ) -1 .
The statement now follows from the fact that, on a dense open set, the spectra of h ± ip do not intersect the real axis.

Lemma 4.4 is interesting when we apply the Roothaan or the ODA, because it means that, as soon as p = 0, most often we will have no choice for µ n+1 and we will take δ n+1 = 0. Saying differently, it is really reasonable to assume that the sequence generated by the Roothaan and the ODA are uniformly well posed (of course when the final state is believed to have a non vanishing pairing), as we did in Theorem 4.1 and 4.2.

Even if it is in general smooth, the function ν can still vary quickly and this will be the case when the pairing term p is small. The appropriate method to find the solution of ν F (µ) = N/2 then depends on the properties of ν F . If the Hamiltonian F has a large enough pairing matrix p, then ν F is smooth and we can use a Newton-like method to solve the equation ν F (µ) = N/2. A trial chemical potential µ 0 being given, we compute the derivative ∂ν F /∂µ(µ 0 ) using Formula (4.24) and then let

µ 1 := µ 0 + N/2 -ν F (µ 0 ) ∂ν F ∂µ (µ 0 ) -1
.

The method can be iterated until convergence of µ n towards the desired µ. The convergence is very fast, as soon as ν F is smooth.

If the Hamiltonian F has a small pairing matrix p, the function ν F will be smooth but close to a step function. Its derivative varies very quickly and the previous Newton method is not appropriate. In this case we can use a simple bisection method. The bounds on ν F (µ) for large |µ| can be used to find a good starting interval [µ l , µ r ] such that ν F (µ l ) < N/2 and ν F (µ r ) > N/2.

We have to find a new µ n+1 at each step of the Roothaan or ODA. It is of course not efficient to find µ n+1 with a very high precision all along the algorithm. Dechargé and Gogny advice in Section II.E of [START_REF] Dechargé | Hartree-Fock-Bogolyubov calculations with the D1 effective interaction on spherical nuclei[END_REF] to apply the Newton scheme only once at each step. This means that

µ n+1 = µ n + N/2 -ν n (µ n ) ∂ν n ∂µ (µ n ) -1
where ν n is the function ν corresponding to F Υ = F Υn . This is then the same as doing perturbation theory on first order. We use a slightly different strategy which we explain in the next section.

Numerical results

In this section, we present some numerical results for two very simple interactions W . We start by considering in Section 5.2 a purely 3-dimensional gravitational model in which

W (x) = - g |x| , g > 0.
Then we consider in Section 5.3 a (repulsive) Coulomb potential which is perturbed at intermediate distances by an attractive effective potential, as is usually employed in nuclear physics:

W (x) = κ |x| -a 1 exp -b 1 |x| 2 + a 2 exp -b 2 |x| 2 , κ > 0.
In the next section we quickly explain our numerical technique to treat these two simple systems.

Method

To simulate our physical systems, we have used the open source software Scilab [START_REF]Scilab: The free software for numerical computation[END_REF]. Our potential W is always real, radial and spin-independent. To reduce the numerical cost we have therefore always imposed the spin, time-reversal and spherical symmetry. This means that we have to cope with ℓ max + 1 real and symmetric N b × N b matrices G ℓ and A ℓ (the one-particle density matrix and the pairing density matrix in the ℓth angular momentum sector). The total energy of the system is given by Equation (3.29) and we have to impose the constraints (3.27) and (3.28) which we recall here for convenience:

0 Υ ℓ SΥ ℓ Υ ℓ , with Υ ℓ := G ℓ A ℓ A ℓ S -1 -G ℓ and S = S 0 0 S , (5.1) ℓmax ℓ=0
(2ℓ + 1) Tr(SG ℓ ) = N/2. (5.2)
We choose a simple basis set (χ 1 , ..., χ N b ) of L 2 ([0, ∞), r 2 dr), made of "hat functions" associated with a chosen grid

0 = r 0 < r 1 < • • • < r N b < r N b +1 := r max .
We impose Dirichlet boundary conditions at r max . We have tested different types of grids and there was no important difference between them. The results presented here are all with regular grids. As we will explain later, for a given basis size N b , the results usually depend a lot on the value of the radius r max of the ball in which the system is placed.

Our main goal is to investigate the existence of pairing. We therefore always start by doing a precise Hartree-Fock calculation, for which we use the Optimal Damping Algorithm described in Section 4.2. We take as initial state a simple uniform state

G init = N 2 Tr(S) Id N b (5.3)
and we run HF until convergence. We have observed a global stability of the results with respect to initial states, hence the previous simple choice is appropriate (but more clever choices might decrease the total number of iterations). Then, we use the converged HF state G opt as initial datum for the HFB algorithm. Of course we have to perturb it a little bit since any HF solution is also an HFB solution. We proceed as follows. Assuming that the overlap matrix S = Id N b and that ℓ max = 0 for simplicity, the optimal HF state G opt can be written in the form

G opt = N/2 k=1 v k v T k ,
where v k are the N/2 first eigenvectors of the mean-field matrix h,

hv k = ǫ k v k .
We then choose a number n v of valence orbitals and a mixing parameter θ, and we perturb G opt as follows

G ′ init = N/2-nv k=1 v k v T k + θ N/2 k=N/2-nv +1 v k v T k + (1 -θ) N/2+nv k=N/2+1 v k v T k , A ′ init = θ(1 -θ) N/2+nv k=N/2-nv +1 v k v T k .
In most cases, we have observed that n v = 1 and θ = 0.95 works perfectly well, that is, the algorithm escapes from the HF solution G opt and converges towards an optimal HFB state. But other values of n v and θ seem to work fine also.

When the maximum angular momentum ℓ max is larger than 0, we often first run the algorithm with ℓ max = 0 for a few iterations before switching to the actual value of ℓ max . We stop the algorithm when the commutators [F ℓ n , Υ ℓ n ] are smaller than a prescribed error. We know from (3.30) and (3.31) that these commutators must all vanish for an exact solution of the discretized HFB minimization problem. In terms of the matrix S, the right quantity to look at is

ℓmax ℓ=0 S -1 2 F ℓ n Υ ℓ n S -SΥ ℓ n F ℓ n S -1 2
where ||•|| is the usual operator norm for (2N b ) × (2N b ) matrices. There is a similar formula in the HF case [START_REF] Cancès | Computational quantum chemistry: a primer[END_REF].

As we have explained, in the HFB case, ensuring the constraint (5.2) is not as easy as in the HF case. In the beginning of the algorithm, our state Υ is rather close to an HF state by construction. Therefore, the function ν F (µ) defined in Section 4.3 is close to a step function. We choose an error ε and look for the next states Υ ℓ n+1 having a total number of particles ℓmax ℓ=0 (2ℓ + 1) Tr(SG ℓ n+1 ) close to N/2, within the error ε, using a simple bisection method. We use the bisection for a fixed number of global iterations. Then, when the pairing term is large enough, we switch to a Newton method in order to find the state Υ n+1 . We have observed that even if in the beginning several Newton iterations can be employed at each step, usually only one Newton iteration is necessary after a while. To guarantee a good value of the average number of particles in the end, we decrease the error ε on | ℓmax ℓ=0 (2ℓ + 1) Tr(SG ℓ n+1 ) -N/2| along the algorithm.

Pure Newtonian interaction

Model

Here we consider a system of N spin-1/2 neutral particles, only interacting through the Newtonian interaction

W (x) = - g |x| , g > 0.
(5.4)

This potential is strongly attractive at short distances. Since 1/|x| does not decay too fast at infinity, it is also quite attractive at large distances. The kinetic energy does not scale the same as the potential energy. By a simple scaling argument, we can therefore always assume that

g ≡ 1.
This model can be used to describe neutron stars and white dwarfs when N ≫ 1. It has been particularly studied from a theoretical point of view in the pseudo-relativistic case where the kinetic energy is given by T [START_REF] Lieb | Gravitational collapse in quantum mechanics with relativistic kinetic energy[END_REF][START_REF] Lieb | The Chandrasekhar theory of stellar collapse as the limit of quantum mechanics[END_REF][START_REF] Lenzmann | Minimizers for the Hartree-Fock-Bogoliubov theory of neutron stars and white dwarfs[END_REF]. In our simulations we restrict ourselves to the non-relativistic case of the Laplacian T = -∆/(2m) (in units such that m = 1/2). It would be interesting to take N large but this is of course much too difficult from a numerical point of view.

= √ c 4 m 2 -c 2 ∆ -mc 2 , see
As mentioned before, we always impose the spin and time-reversal symmetries, which is perfectly justified for the ground state since the interaction (5.4) satisfies the assumption of the Bach-Fröhlich-Jonsson Theorem 2.2. We also impose spherical symmetry which, on the contrary, is not known to hold for the true ground state.

One advantage of the Newtonian interaction (5.4) is that the operators J and K ℓℓ ′ can be explicitely computed in the basis of hat functions. We have shown in Section 3.3.2 that the energy can be expressed in terms of

(ij|mn) ℓ,ℓ ′ = ∞ 0 r 2 dr ∞ 0 s 2 ds χ i (r) χ j (r) χ m (s) χ n (s) w ℓ,ℓ ′ (r, s)
where

w ℓ,ℓ ′ (r, s) = 1 2 1 -1 W r 2 + s 2 -2rst P ℓ (t) P ℓ ′ (t) dt = - 1 2 1 -1 P ℓ (t) P ℓ ′ (t) √ r 2 + s 2 -2rst dt.
Using the well-known formula

1 √ r 2 + s 2 -2rst = ∞ n=0
min(r, s) n max(r, s) n+1 P n (t) we deduce that

w ℓ,ℓ ′ = - 1 2 ∞ n=0 1 -1 P n P ℓ P ℓ ′ min(r, s) n max(r, s) n+1 .
The integral over the Legendre polynomials is related to the usual Clebsch-Gordan coefficients as follows

1 2 1 -1 P n (t) P ℓ (t) P ℓ ′ (t) dt = ℓ ℓ ′ n 0 0 0 2
and only a finite number of terms are non zero in the sum over n. The final result can be expressed as

(ij|mn) ℓ,ℓ ′ = - ∞ n=0 ℓ ℓ ′ n 0 0 0 2 ∞ 0 r 2 dr ∞ 0
s 2 ds min(r, s) n max(r, s) n+1 χ i (r) χ j (r) χ m (s) χ n (s). (5.5) These integrals can be explicitely computed for hat functions and 0 ℓ, ℓ ′ ℓ max with ℓ max not too large. In our numerical experiments we have put the explicit formulas in Scilab for ℓ max = 1. The integrals were stored in memory during the whole calculation.

Roothaan vs ODA

In the HF case, we have observed that the Roothaan algorithm very often oscillates between two states, none of them being the solution of the problem (as described in Theorem 4.1 and in [START_REF]On the convergence of SCF algorithms for the Hartree-Fock equations[END_REF]). The Roothaan algorithm seems more well behaved in the HFB case. With the model presented in this section, we never got real oscillations for HFB. Sometimes the convergence is improved by using the ODA, but in most cases the Roothaan algorithm always converges towards the same state as the ODA in the end. As we will see later, the situation is very different for the model studied in Section 5.3, which is inspired of nuclear physics.

We start by comparing Roothaan and ODA in the HF case. There, oscillations seem to be related to the size of the gap between the largest filled eigenvalue and the smallest unfilled one. Indeed, oscillations in HF seem to only occur when there is pairing in HFB, an effect which is also well-known to be related to the size of the gap (see, e.g., Theorem 5 in [START_REF] Bach | Bogolubov-Hartree-Fock mean field theory for neutron stars and other systems with attractive interactions[END_REF]). When there is no pairing, the HF Roothaan algorithm always behaves like the ODA. However, the situation is complex and there is no exact rule. Sometimes the Roothaan algorithm does not oscillate even when the gap is rather small and there is pairing.

In Figure 4 we display the value of the energy obtained along the algorithm for the Roothaan and the ODA, for the following choice of parameters: N = 6, N b = 200, ℓ max = 0 and r max = 30. The ODA converges in about 17 iterations, whereas the Roothaan algorithm oscillates. We also show the value of the norms ||G n -G n-1 || and ||G n -G n-2 || along the Roothaan algorithm. The oscillation between two points is clearly demonstrated.

When we decrease the parameter r max but keep N b = 200 constant, the gap is seen to increase slightly and the Roothaan algorithm behaves better. In Table 1, we give the numerical value of the last filled eigenvalue and the corresponding gap. The Roothaan algorithm slowly converges for r max = 25 and it coincides with the ODA when r max = 20. The gap for r max = 20 is 2.5 times the one for r max = 30. We will discuss the occurence of pairing in terms of the parameter r max in the next section. As we have mentioned the Roothaan algorithm is usually much more well behaved in the HFB case. However, sometimes the convergence can be improved dramatically by using the ODA. In Figure 5 we display the energy along the iterations of the algorithm in both the Roothaan and ODA cases, for N b = 500, N = 16, ℓ max = 1 and r max = 10. In this case the Roothaan algorithm is very badly behaved. It passes very close to the HF ground state and it takes it a very long time to escape from it. On the other hand, the ODA does not suffer from this problem and it converges much more rapidly.

r max ǫ N/2 ǫ N/2+1 -ǫ N/

Numerical evidence of pairing

Pairing effects in a finite discretization basis might depend on the properties of the basis. As we have explained, the occurence of pairing is related to the size of the gap in HF theory and this gap varies with the radius r max in which the system is confined. If r max is decreased the system is more condensed and the HF gap increases.

In Figure 6 we display the HF and HFB ground state energies computed for N = 16 in a basis set of size N b = 200, in terms of r max . The HF and HFB curves are distinct for r max large enough and they merge at r max = 6 approximately. This observation is confirmed by the value of the norm of A plotted on the right of the same figure. The minima of the HF and HFB ground state energies are attained at about r max ≃ 10 in the HF case and r max ≃ 10.5 in the HFB case, which is sufficiently far from the merging point. The minima of these curves correspond to the best possible approximation for a given basis size N b (here N b = 200) and a given type of grid (here regular). The difference between the corresponding HF and HFB energies is significant. The HF ground state energy at r max = 10 is -19.232176 (in our units in which m = 1/2 and e = 1), whereas the HFB ground state energy at r max = 10.5 is -19.240176. The norm of the pairing matrix A is rather large at this point:

||A|| = Tr(SA 0 SA 0 ) + 3 Tr(SA 1 SA 1 ) ≃ 0.462129.
This goes in favour of the conclusion that pairing really occurs for N = 16 in this model. This intuition is confirmed by a more precise calculations with N b = 500 which we discuss below.

The observation of pairing requires to have an appropriate r max but it does not require to have a very large basis set. Even for N b = 30 and r max = 10.5, we already find that the HFB energy is approximately -19.078416 whereas the HF energy is about -19.072954. The corresponding norm of the pairing matrix A is ||A|| ≃ 0.424124.

Pairing is a subtle effect which decreases the energy by a small amount (much less than one percent here). Catching this effect requires to be very careful when choosing the radius r max . Taking r max too small might lead to the conclusion that there is no pairing. In our simulations we have always observed the occurence of pairing, but provided we choose r max appropriately. The values of r max at which the HF and HFB energies attain their minimum were always found on the right of the merging point of the two curves. In Table 2 below we give our results for N b = 200 and N = 6, 10, 16 and 20. The HFB ground state energy is always smaller than the HF energy.

In the paper [START_REF] Lenzmann | Minimizers for the Hartree-Fock-Bogoliubov theory of neutron stars and white dwarfs[END_REF], Lenzmann and Lewin have rigorously studied the gravitational model of this section. They showed the existence of a ground state in both the HF and HFB cases. But, so far, no proof that pairing occurs has been provided. The numerical results of this section tend to show that there is actually always pairing, at least for N not too large.

Properties of the HFB ground state

In Table 2 below we give our results for N b = 200 and N = 6, 10, 16 and 20, for the optimal values of r max . With ℓ max = 1 we have observed that the shells are filled alternatively. In HF theory, the cases N = 10 and N = 16 correspond to closed shells, whereas for N = 6 and N = 20 the last shell is only partially filled. This is a simple explanation for the fact that the pairing matrix is much bigger in these cases.

In Table 3 we display the occupation numbers for the optimal HFB ground state in the

Quality of the approximation in terms of the number N b of points

In Table 4 we display the HF and HFB ground state energies for N = 16, ℓ max = 1 for the the optimal value of r max , in terms of the number of discretization points N b of the regular grid. The convergence is not very fast, but we see that the difference between the HF and the HFB energy, as well as the norm of the pairing matrix are of the same order for small N b as they are for larger N b 's. From this observation we can conclude that it is probably not necessary to take N b very large in order to decide whether pairing occurs or not.

A simplified model for protons and neutrons

In this section we report on our numerical results concerning a simple model inspired of nuclear physics. The interaction between protons and neutrons is not a fundamental law of nature because these are composite particles made of quarks, which interact through weak, strong and electrostatic forces. A common procedure used in nuclear physics is to use empirical forces [START_REF] Ring | The nuclear many-body problem[END_REF] which involve a small number of parameters which are fitted to experiment or to the known behavior of the model in some limits. The most common forces used in practice are the so-called Skyrme [START_REF] Skyrme | The effective nuclear potential[END_REF] and Gogny [START_REF] Gogny | Proceedings of the International Conference on Nuclear Physics[END_REF][START_REF]Proceedings of the International Conference on Nuclear Self-Consistent Fields[END_REF][START_REF] Dechargé | Hartree-Fock-Bogolyubov calculations with the D1 effective interaction on spherical nuclei[END_REF] forces and they depend nonlinearly on the state itself. Here we consider an effective force which is fixed and does not depend on the quantum state. We also take it spin-independent and isospin-independent.

Our goal is to test some simple ideas and not to do a real nuclear physics calculation.

Model

The nucleon-nucleon potential has been observed to be repulsive at short distances and only attractive at medium distances. It decays very fast at infinity. A simple choice to describe this is to take W (x) = κ |x| a 1 e -b1|x| 2 + a 2 e -b2|x| 2 , (5.6) with a 2 , a 1 > 0, b 1 < b 2 . The constant κ is 1 for the proton-proton interaction and 0 for the proton-neutron and the neutron-neutron interaction. The other constants usually also depend on the isospin (the quantum variable which determines whether a nucleon is a neutron or a proton). For simplicity we work here with particles having a definite isospin. This means that we assume to have either only protons or only neutrons. In particular we want to ask for which strength of the effective force it becomes possible for the protons to overcome their Coulomb repulsion and form a bound state. In reality a nucleus is made of a certain number of protons and neutrons and one has to use a different HFB state for each species.

In our applications we have chosen for simplicity b 1 = 1, b 2 = 4, a 1 = a = 2 a 2 /3. This means that the effective force takes the form

W (x) = κ |x| + a 3 2 e -4|x| 2 -e -|x| 2 .
(5.7)

When κ = 1, this force is purely repulsive for a 2.87 and it becomes attractive at intermediate distances for larger a's. The corresponding force is displayed in Figure 7 for a = 1 and κ = 0. One can ask several questions concerning the model considered in this section:

(1) For which value of a does a system of N identical nucleons bind in Hartree-Fock theory? (2) Is there always pairing when there is binding?

(3) Can pairing effects allow for binding with a smaller a than in HF theory?

These questions are mostly of academic nature for the very simplified model considered in this section. But investigating the same problems with more realistic forces is very important from a physical point of view. From a mathematical point of view, nothing seems to be known for simple models of the same form as in this section. It is not even known that binding always occurs for a large enough with the previous interaction. We hope that our calculations will stimulate some further mathematical studies.

Some computational details

We always minimize over states having the spin, time-reversal and rotation symmetries. The Bach-Fröhlich-Jonsson Theorem 2.2 does not apply to the model of this section, hence we are making a further approximation here.

For such symmetric states we have shown in Section 3.3.2 that the energy can be expressed in terms of (ij|mn) ℓ,ℓ ′ = ∞ 0 r 2 dr ∞ 0 s 2 ds χ i (r) χ j (r) χ m (s) χ n (s) w ℓ,ℓ ′ (r, s)

where, for the model considered in this section,

w ℓ,ℓ ′ (r, s) = 1 2 1 -1 W r 2 + s 2 -2rst P ℓ (t) P ℓ ′ (t) dt = 1 2 1 -1 P ℓ (t) P ℓ ′ (t) κ √ r 2 + s 2 -2rst
a 1 e -b1(r 2 +s 2 -2rst) + a 2 e -b2(r 2 +s 2 -2rst) dt.

For 0 ℓ, ℓ ′ ℓ max with ℓ max not too large, the Gaussian integrals can be computed exactly and it is possible to find the exact expression of w ℓ,ℓ ′ (r, s).

The computation of the integral (ij|mn) ℓ,ℓ ′ against hat functions is much more tedious, however. It is easy to find an exact expression for the Coulomb part, but not so simple for the Gaussian part. So we have performed a numerical calculation of these integrals. Since we have of the order of (N b ) 4 integrals, we could not take N b too large. The results of the previous section indicated that the existence of pairing effects does not depend very much on the size of the basis.

Slow convergence and oscillations of Roothaan

We have observed that the Roothaan algorithm almost always oscillates, even in the HFB case (see some examples in Figures 8,9 and 10). This is in stark contrast with the results of the previous section where the Roothaan algorithm was almost always converging. Sometimes it very slowly converges in the HF case (see, e.g., Figure 9). However we have always obtained convergence for the HF Roothaan algorithm when a is small enough, that is, when it is expected that there is actually no binding. For the case displayed in Figure 9 we have a = 20 but the critical a is about ≃ 24 (see the next section).

We conclude that using the ODA is very important for such attractive potentials. The same might be true with the more involved forces used in nuclear physics.

The critical strength

In finite dimension there is always a minimizer. Saying differently, since the particles are trapped in a ball, they always bind. Furthermore we work with rotation-invariant states. So, for the true model in infinite dimension, the particles escaping to infinity cannot form a bound state of the same kind because they are too far from the (fixed) center of symmetry. In this special case they will spread out and have a vanishing energy.

In Hartree-Fock theory, we conclude that we can detect the loss of binding by looking at the last filled HF eigenvalue. When it crosses 0, this corresponds to the last particle becoming a scattering state. We can therefore choose as definition for the critical strength a, the value at which this eigenvalue is 0. In finite dimensional Hartree-Fock-Bogoliubov theory things are less clear and we will not discuss the problem of binding. In our simulations we have observed that the HFB ground state density was always rather close to the HF ground state density, which suggests that there is binding in HFB as well.

We have made some calculations for N b = 50 and N = 4. We found that the critical strength is about a c ≃ 23.5 in the proton-proton case and a c ≃ 17.5 in the neutron-neutron case. In Figure 11 we display the HF and HFB energies as functions of the parameter a, for N = 4 and κ = 1 (proton-proton case). Figure 12 is the equivalent result for κ = 0 (neutron-neutron case). For these calculations we have chosen r max = 3 which is the optimal choice for a in a neighborhood of the critical value. Like in the previous section the results depend on the radius of the ball in which the system is confined. We see that there is always pairing, in the sense that the HFB curve is below the HF curve. This is even more manifest in the neutron-neutron case for which the potential is much more attractive than for protons, which repel with the Coulomb potential. Also, the norm of the pairing matrix A does not vary too much with a, it stays between 0.80 and 0.95 for a in the range 15 a 30, for both κ = 0 and κ = 1.

From these numerical results we can conclude that pairing seems to happen in this model, for any strength a for which there is binding in Hartree-Fock theory. It is an interesting problem to actually prove that pairing always occurs, for instance for a large enough. We are not aware of any result of this kind.

0 0 1 can

 1 be written for the matrices G and A in the form S are defined by Υ := G A A * S -1 -G and S := S 0 0 S . (3.10)

  and a = a T = a (3.21) (the 2 × 2 matrix refers to the spin variables), if and only if Σ k ΓΣ * k = Γ, for k = 1, 2, and Γ = Γ. (3.22)

r 2 dr and ℓmax ℓ=0 (2ℓ + 1 )

 ℓ=01 Tr(SG ℓ ) = N/2. (3.28)

Lemma 4 . 1 (

 41 The Roothaan algorithm is well defined). The function λ ∈ [0, N b ] → I Υ (λ) is convex, hence left and right differentiable. For any λ ∈ [0, N b ], the minimizers Υ ′ of I Υ (λ) are exactly the states of the form

Theorem 4 . 1 (

 41 Convergence of the Roothaan algorithm). Assume that 0 < N/2 < N b . Let Υ 0 be an initial HFB state such that the sequence (Υ n ) generated by the Roothaan algorithm is uniformly well posed. Then• The sequence Ẽ(Υ 2n , Υ 2n+1 ) decreases towards a critical value of Ẽ;

Fig. 1 .

 1 Fig. 1. The Optimal Damping Algorithm of Cancès & Le Bris in the HFB case.

Fig. 2 .

 2 Fig.2. The function ν F (µ) which gives the average number of particles in the state 1 (-∞,0) (F -µN), in terms of the chemical potential µ. The pairing term in F is equal to tp with t = 0 (Hartree-Fock case, red curve), t = 0.1 (blue curve) and t = 1 (green curve).

Fig. 3 .

 3 Fig. 3. The eigenvalues of F -µN in terms of µ for t = 0.1.

Fig. 4 .

 4 Fig. 4. Left: HF energy along the iterations for the Roothaan Algorithm (blue) and the ODA (red). Right: Values of ||Gn -G n-1 || (red) and ||Gn -G n-2 || (blue) along the Roothaan algorithm, showing the oscillations between two states. Here N = 6, N b = 200, ℓmax = 0 and rmax = 30.

Fig. 5 .

 5 Fig. 5. HFB energy along the iterations of the Roothaan (blue) and the ODA (red) for N b = 500, N = 16, ℓmax = 1 and rmax = 10. The optimal HF energy is also displayed (green).

Fig. 6 .

 6 Fig. 6. Value of the ground state HF and HFB energies (top) and of the norm of the pairing matrix A (bottom), as functions of the radius rmax in which the system is confined, for N = 16, N b = 200 and ℓmax = 1.

Fig. 7 .

 7 Fig. 7. The effective force -e -|x| 2 + 3 e -4|x| 2 /2 used in our calculation of Section 5.3. The (repulsive) Coulomb potential must be added for protons.

Fig. 8 .

 8 Fig. 8. Energy along the iterations in the HF (left) and HFB (right) cases, for the Roothaan (blue) and the ODA (red), with N = 4, N b = 20, ℓmax = 1, rmax = 3, a = 35 and κ = 1 (proton-proton case).

Fig. 9 .

 9 Fig.9. Same calculation with a = 20 and κ = 1 (proton-proton case). The Roothaan algorithm slowly converges in the HF case and it oscillates in the HFB case but the two values are very close.

Fig. 10 .

 10 Fig. 10. Same calculation with a = 20 and κ = 0 (neutron-neutron case). The Roothaan algorithm oscillates in the HFB case, but the two values are very close.

Fig. 11 .

 11 Fig. 11. Left: Values of the HF (blue) and HFB (red) ground state energies as functions of a, with N = 4, N b = 50, ℓmax = 1, rmax = 3 and κ = 1 (proton-proton case). The vertical line is the value of a for which the last filled eigenvalue vanishes. Right: Values of the two filled HF eigenvalues for the same a.

Fig. 12 .

 12 Fig. 12. Same calculations for κ = 0 (neutron-neutron case).

Table 1 .

 1 Value of the last filled eigenvalue ǫ N/2 and the corresponding gap ǫ N/2+1 -ǫ N/2 in HF, for N = 6, N b = 200 and ℓmax = 0. For rmax = 30 the Roothaan algorithm oscillates and we display the last filled eigenvalue and the gap for the two states.

				2 behavior of HF Roothaan
	20	-0.532430	0.159430	fast convergence
	25	-0.536706	0.081016	slow convergence
	30	-0.529200	0.061928	oscillations
		-0.548554	0.067422

Table 2 .

 2 Results for N b = 200 and ℓmax = 1.

	). This

Table 3 .

 3 Occupation numbers of the HFB minimizer, for N b = 200 and ℓmax = 1.

Table 4 .

 4 Value of the HF and HFB energies for N = 16 and ℓmax = 1 and the (approximate) optimal rmax.

	N b 30	r max HF energy HFB energy difference ||A|| 9 -19.112314 -19.117948 0.005634 0.425604
	50	9	-19.189066 -19.196012	0.006946	0.445728
	100 9	-19.222300 -19.229872	0.007572	0.454173
	150 10	-19.229494 -19.237574	0.008080	0.461725
	200 10	-19.232176 -19.240308	0.008132	0.462363
	250 10	-19.233420 -19.241576	0.008156	0.462659
	300 10	-19.234094 -19.242264	0.008170	0.462821
	400 11	-19.234826 -19.243068	0.008242	0.463905
	500 11	-19.235206 -19.243456	0.008250	0.463985

which commutes with N. The solution of (4.22) is then given by the aufbau principle,

where µ is the (N/2)th eigenvalue of h, counted with multiplicity and δ lives in the corresponding eigenspace. Equivalently,

where the v i 's solve the eigenvalue equation

is the dimension of the direct sum of all the eigenspaces corresponding to the eigenvalues < ǫ N/2 and K ′ = Tr 1 (-∞,ǫ N/2 ] (h) is the dimension of the direct sum of all the eigenspaces corresponding to the eigenvalues ǫ N/2 . The n i 's are chosen such that

Therefore, finding Υ, µ and δ in the Hartree-Fock case only requires to diagonalize h once.

In the Hartree-Fock-Bogoliubov case (p = 0), the situation is more complicated since N does not commute with F. Let us consider the real function

We are interested in solving the equation

In the Hartree-Fock case, ν F is a non-decreasing piecewise constant function. There is a solution µ to ν F (µ ′ ) = N/2 when N/2 belong to the range of ν F . Otherwise, one has to partially fill a shell using the matrix δ.

In the Hartree-Fock-Bogoliubov case, ν F is also non-decreasing and in general it is much smoother when p = 0. The following lemma summarizes some important properties of ν F in both the HF and HFB cases.

where (F -µN)v i = ǫ i v i .

Proof. The behavior of ν F for |µ| ≫ 1 was already studied in Lemma 4.2.