A sparse approximation of the Lieb functional with moment constraints

Abstract

The aim of this paper is to present new sparsity results about the so-called Lieb functional, which is a key quantity in Density Functional Theory for electronic structure calculations of molecules. The Lieb functional was actually shown by Lieb to be a convexification of the so-called Lévy-Lieb functional. Given an electronic density for a system of $N$ electrons, which may be seen as a probability density on ${\mathbb{R}}^{3}$, the value of the Lieb functional for this density is defined as the solution of a quantum multi-marginal optimal transport problem, which reads as a minimization problem defined on the set of trace-class operators acting on the space of electronic wave-functions that are anti-symmetric $L^{2}$ functions of ${\mathbb{R}}^{3N}$, with partial trace equal to the prescribed electronic density. We introduce a relaxation of this quantum optimal transport problem where the full partial trace constraint is replaced by a finite number of moment constraints on the partial trace of the set of operators. We show that, under mild assumptions on the electronic density, there exist sparse minimizers to the resulting moment constrained approximation of the Lieb (MCAL) functional that read as operators with rank at most equal to the number of moment constraints. We also prove under appropriate assumptions on the set of moment functions that the value of the MCAL functional converges to the value of the exact Lieb functional as the number of moments go to infinity. We also prove some rates of convergence on the associated approximation of the ground state energy. We finally study the mathematical properties of the associated dual problem and introduce a suitable numerical algorithm in order to solve some simple toy models.

Conflict of interest statement:

The authors have no competing interests to declare that are relevant to the content of this article.

Data availability statement:

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

Contents

1. 1 Introduction
2. 2 The exact Lieb functional
3. 3 Moment-constrained approximation and sparsity result
   1. 3.1 Tchakaloff’s theorem on Hilbert spaces
   2. 3.2 Existence of sparse minimizers for Moment Constrained Approximation of Lieb (MCAL) functional
4. 4 Some convergence results
   1. 4.1 Convergence of the MCAL functional to the exact Lieb functional
   2. 4.2 Convergence rate of the ground state energy in the bounded domain case
5. 5 Duality results for the MCAL functional
   1. 5.1 Semi-definite positive programming problems
   2. 5.2 Dual MCAL problem
6. 6 Numerical scheme
   1. 6.1 MCAL iterative scheme
      1. 6.1.1 Initialization step
      2. 6.1.2 Iteration $n\geq 1$
   2. 6.2 Property of the MCAL iterative scheme
   3. 6.3 Numerical results
7. 7 Proofs
   1. 7.1 Proof of Theorem 6
   2. 7.2 Proof of Theorem 17
   3. 7.3 Proof of Theorem 24

1 Introduction

The so-called Hohenberg-Kohn or Lévy-Lieb functional plays a fundamental role in Density Functional Theory for electronic structure calculations. For the sake of simplicity, we use here atomic units and neglect the effect of spin in this work. For a given electronic density $\rho\in L^{1}(\mathbb{R}^{3})$, which we assume here to be of integral equal to $1$ for the sake of simplicity, and a given number of electrons $N\in\mathbb{N}^{*}$, the Lévy-Lieb functional $F_{LL}(\rho)$ reads as the solution of the following a minimization problem of the form:

$$\boxed{F_{LL}[\rho]:=\mathop{\inf}_{\begin{subarray}{c}\Psi\in{\cal H}_{1}^{N}\\ \rho_{\Psi}=\rho\end{subarray}}\frac{1}{2}\int_{{{\mathbb{R}}}^{3N}}|\nabla\Psi|^{2}+\int_{\mathbb{R}^{3N}}V|\Psi|^{2},}$$

where

• $\mathcal{H}_{1}^{N}:=\bigwedge_{i=1}^{N}H^{1}(\mathbb{R}^{3})$ is the set of admissible electronic wavefunctions for a system of $N$ electrons with finite kinetic energy, that is the set of antisymmetric functions of $H^{1}(\mathbb{R}^{3N})$;

• for any $\Psi\in\mathcal{H}_{1}^{N}$ and $x\in\mathbb{R}^{3}$, $\rho_{\Psi}$ is the electronic density associated to the wavefunction $\psi$;

• the function $V:(\mathbb{R}^{3})^{N}\to\mathbb{R}_{+}\cup\{+\infty\}$ is the electron-electron Coulomb interaction potential.

There is a wide zoology of electronic structure calculation models which rely on various types of approximations of this Lévy-Lieb functional. Recently, Strictly Correlated Electrons (SCE) based approximation of this functional have drawn an increasing interest from mathematicians because it gives rise to a symmetric multi-marginal optimal transport problem with Coulomb cost, with the number of marginal constraints being equal to the number of electrons in the system. The literature about the SCE approximation (namely the multi-marginal optimal transport with Coulomb cost) is growing considerably. Recent developments include results on the existence and non-existence of Monge-type solutions (e.g., [CD15, CDD15, CFK13, Fri19, BDGG12, CS16, BDPK20]), structural properties of Kantorovich potentials (e.g., [CDMS19, DGN17, GKR19, BCD17]), grand-canonical optimal transport [DMLN22], efficient computational algorithms (e.g., [BCN17, FSV22, CEL⁺19, MG19, KLLY19]) and the design of new density functionals (e.g., [GGGG19, CF15, MUMIGG14, LDMG⁺16]).

Moreover, recent works indicate that the solution of this symmetric Coulomb cost multi-marginal problem (MMOT), which is a probability measure on $\mathbb{R}^{3N}$, is actually a sparse object at least in discrete settings. Two types of discrete settings have been considered so far where such sparsity results have been obtained. On the one hand, the most classical discrete approximation consists in introducing a discrete grid $\mathcal{X}$ of $\mathbb{R}^{3}$. The discrete optimal transport plan is then defined as a discrete probability measure defined on the cartesian product grid $\mathcal{X}^{N}$. Actually, it was proved in [FV18, Vög21] that the discrete optimal transport plan does not charge all the points of the discrete cartesian product grid (of cardinality $|\mathcal{X}|^{N}$) but only a number of points in this grid which scales at most linearly with $M$. Finding the few points of $\mathcal{X}^{N}$ which are actually charged by the discrete optimal transport plan is not a trivial task though, and the GenCol algorithm is a numerical procedure which aims at achieving this task. It has been first proposed in [FSV22], then extended in [FP22] and its convergence has been analyzed for two-marginal problems in [FP23]. On the other hand, an alternative approach which was first considered in [ACEL21] consists in introducing an approximation of the exact multi-marginal transport problems where the marginal constraints are replaced by a finite number of moment constraints associated to a finite number $M$ of ”moment functions” which are real-valued functions defined on $\mathbb{R}^{3}$. Under some natural assumptions, this approximate problem is then equivalent to approximating the solution of the dual problem associated to the exact optimal transport problem, namely the so-called Kantorovich potential, as a linear combination of these moment functions. The solution of this moment-contrained optimal transport problem is still a probability measure defined on $\mathbb{R}^{3N}$ but is also a sparse object in the sense that it can be written as a discrete measure charging a number of points belonging to $\mathbb{R}^{3N}$ which scales at most linearly with the number of moment constraints. Finding the location of these points then reads as a non-convex optimization problem defined on a continuous (and not a discrete set) set, and stochastic gradient algorithms have been proposed in [ACE22] in order to find such optimal points, and numerically tested on three-dimensional settings involving $N=100$ electrons. We also refer the reader to the works [CFM14, BCN17, NP22, Lel22, HCL23] where alternative numerical methods have been proposed for the computation of the SCE limit of the Lévy-Lieb functional, which do not rely on sparsity arguments.

The objective of this work is to prove similar sparsity results for the so-called Lieb functional, which is a convex relaxation of the Lévy-Lieb functional, the expression of which is given under the following form:

$$\boxed{F_{L}[\rho]:=\mathop{\inf}_{\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}),\;\rho_{\Gamma}=\rho}{\rm Tr}\left[\bigg{(}-\frac{1}{2}\Delta+V)\Gamma\bigg{)}\right],}$$ (1)

where $\displaystyle\mathcal{H}_{0}^{N}:=\bigwedge_{i=1}^{N}L^{2}(\mathbb{R}^{3})$, $\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})$ denotes the set of non-negative trace-class self-adjoint operators on $\mathcal{H}_{0}^{N}$ and where $\rho_{\Gamma}$ is the electronic density associated to $\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})$, the precise definition of which will be given below. Actually, problem (1) is a particular instance of quantum optimal transport problem. We refer the reader to [GMP16, GP17] for references on earlier works on closely related types of problems. Notice that problem (1) can be understood as a quantum version of a multi-marginal optimal transport problem. Moreover, it still enjoys the nice property, as the original problem, of being a linear programming problem. Our aim here is to prove that solutions of approximations of problems (1) where the partial trace constraint is relaxed by a finite number of moment constraints enjoy similar sparsity properties than solutions of moment constrained multi-marginal symmetric classical optimal transport problems, such as those which were established in [ACEL21]. More precisely, we prove, using the so-called Tchakaloff’s theorem (notice that for the usual entropic regularization of MMOT we cannot use this kind of approach), that the solutions of moment constrained approximations of (1) can be written under the form $\Gamma=\sum_{k=1}^{K}\alpha_{k}|\Psi_{k}\rangle\langle\Psi_{k}|$, where $K\in\mathbb{N}^{*}$ scales at most linearly with the number of moment constraints, and where for all $1\leq k\leq K$, $\alpha_{k}\in[0,1]$, $\Psi_{k}\in\mathcal{H}_{1}^{N}$ and $|\Psi_{k}\rangle\langle\Psi_{k}|$ is the orthogonal projector of $\mathcal{H}_{0}^{N}$ onto the vectorial space spanned by $\Psi_{k}$ (using bra-ket notation). We will, finally, exploit this sparsity structure in order to propose some numerical scheme in order to approximate the solution of (1). Notice that solving (5) for a small systems can be exploited (as done in some recent works for the Levy-Lieb functional [SPTP23, BVMTG22]) in order to build approximations of the Lieb functional for larger systems by means of machine learning techniques. Let us finally mention here that particular moment-constrained approximations of the Lieb functional have already been considered in [Gar22] for the construction of Kohn-Sham potentials. The novel results brought by the present contribution in comparison to the latter work is (i) the extension of existence and convergence results to more general moment constraints that the one considered in [Gar22]; (ii) the results on the sparsity structure of associated minimizers; (iii) convergence rate of the approximate ground state energy; and (iv) study of the mathematical properties of the associated dual problem.

The outline of the article is the following. In Section 2, we recall some fundamental results about the exact Lieb functional. The moment-constrained approximation we consider here and the associated sparsity result on their minimizers is presented in Section 3. Convergence results of the moment-constrained approximation towards the exact Lieb functional are presented in Section 4.1. In Section 4.2, we also prove some rates of convergence of the associated approximation of the ground state energy to the exact one. Section 5 is devoted to present some results about the dual formulation of the moment-constrained problem in the case of electronic density with support included in bounded domains. We, finally, introduce a numerical method in Section  6 exploiting the sparsity result and the convenient dual formulation as an SDP problem. Some numerical experiments for small systems are then predented.

2 The exact Lieb functional

Let us first introduce some notation together with the problem we consider in this work. We use here atomic units and neglect the influence of spin for the sake of simplicity.

Let $N\in{\mathbb{N}}^{*}$ denote the number of electrons in the molecule of interest. Let us assume that there are $N_{\rm nu}\in{\mathbb{N}}^{*}$ nuclei in the molecule, the positions and electric charges of which are denoted by $R_{1},\ldots,R_{N_{\rm nu}}\in{\mathbb{R}}^{3}$ and $Z_{1},\ldots,Z_{N_{\rm nu}}\in{\mathbb{N}}^{*}$. For all $x\in{\mathbb{R}}^{3}$, let us denote by

$$v_{\rm nu}(x):=-\sum_{n=1}^{N_{\rm nu}}\frac{Z_{n}}{|x-R_{n}|}$$

the Coulomb electric potential generated at $x\in{\mathbb{R}}^{3}$ by the $N_{\rm nu}$ nuclei.

Let ${\cal H}:=H^{1}({\mathbb{R}}^{3})$ and ${\cal H}^{N}:=\bigwedge_{i=1}^{N}H^{1}({\mathbb{R}}^{3})$. For any $\Psi\in{\cal H}^{N}$, we denote by $\|\Psi\|$ its $L^{2}(\mathbb{R}^{3N})$ norm and by $\rho_{\Psi}$ the electronic density associated to the wavefunction $\Psi$, namely the real-valued function defined over $\mathbb{R}^{3}$ as follows:

$$\forall x\in\mathbb{R}^{3},\quad\rho_{\Psi}(x):=N\int_{(\mathbb{R}^{3})^{N-1}}|\Psi(x,x_{2},\ldots,x_{N})|^{2}\,dx_{2}\ldots\,dx_{N}.$$

For a given set of nuclei positions ${\bm{R}}:=(R_{1},\ldots,R_{N_{\rm nu}})$ and charges ${\bm{Z}}:=(Z_{1},\ldots,Z_{N_{\rm nu}})$, one can compute the ground state energy as a minimization over a density $\rho$, that is

$$E[{\bm{R}},{\bm{Z}}]=\inf_{\begin{subarray}{c}\rho\in{\cal I}^{N}\end{subarray}}\left\{F_{LL}[\rho]+\int_{{\mathbb{R}}^{3}}v_{\rm nu}\rho\right\},$$ (2)

where ${\cal I}^{N}:=\{\rho\in L^{1}(\mathbb{R}^{3}),\;\rho\geq 0,\;\sqrt{\rho}\in H^{1}({\mathbb{R}}^{3}),\;\int_{\mathbb{R}^{3}}\rho=N\}$ and

$$F_{LL}[\rho]:=\inf_{\begin{subarray}{c}\Psi\in{\cal H}_{1}^{N}\\ \rho_{\Psi}=\rho\end{subarray}}\left\{\frac{1}{2}\int_{{\mathbb{R}}^{3N}}|\nabla\Psi|^{2}+\int_{{\mathbb{R}}^{3N}}V|\Psi|^{2}\right\}$$ (3)

is called the Levy-Lieb functional. In ((3)), the function $V:(\mathbb{R}^{3})^{N}\to\mathbb{R}_{+}\cup\{+\infty\}$ is defined as follows: for all $(x_{1},\ldots,x_{N})\in(\mathbb{R}^{3})^{N}$,

$$V(x_{1},\ldots,x_{N})=\sum_{1\leq i<j\leq N}\frac{1}{|x_{i}-x_{j}|}.$$ (4)

The Levy-Lieb functional is the central object in Density Functional Theory and its knowledge would allow the computation the electronic ground state energy of any molecule. However, it turns out that $F_{LL}$ is not convex, it is therefore convenient to look at a convexification proposed by Lieb [Lie83a] where the minimization is performed over the set of mixed states instead of the set of pure ones as in (3). More precisely, we consider here the alternative minimization problem

$$F_{L}[\rho]:=\inf_{\begin{subarray}{c}\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})\\ \rho_{\Gamma}=\rho\end{subarray}}{\rm Tr\,}(H_{N}\Gamma),$$ (5)

where $H_{N}:=-\frac{1}{2}\Delta+V$ is a self-adjoint operator on $\mathcal{H}_{0}^{N}$ with domain $D(H_{N})=\mathcal{H}_{2}^{N}:=\bigwedge_{i=1}^{N}H^{2}(\mathbb{R}^{3})$, $\mathfrak{S}_{1}^{+}({\cal H}_{0}^{N})$ denotes the set of trace-class self-adjoint non-negative operators on ${\cal H}_{0}^{N}$. For all $\Gamma\in\mathfrak{S}_{1}^{+}({\cal H}_{0}^{N})$, there exists an orthonormal basis $(\Psi_{i})_{i\in\mathbb{N}^{*}}$ of ${\cal H}_{0}^{N}$ and a non-increasing sequence $(\alpha_{i})_{i\in\mathbb{N}^{*}}$ of non-negative numbers such that

$$\Gamma=\sum_{i=1}^{+\infty}\alpha_{i}|\Psi_{i}\rangle\langle\Psi_{i}|,$$ (6)

using so-called bra-ket notation. Then, the associated electronic density $\rho_{\Gamma}$ is defined as follows: for all $x\in\mathbb{R}^{3}$,

$$\rho_{\Gamma}(x):=N\sum_{i=1}^{+\infty}\alpha_{i}\int_{(\mathbb{R}^{3})^{N-1}}|\Psi_{i}(x,x_{2},\ldots,x_{N})|^{2}\,dx_{2}\ldots\,dx_{N}=\sum_{i=1}^{+\infty}\alpha_{i}\rho_{\Psi_{i}}(x).$$

We know that there exist positive constants $\varepsilon,D>0$ such that $H_{N}+D\geq\varepsilon(-\Delta+{\rm Id})$ (in the sense of self- adjoint operators on $\mathcal{H}_{0}^{N}$). We also denote by $\mathfrak{S}_{1,1}(\mathcal{H}_{0}^{N})$ the set of self-adjoint operators $\Gamma$ on $\mathcal{H}_{0}^{N}$ with finite kinetic energy, i.e. such that ${\rm Tr}\left(|H_{N}+D|^{1/2}\Gamma|H_{N}+D|^{1/2}\right)<+\infty$.

Remark 1.

It can then be easily checked that, $\Gamma\in\mathfrak{S}_{1,1}(\mathcal{H}_{0}^{N})$ if and only if $\Gamma\in\mathfrak{S}_{1}(\mathcal{H}_{0}^{N})$ and ${\rm Tr}(H_{N}\Gamma)<+\infty$. Then, if $\Gamma$ admits an eigendecomposition of the form (6), necessarily $\Psi_{i}\in\mathcal{H}_{1}^{N}$ as soon as $\alpha_{i}>0$.

It is well-known then that the infimum in (3) and (5) is attained.

Remark 2 (Convexification).

It is worth highlighting that $F_{L}$ is indeed the convexification of $F_{LL}$ in the sense that

$$F_{L}[\rho]=\inf_{\begin{subarray}{c}\forall i\geq 1,\;\alpha_{i}\geq 0,\;\rho_{i}\in\mathcal{I}^{N}\\ \sum_{i=1}^{+\infty}\alpha_{i}=1\\ \sum_{i=1}^{+\infty}\alpha_{i}\rho_{i}=\rho\\ \end{subarray}}\sum_{i=1}^{+\infty}\alpha_{i}F_{LL}[\rho_{i}]$$

It is useful noticing that $F_{L}$ admits a dual problem.

Theorem 3 ([Lie83a]).

Duality holds in the sense that

$$F_{L}[\rho]=\sup_{\begin{subarray}{c}v\in L^{3/2}({\mathbb{R}}^{3})+L^{\infty}({\mathbb{R}}^{3})\\ H^{v}_{N}\geq 0\end{subarray}}\left\{\int_{{\mathbb{R}}^{3}}v(x)\rho(x)\mathrm{d}x\right\},$$ (7)

where

$$H^{v}_{N}=H_{N}-\sum_{i=1}^{N}v(x_{i}).$$

The constraint in (7) has to be understood in the sense of self-adjoint operators, namely for all $\Psi\in\mathcal{H}_{1}^{N}$, $\langle\Psi|H^{v}_{N}|\Psi\rangle\geq 0$.

Remark 4.

It is important to notice, for the following, that it can be easily proved that the infimum in (3) and (5) is attained. However, it happens that the supremum in (7) is not attained for most densities $\rho$ (we refer the reader to [LLS19]).

3 Moment-constrained approximation and sparsity result

We focus now on a first approximation of (5) by using the moment constraint approach which has previously been studied in the framework of classical optimal transport [ACEL21, ACE22]. We also refer to [Gar22] where a particular instance of moment-constrained approximation of the Lieb functional has been considered for the computation of Kohn-Sham potentials.

We begin by introducing here some notation. From now on, we fix an electronic density $\rho\in\mathcal{I}_{N}$. Let us recall that we have $\mathcal{F}:=L^{3/2}(\mathbb{R}^{3})+L^{\infty}(\mathbb{R}^{3})\subset L^{1}_{\rho}(\mathbb{R}^{3})$. For any $f\in\mathcal{F}$, we denote by

$$\|f\|_{\mathcal{F}}:=\mathop{\inf}_{\begin{array}[]{c}f_{3/2}\in L^{3/2}(\mathbb{R}^{3}),\;f_{\infty}\in L^{\infty}(\mathbb{R}^{3}),\\ f_{3/2}+f_{\infty}=f\\ \end{array}}\|f_{3/2}\|_{L^{3/2}(\mathbb{R}^{3})}+\|f_{\infty}\|_{L^{\infty}(\mathbb{R}^{3})}.$$

Let $M\in\mathbb{N}^{*}$, given a collection of $M$ functions $\Phi:=(\varphi_{1},\ldots,\varphi_{M})\in\mathcal{F}^{M}$, the main idea of the moment-constrained approximation consists in replacing the density constraint in (5) with the $M$ scalar moment constraints associated to the functions $\varphi_{1},\ldots,\varphi_{M}$, that is

$$\int_{\mathbb{R}^{3}}\varphi_{m}\rho_{\Gamma}=\int_{\mathbb{R}^{3}}\varphi_{m}\rho,\quad\forall m=1,\cdots,M.$$ (8)

Notice that (8) is equivalent to

$$\int_{\mathbb{R}^{3}}\varphi\rho_{\Gamma}=\int_{\mathbb{R}^{3}}\varphi\rho,\quad\forall\varphi\in{\rm Span}\{\Phi\}.$$ (9)

We denote by $\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N},\Phi,\rho)$ the set of $\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})$ satisfying constraints (8) (or equivalently (9)).

In the following, we show that there exists at least one solution to the corresponding moment-constrained Lieb optimization problem admits a sparse solution $\Gamma_{\rm opt}^{\Phi}$, such that there exists an integer $K\leq M+2$, weights $\omega_{1},\cdots\omega_{K}\geq 0$ and wavefunctions $\Psi_{1},\cdots,\Psi_{K}\in{\cal H}_{1}^{N}$ such that

$$\sum_{k=1}^{K}\omega_{k}=1\quad\text{and}\quad\Gamma_{\rm opt}^{\Phi}=\sum_{k=1}^{K}\omega_{k}|\Psi_{k}\rangle\langle\Psi_{k}|.$$ (10)

In other words, we will show that there exists a finite-rank minimizer $\Gamma_{\rm opt}^{\Phi}$ the rank of which is at most $K\leq M+2$.

3.1 Tchakaloff’s theorem on Hilbert spaces

Let us first recall the following proposition which is an immediate consequence of Tchakaloff’s theorem, see [BT06]. For any Hilbert space $\mathcal{H}$, we denote by $\mathcal{B}(\mathcal{H})$ the Borel $\sigma$-algebra of $\mathcal{H}$.

Proposition 5.

Let $\mu$ be a Borelian measure on a Hilbert space $\mathcal{H}$ concentrated on a Borel set $\mathcal{A}\in\mathcal{B}(\mathcal{H})$. Let $J_{0}\in{\mathbb{N}}^{*}$ and $\Lambda:\mathcal{H}\to{\mathbb{R}}^{J_{0}}$ a Borel measurable map. Assume that the first moments of $\Lambda_{\sharp}\mu$ exists, that is

$$\int_{{\mathbb{R}}^{J_{0}}}\|x\|\mathrm{d}\Lambda_{\sharp}\mu(x)=\int_{\mathcal{H}}\|\Lambda(\Psi)\|\mathrm{d}\mu(\Psi)<+\infty,$$

where $\|\cdot\|$ denotes the Euclidean norm of ${\mathbb{R}}^{J_{0}}$. Then there exists an integer $1\leq K\leq J_{0}$, elements $\Psi_{1},\cdots,\Psi_{K}\in\mathcal{A}$ and weights $\omega_{1},\cdots,\omega_{K}>0$ such that

$$\forall j=1,\cdots,J_{0},\;\int_{\mathcal{H}}\Lambda_{j}(\Psi)\mathrm{d}\mu(\Psi)=\sum_{k=1}^{K}\omega_{k}\Lambda_{j}(\Psi_{k})=\int_{\mathcal{H}}\Lambda_{i}(\Psi)\,d\mu_{d}(\Psi),$$

where $\Lambda_{j}$ is the $j-$th component of $\Lambda$, and $\displaystyle\mu_{d}=\sum_{k=1}^{K}\omega_{k}\delta_{\Psi_{k}}$.

The main idea of the proof of the sparsity result announced above is to define a measure associated to an operator $\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})$. Assume that the operator $\Gamma$ can be written as

$$\Gamma=\sum_{i=1}^{+\infty}\alpha_{i}|\Psi_{i}\rangle\langle\Psi_{i}|$$ (11)

for some sequence $(\Psi_{i})_{i\in\mathbb{N}^{*}}$ of normalized functions of $\mathcal{H}_{0}^{N}$ and non-negative real numbers $(\alpha_{i})_{i\in\mathbb{N}^{*}}$ such that $\sum_{i\in\mathbb{N}^{*}}\alpha_{i}=N$. Then we can define a Borelian measure $\mu_{\Gamma}:\mathcal{B}(\mathcal{H}_{0}^{N})\to{\mathbb{R}}_{+}$ associated to the decomposition (11) of the operator $\Gamma$ as

$$\mu_{\Gamma}=\sum_{i=1}^{+\infty}\alpha_{i}\delta_{\Psi_{i}}.$$

Naturally, there is no unique such measure $\mu_{\Gamma}$ associated with an operator $\Gamma$ since it heavily depends on the decomposition (11). However, we will see in the following that this is not a problem for our purpose here.

3.2 Existence of sparse minimizers for Moment Constrained Approximation of Lieb (MCAL) functional

In the following, we denote by $\mathds{1}$ the function defined over $\mathbb{R}^{3}$ which is identically equal to $1$.

We then have the following theorem, the proof of which is postponed to Section 7.1.

Theorem 6.

Let $\rho\in\mathcal{I}_{N}$, $M\in\mathbb{N}^{*}$ and $\Phi:=(\varphi_{1},\ldots,\varphi_{M})\in\mathcal{F}^{M}$ such that $\mathds{1}\in{\rm Span}\{\Phi\}$. Let us assume in addition that

• there exists a non-negative non-decreasing continuous function $\theta:{\mathbb{R}}_{+}\to{\mathbb{R}}_{+}$ such that $\displaystyle\theta(r)\mathop{\longrightarrow}_{r\to+\infty}+\infty$ and $C_{\rho}:=\int_{\mathbb{R}^{3}}\theta(|x|)\rho(x)\,dx<+\infty$.

For all $C>0$, let us introduce the Moment-Constrained Approximation of the Lieb functional (MCAL)

$$\boxed{F_{L,\theta}^{\Phi,C}[\rho]:=\inf_{\begin{subarray}{c}\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N},\Phi,\rho)\\ {\rm Tr\,}(\Theta\Gamma)\leq C\end{subarray}}{\rm Tr\,}(H_{N}\Gamma),}$$ (12)

where $\Theta(x_{1},\ldots,x_{N}):=\frac{1}{N}\sum_{i=1}^{N}\theta(|x_{i}|)$ for all $x_{1},\ldots,x_{N}\in\mathbb{R}^{3}$. Then, for all $C\geq C_{\rho}$, $F_{L,\theta}^{\Phi,C}[\rho]$ is finite and a minimum. Moreover, for all $C\geq C_{\rho}$, there exists a minimizer $\Gamma^{\Phi,C}_{{\rm opt},\theta}$ to (12) such that $\Gamma^{\Phi,C}_{{\rm opt},\theta}=\sum_{k=1}^{K}\omega_{k}|\Psi_{k}\rangle\langle\Psi_{k}|$, for some $1\leq K\leq M+1$, with $\omega_{k}\geq 0$ and $\Psi_{k}\in\mathcal{H}_{1}^{N}$ for all $1\leq k\leq K$.

Remark 7.

Let us remark that the existence of a minimizer to a moment-constraint approximation of the Lieb functional has been investigated in [Gar22][Theorem 3.1]. More precisely, in the latter work, the author considers moment functions $(\varphi_{m})_{m\in\mathcal{M}}\subset L^{\infty}(\mathbb{R}^{3},\mathbb{R}_{+})$, where $\mathcal{M}$ is a countable subset of $\mathbb{N}^{*}$, which forms a partition of unity of $\mathbb{R}^{3}$ i.e. such that

$$\sum_{m\in\mathcal{M}}\varphi_{m}=\mathds{1}.$$

In particular, $\mathds{1}\in{\rm Span}\{\varphi_{m},\;m\in\mathcal{M}\}$. Note that in Theorem 6, assumption (A$\theta$) can be seen as an additional condition on $\rho$ which enables to obtain tightness of minimizing sequences. Instead, the author of [Gar22] does not require additional conditions on $\rho$ but considers a tightness condition on the set $(\varphi_{m})_{m\in\mathcal{M}}$ which reads as

$$\mathop{\lim}_{R\to+\infty}\sum_{\begin{array}[]{c}m\in\mathcal{M}\\ ({\rm Supp}\;\varphi_{m})\cap B_{R}^{c}\neq\emptyset\\ \end{array}}\int_{\mathbb{R}^{3}}\rho\varphi_{m}=0,$$ (13)

where for all $R>0$, $B_{R}$ denotes the open ball of $\mathbb{R}^{3}$ of radius $R$ centered at $0$. Note that our existence result, up to the cost of assuming that $\rho$ satisfies (A$\theta$), allows to treat moment constraints for which the tightness condition (13) does not hold. For instance, one can consider a family of moment functions $(\varphi_{m})_{1\leq m\leq M}$ where $(\varphi_{m})_{1\leq m\leq M-1}$ are the characteristic functions of cells of a mesh associated to a bounded subdomain $\Omega\subset\mathbb{R}^{3}$ and $\varphi_{M}=\mathds{1}_{\Omega^{c}}$). It can then be easily checked that such a family does not satisfy condition (13).

Proposition 8 (Lower semi-continuity).

Suppose $\rho_{n}\in\mathcal{I}_{N}$ such that $\rho_{n}\rightharpoonup\rho\in\mathcal{I}_{N}$ in $L^{1}$ then $\liminf F_{L,\theta}^{\Phi,C}[\rho_{n}]=F_{L,\theta}^{\Phi,C}[\rho]$.

Proof.

The proof is a straightforward adaptation of the proof of Theorem 6. Assume that $a_{n}=F_{L,\theta}^{\Phi,C}[\rho_{n}]\to a$ exists then up to the extraction of a subsequence, there exists a trace-class operator $\Gamma_{\infty}\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})$ such that

$$\left((H_{N}+D)^{1/2}\Gamma_{n}(H_{N}+D)^{1/2}\right)_{n\in\mathbb{N}}\leq a_{n}+1/n$$

weakly converges in the sense of trace-class operators to $(H_{N}+D)^{1/2}\Gamma_{\infty}(H_{N}+D)^{1/2}$ as $n$ goes to infinity. Moreover, we have that

$$\liminf{\rm Tr\,}(H_{N}\Gamma_{n})\geq{\rm Tr\,}(H_{N}\Gamma_{\infty}).$$

In particular $\Gamma_{n}$ satisfies the right moment constraints associated to $\rho_{n}$ as well as ${\rm Tr\,}(\Theta\Gamma_{n})\leq C$. Then by using the same arguments as in step 2 of the proof above we deduce that $\Gamma_{\infty}$ is admissible for $F_{L,\theta}^{\Phi,C}[\rho]$. It follows then

$$F_{L,\theta}^{\Phi,C}[\rho]\leq{\rm Tr\,}(H_{N}\Gamma_{\infty})\leq\liminf F_{L,\theta}^{\Phi,C}[\rho_{n}].$$

∎

Remark 9.

We see from the proof of Theorem 6 that assumption (A$\theta$) is needed in order to obtain tightness of the sequence of kernel functions $(\gamma_{n})_{n\in\mathbb{N}}$. This is needed because we are considering operators defined on the space $\mathcal{H}_{0}^{N}=\bigwedge_{i=1}^{N}L^{2}(\mathbb{R}^{3})$. Notice that such a technical assumption is not needed in the case when one considers operators acting on functions acting on a finite domain with Dirichlet boundary conditions. We state such a result below without giving its proof since it follows exactly the same lines as the proof of Theorem 6.

Let $\Omega\subset\mathbb{R}^{3}$ be a bounded subdomain of $\mathbb{R}^{3}$. We then denote by $\mathcal{H}_{0}^{N}(\Omega):=\bigwedge_{i=1}^{N}L^{2}(\Omega)$, $\mathcal{H}_{1}^{N}(\Omega):=\bigwedge_{i=1}^{N}H^{1}_{0}(\Omega)$, $\mathcal{H}_{2}^{N}(\Omega):=\bigwedge_{i=1}^{N}(H^{2}(\Omega)\cap H^{1}_{0}(\Omega))$ and $\mathcal{F}(\Omega):=L^{\infty}(\Omega)+L^{3/2}(\Omega)$. The operator $H_{N,\Omega}:=-\frac{1}{2}\Delta+V$ is then a self-adjoint bounded from below operator acting on $\mathcal{H}_{0}^{N}(\Omega)$ with domain $D(H_{N,\Omega}):=\mathcal{H}_{2}^{N}(\Omega)$. We also denote by $\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega))$ the set of non-negative self-adjoint trace-class operators on $\mathcal{H}_{0}^{N}(\Omega)$. We also define $\mathcal{I}_{N}(\Omega)$ the set of function $\rho\in\mathcal{I}_{N}$ with support included in $\Omega$. For any $M\in\mathbb{N}^{*}$ and any $\Phi:=(\varphi_{m})_{1\leq m\leq M}\subset\mathcal{F}(\Omega)$ and $\rho\in\mathcal{I}_{N}(\Omega)$, we introduce $\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega),\Phi,\rho)$ the set of $\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega))$ such that

$$\int_{\Omega}\rho_{\Gamma}\varphi_{m}=\int_{\Omega}\rho\varphi_{m},\quad\forall 1\leq m\leq M.$$

Then, the following theorem holds:

Theorem 10.

Let $\rho\in\mathcal{I}_{N}(\Omega)$, $M\in\mathbb{N}^{*}$ and $\Phi:=(\varphi_{1},\ldots,\varphi_{M})\in\mathcal{(}\mathcal{F}(\Omega))^{M}$ such that $\mathds{1}|_{\Omega}\in{\rm Span}\{\Phi\}$. Let us introduce

$$\boxed{F_{L,\Omega}^{\Phi}[\rho]:=\inf_{\begin{subarray}{c}\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega),\Phi,\rho)\\ \end{subarray}}{\rm Tr\,}(H_{N,\Omega}\Gamma).}$$ (14)

Then, $F_{L,\Omega}^{\Phi}[\rho]$ is finite and there exists a minimizer $\Gamma^{\Phi}_{{\rm opt},\Omega}$ to (14) such that $\Gamma^{\Phi}_{{\rm opt},\Omega}=\sum_{k=1}^{K}\omega_{k}|\Psi_{k}\rangle\langle\Psi_{k}|$, for some $1\leq K\leq M+1$, with $\omega_{k}>0$ and $\Psi_{k}\in\mathcal{H}_{1}^{N}(\Omega)$ for all $1\leq k\leq K$. Moreover, suppose $\rho_{n}\in\mathcal{I}_{N}$ such that $\rho_{n}\rightharpoonup\rho\in\mathcal{I}_{N}$ in $L^{1}$ then $\liminf F_{L,\Omega}^{\Phi}[\rho_{n}]=F_{L,\Omega}^{\Phi}[\rho]$.

In view of the sparsity results we have just proved, it is natural to consider an approximate MCAL problem, where the set of minimizers is restricted to the set of finite-rank operators satisfying moment constraints. More precisely, for a given $K\in\mathbb{N}^{*}$, we consider the following set

$$\mathcal{O}^{C,\Phi_{,}K}_{\theta}:=\left\{\begin{array}[]{c}({\bm{\omega}},{\bm{\Psi}})\in\mathbb{R}_{+}^{K}\times(\mathcal{H}_{1}^{N})^{K},\quad{\bm{\Psi}}:=(\Psi_{1},\ldots\Psi_{K})\in(\mathcal{H}_{1}^{N})^{K},\\ {\bm{\omega}}:=(\omega_{1},\ldots,\omega_{K})\in\mathbb{R}_{+}^{K},\\ \widetilde{\rho}:=\sum_{k=1}^{K}\omega_{k}\rho_{\Psi_{k}},\quad\int_{\mathbb{R}^{3}}\widetilde{\rho}(x)\theta(|x|)\,dx\leq C,\\ \forall 1\leq m\leq M,\;\int_{\mathbb{R}^{3}}\varphi_{m}\widetilde{\rho}=\int_{\mathbb{R}^{3}}\varphi_{m}\rho\\ \end{array}\right\}.$$

The approximate MCAL functional then reads as follows

$$\boxed{F_{L,\theta}^{\Phi,C,K}[\rho]:=\inf_{({\bm{\Psi}},{\bm{\omega}})\in\mathcal{O}^{C,\Phi_{,}K}_{\theta}}\mathcal{J}({\bm{\Psi}},{\bm{\omega}}),}$$ (15)

where

$$\mathcal{J}({\bm{\Psi}},{\bm{\omega}}):=\sum_{k=1}^{K}\omega_{k}\langle\Psi_{k}|H_{N}|\Psi_{k}\rangle.$$

Remark 11.

Notice that as soon as $K\geq M+1$ then we have that $F_{L,\theta}^{\Phi,C,K}[\rho]=F_{L,\theta}^{\Phi,C}[\rho]$.

Remark 12.

Since $\rho\in\mathcal{I}_{N}$ then the set $\mathcal{O}^{C,\Phi_{,}K}_{\theta}$ is no empty. Moreover it can be shown, by standard arguments, that there exists a minimizer to (15).

As in the case of moment constrained optimal transport [ACE22] we can state some interesting mathematical properties on the set of minimizers of the approximate problem (15). First, consider two elements of $\mathcal{O}^{C,\Phi_{,}K}_{\theta}$, then there exists a continuous path in $\mathcal{O}^{C,\Phi_{,}K}_{\theta}$ connecting these two elements and such that $\mathcal{J}$ varies monotonically along it.

Theorem 13.

Let us assume that $K\geq 2M+2$. Let $({\bm{\Psi}}_{0},{\bm{\omega}}_{0}),({\bm{\Psi}}_{1},{\bm{\omega}}_{1})\in\mathcal{O}^{C,\Phi_{,}K}_{\theta}$. Then, there exists a continuous application $\eta:[0,1]\to\mathcal{O}^{C,\Phi_{,}K}_{\theta}$ made of polygonal chain such that $\eta(0)=({\bm{\Psi}}_{0},{\bm{\omega}}_{0})$, $\eta(1)=({\bm{\Psi}}_{1},{\bm{\omega}}_{1})$ and such that the application $t\mapsto\mathcal{J}(\eta(t))$ is monotone.

Since the proof is a straightforward adaptation of the one for [ACE22][Theorem 1], we refer the reader to it. We only highlight that, as we did in the previous sections, given a couple $({\bm{\Psi}},{\bm{\omega}})$ one can always associate a measure $\mu=\sum_{i}^{K}\omega_{i}\delta_{\psi_{i}}$, then by Thchakaloff’s theorem the result follows. An interesting consequence of theorem 13 concerns the minimizers of MCAL: first, as soon as $K\geq 2M+2$ any local minimizer of MCAL (or of problem (15)) is a global minimizer. Secondly, the set of minimizers forms a polygonally connected set.

Corollary 14.

Assume that $K\geq 2M+2$. Then, any local minimizer of (15) is a global minimizer. Moreover, the set of minimizers of (15) is a polygonally connected subset of $\mathcal{O}^{C,\Phi_{,}K}_{\theta}$.

4 Some convergence results

The aim of this section is to gather some convergence results on the MCAL approximation towards solutions of the exact problem.

4.1 Convergence of the MCAL functional to the exact Lieb functional

The aim of this section is to prove that, under some appropriate assumptions, the MCAL functional converges to the exact Lieb functional as the number of moment constraints go to infinity. Let us denote here by $\mathcal{D}(\mathbb{R}^{3})$ the set of $\mathcal{C}^{\infty}$ real-valued functions defined on $\mathbb{R}^{3}$ with compact support.

More precisely, let $\rho\in\mathcal{I}_{N}$ such that there exists a function $\theta:\mathbb{R}_{+}\to\mathbb{R}_{+}$ satisfying assumption (A$\theta$). Let $C_{\rho}:=\int_{\mathbb{R}^{3}}\theta(|x|)\rho(x)\,dx$ and let $C>C_{\rho}$.

For all $n\in\mathbb{N}^{*}$, let $M_{n}\in\mathbb{N}^{*}$ and $\Phi^{n}:=(\varphi^{n}_{m})_{1\leq m\leq M_{n}}\subset\mathcal{F}$ be a sequence of functions belonging to $\mathcal{F}$ and which satisfies $\mathds{1}\in{\rm Span}\{\Phi^{n}\}$ for all $n\in\mathbb{N}^{*}$ together with the following density conditions:

• for all $f\in\mathcal{D}(\mathbb{R}^{3})$,

  $$\mathop{\inf}_{g_{n}\in{\rm Span}\{\Phi^{n}\}}\|f-g_{n}\|_{\mathcal{F}}\mathop{\longrightarrow}_{n\to+\infty}0.$$

Then, we have the following useful lemma that we will use in the sequel.

Lemma 15.

Let $(\widetilde{\rho}_{n})_{n\in\mathbb{N}^{*}}\subset\mathcal{I}_{N}$ such that $\mathop{\sup}_{n\in\mathbb{N}^{*}}\|\sqrt{\widetilde{\rho}_{n}}\|_{H^{1}(\mathbb{R}^{3})}<+\infty$ and such that for all $n\in\mathbb{N}^{*}$,

$$\forall g_{n}\in{\rm Span}\{\Phi^{n}\},\quad\int_{\mathbb{R}^{3}}\widetilde{\rho}_{n}g_{n}=\int_{\mathbb{R}^{3}}\rho g_{n}.$$

Then, $(\widetilde{\rho}_{n})_{n\in\mathbb{N}^{*}}$ converges in the sense of distributions to $\rho$ as $n$ goes to infinity.

Proof.

The proof uses the same lines as the proof of [Gar22][Theorem 3.2]. We rewrite it here for the sake of completeness. Let $f\in\mathcal{D}(\mathbb{R}^{3})$ and let $(f_{n})_{n\in\mathbb{N}^{*}}$ be a sequence of functions such that $f_{n}\in{\rm Span}\{\Phi^{n}\}$ for all $n\in\mathbb{N}^{*}$ and $\displaystyle\|f-f_{n}\|_{\mathcal{F}}\mathop{\longrightarrow}_{n\to+\infty}0$. Then, it holds that

$$\left|\int_{\mathbb{R}^{3}}f(\widetilde{\rho}_{n}-\rho)\right| =\left|\int_{\mathbb{R}^{3}}(f-f_{n})(\widetilde{\rho}_{n}-\rho)\right|$$

$$\leq C\left(\|\sqrt{\rho}\|_{H^{1}(\mathbb{R}^{3})}^{2}+\mathop{\sup}_{n\in\mathbb{N}^{*}}\|\sqrt{\widetilde{\rho}_{n}}\|_{H^{1}(\mathbb{R}^{3})}^{2}\right)\|f-f_{n}\|_{\mathcal{F}},$$

$$\mathop{\longrightarrow}_{n\to+\infty}0.$$

Hence the desired result. ∎

Remark 16.

One example of sequence $(\Phi_{n})_{n\in\mathbb{N}^{*}}$ satisfying (A$\Phi$) is the following: for all $n\in\mathbb{N}^{*}$, let $\Omega_{n}:=(-n,n)^{3}$ and let $\mathcal{T}_{n}:=\{T_{1}^{n},\ldots,T_{N_{n}}\}$ (with $N_{n}:=\#\mathcal{T}_{n}$) be a regular conforming triangular mesh of $\Omega_{n}$, the elements of which have a maximal diameter size $h_{n}$ such that $h_{n}\leq\frac{1}{n}$. Let $M_{n}:=\#\mathcal{T}_{n}+1=N_{n}+1$. Denoting by $\varphi_{m}^{n}:=\mathds{1}|_{T_{m}^{n}}$ for $1\leq m\leq M_{n}-1$ and by $\varphi_{M_{n}}^{n}:=\mathds{1}|_{\Omega_{n}^{c}}$ and by $\Phi^{n}=(\varphi_{m}^{n})_{1\leq m\leq M_{n}}$ for all $n\in\mathbb{N}^{*}$, one can easily check that the sequence $(\Phi^{n})_{n\in\mathbb{N}^{*}}$ satisfies (A$\Phi$).

We then have the following convergence result, which may be seen as an extension of [Gar22][Theorem 3.2] to more general set of moment functions, up to the additional tightness assumption (A$\theta$), the proof of which is postponed to Section 7.2.

Theorem 17.

Let $\rho\in\mathcal{I}_{N}$ such that there exists a function $\theta:\mathbb{R}_{+}\to\mathbb{R}_{+}$ satisfying assumption (A$\theta$). Let $C_{\rho}:=\int_{\mathbb{R}^{3}}\theta(|x|)\rho(x)\,dx$ and $C\geq C_{\rho}$. For all $n\in\mathbb{N}^{*}$, let $M_{n}\in\mathbb{N}^{*}$ and $\Phi^{n}:=(\varphi_{m}^{n})_{1\leq m\leq M_{n}}\subset\mathcal{F}$ such that assumption (A$\Phi$) holds. We assume in addition that there exists $n_{0}\in\mathbb{N}^{*}$ such that $\mathds{1}\in{\rm Span}\{\Phi^{n}\}$ for all $n\geq n_{0}$. Then, for all $n\geq n_{0}$, there exists at least one sparse minimizer to (12) with $\Phi=\Phi^{n}$ in the sense of Theorem 6. Besides, it holds that

$$\mathop{\lim}_{n\to+\infty}F_{L,\theta}^{\Phi^{n},C}[\rho]=F_{L}[\rho].$$ (16)

Moreover, from any sequence $(\Gamma_{n})_{n\geq n_{0}}$ such that $\Gamma_{n}$ is a minimizer for (12) with $\Phi=\Phi^{n}$, one can extract a subsequence which strongly converges in $\mathfrak{S}_{1,1}(\mathcal{H}_{0}^{N})$ to $\Gamma_{\infty}$, where $\Gamma_{\infty}$ is a minimizer of (5).

Like in Section 3.2, we can state a similar result with less technical assumptions in the case when we consider operators acting on functions defined on a bounded subdomain $\Omega\subset\mathbb{R}^{3}$ with Dirichlet boundary conditions. We state such a result here, using the same notation as in Section 3.2, since it follows exactly the same lines of proof as Theorem 17. To this aim, for all $\rho\in\mathcal{I}_{N}(\Omega)$, we introduce the exact Lieb functional defined on the domain $\Omega$ as

$$F_{L,\Omega}[\rho]:=\inf_{\begin{subarray}{c}\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega))\\ \rho_{\Gamma}=\rho\end{subarray}}{\rm Tr\,}(H_{N,\Omega}\Gamma).$$ (17)

Let us point out here that there exists also $\epsilon_{\Omega},D_{\Omega}>0$ such that

$$H_{N,\Omega}+D_{\Omega}\geq\varepsilon_{\Omega}(-\Delta_{\Omega}+1)$$

where $-\Delta_{\Omega}$ refers here to the self-adjoint bounded from below operator on $\mathcal{H}_{0}^{N}(\Omega)$ with domain $\mathcal{H}_{N}^{2}(\Omega)$ (Laplacian with Dirichlet boundary conditions in $\Omega$). We also denote by $\mathfrak{S}_{1,1}\left(\mathcal{H}_{0}^{N}(\Omega)\right)$ the set of operators $\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega))$ such that ${\rm Tr}(-\Delta_{\Omega}\Gamma)<+\infty$.

Theorem 18.

Let $\rho\in\mathcal{I}_{N}(\Omega)$. For all $n\in\mathbb{N}^{*}$, let $M_{n}\in\mathbb{N}^{*}$ and $\Phi^{n}:=(\varphi_{m}^{n})_{1\leq m\leq M_{n}}\subset\mathcal{F}(\Omega)$ such that for all $f\in\mathcal{D}(\Omega)$,

$$\lim_{n\to+\infty}\mathop{\inf}_{g_{n}\in{\rm Span}\{\Phi^{n}\}}\|f-g_{n}\|_{\mathcal{F}(\Omega)}=0.$$

We assume in addition that there exists $n_{0}\in\mathbb{N}^{*}$ such that $\mathds{1}\in{\rm Span}\{\Phi^{n}\}$ for all $n\geq n_{0}$. Then, for all $n\geq n_{0}$, there exists at least one sparse minimizer to (12) with $\Phi=\Phi^{n}$ in the sense of Theorem 6. Besides, it holds that

$$\mathop{\lim}_{n\to+\infty}F_{L,\Omega}^{\Phi^{n}}[\rho]=F_{L,\Omega}[\rho].$$ (18)

Moreover, from any sequence $(\Gamma_{n})_{n\geq n_{0}}$ such that $\Gamma_{n}$ is a minimizer for (12) with $\Phi=\Phi^{n}$, one can extract a subsequence which strongly converges in $\mathfrak{S}_{1,1}(\mathcal{H}_{0}^{N}(\Omega))$ to $\Gamma_{\infty}$, where $\Gamma_{\infty}$ is a minimizer of (5).

4.2 Convergence rate of the ground state energy in the bounded domain case

In this section, we restrict ourselves to the case of a bounded subdomain $\Omega\subset\mathbb{R}^{3}$. Let $M\in\mathbb{N}^{*}$, $\Phi:=(\varphi_{m})_{1\leq m\leq M}\subset\mathcal{F}(\Omega)$ be a set of moment functions. For all $v\in\mathcal{F}(\Omega)$, let us introduce the ground state energy associated to the potential $v$:

$$E[v]:=\inf_{\Psi\in\mathcal{H}_{1}^{N}(\Omega)}\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle=\inf_{\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega))}{\rm Tr\,}(H_{N,\Omega}^{v}\Gamma),$$

where

$$H_{N,\Omega}^{v}:=H_{N,\Omega}-\sum_{i=1}^{N}v(x_{i}).$$

Rewriting the minimization over $\Gamma$ as an external minimization over $\rho\in\mathcal{I}_{N}(\Omega)$ and then as an internal one over all $\Gamma$ such that ${\rm Tr}\;\Gamma=\rho$, it can easily be checked that

$$E[v]=\mathop{\inf}_{\rho\in\mathcal{I}_{N}(\Omega)}\left\{F_{L}[\rho]-\int_{\Omega}v\,d\rho\right\}.$$ (19)

Let us also define by

$$E^{\Phi}[v]:=\mathop{\inf}_{\rho\in\mathcal{I}_{N}(\Omega)}\left\{F^{\Phi}_{L}[\rho]-\int_{\Omega}v\,d\rho\right\}.$$ (20)

Similarly, let us point out that, if $v\in{\rm Span}\{\Phi\}$, rewriting the minimization over $\Gamma$ as an external minimization over $\rho\in\mathcal{I}_{N}(\Omega)$ and then as an internal one over all $\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega),\Phi,\rho)$, it holds that

$$E[v]=E^{\Phi}[v],\quad\forall v\in{\rm Span}\{\Phi\}.$$

We then prove the following approximation result.

Proposition 19.

Let us assume that $v\in L^{\infty}(\Omega)$ and that $\Phi=(\varphi_{m})_{1\leq m\leq M}\subset L^{\infty}(\Omega)$. Then, it holds that

$$|E[v]-E^{\Phi}[v]|\leq 2N\mathop{\min}_{w\in{\rm Span}\{\Phi\}}\|v-w\|_{L^{\infty}(\Omega)}.$$ (21)

Proof.

Let $\displaystyle v^{\Phi}=\mathop{\rm argmin}_{w\in{\rm Span}\{\Phi\}}\|v-w\|_{L^{\infty}(\Omega)}$. Let $\varepsilon>0$ arbitrarily small. Let $\rho$, $\rho^{\Phi}$, $\widetilde{\rho}^{\Phi}$ and $\overline{\rho}^{\Phi}$ be $\varepsilon$-minimizers of $E[v]$, $E[v^{\Phi}]$, $E^{\Phi}[v]$ and $E^{\Phi}[v^{\Phi}]$ respectively. It then holds that

$$E[v^{\Phi}] \leq F_{L}[\rho^{\Phi}]-\int_{\Omega}v^{\Phi}\,d\rho^{\Phi}$$

$$\leq E[v_{\Phi}]+\varepsilon$$

$$\leq F_{L}[\rho]-\int_{\Omega}v^{\Phi}\,d\rho+\varepsilon$$

$$=F_{L}[\rho]-\int_{\Omega}v\,d\rho+\int_{\Omega}(v^{\Phi}-v)\,d\rho+\varepsilon$$

$$\leq E[v]+\int_{\Omega}(v^{\Phi}-v)\,d\rho+2\varepsilon.$$

Using similar calculations, we obtain that

$$E[v]\leq E[v^{\Phi}]+\int_{\Omega}(v-v^{\Phi})\,d\rho^{\Phi}+2\varepsilon.$$

As a consequence, we obtain that

$$|E[v]-E[v^{\Phi}]|\leq\max\left(\int_{\Omega}|v-v^{\Phi}|\,d\rho,\int_{\Omega}|v-v^{\Phi}|\,d\rho^{\Phi}\right)+2\varepsilon\leq N\|v-v^{\Phi}\|_{L^{\infty}(\Omega)}+2\varepsilon.$$

Since $\varepsilon$ can be chosen arbitrarily small, it actually holds that

$$|E[v]-E[v^{\Phi}]|\leq N\|v-v^{\Phi}\|_{L^{\infty}(\Omega)}.$$ (22)

Using similar arguments, we also obtain that

$$|E^{\Phi}[v]-E^{\Phi}[v^{\Phi}]|\leq N\|v-v^{\Phi}\|_{L^{\infty}(\Omega)}.$$ (23)

Collecting (22) and (23) and using the fact that $E[v^{\Phi}]=E^{\Phi}[v^{\Phi}]$ yields the desired result. ∎

Proposition 19 then enables to quantify the rate of convergence of $|E[v]-E^{\Phi^{n}}[v]|$ as $n$ goes to infinity for some particular sequences of moment functions $(\Phi^{n})_{n\in\mathbb{N}}$ provided that $v$ is regular enough. As an illustration, we analyze here the rate of convergence of a numerical method inspired from the external dual charge approach recently proposed in [Lel22].

Corollary 20.

Let $l\geq 0$ and $\Omega$ be a bounded regular subdomain of $\mathbb{R}^{3}$. Let $\mu\in H^{l+1}(\Omega)$ be an external density of charge and define $v\in H^{1}_{0}(\Omega)\cap H^{l+3}(\Omega)$ as the unique solution to

$$\left\{\begin{array}[]{ll}-\Delta v=\mu&\quad\mbox{in }\Omega,\\ v=0&\quad\mbox{on }\partial\Omega.\\ \end{array}\right.$$

Let $(\mathcal{T}_{h})_{h>0}$ be a sequence of triangular regular meshes of $\Omega$ such that

$$h:=\max_{K\in\mathcal{T}_{h}}{\rm diam}(K).$$

Let $k\in\mathbb{N}$ and $P_{h}^{k}\subset L^{\infty}(\Omega)$ be the subspace of continuous $\mathbb{P}_{k}$ finite element functions associated to the mesh $\mathcal{T}_{h}$. We denote by $V_{h,k}$ the subspace of $H^{1}_{0}(\Omega)\cap H^{2}(\Omega)$ containing all functions $v_{h,k}\in H^{1}_{0}(\Omega)\cap H^{2}(\Omega)$ solution to

$$\left\{\begin{array}[]{ll}-\Delta v_{h,k}=\mu_{h,k}&\quad\mbox{in }\Omega,\\ v_{h,k}=0&\quad\mbox{on }\partial\Omega,\\ \end{array}\right.$$

for some $\mu_{h,k}\in P_{h}^{k}$. Let $\Phi_{h,k}$ be a basis of $V_{h,k}$. Then, asuming that $l\leq k$, there exists a constant $C>0$ such that for all $h>0$,

$$\boxed{|E[v]-E^{\Phi_{h,k}}[v]|\leq CNh^{l+1}\|v\|_{H^{l+3}(\Omega)}.}$$

Proof.

Corollary 20 easily follows for the compact embedding $H^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$ and standard interpolation error results associated with finite element approximations. ∎

Remark 21.

Denoting by $M_{h,k}$ the dimension of $V_{h,k}$, it holds that $M_{h,k}=\mathcal{O}\left(\frac{k}{h^{3}}\right)$. As a consequence, the above result implies that the rate of convergence of $E^{\Phi_{h,k}}[v]$ to $E[v]$ decays like $\mathcal{O}\left(\frac{N}{M_{h,k}^{(l+1)/3}}\right)$ where $M_{h,k}$ is the number of moment constraints in the MCAL approximation.

5 Duality results for the MCAL functional

Let us begin by recalling some classical results about semi-definite programming problems and introduce some notation.

5.1 Semi-definite positive programming problems

Let $n\in\mathbb{N}^{*}$. We denote by $\mathcal{S}^{n}$ the set of symmetric matrices of $\mathbb{R}^{n}$. For any $M\in\mathcal{S}^{n}$, the notation $M\succcurlyeq 0$ (respectively $M\succ 0$) is used to mean that $M$ is a semi-definite non-negative (respectively definite positive) matrix. We also denote by $\mathcal{S}^{n}_{+}:=\{M\in\mathcal{S}^{n},\;M\succcurlyeq 0\}$ and by $\mathcal{S}^{n}_{+,*}:=\{M\in\mathcal{S}^{n},\;M\succ 0\}$. For all $M,N\in\mathcal{S}^{n}$, we denote by $\langle M,N\rangle={\rm Tr}(M^{T}N)$ the Frobenius scalar product between $M$ and $N$.

Let $m\in\mathbb{N}^{*}$, $C\in\mathcal{S}^{n}$, $A:\mathcal{S}^{n}\to\mathbb{R}^{m}$ a linear application and $b\in\mathbb{R}^{m}$. We consider here the following (primal) semi-definite positive programming problem:

$$\boxed{P:=\mathop{\inf}_{\begin{array}[]{c}X\in\mathcal{S}^{n}\\ A(X)=b\\ X\succcurlyeq 0\\ \end{array}}\langle C,X\rangle.}$$ (24)

The dual problem associated to (24) then reads as follows:

$$\boxed{D:=\mathop{\sup}_{\begin{array}[]{c}(y,S)\in\mathbb{R}^{m}\times\mathcal{S}^{n}\\ A^{*}(y)+S=C\\ S\succcurlyeq 0\\ \end{array}}\langle b,y\rangle}$$ (25)

where $A^{*}:\mathbb{R}^{m}\to\mathcal{S}^{n}$ is the adjoint of $A$.

We introduce the following sets:

$$\mathcal{A}_{P} :=\left\{X\in\mathcal{S}^{n},\;A(X)=b,\;X\succcurlyeq 0\right\},$$

$$\mathcal{A}^{s}_{P} :=\left\{X\in\mathcal{S}^{n},\;A(X)=b,\;X\succ 0\right\},$$

$$\mathcal{A}_{D} :=\left\{(y,S)\in\mathbb{R}^{m}\times\mathcal{S}^{n},\;A^{*}(y)+S=C,\;S\succcurlyeq 0\right\},$$

$$\mathcal{A}^{s}_{D} :=\left\{(y,S)\in\mathbb{R}^{m}\times\mathcal{S}^{n},\;A^{*}(y)+S=C,\;S\succ 0\right\}.$$

We also denote by ${\rm Sol}_{P}$ and ${\rm Sol}_{D}$ the set of solutions to (24) and (25). Then, we recall the following classical result [AL11, WSV12]:

Theorem 22.

• If $\mathcal{A}_{P}\times\mathcal{A}^{s}_{D}\neq\emptyset$, ${\rm Sol}_{P}$ is non-empty and bounded and $P=D$;

• If $\mathcal{A}^{s}_{P}\times\mathcal{A}_{D}\neq\emptyset$ and $A$ surjective, then ${\rm Sol}_{D}$ is non-empty and bounded and $P=D$;

• If $\mathcal{A}^{s}_{P}\times\mathcal{A}^{s}_{D}\neq\emptyset$ and $A$ surjective, then ${\rm Sol}_{P}$ and ${\rm Sol}_{D}$ are non-empty and bounded and $P=D$.

5.2 Dual MCAL problem

In this section we study the dual problem in the bounded domain case. We know that the dual variable associated to the density $\rho\in\mathcal{I}_{N}(\Omega)$ is a one-body interaction potential of the form $W^{v}(x_{1},\ldots,x_{N}):=\sum_{i=1}^{N}v(x_{i})$ for a given $v\in\mathcal{F}(\Omega)$.

We then consider the following natural dual problem

$$\boxed{D_{L,\Omega}^{\Phi}[\rho]=\sup_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\},\\ \forall\Psi\in\mathcal{H}_{1}^{N}(\Omega),\;\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0\\ \end{array}}\int_{\Omega}vd\rho,}$$ (26)

where

$$H_{N,\Omega}^{v}:=H_{N,\Omega}-\sum_{i=1}^{N}v(x_{i})=H_{N,\Omega}-W^{v}.$$

If we take any $v:=\sum_{m=1}^{M}\alpha_{m}\varphi_{m}\in{\rm Span}\{\Phi\}$ satisfying the above constraints and any $\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega),\Phi,\rho)$ then we have

$$\begin{split}{\rm Tr\,}(H_{N,\Omega}\Gamma)\geq{\rm Tr\,}(W^{v}\Gamma)&=\int_{\Omega}vd\rho_{\Gamma}=\int_{\Omega}\bigg{(}\sum_{m=1}^{M}\alpha_{m}\varphi_{m}\bigg{)}d\rho_{\Gamma}\\ &\geq\int_{\Omega}\bigg{(}\sum_{m=1}^{M}\alpha_{m}\varphi_{m})\bigg{)}d\rho=\int_{\Omega}vd\rho\end{split}$$

which proves that $F_{L,\Omega}^{\Phi}[\rho]\geq D_{L,\Omega}^{\Phi,C}[\rho]$. We would like to prove that this inequality is actually an equality. Let us introduce the ground state energy associated to the potential $v$:

$$E[v]=\inf_{\Psi\in\mathcal{H}_{1}^{N}(\Omega)}\langle\Psi|H_{N}^{v}|\Psi\rangle=\inf_{\Gamma\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega))}{\rm Tr\,}(H_{N}^{v}\Gamma).$$

We rewrite now the minimization over $\Gamma$ as an external minimization over $\rho\in\mathcal{I}_{N}(\Omega)$ and then as an internal one over all $\Gamma$ in $\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N}(\Omega),\Phi,\rho)$ (we are considering the ground state for a potential $v\in{\rm Span}\{\Phi\}$):

$$E[v]=\inf_{\rho\in\mathcal{I}_{N}(\Omega)}\left\{F_{L,\Omega}^{\Phi}[\rho]-\int_{\Omega}v\mathrm{d}\rho\right\}.$$

Notice that $E$ is nothing but the Legendre-Fenchel transform of $F_{L,\Omega}^{\Phi}[\rho]$. On the other hand, we rewrite (26) in the form

$$D_{L,\Omega}^{\Phi}[\rho]=\sup_{v\in{\rm Span}\{\Phi\}}\left\{\int_{\Omega}v\mathrm{d}\rho-E[v]\right\}.$$ (27)

Thus, $D_{L,\Omega}^{\Phi}[\rho]$ is the Legendre transform of $E$. From Proposition 8 and Fenchel duality theorem for convex lower semi-continuous functions we conclude the following

Theorem 23.

Under the assumptions of Theorem 10, we have $F_{L,\Omega}^{\Phi}[\rho]=D_{L,\Omega}^{\Phi}[\rho]$.

We now have the following result which, taking into account the sparsity result of Theorem 10, gives a more convenient formulation of $D_{L,\Omega}^{\Phi}[\rho]$.

Theorem 24.

Under the assumptions of Theorem 10, there exists at least one maximizer to (26), and it holds that

$$D_{L,\Omega}^{\Phi}[\rho] =\mathop{\max}_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\},\\ \forall\Psi\in\mathcal{H}_{1}^{N}(\Omega),\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0\\ \end{array}}\int_{\Omega}v\rho$$

$$=\mathop{\max}_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\},\\ \forall\Psi\in{\rm Span}\{\Psi_{1},\ldots,\Psi_{K}\},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0\\ \end{array}}\int_{\Omega}v\rho,$$

where

$$\Gamma^{\Phi}_{{\rm opt},\Omega}=\sum_{k=1}^{K}\omega_{k}|\Psi_{k}\rangle\langle\Psi_{k}|$$

for some $1\leq K\leq M+1$, with $\omega_{k}>0$ and $\Psi_{k}\in\mathcal{H}_{1}^{N}(\Omega)$ for all $1\leq k\leq K$ is a minimizer of (14).

6 Numerical scheme

The aim of this section is to propose a new numerical scheme using the sparsity of minimizers of the MCAL functional to compute approximations of the Lieb functional. The scheme proposed here requires the resolution of eigenvalue problems for operators acting on $\mathcal{H}_{0}^{N}(\Omega)$, which leads to high-dimensional problems when the number of electrons is large. The combination of the algorithm proposed here with numerical methods dedicated to overcome the curse of dimensionality will be the object of a future work.

We propose here an iterative scheme which shares some common features with the well-known Column Algorithm used for classical optimal transport problems (see for instance [FP22, FSV22]). The aim is to construct at each iteration $n\in\mathbb{N}^{*}$ a finite set of $L^{2}$-normalized wavefunctions $\mathfrak{P}_{n}\subset\mathcal{H}^{2}_{N}(\Omega)$ which will be used to enforce the inequality constraints in the resolution of the MCAL dual problems. More precisely, inequality constraints in small-dimensional dual problems are enforced to hold on the space spanned by the wavefunctions belonging to the set $\mathfrak{P}_{n}$. As a consequence, in our present quantum optimal transprt framework, semi-definite programming problems have to be solved at each iteration instead of linear programming problems for classical optimal transport problems.

6.1 MCAL iterative scheme

We describe the MCAL algorithm in this section. The algorithm takes as input data:

• $\Phi=(\varphi_{1},\ldots,\varphi_{M})\subset L^{\infty}(\Omega)$ set of moment functions;

• $\rho\in\mathcal{I}_{N}$ with support included in $\Omega$;

• $\rho^{m}:=\int_{\Omega}\varphi_{m}\rho$ for all $1\leq m\leq M$.

• $\widetilde{\mathfrak{P}}_{0}\subset\mathcal{H}^{2}_{N}(\Omega)$ initial finite set of $L^{2}$-normalized wavefunctions.

As an output, after $n$ iterations, the algorithm yields $F^{n}$ which is an approximation of the quantity $F_{L,\Omega}^{\Phi}[\rho]$.

We make here the following assumption on the initial set $\widetilde{\mathfrak{P}}_{0}$.

Assumption (A0): let $\widetilde{K}_{0}:={\rm dim}\;{\rm Span}\{\widetilde{\mathfrak{P}}_{0}\}$ and $(\widetilde{\Psi}_{1}^{0},\ldots,\widetilde{\Psi}_{K_{0}}^{0})$ an orthonormal basis of ${\rm Span}\{\widetilde{\mathfrak{P}}_{0}\}$. We assume that there exists $\widetilde{S}:=(\widetilde{S}_{kl})_{1\leq k,l\leq K_{0}}\in\mathcal{S}_{+}^{K_{0}}$ such that for all $1\leq m\leq M$,

$$\sum_{k,l=1}^{K_{0}}\widetilde{S}_{kl}\int_{\Omega}\varphi_{m}\overline{\widetilde{\Psi}_{k}^{0}}\widetilde{\Psi}_{l}^{0}=\rho^{m}.$$

In the case when $d=3$, a way to find such an initial set $\widetilde{\mathfrak{P}}_{0}$ is given in [Lie83b]. In this case, $\widetilde{K}_{0}$ can be chosen to be equal to $1$ and $\widetilde{\Psi}_{1}^{0}$ can be chosen as follows: for all $1\leq k\leq N$ and $x=(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}$, define

$$\phi^{k}(x)=\sqrt{\frac{\rho(x)}{N}}e^{ikf(x_{1})},$$

where for all $x_{1}\in\mathbb{R}$,

$$f(x_{1})=\left(\frac{2\pi}{N}\right)\int_{-\infty}^{x_{1}}\,ds\int_{-\infty}^{+\infty}\,dt\int_{-\infty}^{\infty}\,du\rho(s,t,u).$$

The family $(\phi^{k})_{1\leq k\leq N}$ then forms an orthonormal family of $L^{2}(\mathbb{R}^{3})$ and one may define $\widetilde{\Psi}_{1}^{0}$ as the normalized Slater determinant associated to the family $(\phi^{k})_{1\leq k\leq N}$.

6.1.1 Initialization step

Compute $\widetilde{S}^{0}:=(\widetilde{S}^{0}_{kl})_{1\leq k,l\leq K_{0}}\in\mathcal{S}_{+}^{K_{0}}$ solution to

$$\widetilde{F}^{0}:=\mathop{\min}_{\begin{array}[]{c}(S_{kl})_{1\leq k,l\leq K_{n}}\in\mathcal{S}_{+}^{K_{0}},\\ \forall 1\leq m\leq M,\\ \sum_{k,l=1}^{K_{0}}S_{kl}\int_{\Omega}\varphi_{m}\overline{\widetilde{\Psi}_{k}^{0}}\widetilde{\Psi}_{l}^{0}=\rho^{m}\\ \end{array}}\sum_{k,l=1}^{K_{0}}S_{kl}\langle\widetilde{\Psi}_{k}^{0}|H_{N,\Omega}|\widetilde{\Psi}_{l}^{0}\rangle$$ (28)

Then, it holds that $\widetilde{S}^{0}=\sum_{k=1}^{K_{0}}\omega_{k}^{0}(U_{k}^{0})(U_{k}^{0})^{T}$ where $(\omega_{k}^{0})_{1\leq k\leq K_{0}}\in\mathbb{R_{+}}^{K_{0}}$ are the eigenvalues of $\widetilde{S}^{0}$ (assumed to be ranked in non-increasing order) and for all $1\leq k\leq K_{0}$, $U_{k}^{0}:=(U_{kl}^{0})_{1\leq l\leq K_{0}}\in\mathbb{R}^{K_{0}}$ is a normalized eigenvector associated with $\omega_{k}^{0}$ so that $(U_{1}^{0},\ldots,U_{K_{0}}^{0})$ forms an orthonormal basis of $\mathbb{R}^{K_{0}}$.

Let $\mathcal{K}_{0}:=\max\left\{k\in\{1,\ldots,K_{0}\},\omega_{k}^{0}>0\right\}$. For all $1\leq k\leq\mathcal{K}_{0}$, let $\Psi_{k}^{0}:=\sum_{l=1}^{K_{0}}U_{kl}^{0}\widetilde{\Psi}_{l}^{0},$ and $S^{0}:={\rm diag}(\omega_{1}^{0},\ldots,\omega_{\mathcal{K}_{0}})\in\mathcal{S}_{+,*}^{\mathcal{K}_{0}}$. We also denote by $\mathfrak{P}_{0}:=\bigcup_{1\leq k\leq\mathcal{K}_{0}}\left\{\Psi_{k}^{0}\right\}$.

Remark 25.

It is easy to see that, by construction, it holds that

$$\widetilde{F}^{0}:=\mathop{\min}_{\begin{array}[]{c}(S_{kl})_{1\leq k,l\leq\mathcal{K}_{0}}\in\mathcal{S}_{+}^{\mathcal{K}_{0}},\\ \forall 1\leq m\leq M,\\ \sum_{k,l=1}^{\mathcal{K}_{0}}S_{kl}\int_{\Omega}\varphi_{m}\overline{\Psi_{k}^{0}}\Psi_{l}^{0}=\rho^{m}\\ \end{array}}\sum_{k,l=1}^{K_{0}}S_{kl}\langle\Psi_{k}^{0}|H_{N,\Omega}|\Psi_{l}^{0}\rangle$$ (29)

and that $S^{0}$ is a minimizer to (29).

Remark 26.

Notice also that this initialization step is useless in the case when $\widetilde{K}^{0}=1$.

6.1.2 Iteration $n\geq 1$

Step 1: Let $\mathcal{K}^{n-1}:={\rm dim}\;{\rm Span}\left\{\mathfrak{P}_{n-1}\right\}$ and $(\Psi^{n-1}_{1},\ldots,\Psi^{n-1}_{K_{n-1}})$ be an orthonormal basis of ${\rm Span}\left\{\mathfrak{P}_{n-1}\right\}$. Let $A^{n-1}:=(A^{n-1}_{kl,m})_{1\leq m\leq M,1\leq k,l\leq K_{n-1}}\in\mathbb{R}^{K_{n-1}^{2}\times M}$ be defined by

$$A^{n-1}_{kl,m}:=\int_{\Omega}\varphi_{m}\overline{\Psi^{n-1}_{k}}\Psi^{n-1}_{l}.$$

Let $C^{n-1}:={\rm Ker}(A^{n-1})^{\perp}\subset\mathbb{R}^{M}$ and

$$V^{n-1}:=\left\{v=\sum_{m=1}^{M}c_{m}\varphi_{m},\;c:=(c_{m})_{1\leq m\leq M}\in C^{n-1}\right\}\subset{\rm Span}\{\Phi\}.$$

Compute $v^{n}\in V^{n-1}$ solution to

$$F^{n}=\mathop{\max}_{\begin{array}[]{c}v\in V^{n-1}\\ \forall\Psi\in{\rm Span}\{\mathfrak{P}_{n-1}\},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0.\end{array}}\int_{\mathbb{R}^{3}}v\rho$$ (30)

Remark 27.

Using the results of semi-definite positive programming and using similar arguments as in the proof of Theorem 24, it can be easily checked that there exists at least one maximizer to (30). In addition, any maximizer to (30) is also a maximizer to

$$F^{n}=\mathop{\max}_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\}\\ \forall\Psi\in{\rm Span}\{\mathfrak{P}_{n-1}\},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0.\end{array}}\int_{\Omega}v\rho,$$

since $\int_{\Omega}v\rho=0$ for any $v=\sum_{m=1}^{M}c_{m}\varphi_{m}$ with $c:=(c_{m})_{1\leq m\leq M}\in{\rm Ker}(A^{n-1})$.

Step 2: Compute $\Psi_{0}^{v_{n}}\in\mathcal{H}_{2}^{N}(\Omega)$ a $L^{2}$-normalized solution to

$$H^{v_{n}}_{N,\Omega}\Psi_{0}^{v_{n}}=E(v_{n})\Psi_{0}^{v_{n}},$$ (31)

where $E(v_{n})$ is the smallest eigenvalue of $H^{v_{n}}_{N,\Omega}$.

Step 3: We now distinguish two different cases.

• Case 1: $E(v_{n})<0$

  Define $\widetilde{\mathfrak{P}}_{n}:=\mathfrak{P}_{n-1}\cup\{\Psi_{0}^{v_{n}}\}$. Let $K_{n}:={\rm dim}\;{\rm Span}\{\widetilde{{\mathfrak{P}}}_{n}\}$ and let $\widetilde{\Psi}_{1}^{n},\ldots,\widetilde{\Psi}_{K_{n}}^{n}$ be an orthonormal basis of ${\rm Span}\{\widetilde{{\mathfrak{P}}}_{n}\}$.

  Compute $\widetilde{S}^{n}:=(\widetilde{S}^{n}_{kl})_{1\leq k,l\leq K_{N}}\in\mathcal{S}_{+}^{K_{n}}$ solution to

  $$\widetilde{F}^{n}:=\mathop{\min}_{\begin{array}[]{c}(S_{kl})_{1\leq k,l\leq K_{n}}\in\mathcal{S}_{+}^{K_{n}},\\ \forall 1\leq m\leq M,\\ \sum_{k,l=1}^{K_{n}}S_{kl}\int_{\Omega}\varphi_{m}\overline{\widetilde{\Psi}_{k}^{n}}\widetilde{\Psi}_{l}^{n}=\rho^{m}\\ \end{array}}\sum_{k,l=1}^{K_{n}}S_{kl}\langle\widetilde{\Psi}_{k}^{n}|H_{N,\Omega}|\widetilde{\Psi}_{l}^{n}\rangle$$ (32)

  Then, it holds that $\widetilde{S}^{n}=\sum_{k=1}^{K_{n}}\omega_{k}^{n}(U_{k}^{n})(U_{k}^{n})^{T}$ where $(\omega_{k}^{n})_{1\leq k\leq K_{n}}\in\mathbb{R_{+}}^{K_{n}}$ are the eigenvalues of $\widetilde{S}^{n}$ (assumed to be ranked in non-increasing order) and for all $1\leq k\leq K_{n}$, $U_{k}^{n}:=(U_{kl}^{n})_{1\leq l\leq K_{n}}\in\mathbb{R}^{K_{n}}$ is a normalized eigenvector associated with $\omega_{k}^{n}$ so that $(U_{1}^{n},\ldots,U_{K_{n}}^{n})$ forms an orthonormal basis of $\mathbb{R}^{K_{n}}$.

  Let $\mathcal{K}_{n}:=\max\left\{k\in\{1,\ldots,K_{n}\},\omega_{k}^{n}>0\right\}$. For all $1\leq k\leq\mathcal{K}_{n}$, let $\Psi_{k}^{n}:=\sum_{l=1}^{K_{n}}U_{kl}^{n}\widetilde{\Psi}_{l}^{n},$ and $S^{n}:={\rm diag}(\omega_{1}^{n},\ldots,\omega_{\mathcal{K}_{n}}^{n})\in\mathcal{S}_{+,*}^{\mathcal{K}_{n}}$. We then denote by $\mathfrak{P}_{n}:=\left\{\Psi_{1}^{n},\ldots,\Psi_{{\mathcal{K}}_{n}}^{n}\right\}$.

  Define $n:=n+1$ and proceed with the next iteration.

• Case 2: $E(v_{n})\geq 0$

  Stop the algorithm.

6.2 Property of the MCAL iterative scheme

We prove the following lemma, which states that the sequence of approximations yielded by the MCAL algorithm is non-increasing. Note however that we do not prove here that the sequence converges indeed to $F_{L,\Omega}^{\Phi}[\rho]$.

Lemma 28.

For all $n\geq 1$, it holds that

$$F^{n}\geq\widetilde{F}^{n}=F^{n+1}\geq F_{L,\Omega}^{\Phi}[\rho].$$

Proof.

The first inequality is simple to see since the dual problem associated to (32) is

$$\widetilde{F}_{n} =\mathop{\max}_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\}\\ \forall\Psi\in{\rm Span}\left\{\widetilde{\mathfrak{P}}_{n}\right\},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0.\end{array}}\int_{\mathbb{R}^{3}}v\rho$$

$$\leq\mathop{\max}_{\begin{array}[]{c}\scriptsize v\in{\rm Span}\{\Phi\}\\ \scriptsize\forall\Psi\in{\rm Span}\left\{\mathfrak{P}_{n-1}\right\},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0.\end{array}}\int_{\mathbb{R}^{3}}v\rho$$

$$=F^{n},$$

since $\mathfrak{P}_{n-1}\subset\widetilde{\mathfrak{P}}_{n}$.

The second equality comes from the fact that

$$\widetilde{F}_{n} =\mathop{\min}_{\begin{array}[]{c}(S_{kl})_{1\leq k,l\leq K_{n}}\in\mathcal{S}_{+}^{K_{n}},\\ \forall 1\leq m\leq M,\\ \sum_{k,l=1}^{K_{n}}S_{kl}\int_{\Omega}\varphi_{m}\overline{\widetilde{\Psi}_{k}^{n}}\widetilde{\Psi}_{l}^{n}=\rho^{m}\\ \end{array}}\sum_{k,l=1}^{K_{n}}S_{kl}\langle\widetilde{\Psi}_{k}^{n}|H_{N,\Omega}|\widetilde{\Psi}_{l}^{n}\rangle$$

$$=\mathop{\min}_{\begin{array}[]{c}(S_{kl})_{1\leq k,l\leq{\mathcal{K}}_{n}}\in\mathcal{S}_{+}^{{\mathcal{K}}_{n}},\\ \forall 1\leq m\leq M,\\ \sum_{k,l=1}^{{\mathcal{K}}_{n}}S_{kl}\int_{\Omega}\varphi_{m}\overline{\Psi_{k}^{n}}\Psi_{l}^{n}=\rho^{m}\\ \end{array}}\sum_{k,l=1}^{{\mathcal{K}}_{n}}S_{kl}\langle\Psi_{k}^{n}|H_{N,\Omega}|\Psi_{l}^{n}\rangle.$$

In addition, we know, by definition of ${\mathcal{K}}_{n}$ and of $\Psi_{1}^{n}$, …, $\Psi_{{\mathcal{K}}_{n}}^{n}$ that there exists at least one minimizer to the second minimization problem which is a positive definite matrix, that is the diagonal matrix with entries $\omega_{1}^{n},\ldots,\omega_{{\mathcal{K}}_{n}}$. Using standard results of semi-definite positive programming, it holds that the dual problem associated to the second minimization problem introduced in the last line of the calculations above is precisely

$$F^{n+1}=\mathop{\max}_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\}\\ \forall\Psi\in{\rm Span}\left\{\mathfrak{P}_{n}\right\},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0.\end{array}}\int_{\mathbb{R}^{3}}v\rho=\widetilde{F}_{n}.$$

Hence the desired result. ∎

6.3 Numerical results

The numerical tests presented in this section were performed using Julia. In particular, the finite element code developped in [QC23] was used to solve the eigenvalue problems (31), and the ProxSDP library was used for the resolution of the semi-definite programming problems. The associated code can be found on ZENODO with the DOI 10.5281/zenodo.11669900.

We present in this section some preliminary numerical results on a toy numerical test case with $N=2$, $d=1$ and $\Omega=(-L,L)$ with $L=10$. More precisely, for a given value $D\in\mathbb{N}^{*}$, the solution of problems (31) is approximated using a Galerkin approximation in the finite element ($\mathbb{P}_{1}$) discretization space

$$W^{D}={\rm Span}\left\{\phi^{1}\wedge\phi^{2},\phi^{1},\phi^{2}\in V^{D}\right\}$$

where

$$V^{D}=\left\{\phi\in\mathcal{C}(\Omega)|\quad\phi(-L)=\phi(L)=0,\;\phi|_{\left(-L+\frac{(i-1)2L}{D},-L+\frac{i2L}{D}\right)}\in\mathbb{P}_{1},\;\forall 0\leq i\leq D\right\},$$

and

$$\forall\phi^{1},\phi^{2}\in V^{D},\;\forall x,y\in\Omega,\quad\phi^{1}\wedge\phi^{2}(x,y)=\frac{1}{\sqrt{2}}(\phi^{1}(x)\phi^{2}(y)-\phi^{2}(x)\phi^{1}(y)).$$

The moment functions $(\varphi_{m})_{1\leq m\leq M}$ are chosen to be $\mathbb{P}_{1}$ hat functions associated to a uniform discretization of $\Omega$ so that

$$Z^{M}:={\rm Span}\left\{\varphi_{m},\;1\leq m\leq M\right\}=\left\{v\in\mathcal{C}(\Omega)|\quad\;v|_{\left(-L+\frac{(j-1)2L}{M-1},-L+\frac{j2L}{M-1}\right)}\in\mathbb{P}_{1},\;\forall 0\leq j\leq M-1\right\}.$$

The electronic density $\rho$ of choice is constructed as follows: we define, for all $x\in\Omega$,

$$\phi_{\rm even}(x)=1-\frac{|x|}{L}\;\mbox{ and }\;\phi_{\rm odd}(x)=\begin{cases}1-|2x+L|/L&\mbox{ if }x\leq 0,\\ \frac{|2x-L|}{L}-1&\mbox{ otherwise}.\\ \end{cases}$$

Then, we define $\widetilde{\Psi}^{0}_{1}=\frac{\phi_{\rm even}}{\|\phi_{\rm even}\|_{L^{2}(\Omega)}}\wedge\frac{\phi_{\rm odd}}{\|\phi_{\rm odd}\|_{L^{2}(\Omega)}}$ and $\rho=\frac{1}{2}\left(\frac{|\phi_{\rm even}|^{2}}{\|\phi_{\rm even}\|_{L^{2}(\Omega)}^{2}}+\frac{|\phi_{\rm odd}|^{2}}{\|\phi_{\rm odd}\|_{L^{2}(\Omega)}^{2}}\right).$

We then apply the MCAL algorithm starting from $\widetilde{\mathfrak{P_{0}}}=\left\{\widetilde{\Psi}^{0}_{1}\right\}$.

Let us first highlight the influence of the parameter $q_{\rm vec}$ on the performance of the algorithm in terms of the number of iterations required to achieve numerical convergence. We first conduct a first series of tests with $M=20$ and $D=100$.

Figure 1 highlights the behaviour of the numerical scheme with respect to $q_{\rm vec}$.

Figure 1: Evolution of $F^{n}$ (above) and $|E(v_{n})|$ (below) as a function of $n$ for different values of $M$.

The upper (respectively lower) figure shows the values of $F^{n}$ (respectively $|E(v_{n})|$) as a function of $n$ for different values of $M$. As predicted by our theoretical results, for any value of $q_{\rm vec}$, the sequence $(F^{n})_{n\in\mathbb{N}}$ is non-increasing and we also checked numerically that $\widetilde{F}^{n}=F^{n+1}$ for all $n\in\mathbb{N}$. In constrast, the sequence $(E(v_{n}))_{n\in\mathbb{N}^{*}}$ is not monotonous. We also observe that for any tested value of $q_{\rm vec}$, the sequence $(F^{n})_{n\in\mathbb{N}}$ converges to the same limit value. It seems that for greater values of $q_{\rm vec}$, the number of iterations $n$ needed for the algorithm to converge is lower.

Figure 2 highlights the behaviour of the numerical scheme with respect to the number $M$ of moment constraints. In these tests, $D=100$ and $q_{\rm vec}=4$.

Again, the upper (respectively lower) figure shows the values of $F^{n}$ (respectively $|E(v_{n})|$) as a function of $n$ for different values of $M$. As before, we observe that the sequence $(F^{n})_{n\in\mathbb{N}}$ is non-increasing and we also checked numerically that $\widetilde{F}^{n}=F^{n+1}$ for all $n\in\mathbb{N}$. In constrast, the sequence $(E(v_{n}))_{n\in\mathbb{N}^{*}}$ is not monotonous. We also observe that for any tested value of $M$, the sequence $(F^{n})_{n\in\mathbb{N}}$ converges to some limit value denoted here by $F^{\infty}(M)$ which depends on $M$. We observe again that the value of $F^{\infty}(M)$ does not increase monotonically with $M$, which stems from the fact that the spaces $Z^{M}$ do not form an increasing family of vector spaces for the inclusion. However, is still holds that $Z^{10}\subset Z^{20}\subset Z^{40}$, and we indeed observe that $F^{\infty}(10)\leq F^{\infty}(20)\leq F^{\infty}(40)$, which is coherent with the variational structure of the moment constraint approach studied here. We also observe that the value of $E(v_{n})$ seems to stagnate in most of the numerical tests (except the one corresponding to $M=10$) to a value close to $-10^{-5}$.

Figure 2: Evolution of $F^{n}$ (above) and $|E(v_{n})|$ (below) as a function of $n$ for different values of $M$.

Lastly, Figure 3 shows the plots of the potential $v_{n}$ obtained after running $n=80$ iterations of the MCAL algorithm for various values of $M$ ($M=10,20,30,40$). We observe that the potential value seems to converge to some limit value of $M$ increases. However, the number of moment constraints should definitely be higher to obtain a better accuracy, which was not possible with our current implementation. More evolved versions of the present MCAL algorothm should be designed to alleviate this bottleneck, which will be the object of a future work.

Figure 3: $v_{n}(x)$ as a function of $x$ for $M=10,20,30,40$ ($n=80$)

7 Proofs

We gather in this section the proofs of our main theoretical results.

7.1 Proof of Theorem 6

Proof of Theorem 6.

Step 1: (Finiteness) Since $\rho\in\mathcal{I}_{N}$, there exists at least one element $\Psi_{0}\in\mathcal{H}_{1}^{N}$ such that $\rho_{\Psi_{0}}=\rho$. Denoting by $\Gamma_{0}:=|\Psi_{0}\rangle\langle\Psi_{0}|$, it can then be easily seen that $\Gamma_{0}\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N},\Phi,\rho)$ and that ${\rm Tr\,}(\Theta\Gamma_{0})=\int_{\mathbb{R}^{3}}\theta(|x|)\rho(x)\,dx=C_{\rho}$. Thus, we immediately obtain that for all $C\geq C_{\rho}$, $F_{L,\theta}^{\Phi,C}[\rho]>-\infty$.

Step 2: (Existence of minimizer) Let $(\Gamma_{n})_{n\in\mathbb{N}}$ be a minimizing sequence associated to (12). Then, we know from the proof of Theorem 4.4 of [Lie83a] that, up to the extraction of a subsequence, there exists a trace-class operator $\Gamma_{\infty}\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})$ such that $\left((H_{N}+D)^{1/2}\Gamma_{n}(H_{N}+D)^{1/2}\right)_{n\in\mathbb{N}}$ weakly converges in the sense of trace-class operators to $(H_{N}+D)^{1/2}\Gamma_{\infty}(H_{N}+D)^{1/2}$ as $n$ goes to infinity. To prove that $\Gamma_{\infty}$ is a minimizer to (12), it is sufficient to prove that $\rho_{\Gamma_{\infty}}$ satisfies

$$\forall 1\leq m\leq M,\;\int_{\mathbb{R}^{3}}\rho_{\Gamma_{\infty}}\varphi_{m}=\int_{\mathbb{R}^{3}}\rho\varphi_{m}\;\mbox{ and }\;\int_{\mathbb{R}^{3}}\rho_{\Gamma_{\infty}}(x)\theta(|x|)\,dx={\rm Tr\,}(\Theta\Gamma_{\infty})\leq C.$$

For all $n\in\mathbb{N}$, let us denote by $\tau_{n}\in L^{2}(\mathbb{R}^{3N}\times\mathbb{R}^{3N})$ the kernel of $\Gamma_{n}$ and by $\tau_{\infty}\in L^{2}(\mathbb{R}^{3N}\times\mathbb{R}^{3N})$ the kernel of $\Gamma_{\infty}$. Let us also denote for all $n\in\mathbb{N}$,

$$\gamma_{n}(x_{1},\ldots,x_{N}):=\tau_{n}(x_{1},\ldots,x_{N};x_{1},\ldots,x_{N})$$

and by

$$\gamma_{\infty}(x_{1},\ldots,x_{N}):=\tau_{\infty}(x_{1},\ldots,x_{N};x_{1},\ldots,x_{N})$$

for all $x_{1},\ldots,x_{N}\in\mathbb{R}^{3}$. Let us prove that $(\gamma_{n})_{n\in\mathbb{N}}$ is a tight sequence. Indeed, let $R>0$ and $B_{R}$ be the ball of radius $R$ of $\mathbb{R}^{3N}$. Then, denoting by $\mathds{1}_{B^{c}_{R}}$ the characteristic function of the set $B_{R}^{c}$, it holds that for all $n\in\mathbb{N}$,

$$\int_{B^{c}_{R}}\gamma_{n} =\int_{\mathbb{R}^{3N}}\mathds{1}_{B^{c}_{R}}\gamma_{n}$$

$$\leq\int_{\mathbb{R}^{3N}}\left(\frac{1}{N}\sum_{i=1}^{N}\frac{\theta(|x_{i}|)}{\theta(R)}\right)\gamma_{n}(x_{1},\ldots,x_{N})\,dx_{1}\ldots\,dx_{N}$$

$$=\frac{1}{\theta(R)}{\rm Tr\,}(\Theta\Gamma_{n})\leq\frac{C}{\theta(R)}.$$

Let us denote by $M_{P}$ the multiplication operator by any function $P$ bounded with compact support on ${\mathbb{R}}^{3N}$. We then know from the proof of Theorem 4.4 of [Lie83a] that

$${\rm Tr}(M_{P}\Gamma^{\infty})=\mathop{\lim}_{n\to+\infty}{\rm Tr}(M_{P}\Gamma^{n}).$$

This, together with the tightness result above, yields that $(\rho_{\Gamma_{n}})_{n\in\mathbb{N}}$ weakly converges to $\rho_{\Gamma_{\infty}}$ in $L^{1}(\mathbb{R}^{3})$. It then easily follows that for all $m=1,\cdots,M,$

$$\int_{{\mathbb{R}}^{3}}\varphi_{m}\rho_{\Gamma_{\infty}}=\lim_{n\to+\infty}\int_{{\mathbb{R}}^{3}}\varphi_{m}\rho_{\Gamma_{n}}=\int_{{\mathbb{R}}^{3}}\varphi_{m}\rho$$

and that

$$\int_{{\mathbb{R}}^{3}}\theta(|x|)\rho_{\Gamma_{\infty}}(x)\,dx={\rm Tr\,}(\Theta\Gamma_{\infty})\leq C.$$

The operator $\Gamma_{\infty}$ is thus a minimizer of (12). In particular, since $\mathds{1}\in{\rm Span}\{\Phi\}$, it holds that ${\rm Tr\,}(\Gamma_{\infty})=N$.

Step 3: (Existence of a sparse minimizer)

Let us now introduce the function $\Lambda:\mathcal{H}_{1}^{N}\to{\mathbb{R}}^{M+1}$ such that for all $m=1,\cdots,M,$

$$\Lambda_{m}(\Psi)=\int_{{\mathbb{R}}^{3}}\varphi_{m}(x)\rho_{\Psi}(x)\mathrm{d}x=\int_{{\mathbb{R}}^{dN}}\varphi_{m}(x)|\Psi(x,x_{2},...,x_{N})|^{2}\mathrm{d}x\mathrm{d}x_{2}...\mathrm{d}x_{N},$$

and

$$\Lambda_{M+1}(\Psi)=\langle\Psi|H_{N}|\Psi\rangle.$$

It can then be easily seen that $\Lambda$ is a continuous map on $\mathcal{H}_{1}^{N}$.

Let $\Gamma_{\rm min}$ be a minimizer of (12). Then, there exists a countable index set $\mathcal{J}\subset\mathbb{N}$, an orthonormal family $(\Psi_{j})_{j\in\mathcal{J}}$ of $\mathcal{H}_{0}^{N}$ and a family of positive numbers $(\alpha_{j})_{j\in\mathcal{J}}$ such that $\sum_{j\in\mathcal{J}}\alpha_{j}=N$ (this comes from the fact that $\mathds{1}\in{\rm Span}\{\Phi\}$) and

$$\Gamma_{\rm min}=\sum_{j\in\mathcal{J}}\alpha_{j}|\Psi_{j}\rangle\langle\Psi_{j}|.$$

In addition, it can be easily checked that $\Psi_{j}\in\mathcal{H}_{1}^{N}$ for all $j\in\mathcal{J}$. We then define $\mu_{\rm min}:=\sum_{j\in\mathcal{J}}\alpha_{j}\delta_{\Psi_{j}}$ which is a Borel measure on $\mathcal{B}(\mathcal{H}_{1}^{N})$ since ${\rm Tr\,}(H_{N}\Gamma_{\rm min})$ is finite and ${\rm Tr\,}\Gamma_{\rm min}=N$. It can then be easily checked that

$$\int_{\mathcal{H}_{1}^{N}}\|\Lambda(\Psi)\|\,d\mu_{\rm min}(\Psi)<+\infty.$$

Thus, by Proposition 5, there exist $1\leq K\leq M+1$, $\overline{\Psi}_{1},\cdots,\overline{\Psi}_{K}\in\mathcal{H}_{1}^{N}$ and $\omega_{1},\cdots,\omega_{K}>0$ such that

$$\int_{\mathcal{H}_{1}^{N}}\Lambda(\Psi)\mathrm{d}\mu_{\rm min}(\Psi)=\sum_{k=1}^{K}\omega_{k}\Lambda(\overline{\Psi}_{k}).$$

Denoting by $\Gamma_{K}=\sum_{k=1}^{K}\omega_{k}|\overline{\Psi}_{k}\rangle\langle\overline{\Psi}_{k}|$, it can then be easily checked that $\Gamma_{K}$ is also a minimizer to (12). Hence the desired result.

∎

7.2 Proof of Theorem 17

Proof of Theorem 17.

The first assertion of the theorem is a direct consequence of Theorem 6. Using the same arguments as in the proof of Theorem 6, one can easily obtain that the sequence $\left((H_{N}+D)^{1/2}\Gamma_{n}(H_{N}+D)^{1/2}\right)_{n\geq n_{0}}$ is compact in $\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})$. Thus, up to the extraction of a subsequence there exists $\Gamma_{\infty}\in\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})$ such that ${\rm Tr\,}(H_{N}\Gamma_{\infty})<+\infty$ and such that $\left((H_{N}+D)^{1/2}\Gamma_{n}(H_{N}+D)^{1/2}\right)_{n\geq n_{0}}$ weakly converges to $\left((H_{N}+D)^{1/2}\Gamma_{\infty}(H_{N}+D)^{1/2}\right)_{n\geq n_{0}}$ in the sense of trace-class operators of $\mathfrak{S}_{1}^{+}(\mathcal{H}_{0}^{N})$.

Moreover, following again the same lines of proof, we obtain that the sequence $(\rho_{\Gamma_{n}})_{n\geq n_{0}}$ weakly converges in $L^{1}(\mathbb{R}^{3})$ to $\rho_{\Gamma_{\infty}}$. As a consequence, it holds that $\int_{\mathbb{R}^{3}}\theta(|x|)\rho_{\Gamma_{\infty}}(x)\,dx\leq C$. Moreover, since for all $n\in{\mathbb{N}}^{*}$, $\int_{\mathbb{R}^{3}}\varphi^{n}_{m}\rho_{\Gamma_{n}}=\int_{\mathbb{R}^{3}}\varphi^{n}_{m}\rho$, using Lemma 15, we then obtain that, necessarily, $\rho_{\Gamma_{\infty}}=\rho$. This makes $\Gamma_{\infty}$ admissible for (5) so that we have that ${\rm Tr\,}(H_{N}\Gamma_{\infty})\geq F_{L}[\rho]$. Notice now that for all $n\geq n_{0}$, $-\infty<F^{\Phi^{n},C}_{L,\theta}[\rho]\leq F_{L}[\rho]$. Thus for any converging subsequence of $(F^{\Phi^{n},C}_{L,\theta}[\rho])_{n\geq n_{0}}$ to some limit $F_{L}^{\infty}$, it holds that $-\infty<F_{L}^{\infty}\leq F_{L}[\rho]$. For this subsequence, still denoted by $(F^{\Phi^{n},C}_{L,\theta}[\rho])_{n\geq n_{0}}$ for the sake of simplicity, it holds that $\mathop{\lim}_{n\to\infty}{\rm Tr\,}(H_{N}\Gamma_{n})=F_{L}^{\infty}$, and we then have that

$$F_{L}[\rho]\geq F^{\infty}_{L}\geq{\rm Tr\,}(H_{N}\Gamma_{\infty}).$$

Thus, necessarily, $\Gamma_{\infty}$ is a minimizer of (5). Moreover, $F_{L}^{\infty}=F_{L}[\rho]$ for any extracted subsequence so that $\displaystyle\mathop{\lim}_{n\to+\infty}{\rm Tr}(H_{N}\Gamma_{n})={\rm Tr}(H_{N}\Gamma_{\infty})$. Using the compactness of the Fock space of bounded particle number for the geometric convergence [Lew11][Lemma 2.2, Lemma 2.3], we thus obtain the desired result. ∎

7.3 Proof of Theorem 24

Proof of Theorem 24.

Step 1: Let us first prove that there exists a maximizer to the optimization problem

$$\mathop{\sup}_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\},\\ \forall\Psi\in{\rm Span}\{\Psi_{1},\ldots,\Psi_{K}\},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0\\ \end{array}}\int_{\Omega}v\rho.$$ (33)

We denote here by $\mathcal{S}^{K}$ the set of symmetric matrices of $\mathbb{R}^{K\times K}$. For any $\varphi\in{\rm Span}\{\Phi\}$, let us consider the linear form $l_{\varphi}:\mathcal{S}^{K}\to\mathbb{R}$ defined as follows:

$$\forall S:=(S_{kl})_{1\leq k,l\leq K}\in\mathcal{S}^{K},\quad l_{\varphi}(S):=\int_{\Omega}\varphi(x)\sum_{k,l=1}^{K}S_{kl}\overline{\Psi_{k}}(x)\Psi_{l}(x)\,dx={\rm Tr}(\varphi\Gamma_{S}),$$

where

$$\Gamma_{S}:=\sum_{k,l=1}^{K}S_{kl}|\Psi_{k}\rangle\langle\Psi_{l}|.$$

Let us now consider the vectorial space

$$L:=\{l_{\varphi},\;\varphi\in{\rm Span}\{\Phi\}\}.$$

The space $L$ is a finite-dimensional subspace of the set of linear forms on $\mathcal{S}^{K}$, and its dimension $J$ is lower or equal to the dimension of ${\rm Span}\{\Phi\}$. Let $(\widetilde{l}_{1},\ldots,\widetilde{l}_{J})$ be a basis of $L$. By construction, there exists $\widetilde{\varphi}_{1},\ldots,\widetilde{\varphi}_{J}\in{\rm Span}\{\Phi\}$ such that $\widetilde{l}_{j}=l_{\widetilde{\varphi}_{j}}$ for all $1\leq j\leq J$. Let us then denote by $\widetilde{\Phi}:=\{\widetilde{\varphi}_{1},\ldots,\widetilde{\varphi}_{J}\}$. It can then be easily checked that any element $\varphi$ of ${\rm Span}\{\Phi\}$ can be rewritten as

$$\varphi=\widetilde{\varphi}+\varphi_{0},$$

where $\widetilde{\varphi}\in{\rm Span}\{\widetilde{\Phi}\}$ and $\varphi_{0}\in{\rm Span}\{\Phi\}$ such that $l_{\varphi_{0}}=0$. In particular, this implies that $\int_{\Omega}\varphi_{0}\rho=0$ since for all $\varphi\in{\rm Span}\{\Phi\}$,

$$\int_{\Omega}\varphi\rho=\int_{\Omega}\varphi(x)\sum_{k=1}^{K}\omega_{k}|\Psi_{k}(x)|^{2}\,dx=l_{\varphi}({\rm diag}(\omega_{1},\ldots,\omega_{K})).$$

Thus, proving that there exists a maximizer to (33) is equivalent to proving that there exists a maximizer to

$$\mathop{\sup}_{\begin{array}[]{c}v\in{\rm Span}\{\widetilde{\Phi}\},\\ \forall\Psi\in{\rm Span}\{\Psi_{1},\ldots,\Psi_{K}\},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0\\ \end{array}}\int_{\Omega}v\rho.$$ (34)

Now, by definition of $\widetilde{\varphi}_{1},\ldots,\widetilde{\varphi}_{J}$, it holds that the application $A:\mathcal{S}^{K}\to\mathbb{R}^{J}$ defined so that for all $1\leq j\leq J$ and all $S=(S_{kl})_{1\leq k,l\leq K}$,

$$A(S)_{j}:=\int_{\Omega}\widetilde{\varphi}_{j}\sum_{k,l=1}^{K}S_{kl}\overline{\Psi_{k}}\Psi_{l}$$

is surjective. Indeed, this comes from the fact that ${\rm dim}\;{\rm Rank}(A)={\rm dim}\;L=J$. It can then be easily checked that (34) is then equivalent to the dual semi-definite programming problem:

$$\mathop{\sup}_{\begin{array}[]{c}(y,S)\in\mathbb{R}^{J}\times\mathcal{S}^{K}\\ A^{*}(y)+S=C\\ S\succcurlyeq 0\\ \end{array}}\langle b,y\rangle,$$ (35)

where $b=(b_{j})_{1\leq j\leq J}$ is such that $b_{j}=\int_{\Omega}\widetilde{\varphi}_{j}\rho$ for all $1\leq j\leq J$ and $C=(C_{kl})_{1\leq k,l\leq K}\in\mathcal{S}^{K}$ with

$$C_{kl}:=\langle\Psi_{k}|H_{N,\Omega}|\Psi_{l}\rangle\quad\forall 1\leq k,l\leq K.$$

Indeed, if $(y,S)\in\mathbb{R}^{J}\times\mathcal{S}^{K}$ is a maximizer to (35), it holds that $v=\sum_{j=1}^{J}y_{j}\widetilde{\varphi}_{j}$ is a maximizer to (34), and thus to (33).

The primal problem associated to (35) reads as

$$\mathop{\inf}_{\begin{array}[]{c}X\in\mathcal{S}^{K}\\ A(X)=b\\ X\succcurlyeq 0\\ \end{array}}\langle C,X\rangle,$$ (36)

Let us also remark that $\int_{\Omega}\rho\varphi=\int_{\Omega}\rho_{\Gamma_{S}}\varphi$ for all $\varphi\in{\rm Span}\{\Phi\}$ if and only if $A(S)=b$. Thus, this implies that there exists at least one minimizer $X$ to (36) which is given by $X={\rm diag}(\omega_{1},\ldots,\omega_{K})$ and is positive definite. Using Theorem 22, we then obtain the existence of at least one maximizer to (35), and hence to (33) and (34).

Step 2: To conclude the proof of the desired result, it only remains to show that

$$D_{L,\Omega}^{\Phi}[\rho] :=\mathop{\sup}_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\},\\ \forall\Psi\in\mathcal{H}_{1}^{N}(\Omega),\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0\\ \end{array}}\int_{\Omega}v\rho$$

$$=\mathop{\sup}_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\},\\ \forall\Psi\in{\rm Span}\{\Psi_{1},\ldots,\Psi_{K}\},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0\\ \end{array}}\int_{\Omega}v\rho.$$

On the one hand, it holds from Theorem 23, that $D_{L,\Omega}^{\Phi}[\rho]=F_{L,\Omega}^{\Phi}[\rho]$. On the other hand, using similar arguments as in the proof of Theorem 23, it holds that

$$\mathop{\sup}_{\begin{array}[]{c}v\in{\rm Span}\{\Phi\},\\ \forall\Psi\in\mathcal{W},\quad\langle\Psi|H_{N,\Omega}^{v}|\Psi\rangle\geq 0\\ \end{array}}\int_{\Omega}v\rho=\mathop{\inf}_{\Gamma\in{\mathfrak{S}}_{1}^{+}(\mathcal{W},\Phi,\rho)}{\rm Tr}(H_{N,\Omega}\Gamma),$$

where $\mathcal{W}:={\rm Span}\{\Psi_{1},\ldots,\Psi_{K}\}$. Since, by definition of $\Psi_{1}$, …, $\Psi_{K}$, it holds that $\displaystyle F_{L,\Omega}^{\Phi}[\rho]=\mathop{\inf}_{\Gamma\in{\mathfrak{S}}_{1}^{+}(\mathcal{W},\Phi,\rho)}{\rm Tr}(H_{N,\Omega}\Gamma)$, we obtain the desired result. ∎