Rigorous Derivation of the Gross-Pitaevskii Equation

The time dependent Gross-Pitaevskii equation describes the dynamics of initially trapped Bose-Einstein condensates. We present a rigorous proof of this fact starting from a many-body bosonic Schroedinger equation with a short scale repulsive interaction in the dilute limit. Our proof shows the persistence of an explicit short scale correlation structure in the condensate.

with the normalization |ϕ(r)| 2 d 3 r = 1.The coupling constant is σ = 8πN a where N is the number of particles and a is the scattering length of the interaction potential.A single orbital theory might indicate that the full Nbody wave function essentially factorizes.This, however, is usually not the case: very short scale interactions may introduce strong pair-correlations that substantially influence the energy of the system and the dynamics even on larger scales.In particular, the emergence of the scattering length in σ is a correlation effect and it is remarkable that this correlation can be consistently modeled by the same coupling constant along the whole evolution.
The GP equation thus implicitly postulates a persistent two-scale structure of the state.On short scales, the state exhibits a time independent pair-correlation, while on large scales it is given by the product of N copies of a time dependent orbital.It is a fundamental problem in many-body theory to show that the GP postulation of a two-scale structure is a rigorous consequence of the first principle many-body Schrödinger equation.While many theoretical work addressed this question for low energy states, there is no rigorous result from a dynamical point of view.In this letter, we shall show that not only the two-scale structure is preserved, but it even emerges dynamically for a class of initial data.
On the mathematically rigorous level, the GP theory has only been verified for the ground state of the trapped Bose gas.Lieb, Seiringer and Yngvason [6] considered a model where the scattering length a of the interaction varies with N , so that N a remains constant.
For repulsive interactions they proved that the ground state energy per particle in the GP limit (N → ∞, a 0 = N a = fixed) is given by a variational principle, min{E GP (ϕ) : ϕ = 1}, where is the GP energy functional.This result was extended to rotating Bose systems that model superfluidity in [7].In [8], Lieb and Seiringer have also rigorously established that the ground state of the trapped Bose gas exhibits a complete Bose-Einstein condensation (BEC) by showing that the one particle density matrix of the ground state converges to the one dimensional projection |ϕ GP ϕ GP |, where ϕ GP is the minimizer of the functional E GP .
The ground state has been known to have a non-trivial pair correlation.Mathematically, this was first exploited for the homogeneous dilute Bose gas, when Dyson [9] used a specific nearest neighbor correlation in his trial state to verify the upper bound on the ground state energy obtained from Bogoliubov theory [10].The proof of the lower bound, achieved by Lieb and Yngvason forty years later [11], demonstrates that low energy states necessarily contain a characteristic short scale pair correlation.In the GP limit, as shown in [6], the trapping potential introduces a second length scale but it leaves the short scale structure intact.
The works mentioned above justify the two-scale hypothesis of the GP theory based on the energy minimizing principle.In this letter, we outline a mathematical proof (see [12,13] for details) of the fact that the GP theory also describes the dynamics of trapped Bose-Einstein condensates.Assuming that at time t = 0 we have a condensate, we show, under some conditions on the interaction potential and on the initial state, that the system still exhibits complete condensation at times t = 0, and that the condensate wave function evolves according to the GP equation (1).Our results cover two types of initial data: type I states with the physical two-scale structure as in [6,8] and type II states with no short scale structure.For type I initial data, the main difficulty is to demonstrate the persistence of the two-scale structure even for states that are far from the ground state.Energy conservation alone does not prevent the destruction of the short scale structure.The key observation is that higher moments of the energy are also conserved and they impose stringent restrictions on the state.In Theorem 2 below we show that a finite second moment of the energy per particle forces a pair correlation described by the zero energy scattering function of the interaction potential.For type II initial data, our results imply that the short scale structure emerges dynamically.See the discussion after Theorem 1.
We now define the model precisely.The Hamiltonian for N bosons in three dimensions is given by , where V is a positive, spherically symmetric, compactly supported, smooth potential with scattering length a 0 .By scaling, the scattering length of V N is a 0 /N .We assume that the strength of the interaction, measured by the dimensionless quantity α = V (r)|r| −1 dr + sup r∈R 3 r 2 V (r), is sufficiently small.
Let ψ N be the ground state of (3) and let γ for any fixed k.In particular, (4) for k = 1 means that ψ N exhibits complete condensation.After instantaneously removing the trap, the evolution of the system is generated by the Hamiltonian Theorem 1.Let ψ N,t be the solution to the Schrödinger equation i∂ t ψ N,t = H N ψ N,t with the initial condition ψ N,0 = ψ N , and let γ (1) N,t be its one particle marginal density.Then, for any fixed time t ∈ R, ψ N,t exhibits complete Bose-Einstein condensation, that is where ϕ t solves the Gross-Pitaevskii equation with initial data ϕ t=0 = ϕ GP .The convergence in (6) is in the sense that Tr K (γ Our proof, in fact, allows for more general initial data.Let H 1 (R 3 ) := {ϕ(r) : |ϕ| 2 + |∇ϕ| 2 < ∞}.The technical condition we need to conclude (6), (7) is that the initial state ψ N,0 asymptotically factorizes in the sense that there exists ϕ ∈ H 1 (R 3 ) and for every k ≥ 1 there exists a family {ξ The function ϕ is then the initial data for (7).As we prove in Appendix C of [13], the ground state ψ N of the Hamiltonian (3) satisfies (8).This shows that the pair correlations of ψ N , which affect the dynamics in a non-trivial way, are not in conflict with the asymptotic factorization as long as the convergence in ( 8) is in a weaker topology than the energy norm.The local correlation structure of ψ N lives on a length scale 1/N and its effect, measured in L 2 -sense in (8), is negligible.However, the energy of ψ N does not converge to the energy of the factorized state, and, similarly, the convergence in (4) does not hold in the H 1 space or in energy sense.
A typical type II initial data is given by a product wave function ψ N,0 = ϕ ⊗N with ϕ ∈ H 1 (R 3 ).For ψ N,0 the condition (8) holds trivially.Our result thus shows that a product initial state also evolves according to the GP equation ( 7) (with initial data ϕ t=0 = ϕ) despite that its total energy is larger than E GP (ϕ).In fact, it is easy to see that the energy per particle in the factorized state, ϕ ⊗N , ( H N /N )ϕ ⊗N , is asymptotically given by a functional similar to (2) but with 8πa 0 replaced with its Born approximation, b 0 = V (r)d 3 r.Note that b 0 > 8πa 0 .This indicates that the excess energy originating from the emergence of the scattering length is dispersed into modes that do not influence the evolution of the condensate.Hence, for product initial data, the scattering length emerges dynamically and the energy functional of the single particle orbitals does not predict the evolution equation.
Outline of the proof.For k = 1, . . ., N , let γ where Tr k+1 denotes the partial trace over r k+1 .By the Banach-Alaouglu theorem, γ  ∞,t to be described by the infinite hierarchy with σ = b 0 .However, the coupling constant in the correct equations for γ develops a short scale structure which lives on the same length scale ℓ ≃ 1/N as the potential V N .Since γ is not constant on this length scale, one cannot apply the formal limit N V N → b 0 δ in (9).
In the following Theorem 2 we show that for r j close to r k+1 , the short scale structure of γ (k+1) N,t is described by the function f N (r j − r k+1 ), where f N is the solution to the zero energy scattering equation with boundary condition f N (r) → 1 as |r| → ∞.By scaling, f N (r) = f (N r), where f is a solution of −∆ + 1 2 V f = 0 with the same boundary condition.Therefore the correct value of σ in ( 10) is given by Note that the short scale structure disappears from γ (k) ∞,t after taking the weak limit, but it still affects the limiting macroscopic dynamics.
The infinite hierarchy (10) has a factorized solution: the family γ (10) with σ = 8πa 0 if and only if ϕ t is a solution to the GP equation (7).Therefore, to identify the limiting density, it suffices to show that (i) any limit point {γ k=1 is a solution of the infinite hierarchy (10) with σ = 8πa 0 , and (ii) the solution to the infinite hierarchy is unique.The proofs of (i), (ii) rely on estimates on the energy distribution of the initial wave function and on the factorization property (8).
(i) Convergence to the infinite hierarchy.The following key theorem identifies the short scale structure assuming a bound on the second moment of the energy.This is the main ingredient in the proof of (10) Theorem 2. Suppose α is small enough.Then there exists a constant C > 0 such that for all 1 ≤ i < j ≤ N and all ψ ∈ L 2 (R 3N , dR) symmetric with respect to permutations.Here R = (r 1 , . . ., r N ) ∈ R 3N .
From this theorem, we obtain the a-priori bound uniformly in N and t provided the initial wave function 13) identifies the short scale structure of ψ N,t with a precision 1/N .Furthermore, since the left side of ( 12) is a constant of motion, (13) shows that the separation between the singular short scale structure and the regular part of ψ N,t is preserved by the time evolution.
Outline of the proof of Theorem 2: we define h m and, using the permutation symmetry of ψ and V ≥ 0, we obtain, for arbitrary i = j, From the definition (11) of f N (r), we find with L i = −∆ ri + 2∇ ri (log f N (r i − r j ))∇ ri and an analogous identity for L j .Using integration by parts, and with Eq. ( 12) now follows because and therefore the second term on the right hand side of the last equation can be controlled by the first one (for α small enough) using the operator inequality |r| −2 ≤ −C∆ r (Hardy inequality).
(ii) Uniqueness of the infinite hierarchy.The first step is to prove a-priori bounds in a certain Sobolev norm.
∞,t be any weak limit point of γ N,t , then the following estimate holds uniformly in time can be obtained through (8).The difficulty is that the norm ||| • ||| k cannot be directly controlled by TrH k N (•) and in fact |||γ because of the singular short scale structure.It is only after taking the weak limit N → ∞ that the short scale structure disappears and (14) can be proven.
For illustration, consider the case k = 2 discussed in Theorem 2. The estimate (13) implies that ψ N,t (R) ∼ f N (r i − r j )φ t (R) where φ t is a smooth function in the variable r i − r j .Together with the fact that dr i dr j |∇ ri ∇ rj f N (r i − r j )| 2 ∼ N , we have that dR ∇ ri ∇ rj ψ N (R) 2 → ∞ as N → ∞.However, since f N → 1 weakly in L 2 -sense (but not in energy sense), a Sobolev estimate on the limit of the density matrix of ψ N,t can be deduced from (13).
The uniqueness holds in the Sobolev norm ||| • ||| k : (10) with Γ t=0 = Γ and such that |||γ Outline of the proof.Iterating the integral form of (10), we obtain a Dyson series where The error term η     k+m) .The graphical structure of Λ encodes the collision history given in (16).Every line e of Λ carries a regularized free propagator, (α e − p 2 e ± i/t) −1 with a momentum variable p e ∈ R 3 and a frequency α e ∈ R. At each vertex, there is a p-and a α-delta function due to momentum and energy conservation.The kernel of ω (k) m,t (whose variables correspond to the 2k momenta of the roots) is computed by performing 3m momentum integrals and 2k + 3m frequency integrals in each Λ.
Because of the singularity of the interaction at r = 0, each graph is potentially ultraviolet divergent.Power counting suggests, however, that the integrals are finite.Suppose we cutoff all momentum integrals at |p| ≃ β ≫ 1, and all frequency integrals at |α| ≃ β 2 ; then the integration volume scales as β 3(3m)+2(2k+3m) = β 4k+15m .The m momentum and frequency delta functions scale as β −5m and the 2k + 3m propagators scale as β −2(2k+3m) .The a-priori estimate (14) shows that γ (k+m) β −5(k+m) .Since 4k+15m < 5m+2(2k+3m)+ 5(k + m), the integrals should be convergent in the ultraviolet regime β ≫ 1.To make this argument rigorous, we use an integration scheme, dictated by the structure of the graph.We start by integrating the momenta and frequency of the leaves; these integrals are convergent because of the a-priori estimates (14) and they also provide a momentum decay on the lines adjacent to the leaves.We then iterate this procedure, integrating all momenta and frequencies by moving from the right to the left of the graph and transferring a suitable momentum decay.

Conclusion.
We have proven that Bose-Einstein condensates evolve according to the Gross-Pitaevskii equation.This provides a mathematical description of recent experiments on the evolution of initially trapped condensates.On the theoretical level, our result for factorized initial wave functions is more surprising.The emergence of the scattering length in the GP equation is a consequence of the short scale correlations.We show, however, that the GP equation is correct even if the initial state is uncorrelated (product).Our result thus suggests that the many body dynamics builds up correlations on the length scale ℓ ≃ 1/N in a very short time.Since the emergence of correlations reduces the local energy, one may ask what happens with the excess energy.Although we are not able to give a mathematically rigorous answer, we believe that the excess energy is transferred to incoherent excitations living on intermediate length scales ℓ, with N −1 ≪ ℓ ≪ 1.Since the macroscopic dynamics described by the GP equation is affected only by the structure of the wave function on length scales of O(1) and O(1/N ), these mesoscopic excitations have no influence on the evolution of the condensate.
denote the k-particle marginal density of ψ N,t with a consistent normalization Tr γ (k) N,t = 1.The time evolution is governed by a hierarchy of N coupled equations, commonly known as the BBGKY hierarchy has at least one weak* limit point, γ , in the space of trace class operators.Since formally N V N (r) = N 3 V (N r) → b 0 δ(r), one may naively expect from (9) the time evolution of γ (k) } k≥1 of the family of densities {γ

FIG. 1 :
FIG. 1: Diagrammatic expansion for ω(k) m,t as sum of Feynman graphs with m vertices, 2(k + m) leaves and 2k roots has the same form as ω with U (k+n) sn γ(k+n) replaced by the full evolution γ (k+n) sn .To prove the convergence of the expansion (15) as n → ∞, we expand each term into a sum of contributions associated with certain Feynman graphs.A typical graph Λ contributing to ω is drawn in Fig.1.It has m four-valent vertices and 2k +3m lines.The external lines on the left (called roots) correspond to the 2k momenta variables of the operator kernel of ω . The 2(k + m) external lines on the right (called leaves) represent the kernel of γ (