Thierry Jecko
email: jecko@math.cnrs.fr

Camille Noûs
email: camille.nous@cogitamus.fr

On the analyticity of electronic reduced densities for molecules

Keywords: Analytic elliptic regularity, molecular Hamiltonian, electronic reduced densities, electronic reduced density matrices, Coulomb potential, twisted differential calculus

We consider an electronic bound state of the usual, non-relativistic, molecular Hamiltonian with Coulomb interactions and fixed nuclei. Away from appropriate collisions, we prove the real analyticity of all the reduced densities and density matrices, that are associated to this bound state. We provide a similar result for the associated reduced current density.

1 Introduction.

The notions of density and density matrices are quite old tools in the treatment of physical many body problems. For instance, the Thomas-Fermi theory goes back to the year 1927 (cf. [Li]). Nowadays these tools do play a central rôle in one important approach of the molecular problem, namely the Density Functional Theory (DFT) (cf. [E, LiSe]). Regularity properties of the related reduced density and density matrices are useful to develop the DFT. For instance, the real analyticity of a density was recently used for Hohenberg-Kohn Theorems in [G]. The purpose of the present paper is to prove, in a large region of the configuration space, the real analyticity of all reduced densities and density matrices of a pure state for a molecular Hamiltonian with fixed nuclei.

We consider a molecule with N moving electrons, with N > 1, and L fixed nuclei, with L ≥ 1 (Born-Oppenheimer idealization). While the L distinct vectors R 1 , • • • , R L ∈ R 3 denote the positions of the nuclei, the positions of the electrons are given by x 1 , • • • , x N ∈ R 3 . The charges of the nuclei are respectively given by the positive Z 1 , • • • , Z L and the electronic charge is -1. In this picture, the Hamiltonian of the electronic system is

H := N j=1 -∆ x j - L k=1 Z k |x j -R k | -1 + 1≤j<j ≤N |x j -x j | -1 + E 0 ,(1.1) where E 0 = 1≤k<k ≤L Z k Z k |R k -R k | -1
and -∆ x j stands for the Laplacian in the variable x j . Here we denote by | • | the euclidian norm on R 3 . Setting ∆ := N j=1 ∆ x j , we define the potential V of the system as the multiplication operator satisfying H = -∆ + V . Thanks to Hardy's inequality

∃c > 0 ; ∀f ∈ W 1,2 (R 3 ) , R 3 |t| -2 |f (t)| 2 dt ≤ c R 3 |∇f (t)| 2 dt , (1.2)
one can show that V is ∆-bounded with relative bound 0. Therefore the Hamiltonian H is self-adjoint on the domain of the Laplacian ∆, namely W 2,2 (R 3N ) (see Kato's theorem in [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF], p. 166-167).

From now on, we fix an electronic bound state ψ ∈ W 2,2 (R 3N ) \ {0} such that, for some real E, Hψ = Eψ. We point out (cf. [S, Z]) that such a bound state exists at least for appropriate E ≤ E 0 (cf. [FH]) and for

N < 1 + 2 L k=1 Z k .
Associated to that bound state ψ, we consider the following objects. Let k be an integer such that 0 < k < N . Let ρ k : (R 3 ) k → R be the almost everywhere defined, L 1 (R 3k )function given by, for x = (x 1 ;

• • • ; x k ) ∈ R 3k , ρ k (x) = R 3(N -k)
ψ(x; y) 2 dy .

(1.3) Define γ k : (R 3 ) 2k → C as the almost everywhere defined function given by, for x = (x 1 ;

• • • ; x k ) ∈ R 3k and x = (x 1 ; • • • ; x k ) ∈ R 3k , γ k (x; x ) = R 3(N -k)
ψ(x; y)ψ(x ; y) dy .

(1.4)

By [K], ψ is actually a continuous function. Therefore, the objects ρ k and γ k are everywhere defined and, ρ k (x) = γ k (x; x) holds true everywhere. It turns out that γ k may be seen as the kernel ("Green function") of a trace class operator (see [Le, LiSe]). We call the function ρ k the k-particle reduced density and the kernel γ k the k-particle reduced density matrix.

When k = 1, ρ := ρ 1 is often simply called the (electronic) density. In this case, we also introduce the reduced current density, defined as the almost everywhere defined, L 1 (R 3 )function C : R 3 → R 3 given by, for x ∈ R 3 , the following imaginary part

C(x) = R 3(N -1)
∇ x ψ(x; y)ψ(x; y) dy .

(1.5)

From a physical point of view, the previous objects differ from the true physical ones by some prefactor (see [E, Le, LiSe, LSc]).

It is useful to introduce the following subsets of R 3k . Denoting for a positive integer p by [[1; p]] the set of the integers j satisfying 1 ≤ j ≤ k, the closed set

C k := x = (x 1 ; • • • ; x k ) ∈ R 3k ; ∃(j; j ) ∈ [[1;
k]] 2 ; j = j and x j = x j gathers all possible collisions between the first k electrons while the closed set

R k := x = (x 1 ; • • • ; x k ) ∈ R 3k ; ∃j ∈ [[1; k]] , ∃ ∈ [[1; L]] ; x j = R
groups together all possible collisions of these k electrons with the nuclei. We set U

(1)

k := R 3k \ (C k ∪ R k ),
which is an open subset of R 3k . Now, we consider two sets of positions for the first k electrons and introduce the set of all possible collisions between positions of differents sets, namely

C (2) k :=    (x; x ) ∈ (R 3k ) 2 ; x = (x 1 ; • • • ; x k ) , x = (x 1 ; • • • ; x k ) , ∃(j; j ) ∈ [[1; k]] 2 ; x j = x j    .
We introduce the open subset of (R 3k ) 2 defined by

U (2) k = U (1) k × U (1) k \ C
(2) k .

In the k = 1 case, we note that C 1 = ∅, U

(1)

1 = R 3 \ {R 1 , • • • , R L } and U (2) 1 = R 3 \ {R 1 , • • • , R L } 2 \ D ,
where D = {(x; x ) ∈ (R 3 ) 2 ; x = x } is the diagonal of (R 3 ) 2 . We have the following Theorem 1.1. [START_REF] Fournais | The electron density is smooth away from the nuclei[END_REF][START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF][START_REF] Th | A new proof of the analyticity of the electronic density of molecules[END_REF].

The one-particle reduced density ρ := ρ 1 is real analytic on U

(1)

1 = R 3 \ {R 1 , • • • , R L }.
The original proof was given in [START_REF] Fournais | The electron density is smooth away from the nuclei[END_REF][START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF] and relies on a clever decomposition of the integration domain of the integral defining ρ into pieces, on which one can perform an appropriate change of variables. A different proof of Theorem 1.1 is provided in [J].

We refer to [START_REF] Fournais | Non-isotropic cusp conditions and regularity of the electron density of molecules at the nuclei[END_REF][START_REF] Fournais | Estimates on derivatives of Coulombic wave functions and their electron densities[END_REF] for a refined information on ρ. In [START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF], a similar analyticity result on γ 1 was announced, without proof. More precisely, it was claimed there that γ 1 should be real analytic on U

(1)

1 ×U

(1)

1 . Adapting the method of [START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF], a shlighly weaker result on γ 1 was proved in [HS], namely Theorem 1.2. [HS]. The one-particle reduced density matrix γ 1 is real analytic on U

(2) 1 .

Note that

U (2) 1 is indeed a strict subset of U (1) 1 × U (1) 1 .
In [START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF], it was also announced without proof that ρ 2 should be real analytic on

U (1) 2 = R 6 \ (C 2 ∪ R 2 ) = R 3 \ {R 1 , • • • , R L } 2 \ D = U (2)
1 .

In the present paper, we work with the tools and the method, used in [J] to get Theorem 1.1, to prove the following three results.

Theorem 1.3. The reduced current density C is real analytic on U

(1)

1 = R 3 \ {R 1 , • • • , R L }.
Theorem 1.4. For all integer k with 0 < k < N , the k-particle reduced density ρ k is real analytic on U

(1) k .

Theorem 1.5. For all integer k with 0 < k < N , the k-particle reduced density matrix γ k is real analytic on U

(2)

k = (U (1) k × U (1) k ) \ C (2) k .
Theorem 1.4 for k = 1 coincides with Theorem 1.1. The proof of Theorem 1.4 below is, when restricted to the k = 1 case, just a rewriting of the proof of Theorem 1.1 in [J] and provides a direct justification of Theorem 1.3. We note that Theorem 1.5 for k = 1 is exactly Theorem 1.2 and that Theorem 1.4 with k = 2 fits precisely to the announced result on ρ 2 in [START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF].

The set of all possible collisions between particles is C N ∪ R N and the potential V is real analytic precisely on R 3N \ (C N ∪ R N ). Classical elliptic regularity applied to the equation Hψ = Eψ shows that ψ is also real analytic on R 3N \ (C N ∪ R N ). We refer to [START_REF] Ammann | Regularity for eigenfunctions of Schrödinger operators[END_REF][START_REF] Fournais | Non-isotropic cusp conditions and regularity of the electron density of molecules at the nuclei[END_REF][START_REF] Fournais | Positivity of the spherically averaged atomic one-electron density[END_REF][START_REF] Fournais | Analytic structure of many-body coulombic wave functions[END_REF][START_REF] Fournais | Estimates on derivatives of Coulombic wave functions and their electron densities[END_REF][START_REF] Kato | On the eigenfunctions of many-particle systems in Quantum Mechanics[END_REF] for more information on the behaviour of ψ near collisions. It turns out that, at least at some places in C N ∪ R N , ψ is not real analytic (cf. [START_REF] Fournais | Non-isotropic cusp conditions and regularity of the electron density of molecules at the nuclei[END_REF][START_REF] Kato | On the eigenfunctions of many-particle systems in Quantum Mechanics[END_REF]). Thus, in the definition (1.3) of ρ k , for an integer k with 0 < k < N , the domain of integration does contain such singularities. Therefore it is not clear a priori that ρ k is real analytic somewhere. This is however true, by Theorem 1.4, on the complement of the set of all possible collisions of the "external" electrons, namely on U

(1)

k . The proofs of Theorems 1.1 and 1.2 do not directly use the analyticity of ψ. This is also the case of the proofs of Theorems 1.3, 1.4, and 1.5, below.

It turns out that our proofs of Theorems 1.4 and 1.5 have a common structure that we want to describe now. Using an appropriate, x-dependent, unitary operator U x , that acts on L 2 (R 3(N -k) ), we locally transform the equation Hψ = Eψ into an elliptic differential equation (U x HU -1

x )U x ψ = EU x ψ (the "twisted" equation), the coefficients of which have nice analytic properties. Applying elliptic regularity to the latter equation, we obtain the real analyticity of x → U x ψ ∈ L 2 (R 3(N -k) ) that yields the real analyticity of ρ k : x → U x ψ 2 ( • being the norm on L 2 (R 3(N -k) )). For γ k , we follow the same lines but with a (x; x )-dependent, unitary operator U (x;x ) , and two "twisted" differential equations in the variables (x; x ; y). These proofs are performed in Section 4. Let us now compare the two available methods to prove such results, namely the one mentioned just above and the one used in [START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF][START_REF] Hearnshaw | Analyticity of the one-particle density matrix[END_REF] to get Theorems 1.1 and 1.2. We point out that elliptic regularity was also an important tool in [START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF][START_REF] Hearnshaw | Analyticity of the one-particle density matrix[END_REF]. Another common tool is the use of a x-dependent change of variables on the "internal" variables y ∈ R 3(N -k) , since the above unitary operators are both the implementation in L 2 (R 3(N -k) ) of such a change of variables. However, we use here global changes of variables in contrast to the local ones in [START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF][START_REF] Hearnshaw | Analyticity of the one-particle density matrix[END_REF]. While, there, the integrals on R 3(N -1) defining ρ 1 and γ 1 are split in several pieces, on which an appropriate change of variables is performed to make the real analyticity apparent, we view here ρ k (resp.

γ k ) as the L 2 (R 3(N -k) )-norm (resp. the L 2 (R 3(N -k) )-scalar product) of one (resp. two) L 2 (R 3(N -k) )-valued, real analytic function(s).
In view of the behaviour of ψ near some points in C N ∪ R N , that was established in [START_REF] Fournais | Analytic structure of many-body coulombic wave functions[END_REF], it is quite natural to expect that, for an integer k with 0 < k < N , neither ρ k nor γ k is real analytic everywhere. Their exact domain of real analyticity is still an open question. The domains obtained in Theorems 1.4 and 1.5 are the largest ones that can be reached by the method used in this paper, as explained in Remarks 4.2 and 4.4.

As already pointed out in [J], the twisted calculus, that we use here, allows us to treat more general Hamiltonians than H. We introduce, in Section 5, a class of elliptic, second order differential operators with Coulomb singularities. We observe that the conjugation by the twist actually preserves the ellipticity. Assuming that there exists a bound state for such an operator, we show in Theorem 5.1 that our previous, main results hold true for it. Moreover, we expect that this should extend to a class of elliptic, pseudo-differential operators with Coulomb singularities, since the used tools are available for such operators. We refer to [MS] for more information on the twisted (pseudo-)differential calculus.

Maybe it would be useful for the DFT to extend the present results on pure states to general states (see Sections 3.1.4 and 3.1.5 in [LiSe] for details). We expect that Theorem 3.9 below could be a good starting point for this purpose.

The paper is organized as follow. In Section 2, we provide a general notation, basic facts on real analyticity, and a result on elliptic regularity from [J]. In Section 3, we review the twisted calculus based on Hunziker's twist, show how it can be used to "desingularize" a differential equation with Coulomb singularities, and apply elliptic regularity to the corresponding twisted equation. Section 4 is devoted to the proof of our main results, namely Theorems 1.3, 1.4, and 1.5. The extension of these results to a larger class of Hamiltonians is performed in Section 5. Finally, we gathered in Appendix A some elementary computations that are used in the main text and provide in Appendix B basic material on global pseudo-differential operators.

Acknowledgments: The authors thanks S. Fournais and T. Østergaard Sørensen for fruitful discussions related to the subject of this paper. Dedicated to K.

2 Basic tools and elliptic regularity.

To prepare the proof of the main results, we recall basic tools and state a general result on the analytic regularity of solutions of some elliptic equations.

Notation and basic tools.

We start with a general notation. We denote by C the field of complex numbers, and, for z ∈ C, z (resp. z) stands for the real part (resp. imaginary part) of z. Let p be a positive integer. Recall that, for u ∈ R p , we write |u| for the euclidian norm of u. Given such a vector u ∈ R p and a nonnegative real number r, we denote by B(u; r [ (resp. B(u; r]) the open (resp. closed) ball of radius r and centre u, for the euclidian norm | • | in R p . In the one dimensional case, we use the following convention for (possibly empty) intervals:

for (a; b) ∈ R 2 , let [a; b] = {t ∈ R; a ≤ t ≤ b}, [a; b[= {t ∈ R; a ≤ t < b}, ]a; b] = {t ∈ R; a < t ≤ b},
and ]a; b[= {t ∈ R; a < t < b}. We denote by N the set of nonnegative integers and set N * = N \ {0}. If p ≤ q are non negative integers, we set [[p; q]] := [p; q] ∩ N, [[p; q[[ = [p; q[∩N, ]]p; q[[ =]p; q[∩N, and ]]p; q]] :=]p; q] ∩ N. Given an open subset O of R p and n ∈ N, we denote by W n,2 (O) the standard Sobolev space of those L 2 -functions on O such that, for n ∈ [[0; n]], their distributional derivatives of order n belong to L 2 (O). We denote by •, • n (resp. • n ) the right linear scalar product (resp. the norm) on W n,2 (O). In particular, W 0,2 (O) = L 2 (O). Without reference to p and O, we denote by • (resp. •, • ) the L 2 -norm (resp. the right linear scalar product) on L 2 (O). For two Banach spaces (A, • A ) and (B, • B ), the space L(A; B) of continuous linear maps from A to B is also a Banach space for the operator norm • L(A; B) defined by

∀M ∈ L(A; B) , M L(A; B) := sup a∈A\{0} M (a) B a A .
We simply denote L(A; A) by L(A).

Let p be a positive integer and O an open subset of R p . Let (A, • A ) be a Banach space and θ : O

x = (x 1 ; • • • ; x p ) → θ(x) ∈ A. For j ∈ [[1; p]],
we denote by ∂ j θ or ∂ x j θ the j'th first partial derivative of θ. For α ∈ N p and x ∈ R d , we set

D α x := (-i∂ x ) α := (-i∂ x 1 ) α 1 • • • (-i∂ xp ) αp , D x = -i∇ x , x α := x α 1 1 • • • x αp p , |α| := α 1 + • • • + α p , α! := (α 1 !) • • • (α p !), |x| 2 = x 2 1 + • • • + x 2
p , and x := (1 + |x| 2 ) 1/2 . We choose the same notation for the length |α| of a multiindex α ∈ N p and for the euclidian norm |x| of a vector x ∈ R p but the context should avoid any confusion.

We denote by C

∞ (O; A) (resp. C ∞ b (O; A), resp. C ∞ c (O; A), resp. C ω (O; A)
) the vector space of functions from O to A which are smooth (resp. smooth with bounded derivatives, resp. smooth with compact support, resp. real analytic). In the case A = C, we simply write

C ∞ (O) (resp. C ∞ b (O), resp. C ∞ c (O), resp. C ω (O)
). Now we recall basic facts on analytic functions. Take again a Banach space (A, • A ), a positive integer p, and an open subset O of R p . For a smooth function θ : O -→ A, the following properties are equivalent:

θ ∈ C ω (O; A); ∃A > 0 , ∀ compact K ⊂ O , ∀α ∈ N p , sup x∈K (D α x θ)(x) A ≤ A |α|+1 • (α!) ; (2.1) ∃A > 0 , ∀ compact K ⊂ O , ∀α ∈ N p , sup x∈K (D α x θ)(x) A ≤ A |α|+1 • (|α|!) ; (2.2) ∃A > 0 , ∀ compact K ⊂ O , ∀α ∈ N p , sup x∈K (D α x θ)(x) A ≤ A |α|+1 • (1 + |α|) |α| . (2.3)
We refer to [START_REF] Hörmander | An introduction to complex analysis in several variables[END_REF][START_REF] Hearnshaw | Analyticity of the one-particle density matrix[END_REF][START_REF] Th | A new proof of the analyticity of the electronic density of molecules[END_REF] for details. Let (B, • B ) be another Banach space and f : A 2 -→ B be a continuous bilinear map.

Then, for any (

θ 1 ; θ 2 ) ∈ (C ω (O; A)) 2 , the map f (θ 1 ; θ 2 ) : O x → f (θ 1 (x); θ 2 (x)) ∈ B is real analytic.
Consider the case where A = W n,2 (O), for some integer n, and recall that •, • n (resp.

• n ) is the right linear scalar product (resp. the norm) on A. Take (

θ 1 ; θ 2 ) ∈ (C ω (O; A)) 2 . Then the map θ 1 , θ 2 n : O -→ C, defined by θ 1 , θ 2 n (x) = θ 1 (x), θ 2 (x) n , is real analytic. So is also the real valued map θ 1 2 n : O x → θ 1 (x) 2 n . In particular, taking n = 0, the maps θ 1 , θ 2 : O x → θ 1 (x), θ 2 (x) ∈ C and θ 1 2 : O x → θ 1 (x) 2 ∈ R + are also real analytic.
Finally we need some specific notation to describe the structure of the considered Hamiltonian and of our results. For a A-valued function θ : R 3 z = (z (1) ; z (2) ; z (3) ) → θ(z) ∈ A, we denote by dθ or d z θ the total derivative of θ. According to the previous notation, we denote, for j ∈ {1; 2; 3}, by ∂ j θ or ∂ z (j) θ the j'th first partial derivative. For a multiindex α = (α 1 ; α

2 ; α 3 ) ∈ N 3 , we set D α z := (-i∂ z ) α := (-i∂ 1 ) α 1 (-i∂ 2 ) α 2 (-i∂ 3 ) α 3 , D z = -i∇ z , z α = (z (1) ) α 1 (z (2) ) α 2 (z (3) ) α 3 , and |α| = α 1 + α 2 + α 3 . Let m be a positive integer. We write a point x ∈ R 3m as x = (x 1 ; • • • ; x m ). For any multiindex α = (α 1 ; • • • ; α m ) ∈ N 3m , we set D α x := (-i∂ x ) α := (-i∂ x 1 ) α 1 • • • (-i∂ xm ) αm , D x = -i∇ x , x α = x α 1 1 • • • x αm m , and |α| = |α 1 | + • • • + |α m |. Let m and p be positive integers. We write x ∈ R 3m (resp. y ∈ R 3p ) as x = (x 1 ; • • • ; x m ) (resp. y = (y 1 ; • • • ; y p )). For a A-valued function θ : R 3m × R 3p (x; y) → θ(x; y) ∈ A, let d x θ (resp. d y θ
) be the total derivative of θ w.r.t. x (resp. y). For n ∈ N, W n,2 (R 3p ) denotes the standard Sobolev space of L 2 -functions on R 3p admitting, up to order n, L 2 distributional derivatives. In particular, W 0,2 (R 3p ) = L 2 (R 3p ). We set W n = W n,2 (R 3p ) and let B n = L(W n ; W 0 ) be the Banach space of linear continuous maps from W n to W 0 .

Elliptic regularity.

In this subsection, we extend, essentially with the same proof, a regularity result that was proved but not stated in [J].

Let (m; p) ∈ (N * ) 2 . Let Ω be an open subset of R 3m . We consider second order differential operators P on Ω × R 3p given by

P = (α; β)∈N 3m ×N 3p |α|+|β|≤2 c αβ (x; y) D α x D β y , (2.4)
where the coefficients

c αβ belong to C ∞ b (Ω × R 3p ; C). In particular, for x ∈ Ω, the multi- plication operator C αβ (x; •) by c αβ (x; •) belongs to L(W n ), for all n ∈ N.
The principal symbol of those operators P is the smooth, complex-valued function σ P , that is defined on

Ω × R 3p × R 3m × R 3p by σ P (x; y; ξ; η) = (α; β)∈N m ×N p |α|+|β|=2 c αβ (x; y) ξ α η β .
(2.5)

We assume that P is (globally) elliptic in the following sense: there exists C P > 0 such that, for all (x; y; ξ

; η) ∈ Ω × R 3p × R 3m × R 3p , |ξ| 2 + |η| 2 ≥ 1 =⇒ σ P (x; y; ξ; η) ≥ C P |ξ| 2 + |η| 2 .
(2.6) Furthermore, we require that, for (α; β) ∈ N 3m × N 3p with |α| + |β| ≤ 2, the map Ω x → C αβ (x; •) ∈ B 0 is well-defined and real analytic. Under these assumptions on P , we note that, for all α ∈ N 3m with |α| ≤ 2, the function a α : Ω -→ B 2-|α| , which maps each x ∈ Ω to the differential operator

a α (x) = β∈N 3p |β|≤2-|α| c αβ (x; •) D β y (2.7)
on R 3p , is well-defined and real analytic.

Definition 2.1. For (m; p) ∈ (N * ) 2
and Ω an open subset of R 3m . We denote by Diff 2 (Ω) the set of all the differential operators P on Ω × R 3p satisfying the previous requirements.

Theorem 2.2. [J].

Let P be differential operator in the class Diff 2 (Ω). Let W : Ω -→ B 1 be a real analytic map. Take ϕ ∈ W 2,2 (Ω × R 3p ) such that (P + W )ϕ = 0. Then, the map Ω x → ϕ(x; •) ∈ W 2 is real analytic.
Remark 2.3. If, in Theorem 2.2, one replaces 3m by some positive integer n and one sets p = 0 and W = 0 with the convention that Ω × R 0 = Ω and W 2 = C, then one recovers exactly Hörmander's Theorem 7.5.1 in [START_REF] Hörmander | Linear partial differential operators[END_REF]. We assign Theorem 2.2 to [J] since, although it was not stated there, it was proved in a particular case and the corresponding proof directly extends to the present, general case. In [START_REF] Fournais | The electron density is smooth away from the nuclei[END_REF][START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF][START_REF] Hearnshaw | Analyticity of the one-particle density matrix[END_REF], one can find similar results in a slightly narrower framework. Therefore, Theorem 2.2 is more or less known in the literature.

Nevertheless, we provide below, for completeness, a proof of Theorem 2.2. To make it more transparent and facilitate the comparison with known results of this kind, we state intermediate results in some lemmata.

For convenience, we choose to use an appropriate pseudo-differential calculus. Since we are concerned with a local result (in Ω), it is natural to take a local one (see, for instance, in [START_REF] Hörmander | The analysis of linear partial differential operators III[END_REF] on page 83 or in [T], section 7.10). However, we prefer to use the basic global calculus presented in the beginning of Chapter 18 in [START_REF] Hörmander | The analysis of linear partial differential operators III[END_REF], since we think that it is more accessible.

Let P ∈ Diff 2 (Ω). To prove Theorem 2.2, it suffices to show the required real analyticity near each point in Ω. Let x 0 ∈ Ω and let Ω 0 be a bounded neighbourhood of x 0 such that Ω 0 ⊂ Ω. As a first step, we construct an extension P ∈ Diff 2 (R 3m ) of the restriction of P to Ω 0 . In view of (2.6), we observe that the real part σ P of σ P does not vanish, except for (ξ; η) = (0; 0). Since it is continuous, it must keep the same sign σ 0 on

Ω × R 3p × (R 3(m+p) \ {(0; 0)}) .
Furthermore, (2.6) also implies that, for (x; y; ξ

; η) ∈ Ω × R 3p × R 3(m+p) , σ 0 • σ P (x; y; ξ; η) ≥ C P |ξ| 2 + |η| 2 . (2.8) Let χ ∈ C ∞ c (R 3m
) be a cut-off function with support in Ω such that 0 ≤ χ ≤ 1 and, on Ω 0 , χ = 1. Recalling that P is given by (2.4), we set

P = (α; β)∈N 3m ×N 3p |α|+|β|≤2 ĉαβ (x; y) D α x D β y ,
where ĉαβ (x; y) = χ(x)c αβ (x; y), for all (x; y) ∈ R 3(m+p) and |α| + |β| < 2, and where, for all (x; y; ξ; η) ∈ R 6(m+p) , σ P (x; y; ξ; η) :=

(α; β)∈N m ×N p |α|+|β|=2 ĉαβ (x; y) ξ α η β := χ(x) • σ P (x; y; ξ; η) + σ 0 C P • 1 -χ(x) • |ξ| 2 + |η| 2 .
(2.9)

We observe that (2.8) with σ P replaced by σ P holds true on R 6(m+p) . In particular, P is elliptic. We immediately verify that P ∈ Diff 2 (R 3m ). We note that P coincides with P as differential operator on Ω 0 . We observe that the symbol of P belongs to the global class S 2 (R 6N ) of symbols, for N = m + p, and, since it is elliptic, we can construct a parametrix (cf. Appendix B). Thus there exist two pseudo-differential operators Q = q(x; y; D x ; D y ) and R = r(x; y; D x ; D y ) with symbols q, r ∈ S -2 (R 6N ) such that Q P = I -R (where I is the identity differential operator) and

Q , R ∈ L W k,2 (R 3N ); W k+2,2 (R 3N ) .
(2.10)

Proof of Theorem 2.2: We take ϕ ∈ W 2,2 (Ω×R 3p ) such that (P +W )ϕ = 0. Let x 0 ∈ Ω. We use the above construction. Let Ω 0 a neighbourhood of

x 0 such that Ω 0 ⊂ Ω 0 . Let χ 0 ∈ C ∞ c (R 3
) with χ 0 = 1 near Ω 0 and with support in Ω 0 . Thus χχ 0 = χ 0 , where χ is the cut-off function appearing in the definition of P .

Lemma 2.4. We have

χ 0 ϕ ∈ C ∞ c (R 3m ; W 2 ).
Proof: Since Q P = I -R, P χ 0 = χP χ 0 , and (P + W )ϕ = 0, we obtain

χ 0 ϕ = Rχ 0 ϕ + Qχ[P , χ 0 ]ϕ -Qχ 0 W ϕ , (2.11)
where the commutator [P , χ 0 ] is a first order diffential operator in x. Since

ϕ ∈ W 2,2 (Ω × R 3p ), ϕ ∈ W 1,2 (Ω; W 1 ) and χ 0 ϕ ∈ W 1,2 (R 3m ; W 1
). We successively observe that

W ϕ ∈ W 1,2 (Ω; W 0 ) , χ 0 W ϕ ∈ W 1,2 (R 3m ; W 0 ) , and QχW χ 0 ϕ ∈ W 2,2 (R 3m ; W 1 ) .
Similarly, we have

[P , χ 0 ]ϕ ∈ W 0,2 (Ω; W 1 ) , χ[P , χ 0 ]ϕ ∈ W 0,2 (R 3m ; W 1 ) , Qχ[P , χ 0 ]ϕ ∈ W 2,2 (R 3m ; W 1 ) and Rχ 0 ϕ ∈ W 3,2 (R 3m ; W 1 ) .
Thus χ 0 ϕ ∈ W 2,2 (R 3m ; W 1 ), by (2.11). Now, we apply the same argument again but with this new information, yielding χ 0 ϕ ∈ W 2,2 (R 3m ; W 2 ). In this way, we prove by induction that, for all k ∈ N,

χ 0 ϕ ∈ W k,2 (R 3m ; W 2 ), yielding χ 0 ϕ ∈ C ∞ (R 3m ; W 2 ).
Lemma 2.5. There exists

C e > 0 such that, for all v ∈ C ∞ c (R 3m ; W 2 ) with support in Ω 0 , for all r ∈ {0; 1; 2}, for all α ∈ N 3m , |α| + r ≤ 2 =⇒ D α x v L 2 (Ω 0 ;Wr) ≤ C e (P + W )v L 2 (Ω 0 ;W 0 ) .
(2.12)

Proof: Let P 2 := σ P (x; y; D x ; D y ) and P 1 = P -P 2 . We consider the following operators as operators in L 2 (Ω 0 ; W 0 ) with domain W 2,2 (Ω 0 × R 3p ). By ellipticity, ∆ = ∆ x + ∆ y is P 2 -bounded. By assumption, W is (-∆ x ) 1/2 -bounded. Since the operators (-∆ x ) 1/2 and D α x D β y , for |α| + |β| ≤ 1, are ∆-bounded with relative bound less than one, so are also W and, thanks to the boundedness assumption on its coefficients, P 1 . Thus W and P 1 are P 2 -bounded with relative bound less than one. This implies that W is P -bounded with relative bound less than one. For |α|

+ |β| ≤ 2, D α x D β y is ∆-bounded.
The above properties imply that it is also (P + W )-bounded, yielding (2.12).

For ≥ 0, let Ω := {x ∈ Ω 0 ; d(x; R 3m \ Ω 0 ) > } and, for r ∈ N, denote the L 2 (Ω ; W r )- norm of v by N ;r (v). Note that χ 0 ϕ ∈ C ∞ c (R 3m ; W 2 ) with support in Ω 0 , by Lemma 2.4.
Lemma 2.6. There exists B > 0 such that, for all > 0, j ∈ N, r ∈ {0; 1; 2}, and α ∈ N 3m ,

r + |α| < 2 + j =⇒ r+|α| N j ;r (D α x ϕ) ≤ B r+|α|+1 .
(2.13)

Proof: Let D 0 be the diameter of Ω 0 . Since P ∈ Diff 2 (Ω), there exists C a > 0 such that, for all α ∈ N 3 , 0

≤ 1 ≤ D 0 , |α| 1 |β|≤2 sup x∈Ω 1 ∂ α x a β B 2-|β| ≤ C |α|+1 a • (|α|!) .
(2.14) by (2.2). Using that (P + W )ϕ = 0, it suffices to follow the proof of (7.5.9) in [START_REF] Hörmander | Linear partial differential operators[END_REF], p. 179-180, to get the desired result.

To complete the proof of Theorem 2.2, it suffices, thanks to (2.13), to follow the argument on page 180 in [START_REF] Hörmander | Linear partial differential operators[END_REF] and use (2.2).

Remark 2.7. We point out that one can prove Theorem 2.2 without pseudodifferential technics. We used them in the proof of Lemma 2.4. In fact, it turns out that the arguments used in [START_REF] Fournais | Analyticity of the density of electronic wave functions[END_REF] (see the proof of Lemma 3.1) and in [HS] (see Section 3.3) can be adapted to get a different proof of those lemma.

3 Twisted differential calculus.

In this section, we recall the notion of "twist" used in [Hu, KMSW, MS] to "desingularize" a class of Hamiltonians with Coulomb-like singularities.

3.1 Framework.

We need some more notation. Denoting by S 2 the unit euclidian sphere of R 3 , let ṽ0 : S 2 -→ C and η :]0; +∞[-→ R + be two functions such that, for some η 0 ∈]0; 1[ and

B 0 > 0, (1/B 0 ) ≤ |ṽ 0 | ≤ B 0 and ∀t > 0 , sup t(1-η 0 ) ≤ s ≤ t(1+η 0 ) η(s) ≤ B 0 • η(t) . (3.1)
We assume further that the multiplication operator by the function

η(| • |) : R 3 z → η(|z|) ∈ R + belongs to L(W 1,2 (R 3 ); L 2 (R 3 )).
The multiplication operator by the function v 0 : R 3 \{0} -→ C, defined by v 0 (z) = η(|z|)ṽ 0 (z/|z|), also belongs to

L(W 1,2 (R 3 ); L 2 (R 3 )). Definition 3.1. Let V be the class of functions v : R 3 -→ C such that v ∈ C ω (R 3 \ {0})
and there exists C > 0 such that

∀α ∈ N 3 , ∀z ∈ R 3 \ {0} , |z| |α| • D α v(z) ≤ C 1+|α| • |α|! • |v 0 (z)| . (3.2)
We observe that, if v ∈ V, then the multiplication operator by v belongs to the space

L(W 1,2 (R 3 ); L 2 (R 3 )).
We consider a system of m + p electrons, interacting to one another and moving under the influence of L fixed nuclei. Let I e := {1; • • • ; m} be the set of indices for the "external" electronic variables x ∈ R 3m , let I i := {1; • • • ; p} be the one of "internal" electronic variables y ∈ R 3p , and let I n := {1; • • • ; L} be the one of the nuclear variables R ∈ R 3 . It is convenient to introduce the disjoint union of those sets, namely

I e I i I n := (a; j) ; a ∈ {e; i; n} , j ∈ I a .
For c = (a; j) ∈ I e I i I n , let π 1 (c) = a and π 2 (c) = j. We define

X c = x j if π 1 (c) = e, X c = y j if π 1 (c) = i, and X c = R j if π 1 (c) = n. Let P := {c; c } ; (c; c ) ∈ I e I i I n 2 and c = c
be the set of all possible particle pairings. Let Ω be an open subset of R 3m . We consider potentials V where

∀(x; y) ∈ Ω × R 3p , V (x; y) = {c; c }∈P V {c; c } (X c -X c ) , (3.3)
with V {c; c } ∈ V (cf. Definition 3.1), for all {c; c } ∈ P.

Recall that W n = W n,2 (R 3p ), for n ∈ N, and B 1 = L(W 1 ; W 0 ). We note that, if v ∈ V and {c; c } ∈ P, then, for any fixed x, the multiplication by the almost everywhere defined function v(X c -X c ) of y belongs to B 1 , thanks to (3.2) and the assumption on v 0 .

Recall that we introduced in Definition 2.1 a class of differential operators. The class of Hamiltonians, that we want to treat, is defined in the following:

Definition 3.2. For (m; p) ∈ (N * ) 2
and Ω an open subset of R 3m . We denote by Ham 2 (Ω) the set of all the Hamiltonians H on Ω × R 3p such that H = P + V , where P ∈ Diff 2 (Ω) (cf. Definition 2.1) and V satisfies (3.3).

Heuristic.

In this subsection, we want to motivate the notion of "twist", that is going to play a crucial rôle in the proof of the main results. To explain the basic idea, let us consider the following situation.

Take the Hamiltonian H in (1.1) with N = 2. We are not able to apply Theorem 2.2 with x = x 1 and y = x 2 to the equation Hψ = Eψ, since the map

x 1 → (x 2 → |x 1 -x 2 | -1 ) ∈ B 1 = L(W 1,2 (R 3 ); L 2 (R 3
)) is admittedly well defined by (1.2) but is not analytic. The problem lies in the dependency of the Coulomb singularity on x 1 . Assume that, for some open subset Ω of R 3 and some x 0 ∈ R 3 , we have a smooth map

f : Ω × R 3 → R 3 such that, for all x 1 ∈ Ω, f (x 1 ; •) is a smooth diffeomorphism of R 3 satisfying f (x 1 ; x 0 ) = x 1 . For such x 1 , let U x 1 be the unitary implementation of f (x 1 ; •) on L 2 (R 3 ) (the "twist"). For any function ϕ ∈ W 1,2 (R 3 ) and for x 1 ∈ Ω, U x 1 |x 1 -•| -1 ϕ (y) = f (x 1 ; x 0 ) -f (x 1 ; y) -1 • U x 1 ϕ (y) .
Now, since f (x 1 ; •) is a diffeomorphism, the singularity occurs when y = x 0 . Therefore, it does not depends on x 1 anymore. We can write

U x 1 |x 1 -•| -1 U -1 x 1 ϕ (y) = |x 0 -y| • f (x 1 ; x 0 ) -f (x 1 ; y) -1 • |x 0 -y| -1 • ϕ(y) .
Assuming further that the map

Ω x 1 → |x 0 -•| • f (x 1 ; x 0 ) -f (x 1 ; •) -1 ∈ L ∞ (R 3 ) (3.4) is smooth, we see that Ω x 1 → U x 1 |x 1 -•| -1 U -1 x 1 ∈ B 1 = L(W 1,2 (R 3 ); L 2 (R 3 )) is also smooth, thanks to (1.2).
What can we do with the other singular terms in H? Since the terms |x 1 -R | -1 do not depend on x 2 , the previous conjugation has no effect on them. But, if x 1 stays far away from all R , they are not singular. We are left with the terms |x 2 -R | -1 , which are a constant function of x 1 with values in B 1 by (1.2). Since we want to perform the above conjugation to treat the term |x 1 -x 2 | -1 , we need to verify that this conjugation does not distroy the nice property of these terms. If we require that, for all and all x 1 ∈ Ω, f (x 1 ; R ) = R , and the smoothness of the map (3.4) with x 0 replaced by R , the above computation shows that the map Ω

x 1 → U x 1 |x 2 -R | -1 U -1
x 1 ∈ B 1 is also smooth. With some more work, we can prove real analyticity of the previous terms. We just have seen that the "twisted" potential

U x 1 V U -1
x 1 is a real analytic, B 1 -valued function. What about the "twisted" kinetic energy? A simple computation shows that the operators ∇ x 1 and ∇ x 2 are transformed by the twist into differential operators in the (x 1 ; x 2 ) variables with explicit coefficients in terms of the derivatives of f (cf. (3.16) and (3.18)). As we shall see below, the ellipticity of H is preserved by the "twist" and the coefficients of the "twisted" operator have some property of real analyticity. Coming back to the original problem, we do not apply Theorem 2.2 with x = x 1 and y = x 2 to the equation Hψ = Eψ but to the "twisted" equation (U x 1 HU -1

x 1 )U x 1 ψ = EU x 1 ψ. Therefore, we obtain some analytic regularity of the "twisted" bound state, namely U x 1 ψ, not on ψ itself. But, if we are interested in the regularity of ρ 1 , we are done. Indeed, since we can write ρ 1 (x 1 ) = ψ(x 1 ; •) 2 ( • being the L 2 -norm on R 3 ) and since U x 1 is unitary on L 2 (R 3 ), ρ 1 (x 1 ) = U x 1 ψ(x 1 ; •) 2 and ρ 1 is real analytic on Ω by composition. We just have sketched, in the N = 2 case, the proof of the smoothness of ρ 1 , which was performed in [J], provided there exists a map f satisfying the above requirements. Actually it does (see Subsection 3.3 below).

Hunziker's twist.

With the above strategy in mind, we recall Hunziker's twist (cf. [Hu, KMSW]) and provide basic properties of it.

We take a real-valued cut-off function τ ∈ C ∞ c (R 3 ) such that τ (x) = 0 if |x| ≥ 1, and τ (0) = 1. Let U be the open subset of R 3m defined by

U = x ∈ R 3m ; ∀(c; c ) ∈ P; π 1 (c) = e and π 1 (c ) ∈ {e; n} , X c = X c = x ∈ R 3m ; ∀(j; j ) ∈ [[1; m]] 2 , ∀ ∈ [[1; L]] , x j = R (3.5)
and j = j =⇒ x j = x j .

Let x 0 = (x 0 1 ; • • • ; x 0 m ) ∈ U.
We can find some r 0 > 0 such that, for (j;

j ) ∈ [[1; m]] 2 , for ∈ [[1; L]], |x 0 j -R | ≥ r 0 and j = j =⇒ |x 0 j -x 0 j | ≥ r 0 . (3.6)
Let f : R 3m × R 3 -→ R 3 be the smooth function defined by

f (x; z) = z + m j=1 τ r -1 0 (z -x 0 j ) • (x j -x 0 j ) . (3.7)
Notice that, by (3.6), we have, for all x ∈ R 3m ,

∀j ∈ [[1; m]] , ∈ [[1; L]] , f (x; x 0 j ) = x j and f (x; R ) = R . (3.8) For all (x; z) ∈ R 3m × R 3 , the total differential of f w.r.t.
x is given by

(d x f )(x; z) = m j=1 τ r -1 0 (z -x 0 j ) dx j (3.9)
where dx j : R 3m

x → x j ∈ R 3 , for 1 ≤ j ≤ m. In particular, the map

R 3m x → (d x f )(x; •) ∈ C ∞ b (R 3 ; L(R 3m ; R 3
)) is well-defined and constant. Using (3.9) and setting

C(τ ) = sup z∈R 3 dτ (z) L(R 3 ) < +∞ ,
we can write, for all j ∈ [[1; m]], for all (z; z ) ∈ (R 3 ) 2 , and for all x ∈ R 3m , (d

x j f )(x; z) -(d x j f )(x; z ) L(R 3m ;R 3 ) = τ r -1 0 (z -x 0 j ) -τ r -1 0 (z -x 0 j ) ≤ r -1 0 C(τ ) |z -z | .
(3.10) For all (x; z) ∈ R 3m × R 3 , the total differential of f w.r.t. z is given by, for z ∈ R 3 ,

(d z f )(x; z) • z = z + r -1 0 m j=1 (dτ ) r -1 0 (z -x 0 j ) • z (x j -x 0 j ) . (3.11) For δ > 0, set Ω(δ) := B(x 0 1 ; δ[× • • • × B(x 0 m ; δ[.
Recall that the class V (cf. Definition 3.1) depends on a parameter η 0 (see (3.1)). By (3.11), we can find some δ 0 ∈]0; r 0 /2], depending on C(τ ) and r 0 , such that ∀x ∈ Ω(δ 0 ) , sup

z∈R 3 (d z f )(x; z) -I 3 L(R 3 ) ≤ min η 0 ; 1 -η 0 < 1 .
(3.12)

We note that the requirement δ 0 ∈]0; r 0 /2] yields Ω(δ 0 ) ⊂ U, by (3.6) and the triangle inequality. Thanks to (3.12), for each x ∈ Ω(δ 0 ),

f (x; •) is a C ∞ -diffeomorphism of R 3 .
We denote by f -1 (x; •) its inverse. Furthermore, for all x ∈ Ω(δ 0 ), for all (z; z )

∈ (R 3 ) 2 , (1 -η 0 ) |z -z | ≤ f (x; z) -f (x; z ) ≤ (1 + η 0 ) |z -z | . (3.13)
By a direct computation using (3.9) and (3.11), one can express the derivative d z f -1 and d x f -1 in terms of d z f (see Appendix A for details). In particular, one can check

that f ∈ C ∞ b (Ω(δ 0 ) × R 3 ; R 3 ) and f -1 ∈ C ∞ b (Ω(δ 0 ) × R 3 ; R 3 ). Consider the map F : Ω(δ 0 ) × R 3p -→ R 3p , defined by F (x; y) = (f (x; y 1 ); • • • ; f (x; y p )). For each x ∈ Ω(δ 0 ), F (x; •) is a C ∞ -diffeomorphism of R 3p
and we denote its inverse map by F -1 (x; •). Explicit expressions for d x F , d y F , d x F -1 , and d y F -1 , are provided in Appendix A. In particular, we observe that

F ∈ C ∞ b (Ω(δ 0 ) × R 3p ; R 3p ) and F -1 ∈ C ∞ b (Ω(δ 0 ) × R 3p ; R 3p ).
We introduce the map U : Ω(δ 0 ) x → U x ∈ B 0 (with B 0 = L(W 0 ) and W 0 = L 2 (R 3p )), with values in the set of unitary operators on W 0 , defined by, for θ ∈ W 0 ,

U x θ (y) = Det(d y F )(x; y) 1/2 • θ F (x; y) . (3.14)
For x ∈ Ω(δ 0 ), the map U x is the unitary implementation of F (x; •) on L 2 (R 3p ) and

U -1 x θ (y) = Det d y F -1 (x; y) 1/2 • θ F -1 (x; y) . (3.15)
Denoting by A T the transposed of a linear map A and by DetA its determinant, we see by a direct computation that, as differential operators on Ω(δ 0 ) × R 3p ,

U -1 x D x U x = D x + J 1 (F ) D y + J 2 (F ) , (3.16) U x D x U -1 x = D x + J 1 F -1 D y + J 2 F -1 , (3.17) U -1 x D y U x = J 3 (F ) D y + J 4 (F ) , (3.18)
U x D y U -1 x = J 3 F -1 D y + J 4 F -1 , (3.19)
where, for x ∈ Ω(δ 0 ) and y ∈ R 3p , with G = F and

G -1 = F -1 or with G = F -1 and G -1 = F , J 1 (G)(x; y) := d x G (x; y ) T x; y = G -1 (x; y) , J 2 (G)(x; y) := -i Det d y G -1 (x; y) 1/2 D x Det d y G (x; y ) 1/2 y =G -1 (x;y) , J 3 (G)(x; y) := d y G (x; y ) T x; y = G -1 (x; y) , J 4 (G)(x; y) := -i Det d y G -1 (x; y) 1/2 D y Det d y G (x; y ) 1/2 y =G -1 (x;y)
.

Thanks to the computation written in Appendix A, we note that, on one hand, for G ∈ {F ; F -1 },

J 1 (G) ∈ C ∞ b Ω(δ 0 ) × R 3p ; L R 3p ; R 3m , iJ 2 (G) ∈ C ∞ b Ω(δ 0 ) × R 3p ; R 3m , (3.20) J 3 (G) ∈ C ∞ b Ω(δ 0 ) × R 3p ; L R 3p , iJ 4 (G) ∈ C ∞ b Ω(δ 0 ) × R 3p ; R 3p , (3.21)
and, on the other hand, if we denote by J the multiplication operator by any component of J k (G), k ∈ {1; 2; 3; 4} and G ∈ {F ; F -1 },

J ∈ C ω Ω(δ 0 ); B 0 . (3.22)
Furthermore (cf. Appendix A), there exists C > 0 such that, for G ∈ {F ; F -1 }, for all (x; y) ∈ Ω(δ 0 ) × R 3p , J 3 (G)(x; y) is invertible and

C -1 ≥ J 3 (G)(x; y) L(R 3p ) + J 3 (G)(x; y) -1 L(R 3p ) ≥ C . (3.23) Now, using (3.18), D y U x = U x J 3 (F )D y -J 4 (F )
thus, for all x ∈ Ω(δ 0 ), we see, by induction and by (3.21), that U x leaves W n invariant, for all n ∈ N. Similarly, we can show, using (3.19) and (3.21), that, for all x ∈ Ω(δ 0 ), U -1

x also leaves W n invariant, for all n ∈ N. By the same arguments, we can check, starting from (3.16) and (3.18), that U preserves the space W 2,2 (Ω(δ 0 ) × R 3p ).

Potential regularization by a twist.

We first consider the action by conjugation of the twist U on the potential part V of Hamiltonians H in the class Ham 2 (Ω).

Proposition 3.3. Take x 0 ∈ U ⊂ R 3m . We choose r 0 > 0 and δ 0 ∈]0; r 0 /2] such that (3.6) and (3.12) hold true. Recall that the neighbourhood Ω(δ

0 ) := B(x 0 1 ; δ 0 [× • • • × B(x 0 m ; δ 0 [ of x 0 is included in U.
Then, for any v ∈ V and any {c; c } ∈ P, the map

Ω(δ 0 ) x → U x v(X c -X c )U -1
x ∈ B 1 is well-defined and real analytic.

Remark 3.4. We point out that the Coulomb potential

| • | -1 does belong to V. Indeed, (3.2), for v = v 0 = | • | -1 ,
was proved in [J] and the multiplication by this v 0 belongs to L(W 1,2 (R 3 ); L 2 (R 3 )) by Hardy's inequality (1.2).

Remark 3.5. A proof of Proposition 3.3 is essentially available in [Hu, J]. However, the proof in [Hu] is restricted to Coulomb interactions. We provide below a more transparent version of the proof in [J].

Proof of Proposition 3.3: Let x ∈ Ω(δ 0 ). Recall that, for k ∈ {0; 1}, U x and U -1 x leave W k invariant. Recall that, for all {c; c } ∈ P, the multiplication by v(X c -X c ) belongs to B 1 . Therefore the map w : Ω(δ

0 ) x → U x v(X c -X c )U -1
x ∈ B 1 is well-defined. We treat all possible cases for {c; c } ∈ P.

If π 1 (c) = π 1 (c ) = n, then w : Ω(δ 0 ) x → v(X c -X c ) is a constant. Assume that π 1 (c) = e and π 1 (c ) ∈ {e; n} or that π 1 (c ) = e and π 1 (c) ∈ {e; n}. Since Ω(δ 0 ) ⊂ U, X c -X c does not vanish on Ω(δ 0 ). Thus w : Ω(δ 0 ) x → v(X c -X c ) is real analytic, by composition. Assume that π 1 (c) = e and π 1 (c ) = i. Let j = π 2 (c) and k = π 2 (c ). Then, for x ∈ Ω(δ 0 ), U x v(X c -X c )U -1
x is the multiplication operator on W 1 by the function

y → v x j -f (x; y k ) = v f (x; x 0 j ) -f (x; y k ) ,
by (3.8). Let j ∈ I e . Using the fact that d x f does not depend on x (see (3.9)), we have, for α ∈ N 3 and y k = x 0 j ,

∂ α x j v f (x; x 0 j ) -f (x; y k ) (3.24) = (∂ α v) f (x; x 0 j ) -f (x; y k ) • τ r -1 0 (x 0 j -x 0 j ) -τ r -1 0 (y k -x 0 j ) |α| .
Using successively (3.10), (3.13), and (3.2), we obtain

∂ α x j v f (x; x 0 j ) -f (x; y k ) ≤ (∂ α v) f (x; x 0 j ) -f (x; y k ) • r -1 0 C(τ ) |x 0 j -y k | |α| , ≤ (∂ α v) f (x; x 0 j ) -f (x; y k ) • r 0 (1 -η 0 ) -1 C(τ ) f (x; x 0 j ) -f (x; y k ) |α| , ≤ C 1+|α| 1 • |α|! • v 0 f (x; x 0 j ) -f (x; y k ) ,
for some C 1 > 0. Now, by the structure of v 0 and by (3.13) and (3.1), we derive that, for some positive constants C 2 and C 3 ,

∂ α x j v f (x; x 0 j ) -f (x; y k ) ≤ C 1+|α| 2 • |α|! • η f (x; x 0 j ) -f (x; y k ) , ≤ C 1+|α| 3 • |α|! • η |x 0 j -y k | . (3.25)
Now, one can check, by induction on |α| and by the fact that the multiplication by η(| • |) belongs to B 1 , that, on Ω(δ 0 ), the ∂ α x j partial derivative of w exists and maps each x ∈ Ω(δ 0 ) to the multiplication operator on W 1 by the almost everywhere defined function (3.24) of y. Thus, w ∈ C ∞ (Ω(δ 0 ); B 1 ). Using (3.25) again, we obtain that sup

x∈Ω(δ 0 ) (D α x j w)(x) B 1 ≤ C 1+|α| 4
• |α|! , (3.26)

for some C 4 > 0. By (2.2), w is real analytic w.r.t. x j . Since this holds true for all j , w is real analytic. Similarly, we can treat the case π 1 (c) = i and π 1 (c ) = e.

Assume that π 1 (c) = n and π 1 (c ) = i. Let = π 2 (c) and k = π 2 (c ). For x ∈ Ω(δ 0 ), w(x) is the multiplication operator in B 1 by the function

y → v R -f (x; y k ) = v f (x; R ) -f (x; y k ) , (3.27)
by (3.8). We redo the computation (3.24) -(3.25), with x 0 j replaced by R . Furthermore, we also can check that the ∂ α

x j partial derivative of w exists and maps each x ∈ Ω(δ 0 ) to the multiplication operator on W 1 by the almost everywhere defined ∂ α

x j partial derivative of (3.27). Using the new estimate (3.25), we derive a bound like (3.26). As above, we conclude that w is real analytic. Again, we have a similar treatment of the case π 1 (c) = i and π 1 (c ) = n. We are left with the case π 1 (c) = i and π 1 (c ) = i. We set k = π 2 (c) and k = π 2 (c ). By assumption, k = k . Using (3.8) again, we see that, for x ∈ Ω(δ 0 ), w(x) is the multiplication operator in B 1 by the function y → v f (x; y k ) -f (x; y k ) .

(3.28)

Once again, we perform the computation (3.24) -(3.25) and follow the above arguments, using this time that the multiplication by y → η(|y k -y k |) belongs to B 1 . This yields the real analyticity of w.

Twisted regular differential calculus.

Next, we let act the twist U by conjugation on the regular, differential part P ∈ Diff 2 (Ω) of Hamiltonians H in the class Ham 2 (Ω) (cf. Definitions 2.1 and 3.2). It turns out that this action preserves this class Diff 2 (Ω).

Proposition 3.6. Let x 0 ∈ U ∩ Ω ⊂ R 3m . Take r 0 > 0 and δ 0 ∈]0; r 0 /2] such that (3.6) and (3.12) hold true, and such that Ω(δ

0 ) := B(x 0 1 ; δ 0 [× • • • × B(x 0 m ; δ 0 [⊂ Ω. We recall that Ω(δ 0 ) ⊂ U. Let us define U * : Ω(δ 0 ) x → U -1
x ∈ B 0 . Let P ∈ Diff 2 (Ω) (see Definition 2.1).Then the operator composition P = U P U * defines a differential operator on Ω(δ 0 ) × R 3p that belongs to Diff 2 0 )).

Remark 3.7. In [Hu, J], it was proven that the twisted Laplacian U ∆U * is elliptic, when m = 1. Proposition 3.6 shows that the conjugation by the twist of an elliptic operator is also elliptic. This property is actually used in [KMSW, MS] in a semiclassical framework.

Remark 3.8. In the proof of Proposition 3.6 below, we essentially follow [Hu, J]. Observing that U is a Fourier integral operator (cf. [START_REF] Grigis | Microlocal analysis for differential operators. An introduction[END_REF][START_REF] Hörmander | The analysis of linear partial differential operators IV[END_REF]), one can use composition rules for such operators to derive the relationship between the principal symbols of P and P , given in (3.30) below.

Proof of Proposition 3.6: Take P ∈ Diff 2 (Ω), that is of the form (2.4). We first observe that, for (α; β) ∈ N 3m × N 3p with |α| + |β| ≤ 2, U c αβ U * is, as differential operator, the multiplication operator by the function Ω(δ 0 ) × R 3p (x; y) → c αβ x; F (x; y) , which belongs to C ∞ b (Ω(δ 0 ) × R 3p ; C). Using (3.17) and (3.19) componentwise and taking advantage of (3.20) and (3.21) for G = F -1 , we see that P = U P U * has the form (2.4) on Ω(δ 0 ) × R 3p with coefficients in C ∞ b (Ω(δ 0 ) × R 3p ; C). We are left with the proof of the ellipticity of P . To this end, we compute its principal symbol σ P . Let us denote by "•" the scalar product in R d , for any d ∈ N * . It is well-known that the total symbol of P is the map

S P : Ω(δ 0 ) × R 3p × R 3m × R 3p -→ C defined by S P (x; y; ξ; η) = e -ix•ξ-iy•η P e ix•ξ+iy•η . (3.29)
Using this, we obtain in Appendix A that, for (x; y; ξ

; η) ∈ Ω(δ 0 ) × R 3p × R 3m × R 3p , σ P (x; y; ξ; η) = σ P x; F (x; y); ξ + J 1 F -1 (x; y) η; J 3 F -1 (x; y) η .
(3.30)

Since P is elliptic, there exists C > 0 such that (2.6) holds true with P replaced by P and C replaced by C. Take (ξ; η) satisfying |ξ| 2 + |η| 2 ≥ 1. Then

C -1 σ P (x; y; ξ; η) ≥ ξ + J 1 F -1 (x; y)η 2 + J 3 F -1 (x; y)η 2 .
(3.31) By (3.20), we can find M > 0 such that, for all (x; y) ∈ Ω(δ 0 ) × R 3p ,

J 1 F -1 (x; y) L(R 3p ;R 3m ) ≤ M . Now let S = √ 1 + 4M 2 . Consider first the case where S|η| ≤ (|ξ| 2 + |η| 2 ) 1/2 . We have 2M |η| ≤ |ξ|. Thus ξ + J 1 F -1 (x; y)η 2 ≥
|ξ| 2 4 and, using (3.23), we obtain from (3.31) the lower bound

C -1 σ P (x; y; ξ; η) ≥ min(1/4; C -2 )(|ξ| 2 + |η| 2 ) . If, now, S|η| ≥ (|ξ| 2 + |η| 2 ) 1/2 , it follows from (3.23) and (3.31) that C -1 σ P (x; y; ξ; η) ≥ C -2 S -2 (|ξ| 2 + |η| 2 ) .
This yields (2.6). Thus P ∈ Diff 2 (Ω(δ 0 )).

Elliptic regularity for twisted Hamiltonians.

Now, we conjugate by the twist U any Hamiltonians H of the class Ham 2 (Ω) (cf. Definition 3.2) and show that we can apply Theorem 2.2 to the twisted operator.

Theorem 3.9. Let Ω be an open subset of R 3m and x 0 ∈ U ∩ Ω. Then, one can find a neighbourhood Ω 0 of x 0 and a unitary operators valued map U : Ω 0 x → U x ∈ B 0 such that the following holds true: For any Hamiltonian H ∈ Ham 2 (Ω) on Ω × R 3p ⊂ R 3(m+p) , for any Ψ ∈ W 2,2 (Ω × R 3p ) such that HΨ = 0, the map U Ψ : Ω 0

x → U x Ψ(x; •) ∈ W 2 is well-defined and real analytic. Furthermore, for any α ∈ N 3m with 1 ≤ |α| ≤ 2, the map U (D α

x Ψ) : Ω 0 x → U x (D α

x Ψ)(x; •) ∈ W 2-|α| is well-defined and real analytic.

Proof: Given x 0 ∈ U ∩ Ω, take Ω 0 = Ω(δ 0 ) satisfying the assumptions of Proposition 3.6 and consider the twist U : Ω 0 → B 0 defined by (3.14). Recall that

U * : Ω 0 x → U -1 x ∈ B 0 . Take H = P + V ∈ Ham 2 (Ω) and Ψ ∈ W 2,2 (Ω × R 3p
) such that HΨ = 0. Applying Proposition 3.3 to each term of V and Proposition 3.6 to P , we see that P = U P U * and W = U V U * satisfy the assumptions of Theorem 2.2 on Ω 0 , yielding the real analyticity of the map U Ψ : Ω 0 x → U x Ψ(x; •) ∈ W 2 . By (3.17), we have, for x ∈ Ω 0 , setting J 1 := J 1 (F -1 ) and J 2 := J 2 (F -1 ),

U x D x Ψ (x; •) = U x D x U -1 x • U x Ψ(x; •) = D x U x Ψ(x; •) + J 1 (x; •) • D y U x Ψ(x; •) + J 2 (x; •) U x Ψ(x; •) . By (3.20), this implies that, for α ∈ N 3m with |α| = 1, U (D α x Ψ) : Ω 0 x → U x (D α x Ψ)(x; •) ∈ W 1 is well-defined and real analytic. Letting act U x D x U -1
x on the above equality and using (3.17), (3.20), and (3.21), we obtain that, for α ∈ N 3m with |α| = 2, U (D α

x Ψ) : Ω 0 x → U x (D α

x Ψ)(x; •) ∈ W 0 is well-defined and real analytic.

Remark 3.10. We observe that the present framework (the one in Theorem 3.9) is closer to the one in [KMSW] than the one in Hunziker's original paper [Hu]. This comes from the fact that, here and in [KMSW], the considered Hamiltonian contains derivatives w.r.t. the "external" variable x, whereas it is not the case in [Hu]. We further point out that a variant of Hunziker's twist is used in a similar way in [MM].

4 Proof of the main results.

In this section, we provide a proof for our main results, namely Theorem 1.3, Theorem 1.4, and Theorem 1.5. We first prove the result on the current density C.

Proof of Theorem 1.3: , by (3.5). Note that H := H -E ∈ Ham(R 3 ) and Hψ = 0. Applying Theorem 3.9 to Ω = R 3m = R 3 and x 0 = x 0 ∈ U ∩ Ω, there exists some neighbourhood Ω 0 of x 0 and a map U : Ω 0 x → U x ∈ B 0 , the values of which are unitary operators on W 0 , such that, for all α ∈ N 3 with |α| ≤ 1, the map

Let x 0 ∈ R 3 \ {R 1 , • • • , R L }. Notice that, if we set m = 1 and p = N -1, then U = R 3 \ {R 1 , • • • , R L }
U (D α x ψ) : Ω 0 x → U x (D α x ψ)(x; •) ∈ W 2-|α|
is well-defined and real analytic. Let α ∈ N 3 with |α| = 1. For x ∈ Ω 0 , we can write, since U x is unitary,

D α x ψ (x; •) , ψ(x; •) = U x D α x ψ (x; •) , U x ψ(x; •) .
As the scalar product of real analytic, W 0 -valued maps, the map Ω 0 x → (D α x ψ), ψ is also real analytic. This gives the desired result by (1.5).

Next we come to the proof of Theorem 1.4 on the densities ρ k . We also comment on the limitation on the analyticity domain in this result (cf. Remark 4.2).

Proof of Theorem 1.4: Take an integer k such that 0 < k < N . Let

x 0 = (x 0 1 ; • • • ; x 0 k ) ∈ U (1)
k . Setting m = k and p = N -k, we observe that U

(1) k coincide with the set U, defined in (3.5). Since -∆ = -∆ x -∆ y ∈ Diff 2 (R 3m ) and V -E is a potential of the type V in (3.3) on R 3m × R 3p , by Remark 3.4, we may apply Theorem 3.9 to H = H -E, Ω = R 3m , Ψ = ψ, and x 0 = x 0 ∈ U ∩ Ω. Thus there exist some neighbourhood Ω 0 of x 0 and a map U : Ω 0 → B 0 , with values in the set of unitary operators on W 0 , such that the map U ψ : Ω 0 x → U x ψ(x; •) ∈ W 2 is well-defined and real analytic. For x ∈ Ω 0 , the density ρ k (x) is, by definition (see (1.3)), the squared W 0 -norm of ψ(x; •), which is also the one of U x ψ(x; •), since U x is unitary on W 0 . Thus, ρ k is real analytic on Ω 0 .

Remark 4.1. In the k = 1 case, the above proof reproduces the one of [J] on the density ρ = ρ 1 . Thanks to Hunziker's twist, the proof in the general case is essentially identical.

Remark 4.2. Theorem 1.4 ensures the real analyticity of ρ k on U

(1) k . We actually can give two obstacles to the extension of this property to a larger set by the method used in Section 3. If we try to apply this method near a point x 0 ∈ C k ∪ R k , we have to work near the singularity of a term of the form |x j -R | -1 or |x j -x j | -1 with 1 ≤ j = j ≤ k and such a singular term would be unaffected by any twist. Thus, we cannot apply Theorem 2.2 to the twisted equation to get the result of Theorem 3.9. The second obstacle takes place in the construction of the twist near such a point x 0 ∈ C k ∪ R k . To explain this, let us use the notation of Subsection 3.3 and consider a point x 0 such that, either x 0 j = R or x 0 j = x 0 j , for j = j . Assume that we have a function f defined on a neighbourhood of such x 0 times R 3 such that, for x close to x 0 , f (x; •) is a diffeomorphism on R 3 that satisfies (3.8). If x 0 j = x 0 j , then, for x close to x 0 with different x j and x j , we would have x j = f (x; x 0 j ) = f (x; x 0 j ) = x j , which is a contradiction. If x 0 j = R , then, for x close to x 0 with x j = R , we would have x j = f (x; x 0 j ) = f (x; R ) = R , which is again a contradiction.

Finally, we prove Theorem 1.5 concerning the density matrices γ k . Again we provide a comment on the limitation on the analyticity domain in this result (cf. Remark 4.4).

Proof of Theorem 1.5: Let k be an integer between 0 and N . Let

x 0 = (x 0 1 ; • • • ; x 0 k ) and x 0 = (x 0 1 ; • • • ; x 0 k ) two points in (R 3 ) k such that (x 0 ; x 0 ) ∈ U (2)
k . Setting m = 2k and p = N -k, we note that U = U

(2) k . Denoting by y the variable in R 3(N -k) = R 3p , we consider the Hamiltonians H 1 and

H 2 on (R 3k ) 2 × R 3(N -k) = R 3m × R 3p defined by H 1 = -∆ x -∆ x -∆ y + V (x; y) and H 2 = -∆ x -∆ x -∆ y + V (x ; y)
and the almost everywhere defined functions Ψ 1 : R 3m ×R 3p (x; x ; y) → ψ(x; y) and Ψ 2 : R 3m × R 3p (x; x ; y) → ψ(x ; y). We note that H 1 ∈ Ham 2 (R 3m ) and H 2 ∈ Ham 2 (R 3m ). From Hψ = Eψ, we derive that H 1 Ψ 1 = EΨ 1 and H 2 Ψ 2 = EΨ 2 , in the distributional sense. Furthermore, for any bounded subset B of R 3m , we observe that

Ψ 1 ∈ W 2,2 B × R 3p y and Ψ 2 ∈ W 2,2 B × R 3p y .
Let B 0 be a bounded neighbourhood of (x 0 ; x 0 ). Now we apply Theorem 3.9 to Ω = B 0 and x 0 = (x 0 ; x 0 ) ∈ U ∩ Ω. Therefore there exist a bounded neighbourhood Ω 0 of (x 0 ; x 0 ) and a map U : Ω 0 → B 0 , with values in the set of unitary operators on W 0 , such that the maps U Ψ 1 : Ω 0 (x; x ) → U x;x Ψ 1 (x; x ; •) ∈ W 2 and U Ψ 2 : Ω 0 (x; x ) → U x;x Ψ 2 (x; x ; •) ∈ W 2 are well-defined and real analytic. Now, by definition of γ k (see (1.4)), we have, for (x; x ) ∈ Ω 0 ,

γ k (x; x ) = Ψ 1 (x; x ; •) , Ψ 2 (x; x ; •) = U x;x Ψ 1 (x; x ; •) , U x;x Ψ 2 (x; x ; •) , since U x;
x is unitary. As the scalar product of real analytic, W 0 -valued maps on Ω 0 , γ k is also real analytic on Ω 0 .

Remark 4.3. We note that our proofs of Theorems 1.3, 1.4, and 1.5, have a common structure. They all use an appropriate twist and Theorem 3.9. Differences between them occur in the used set of variables.

(2) k = (U (1) k ×U (1) k )\C (2) k .
5 Extensions.

In Section 3, we worked in a larger framework than the one that would be needed to treat the physical operator H in (1.1). It is thus natural to expect that our main results on H extend to a larger class of Hamiltonians. This is the case, as we shall see in this Section.

We consider again a system of N electrons, interacting to one another and moving under the influence of L fixed nuclei. Let I := {1; • • • ; N } be the set of indices for the electronic variables x ∈ R 3N and let I n := {1; • • • ; L} be the one of the nuclear variables R ∈ R 3 . It is convenient to introduce the disjoint union of those sets, namely

I

I n := (a; j) ; a ∈ {e; n} , j ∈ I a .

For c = (a; j) ∈ I e I n , let π 1 (c) = a and π 2 (c) = j. We define

X c = x j if π 1 (c) = e and X c = R j if π 1 (c) = n. Let P := {c; c } ; (c; c ) ∈ I e I n 2 and c = c
be the set of all possible particle pairings. We consider a potential V where

∀x ∈ R 3N , V (x) = {c; c }∈ P V {c; c } (X c -X c ) , (5.1)
with V {c; c } ∈ V (cf. Definition 3.1), for all {c; c } ∈ P. Now, let P = P (x; D x ) a differential operator on R 3N of the form

P = α∈N 3N |α|≤2 cα (x) D α x , (5.2)
where the coefficients cα belong to C ∞ b (R 3N ; C). Let H = P + V and assume that there exists ψ ∈ W 2,2 (R 3N ) such that H ψ = 0. For 0 < k < N , we still denote by ρ k , γ k , and C the objects that are respectively defined by (1.3), (1.4), and (1.5), with ψ replaced by ψ.

Theorem 5.1. Under the above assumptions on H and ψ, the results of Theorems 1.3, 1.4, and 1.5, are valid.

Proof: We can follow the above proofs of Theorems 1.3, 1.4, and 1.5. In the later, we use the elliptic operators

H 1 = -∆ x -P (x; y; D x ; D y ) + V (x; y) and H 2 = -∆ x -P (x ; y; D x ; D y ) + V (x ; y)
and the almost everywhere defined functions Ψ 1 : R 3m × R 3p (x; x ; y) → ψ(x; y) and Ψ 2 : R 3m × R 3p (x; x ; y) → ψ(x ; y). Remark 5.2. We point out that the present class of operators includes Laplace-Beltrami operators associated to appropriate metrics. One can also replace D x by D x -A(x), where A is a real analytic vector potential in C ∞ b (R 3N ; R 3 ), and add an external electric potential V e in C ∞ b (R 3N ; C), that is also real analytic. Notice that the class V in Definition 3.1 contains non-radial pair potentials. Among the admissible functions η used to define V, we have 0 < t → t -1 but we can also choose 0 < t → t -1 (ln(t)) -, with > 0. We observe that this class V is contained in the class of pair potentials used in [FS].

(the term containing z being in the kth position) and, for y ∈ R 3p , by (3.11), each term containing dτ being in the kth position. Since we can write, for x ∈ Ω(δ 0 ) and y ∈ R 3p ,

(d y F )(x; y) • y = p k=1 0 ; • • • ; 0 ; (d z f )(x; y k ) • y k ; 0 ; • • • ; 0 ∈ R 3p = y + r -1 0 p k=1 m j=1 0 ; • • • ; 0 ; (dτ ) r -1 0 (y k -x 0 j ) • y k (x j -x 0 j ) ; 0 ; • • • ; 0 ,
F -1 x; y = f -1 (x; y 1 ) ; • • • ; f -1 (x; y p ) ∈ R 3p , we have, for x ∈ R 3m , d x F -1 x; F (x; y) • x = d x f -1 x; f (x; y 1 ) • x ; • • • ; d x f -1 x; f (x; y p ) • x ∈ R 3p ,
where d x f -1 is given by (A.5). Finally, for y ∈ R 3p , we have

d y F -1 x; F (x; y) • y = p k=1 0 ; • • • ; 0 ; d z f -1 x; f (x; y k ) • y k ; 0 ; • • • ; 0 ∈ R 3p ,
where d z f -1 is explicited in (A.3). We derive from (A.6) that

F ∈ C ∞ b Ω(δ 0 ) × R 3p ; R 3p and F -1 ∈ C ∞ b Ω(δ 0 ) × R 3p ; R 3p . (A.7)
Next, we compute the terms J k (G), for 1 ≤ k ≤ 4 and G ∈ {F ; F -1 }, that appear in (3.16), (3.17 Finally, we compute the principal symbol σ P appearing in the proof of Proposition 3.6.

In the vector space of polynomials in R 3m × R 3p with C ∞ b (Ω(δ 0 ) × R 3p ; C)-coefficients, we denote by ≡ the equality of polynomials modulo polynomials of degree less than 2. Then, for (x; y; ξ; η) ∈ Ω(δ 0 ) × R 3p × R 3m × R 3p , we have, by (3.29), (3.17), and (3.19), σ P (x; y; ξ; η) ≡ S P (x; y; ξ; η) ≡ (α; β)∈N 3m ×N 3p |α|+|β|=2 c αβ x; F (x; y) ξ + J 1 F -1 (x; y) η + J 2 F -1 (x; y) α × J 3 F -1 (x; y) η + J 4 F -1 (x; y) β ≡ (α; β)∈N 3m ×N 3p |α|+|β|=2 c αβ x; F (x; y) ξ + J 1 F -1 (x; y) η α J 3 F -1 (x; y) η β and the last term is actually σ P (x; y; ξ; η), since it is homogeneous of degree 2. Thus σ P (x; y; ξ; η) = σ P x; F (x; y); ξ + J 1 F -1 (x; y) η; J 3 F -1 (x; y) η , yielding (3.30).

B A global pseudo-differential calculus.

In this section, we describe a global pseudo-differential calculus that is suited for Section 2. It is the one developed in the beginning of Chapter 18 in [START_REF] Hörmander | The analysis of linear partial differential operators III[END_REF].

With the notation of Subsection 2.2, we set N = m + p. For k ∈ Z, the class S k in [START_REF] Hörmander | The analysis of linear partial differential operators III[END_REF] (p. 65-75) is the set of smooth, complex valued functions σ on R 6N such that, for all (α; β) ∈ (N 3N ) 2 , there exists a constant C α;β > 0 such that, for all (x; y; ξ; η) ∈ R 6N , R 3N e i(x;y)•(ξ;η) σ(x; y; ξ; η) û(ξ; η) dξ dη , where (x; y) • (ξ; η) is the usual scalar product in R 3N of the vectors (x; y) and (ξ; η) and where û denotes the Fourier transform of u. It turns out that, for any s ∈ R, this operator σ(x; y; D x ; D y ) extends to a bounded operator from W s,2 (R 3N ) to W s-k,2 (R 3N ). If σ is the constant function equal to one on R 3N then σ ∈ S 0 and σ(x; y; D x ; D y ) is the identity operator I. For (k 1 ; k 2 ) ∈ Z 2 , for σ 1 ∈ S k 1 and σ 2 ∈ S k 2 , the composition σ 1 (x; y; D x ; D y )σ 2 (x; y; D x ; D y ) is given by σ(x; y; D x ; D y ) where σ ∈ S k 1 +k 2 and σ = σ 1 σ 2 + r with r ∈ S k 1 +k 2 -1 . Now, we are in a position to construct the parametrix Q needed in Subsection 2.2. We observe that P = σ(x; y; D x ; D y ) where its symbol σ belongs to S 2 . Furthermore, one can write σ = σ P + r with σ P ∈ S 2 and r ∈ S 1 . Let τ P ∈ C ∞ c (R 3N ) such that τ (ξ; η) = 1 if |ξ| 2 + |η| 2 ≤ 1. By the ellipticity assumption of P (see (2.6)), the symbol

q 1 := 1 -τ P • σ -1 P ∈ S -2
and, setting Q 1 = q 1 (x; y; D x ; D y ), there exist symbols r 1 , r 2 , r 3 ∈ S -1 such that Q 1 P = (σq 1 )(x; y; D x ; D y ) + r 1 (x; y; D x ; D y ) = (σ P q 1 )(x; y; D x ; D y ) + r 2 (x; y; D x ; D y ) = I + r 3 (x; y; D x ; D y ) .

Let Q = (I -r 3 (x; y; D x ; D y ))Q 1 . Then, Q = q(x; y; D x ; D y ) with q ∈ S -2 and there exists r 4 ∈ S -2 such that Q P = I -r 4 (x; y; D x ; D y ). Setting R := r 4 (x; y; D x ; D y ), we have Q P = 1 -R and (2.10) holds true.

Remark 4. 4 .

 4 According to the first obstacle mentioned in Remark 4.2, the method of Section 3 requires to work on U , because of the second obstacle mentioned there, we also need to exclude the set C(2) k and, therefore, work on U

ϕϕϕ

  ),(3.18), and (3.19). Let G ∈ {F ; F -1 }. Let us consider the smooth functions ρ ± ∈ C ∞ (Ω(δ 0 ) × R 3p ; R) defined byρ + (x; y) = Det (d y G)(x; y) 1/2 and ρ -(x; y) = Det d y G -1 (x; y) 1/2 .From the identity G -1 (x; G(x; y)) = (x; y), we derive that ρ+ (x; y)ρ -(x; G(x; y)) = 1. Let ϕ ∈ C ∞ (Ω(δ 0 ) × R 3p ; R). We set (U G ϕ)(x; y) = ρ + (x; y) ϕ x; G(x; y) and U -1 G ϕ (x; y) = ρ -(x; y) ϕ x; G -1 (x; y) . Since d x U -1 G ϕ (x; y) = ϕ x; G -1 (x; y) d x ρ -(x; y) + ρ -(x; y) (d x ϕ) x; G -1 (x; y) + ρ -(x; y) (d y ϕ) x; G -1 (x; y) d x G -1 (x; y) , we have, for x ∈ R 3m , (x; y)x = ϕ x; G -1 (x; y) ∇ x ρ -(x; y) • x + ρ -(x; y) (∇ x ϕ) x; G -1 (x; y) • x + ρ -(x; y) (∇ y ϕ) x; G -1 (x; y) • d x G -1 (x; y)x ,with scalar products in R 3m , R 3m , and R 3p , respectively. This yields∇ x U -1 G ϕ (x; y) = ϕ x; G -1 (x; y) ∇ x ρ -(x; y) + ρ -(x; y) (∇ x ϕ) x; G -1 (x; y) + ρ -(x; y) d x G -1 T (x; y)(∇ y ϕ) x; G -1 (x; y) .Thus, using ρ + (x; y)ρ -(x; G(x; y)) = 1,U G ∇ x U -1 G ϕ (x; y) = ρ + (x; y)ϕ(x; y) (∇ x ρ -) x; G(x; y) + (∇ x ϕ)(x; y) + d x G -1 T x; G(x; y) (∇ y ϕ)(x; y) . (x; y) = ϕ x; G -1 (x; y) d y ρ -(x; y) + ρ -(x; y) (d y ϕ) x; G -1 (x; y) d y G -1 (x; y) ,we have, for y ∈ R 3p , (x; y)y = ϕ x;G -1 (x; y) ∇ y ρ -(x; y) • y + ρ -(x; y) (∇ y ϕ) x; G -1 (x; y) • d y G -1 (x; y)y ,with scalar products in R 3p . This yields∇ y U -1 G ϕ (x; y) = ϕ x; G -1 (x; y) ∇ y ρ -(x; y) + ρ -(x; y) d y G -1 T (x; y)(∇ y ϕ) x; G -1 (x; y) . Thus, using ρ + (x; y)ρ -(x; G(x; y)) = 1, U G ∇ y U -1 G ϕ (x; y) = ρ + (x; y)ϕ(x; y) (∇ y ρ -) x; G(x; y) + d x G -1 T x; G(x; y) (∇ y ϕ)(x; y) .This proves (3.16), (3.17),(3.18), and (3.19). Now, taking into account (A.7), we obtain (3.20) and (3.21). Using (A.2) and (A.3), we see that (3.23) holds true.

( 1 +

 1 |ξ| 2 + |η| 2 ) |β|/2 ∂ α x;y ∂ β ξ;η σ(x; y; ξ; η) ≤ C α;β (1 + |ξ| 2 + |η| 2 ) k/2 ,where ∂ α x;y (resp. ∂ β ξ;η ) stands for the α (resp. β) partial derivative w.r.t. (x; y) (resp. (ξ; η)). Any such symbol σ defines an operator σ(x; y; D x ; D y ) from the Schwarz space S on R 3N into itself in the following linear way: for any u ∈ S, for any (x; y) ∈ R 3N , σ(x; y; D x ; D y )u (x; y) = (2π) -3N

Appendix.

A Basic computations.

In this section, we provide explicit formulae for several objects related to the twist U in Subsection 3.3 and the computation of the principal symbol of some differential operator, that is needed in Subsection 3.5.

First of all, let us notice that, from the expressions of d x f and d z f , there exists some τ , r 0 , and x 0 dependent constant M , such that, denoting by I 3 the identity on R 3 ,

Using (3.12), we can find some c > 0 such that

and z ∈ R 3 , we have

By differentiation of this equality and by using (3.9), we obtain, with equalities in L(R 3 ), in L(R 3 ), and in L(R 3m ; R 3 ), respectively,

Now, using (A.1), (A.2), (A.3), (A.4), and (A.5), it is straightforward to check that

From the very definition of F , we note that, for x ∈ Ω(δ 0 ), x ∈ R 3m , and y ∈ R 3p ,

and the k-th component (with k ∈ [[1; p]]) in the latter expression is given by


Lamberto Dell'elce

Alesia Herasimenka

Aaron J Rosengren

Nicola Baresi

EFFICIENT NUMERICAL SOLUTION OF THE LOW-THRUST LAMBERT PROBLEM

An algorithm for the numerical solution of low-thrust Lambert's problem is proposed. After averaging the extremal flow of the optimal control Hamiltonian, a one-parameter family of solutions of a reduced-order two-point boundary value problems is achieved by means of a differential continuation scheme. Sensitivities of the shooting function are then used in conjunction with an ad hoc near-identity transformation between averaged and osculating variables to achieve an accurate solution for all longitudes of the departure and arrival orbits. Hence, a single simplified shooting problem has to be solved to approximate the solution for any combination of departure and arrival dates (i.e., to draw a so-called "pork-chop chart"). Both the averaged flow and the near identity transformation are efficiently evaluated via the fast Fourier transform algorithm, yielding a fully numerical procedure.

INTRODUCTION

Lambert's problem in its classical form consists in finding a Keplerian orbit joining two position vectors in a given transfer time. Solutions of this problem are extensively used for preliminary mission design since they offer the identification of launch opportunities and a rough evaluation of their fuel cost by assuming impulsive maneuvers at the two boundary points. Pork-chop plots are often generated to illustrate total ∆V cost as a function of departure and arrival dates by recursively solving Lambert's problems. Recent interest was brought to the low-thrust counterpart of the problem. Because impulsive maneuvers are not allowed, the entire state vector is imposed at the two boundaries, i.e., both position and velocity of the satellite have to match the ones of departure and arrival bodies. Differently of the original problem, no exact closed form solution is available, so that assumptions on the shape and thrust direction are introduced to achieve efficient approximations [START_REF] Izzo | Lambert's problem for exponential sinusoids[END_REF] . [START_REF] Avanzini | Solution of Low-Thrust Lambert Problem with Perturbative Expansions of Equinoctial Elements[END_REF] Equivalent of pork-chop charts for low-thrust problems are referred to as bacon plots. [START_REF] Woolley | Low-thrust trajectory Bacon plots for Mars mission design[END_REF] In our view, two problems need to be addressed when tackling low-thrust transfers on a fixed maneuvering time. First, minimum thrust magnitude necessary to carry out the maneuver has to be identified, which, in turn, is a reachability problem. Second, once a sufficiently-large thrust is chosen, minimum energy maneuvers can be found. This work focuses on the first problem. A numerical methodology based on the averaging of the extremal flow of the optimal control system 4 is proposed. First, a reduced-order solution of the averaged two-point boundary value problem (TPBVP) parametrized by the costate of the fast variable is solved. This step requires the solution of a single shooting problem followed by a numerical continuation procedure. This problem is independent of the thrust magnitude. Second, sensitivity of the shooting functions are computed. These sensitivities are then used to compute perturbations of the averaged TPBVP associated to short-periodic variations and second-order terms. Finally, the thrust magnitude is identified by reconstructing a first-order approximation of the fast variables from the averaged solution.

Steps 1 and 2 are independent of the boundary conditions on the longitude, which is a major advantage of the methodology since these are the only two computationally demanding steps. In this way, a single one-parameter family of solutions of the averaged problem is sufficient to approximate the solution of the original problem for any departure and arrival dates. The methodology is finally applied to an Earth-Venus transfer. Discontinuities of the minimum thrust magnitude map are discussed and related to the cut locus of the problem. We ignore mass consumption and variable thrust profile in this case study because it is problem dependent. However, their inclusion is straightforward by appending mass to the vector of slow variables.

LOW-THRUST LAMBERT'S PROBLEM

Denote by x ∈ M ⊂ R n a set of integrals of motions of the two body problem (specifically, we use equinoctial elements) and by ϕ ∈ S 1 the mean longitude. Hereafter, x and ϕ are referred to as slow and fast variables, respectively. Let t 0 and t f be the desired departure and arrival dates. We are interested in finding the minimum thrust magnitude, ε, necessary to carry out the transfer between two Keplerian orbits with elements (x 0 , ϕ 0 (t 0 )) and (x f , ϕ f (t f )) in a time t f -t 0 . The notation ϕ 0 (t 0 ) denotes the value of the mean longitude of the departure orbit at epoch t 0 . For example, if an interplanetary transfer is envisaged, ϕ 0 (t) is given by the ephemeris of the departure planet at time t. Analogous notation holds for ϕ f (t f ). The optimal control problem that we tackle is thus min ε subject to:

d x d t = ε m i=1 f i (x, ϕ) u i d ϕ d t = ω(x) + ε m i=1 g i (x, ϕ) u i x(t 0 ) = x 0 , ϕ(t 0 ) = ϕ 0 (t 0 ) x(t f ) = x f , ϕ(t f ) = ϕ f (t f ) u ≤ 1 ∀ t ∈ [t 0 , t f ] (1)
where m = 3, and u is the control variable taking values in the unit sphere of R m . Vector fields f i and g i are periodic with respect to the fast variable, and they stem from Gauss variational equations (GVE) of the chosen set of orbital elements.

The pre-Hamiltonian of Problem (1) is defined as

H = λ ω(I) + ε m i=1 [p • f i (x, ϕ) + λg i (x, ϕ)] u i (2)
where p and λ denote adjoints of slow and fast variables, respectively. Application of Pontryagin maximum principle (PMP) yields the control action of the extremal flow u * (x, ϕ, p, λ), namely

u * = arg max u ≤1 H (x, p, ϕ, λ, ε, u) = K K where the i-th component of K is K i = p • f i (x, ϕ) + λg i (x, ϕ).
Replacing u * into Eq. ( 2) yields the Hamiltonian of the extremal flow:

H = λ ω(x) + εK(x, p, ϕ, λ) (3)
where K := K . Transversality conditions are such that adjoints are free at both t 0 and t f , so that necessary conditions for the optimality of Problem (1) are satisfied by finding a zero of the shooting function S(p(t 0 ), λ(t 0 ), ε) defined as

S(p(t 0 ), λ(t 0 ), ε) =   x * (t f ) -x f ϕ * (t f ) -ϕ f (t f ) p(0) -1   (4)
where x * (t) and ϕ * (t) denote trajectories of x and ϕ emanated by H with initial conditions x 0 , ϕ 0 , p(t 0 ), λ(t 0 ), and small parameter ε. The last equation of S is an arbitrary normalizing condition due to the homogeneity of H with respect to adjoints.

AVERAGED PROBLEM

States and adjoints are decomposed as

x = X + ε [δX + ν X (X, P , φ)] p = P + ε [δP + ν P (X, P , φ)] ϕ = ϕ 0 + Φ ε + δΦ + εν Φ (X, P , φ) λ = ε (Λ + ν Λ (X, P , φ)) + ε 2 δΛ (X, P , φ, Λ, δX, δP ) (5)
Hereafter, X and P are referred to as averaged state and adjoint, respectively, and ν (•) denotes short-periodic variations. Second-order terms * , namely δX, δP , δΦ, and δΛ, are necessary to have a precise enough approximation of the fast variable. The power of ε that multiplies each term of the decomposition (i.e., ε -1 , ε 0 , ε 1 , ε 2 ) is used to make explicit each respective formal order of magnitude, so that X, P , Φ, Λ, δX, δP , δΦ, δΛ, and any ν (•) (once stripped out of their power of ε) are all O(1) throughout the entire trajectory. Specifically, gross behavior of the fast variable is captured by Φ, so that

εϕ(t) = Φ(t) + O(ε) = O(1) when t = O ε -1 .
Similarly, adjoint Λ is well defined when ε approaches zero as discussed in. [START_REF] Dell'elce | On the convergence of time-optimal maneuvers of fastoscillating control systems[END_REF] Short-periodic variations and δΛ are evaluated at the angle:

φ = ϕ 0 + Φ ε + δΦ
The decomposition in Eq. ( 5) is not standard since second-order averaging is generally carried out by including second-order terms in the equations of motion of the averaged state without introducing new variables such as δX and δP . The main advantage of the proposed decomposition is that both the averaged counterpart of Problem (4) and the computation of second-order terms are εindependent. We also note that second-order short periodic variations are neglected since they do not affect the order of magnitude of the error in ϕ, so that Eq. ( 5) is not a full second-order approximation.

The average operator acting on a generic function f (ϕ) periodic with respect to ϕ is denoted with an overline, and it is defined as

f := 1 2π 2π 0 f (ϕ) dϕ
When the operator is applied to a function with other arguments than ϕ, these arguments are kept constant in the integral.

Equations of motion of averaged, short-periodic, and second-order terms are deduced by first substituting Eq. ( 5) into the Hamilton equations associated to H(x, p, ϕ, λ, ε), and then by expanding these equations up to the second-order in ε. Time derivatives of secular and second-order terms are finally achieved by averiging the resulting equations and collecting terms of the appropriate order (hence, they are independent of the fast variable by construction), whereas formal solution of short periodic variations is obtained by compensating first-order oscillatory dynamics. Outcome of this process is detailed in the reminder of the section.

To simplify the notation, omitted arguments of functions and their derivatives are understood as the average value of slow variables, namely X, P , εΛ = 0, and fast variable φ.

For example K = K(X, P , φ, 0), K = K(X, P , 0), ∂ K ∂ X = ∂ K ∂ x x=X,p=P ,ϕ= φ,λ=0 , ω = ω(X), and ∂ ω ∂ X = ∂ ω ∂ x x=X .

First-order terms and averaged shooting problem

Targeting ε-independent equations of motion, we introduce the re-scaling of time variable dτ = εdt Motion of averaged states and adjoints is given by

d X d τ = ∂ K ∂ P , d P d τ = - ∂ ω ∂ X Λ - ∂ K ∂ X , d Φ d τ = ω(X), d Λ d τ = 0 (6)
Equations ( 6) are hamiltonian with the first-order averaged Hamiltonian defined as

H = ω(X)Λ + K(X, P , 0)
Variable Φ is cyclic and its momentum Λ is a constant of motion.

The averaged shooting function, counterpart of Eq. ( 4), is

S(P (0), τ f |x 0 , x f , Λ) = X * (τ f ) -x f P (0) -1 (7)
We note that Problem (7) is independent of both ε and Φ, and it is a minimum time problem. Conversely, maneuvering time is fixed in both Problem (1) and its shooting function, Eq. ( 4). This dualism stems from the relation between minimum-time flow and boundary of the reachable set [5, Chap 10, Theorem 1]. The notation S(P (0), τ f |x 0 , x f , Λ) is introduced to emphasize that P (0) and τ f are shooting variables, whereas x 0 , x f , and Λ are parameters. Problems ( 4) and ( 7) have the same boundary conditions of slow variables, namely x 0 and x f . Corrections due to short-periodic variations and second-order terms are accounted for by computing sensitivities of S, as detailed in the next Section.

Short-periodic variations

Time derivative of short periodic variations is

d ν (•) d τ = d Φ d τ ∂ ν (•) ∂ φ + ε d δΦ d τ ∂ ν (•) ∂ φ + d X d τ • ∂ ν (•) ∂ X + d P d τ • ∂ ν (•) ∂ P = ω ∂ ν (•) ∂ φ + O(ε) (8)
The average of Eq. ( 8) is zero because derivatives of averaged and second-order variables do not depend on the fast variable and ν (•) = 0 by construction.

Short-periodic variations are obtained by solving

ω ∂ ν X ∂ φ = ∂ K ∂ P - ∂ K ∂ P ω ∂ ν P ∂ φ = - ∂ K ∂ X + ∂ K ∂ X - ∂ ω ∂ X ν Λ ω ∂ ν Φ ∂ φ = ∂ K ∂ λ - ∂ K ∂ λ + ∂ ω ∂ X • ν X ω ∂ ν Λ ∂ φ = - ∂ K ∂ φ (9)
with the conditions to resolve integration constants ν X = 0,

ν P = 0, ν Φ = 0, ν Λ = 0 (10)
Denote by f (k) the k-th coefficient of the Fourier series of a smooth periodic function f (ϕ), namely f (ϕ) = k∈Z f (k) e ikϕ , and by i the imaginary unit. Formal solution of Eqs. ( 9) and ( 10) is

ν X = - i ω k∈Z 0 1 k ∂ K (k) ∂ P e ik φ ν P = i ω k∈Z 0 1 k ∂ K (k) ∂ X - K (k) ω ∂ ω ∂ X e ik φ ν Φ = - i ω k∈Z 0 1 k ∂ K (k) ∂ λ - i k ω ∂ ω ∂ X • ∂ K (k) ∂ P e ik φ ν Λ = - K (X, P , φ, 0) -K (X, P , 0) ω(X) (11)
and it can be verified by substituting it into Eq. (9). We note that ν Λ establishes a first-order equivalence between original and averaged Hamiltonians. Short-periodic variations can be efficiently evaluated by means of the fast Fourier transform (FFT) algorithm.

Second-order terms

Equations of motion of δX, δP , and δΦ are obtained by expanding Eq. ( 3) up to the second order in ε, averaging, and discarding first-order terms (which where already used in the previous sub-sections). This yields * d δX

d τ = ∂ 2 K ∂X ∂P (δX + ν X ) + ∂ 2 K ∂P 2 (δP + ν P ) + ∂ 2 K ∂ϕ ∂P ν Φ + ∂ 2 K ∂λ ∂P (Λ + ν Λ ) d δP d τ = - ∂ 2 K ∂X 2 (δX + ν X ) - ∂ 2 K ∂P ∂X (δP + ν P ) - ∂ 2 K ∂ϕ ∂X ν Φ - ∂ 2 K ∂λ ∂X (Λ + ν Λ ) -Λ ∂ 2 ω ∂x δX - ∂ ω ∂ x δΛ d δΦ d τ = ∂ ω ∂ X • δX + ∂ K ∂ λ (12)
Higher-order terms of Eq. ( 8) do not contribute to Eq. ( 12) because their average is zero, as previously mentioned. If the determination second-order short periodic variations were envisaged, they could not be neglected.

Expansion of the hamiltonian up to the second order in ε enables an explicit expression of δΛ as a function of the other variables, namely

δΛ = - 1 ω ∂ K ∂ X + (Λ + ν Λ ) ∂ ω ∂ X • (δX + ν X ) + ∂ K ∂ P • (δP + ν P ) + ∂ K ∂ Φ ν Φ + ∂ K ∂ λ (Λ + ν Λ ) Denoting ν (k) X := - i ω k ∂ K (k) ∂ P , ν (k)
P := i ω k ∂ K (k) ∂ X - K (k) ω ∂ ω ∂ X , ν (k)
Φ := 1 ω ∂ ω ∂ X • ν (k) X + ∂ K (k) ∂ λ
the k-th coefficients of the series in Eq. ( 11), the averaged value of δΛ (which appears in the time derivative of δP ) is given by

δΛ = - 1 ω ∂ K (0) ∂ X + Λ ∂ ω ∂ X • δX + ∂ K (0) ∂ P • δP + ∂ K (0) ∂ λ Λ - 1 ω Z 0 ∂ K (-k) ∂ X + ν (-k) Λ ∂ ω ∂ X • ν (k) X + ∂ K (-k) ∂ P • ν (k) P -ν Φ (-k) + 1 ω ∂ K (-k) ∂ λ K (k)
where we extensively used the integration rule of the average product of two periodic functions f (ϕ) and g(ϕ)

f (ϕ) g(ϕ) = k∈Z f (k) g (-k) (13)
* The convention

∂ 2 K ∂X ∂P [l,m] = ∂ 2 K ∂Xm ∂P l is used for second derivatives. Indeed, ∂ 2 K ∂X ∂P = ∂ 2 K ∂P ∂X T .
Replacing Eq.(11) into Eq. ( 12), integrating by parts terms that include derivatives of ϕ, and again using the rule of Eq. ( 13), yields

d δX d τ = ∂ 2 K (0) ∂X ∂P δX + ∂ 2 K (0) ∂P 2 δP + ∂ 2 K (0) ∂λ ∂P Λ + k∈Z 0 ∂ 2 K (-k) ∂X ∂P ν (k) X + ∂ 2 K (-k) ∂P 2 ν (k) P - ∂ K (-k) ∂ P ν (k) Φ - ∂ 2 K (-k) ∂λ ∂P K (k) ω d δP d τ = - ∂ 2 K (0) ∂X 2 + Λ ∂ 2 ω ∂X 2 δX - ∂ 2 K (0) ∂P ∂X δP - ∂ 2 K (0) ∂λ ∂X Λ - ∂ ω ∂ X δΛ - k∈Z 0 ∂ 2 K (-k) ∂X 2 ν (k) X + ∂ 2 K (-k) ∂P ∂X ν (k) P - ∂ K (-k) ∂ X ν Φ (k) - ∂ 2 K (-k) ∂λ ∂X K (k) ω d δΦ d τ = ∂ ω ∂ X • δX + ∂ K (0) ∂ λ (14)
Equation ( 14) is a linear system in δX, δP , and δΦ. Although these expressions appear cumbersome, they are obtained by elementary combinations of Fourier coefficients of K and their derivatives. As such, this formulation lends itself to a practical implementation by leveraging on automatic differentiation (AD) of the Hamiltonian 6 and FFT.

ALGORITHM FOR THE GENERATION OF BACON PLOTS

We propose an algorithm for the generation of approximated solutions of Problem (1).

1. Find a one-parameter family of solutions of the average problem ∀ Λ ∈ Λ min , Λ max find P Λ 0 , τ Λ f s.t.: S P Λ 0 , τ Λ f |x 0 , x f , Λ = 0 We note that the same initial conditions of the osculating problem, x 0 and x f , are used here. After finding the solution for Λ = 0, a numerical continuation scheme can be used (Λ serves as continuation parameter). Although fold bifurcations may occur and multiple solutions may be found for the same value of Λ, we prefer to formally use Λ as continuation parameter instead of introducing a pseudo-arc length parameter to streamline the notation.

Compute sensitivities

∀ Λ ∈ Λ min , Λ max evaluate ∂ S ∂ τ f , ∂ S ∂ P 0 , ∂ S ∂ x 0 , ∂ S ∂ x f , ∂ Φ f ∂ τ f , ∂ Φ f ∂ P 0 , ∂ Φ f ∂ x 0 where Φ f = Φ τ Λ f | x 0 , P Λ 0 ,
Λ is the final re-scaled fast variable of the average system.

Evaluate perturbations of the bondary conditions

find δτ Λ f , δP Λ 0 , δΦ Λ f s.t.: ∂ S ∂ τ f δτ Λ f + ∂ S ∂ P 0 δP Λ 0 = - ∂ S ∂ x 0 δX 0 - ∂ S ∂ x f δX f δΦ Λ f = ∂ Φ f ∂ τ f δτ f + ∂ Φ f ∂ P 0 δP 0 + ∂ Φ f ∂ x 0 δX 0 + δΦ τ Λ f
where δX 0 and δX f include both short-periodic variations and second-order terms, namely

δX 0 = -ν X x 0 , P Λ 0 , ϕ 0 (t 0 ) δX f = -ν X x f , P τ Λ f , ϕ f (t f ) -δX τ Λ f
Here, δX τ Λ f and δΦ τ Λ f are evaluated by integrating Eq. ( 14) with initial conditions δX = δP = 0 and δΦ = 0.

4.

Evaluate minimum ε and identify the good members of the family of average solutions

Find Λ s.t.: mod ϕ 0 (t 0 ) + Φ Λ f ε Λ + δΦ Λ f -ϕ f (t f ), 2π = 0 where ε Λ = τ Λ f t f -t 0 -δτ Λ f
is such that the duration of the trajectory equals t f -t 0 .

Concerning the complexity of the methodology: Step 1 implies first the solution of a single average problem for Λ = 0 (this is the only shooting problem that has to be solved without a priori knowledge). Then, differential continuation can be used to obtain the solution in the range Λ min , Λ max . Step 2 requires the numerical integration of the average system's linearization. AD can be used to evaluate derivatives. Specifically, we used the software Hampath * described in [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF] to compute derivatives of both H and S which are necessary in Steps 1 and 2, respectively. Hampath only requires in input the maximized Hamiltonian and the shooting function, and it also offers a continuation algorithm based on differential homotopy that can be used to compute the one-parameter family of solutions P Λ 0 , τ Λ f . Step 3 requires the integration of second-order terms, the evaluation of short-periodic corrections (a single FFT is sufficient to evaluate ν (•) for any ϕ and, as such, for any t 0 or t f ), and the solution of n-dimensional linear systems. Step 4 only requires elementary algebraic operations.

We stress that all steps only use the average system and the near-identity transformation. Departure and arrival dates only appear in Steps 3 and 4, so that outcomes of Steps 1 and 2, which involve the only computationally demanding tasks, need to be evaluated only once for any combination of t 0 and t f .

CASE STUDY

The proposed case study consists of an Earth-Venus low-thrust transfer. Table 1 lists simulation parameters necessary for repeatability.

The relevance of second-order terms is disclosed in Figure 1. Contributions of average, secondorder, and short-periodic variations are pointed out in Figure 1(a), where a sample trajectory of the adjoint p 4 is drawn. Drift between p 4 and P 4 is of order ε at the end of the maneuver, and (for this particular trajectory) including δP 4 is more important than short-periodic variations ν P 4 to reduce the error with the original trajectory. In general, contributions of second-order terms and shortperiodic variations are of the same order of magnitude. Necessity to account for second-order terms stems from the need of a first-order estimation of the fast variable. As depicted in Figure 1(b), error between ϕ and Φ converges to δΦ (which is of order one) as ε approaches zero. Indeed, a drawback of the proposed framework is that quality of the approximation is tied to the small magnitude of ε. Hence, the methodology is not adequate to approximate short transfers, which are characterized by large values of ε.

An example of bacon plot is illustrated in Figure 2(a). Solutions along red-dashed lines have the same maneuvering time, which is larger for upper lines. Variations of ε along these lines are small when long transfers are envisaged, i.e., importance of the position on departure and arrival orbits dwindles when several revolutions are necessary to complete the maneuver. Bacon plots are not smooth. Both jump and derivative discontinuities are possible. The first ones arise when new local solutions are generated, but this does not happen in the case study at hand. Instead, discontinuities of the derivative, which appear in Figure 2(a), occur when two existing local solutions have the same ε. This can be better appreciated by inspecting Figure 2(b), which depicts the value of Λ. In fact, each point of the bacon plot is associated to a precise member of the one-parameter family of average solutions. Discontinuities of Λ (occurring at the interface of yellow and blue regions of this figure) are associated to coexistence of two local solutions with the same cost function ε. This mechanism is further emphasized in Figure 3, which shows the ε of local solutions along the red-dashed line of Figure 2(b). Two distinct solutions are identified for some values of the departure date, but only the one with lower ε is of interest. Red dot indicates a cut point, i.e., where the two local solution have the same ε, and one branch stops tracking the global optimum. At this point, a discontinuity occur in the charts of Figure 2 (derivative of ε, and jump of Λ and of all other shooting variables). In conclusion, discontinuities stem from crossing of the cut locus, whose immediate identification is a major asset of the proposed methodology.

CONCLUSION

This paper offers an efficient numerical solution to the low-thrust Lambert problem which may serve as a tool for preliminary mission design. Overall, the methodology requires the solution of a single simplified TPBVP to draw a full chart of minimum thrust magnitude for any combination of departure and arrival dates. A major advantage of the methodology is the immediate identification of cut points enabling the straightforward discarding of non-minimizing local solutions.

Formal solution of short-periodic variations and an implementable expression of time derivatives of second-order terms are provided. These expressions can be efficiently evaluated by means of the FFT algorithm.

Variable thrust profiles were neglected in this study, but their inclusion in the methodology is straightforward. Perspectives of the work encompass solution of minimum-fuel transfers, and the possibility to account for ∆V at the departure and arrival dates, yielding a "hybrid" Lambert problem.

( a )

 a Adjoint p4. Here, Λ = -1 and ε = 0.05 mN/kg. 2 mN/kg = 10 -2.5 mN/kg = 10 -3 mN/kg = 10 -3.5 mN/kg (b) Error between true fast variable and first-order approximation, ε . We note that such error converges to the second-order correction, δΦ, as ε approaches 0. Here, Λ = -1.

Figure 1 .

 1 Figure 1. Importance of second-order terms.

  (a) Minimum thrust magnitude. (b) Continuation parameter, Λ, of best solutions identified in Step 4.

Figure 2 .Figure 3 .

 23 Figure 2. Earth-Venus transfer for various deparutre and arrival dates.

Table 1 .

 1 Earth-Venus transfer. Elements of departure and arrival orbits.

	Earth	Venus

* The reference "second order" may stride with the order of magnitude these terms in the expansion of Eq. (5). Instead, they are referred to as second order terms because time derivatives of εδX and εδP are O(ε 2 ). For the sake of clarity and without arguing further on this notation, also δΦ and δΛ are referred to as "second order terms" for analogy.

* Website: https://www.hampath.org/. New developments of Hampath are available in the control-toolbox library: https://ct.gitlabpages.inria.fr/gallery/