The Schrödinger equation

\displaystyle i \hbar \partial_t |\psi \rangle = H |\psi\rangle

is the fundamental equation of motion for (non-relativistic) quantum mechanics, modeling both one-particle systems and {N}-particle systems for {N>1}. Remarkably, despite being a linear equation, solutions {|\psi\rangle} to this equation can be governed by a non-linear equation in the large particle limit {N \rightarrow \infty}. In particular, when modeling a Bose-Einstein condensate with a suitably scaled interaction potential {V} in the large particle limit, the solution can be governed by the cubic nonlinear Schrödinger equation

\displaystyle i \partial_t \phi = \Delta \phi + \lambda |\phi|^2 \phi. \ \ \ \ \ (1)

I recently attended a talk by Natasa Pavlovic on the rigorous derivation of this type of limiting behaviour, which was initiated by the pioneering work of Hepp and Spohn, and has now attracted a vast recent literature. The rigorous details here are rather sophisticated; but the heuristic explanation of the phenomenon is fairly simple, and actually rather pretty in my opinion, involving the foundational quantum mechanics of {N}-particle systems. I am recording this heuristic derivation here, partly for my own benefit, but perhaps it will be of interest to some readers.

This discussion will be purely formal, in the sense that (important) analytic issues such as differentiability, existence and uniqueness, etc. will be largely ignored.

— 1. A quick review of classical mechanics —

The phenomena discussed here are purely quantum mechanical in nature, but to motivate the quantum mechanical discussion, it is helpful to first quickly review the more familiar (and more conceptually intuitive) classical situation.

Classical mechanics can be formulated in a number of essentially equivalent ways: Newtonian, Hamiltonian, and Lagrangian. The formalism of Hamiltonian mechanics for a given physical system can be summarised briefly as follows:

  • The physical system has a phase space {\Omega} of states {\vec x} (which is often parameterised by position variables {q} and momentum variables {p}). Mathematically, it has the structure of a symplectic manifold, with some symplectic form {\omega} (which would be {\omega = dp \wedge dq} if one had position and momentum coordinates available).
  • The complete state of the system at any given time {t} is given (in the case of pure states) by a point {\vec x(t)} in the phase space {\Omega}.
  • Every physical observable (e.g., energy, momentum, position, etc.) {A} is associated to a function (also called {A}) mapping the phase space {\Omega} to the range of the observable (e.g. for real observables, {A} would be a function from {\Omega} to {{\mathbb R}}). If one measures the observable {A} at time {t}, one will obtain the measurement {A(x(t))}.
  • There is a special observable, the Hamiltonian {H: \Omega \rightarrow {\mathbb R}}, which governs the evolution of the state {\vec x(t)} through time, via Hamilton’s equations of motion. If one has position and momentum coordinates {\vec x(t) = (q_i(t), p_i(t))_{i=1}^n}, these equations are given by the formulae

    \displaystyle \partial_t p_i = - \frac{\partial H}{\partial q_i}; \partial_t q_i = \frac{\partial H}{\partial p_i};

    more abstractly, just from the symplectic form {\omega} on the phase space, the equations of motion can be written as

    \displaystyle \partial_t \vec x(t) = - \nabla_\omega H(\vec x(t)), \ \ \ \ \ (2)

    where {\nabla_\omega H} is the symplectic gradient of {H}.

Hamilton’s equation of motion can also be expressed in a dual form in terms of observables {A}, as Poisson’s equation of motion

\displaystyle \partial_t A(\vec x(t)) = - \{ H, A \}(\vec x(t))

for any observable {A}, where {\{H,A\} := \nabla_\omega H \cdot \nabla A} is the Poisson bracket. One can express Poisson’s equation more abstractly as

\displaystyle \partial_t A = -\{H,A\}. \ \ \ \ \ (3)

In the above formalism, we are assuming that the system is in a pure state at each time {t}, which means that it only occupies a single point {\vec x(t)} in phase space. One can also consider mixed states in which the state of the system at a time {t} is not fully known, but is instead given by a probability distribution {\rho(t,\vec x)\ dx} on phase space. The act of measuring an observable {A} at a time {t} will thus no longer be deterministic, but will itself be a random variable, whose expectation {\langle A \rangle} is given by

\displaystyle \langle A \rangle(t) = \int_\Omega A(\vec x) \rho(t,\vec x)\ d\vec x. \ \ \ \ \ (4)

The equation of motion of a mixed state {\rho} is given by the advection equation

\displaystyle \partial_t \rho = \hbox{div}( \rho \nabla_\omega H )

using the same vector field {-\nabla_\omega H} that appears in (2); this equation can also be derived from (3), (4), and a duality argument.

Pure states can be viewed as the special case of mixed states in which the probability distribution {\rho(t,\vec x)\ d\vec x} is a Dirac mass {\delta_{\vec x(t)}(\vec x)}. (We ignore for now the formal issues of how to perform operations such as derivatives on Dirac masses; this can be accomplished using the theory of distributions (or, equivalently, by working in the dual setting of observables) but this is not our concern here.) One can thus think of mixed states as continuous averages of pure states, or equivalently the space of mixed states is the convex hull of the space of pure states.

Suppose one had a {2}-particle system, in which the joint phase space {\Omega = \Omega_1 \times \Omega_2} is the product of the two one-particle phase spaces. A pure joint state is then a point {x = (\vec x_1,\vec x_2)} in {\Omega}, where {\vec x_1} represents the state of the first particle, and {\vec x_2} is the state of the second particle. If the joint Hamiltonian {H: \Omega \rightarrow {\mathbb R}} split as

\displaystyle H(\vec x_1,\vec x_2) = H_1(\vec x_1) + H_2(\vec x_2)

then the equations of motion for the first and second particles would be completely decoupled, with no interactions between the two particles. However, in practice, the joint Hamiltonian contains coupling terms between {\vec x_1, \vec x_2} that prevents one from totally decoupling the system; for instance, one may have

\displaystyle H(\vec x_1,\vec x_2) = \frac{|p_1|^2}{2m_1} + \frac{|p_2|^2}{2m_2} + V(q_1-q_2),

where {\vec x_1=(q_1,p_1)}, {\vec x_2=(q_2,p_2)} are written using position coordinates {q_i} and momentum coordinates {p_i}, {m_1,m_2 > 0} are constants (representing mass), and {V(q_1-q_2)} is some interaction potential that depends on the spatial separation {q_1-q_2} between the two particles.

In a similar spirit, a mixed joint state is a joint probability distribution {\rho(\vec x_1,\vec x_2)\ d\vec x_1 d\vec x_2} on the product state space. To recover the (mixed) state of an individual particle, one must consider a marginal distribution such as

\displaystyle \rho_1(\vec x_1) := \int_{\Omega_2} \rho(\vec x_1,\vec x_2) \ d\vec x_2

(for the first particle) or

\displaystyle \rho_2(\vec x_2) := \int_{\Omega_1} \rho(\vec x_1,\vec x_2) \ d\vec x_1

(for the second particle). Similarly for {N}-particle systems: if the joint distribution of {N} distinct particles is given by {\rho(\vec x_1,\ldots,\vec x_N)\ \vec dx_1 \ldots \vec dx_N}, then the distribution of the first particle (say) is given by

\displaystyle \rho_1(\vec x_1) = \int_{\Omega_2 \times \ldots \times \Omega_N} \rho(\vec x_1,\vec x_2,\ldots,\vec x_N)\ d\vec x_2 \ldots d\vec x_N,

the distribution of the first two particles is given by

\displaystyle \rho_{12}(\vec x_1,\vec x_2) = \int_{\Omega_3 \times \ldots \times \Omega_N} \rho(\vec x_1,\vec x_2,\ldots,\vec x_N)\ d\vec x_3 \ldots d\vec x_N,

and so forth.

A typical Hamiltonian in this case may take the form

\displaystyle H(\vec x_1,\ldots,\vec x_n) = \sum_{j=1}^N \frac{|p_j|^2}{2m_j} + \sum_{1 \leq j < k \leq N} V_{jk}(q_j-q_k)

which is a combination of single-particle Hamiltonians {H_j} and interaction perturbations. If the momenta {p_j} and masses {m_j} are normalised to be of size {O(1)}, and the potential {V_{jk}} has an average value (i.e. an {L^1} norm) of {O(1)} also, then the former sum has size {O(N)} and the latter sum has size {O(N^2)}, so the latter will dominate. In order to balance the two components and get a more interesting limiting dynamics when {N \rightarrow \infty}, we shall therefore insert a normalising factor of {\frac{1}{N}} on the right-hand side, giving a Hamiltonian

\displaystyle H(\vec x_1,\ldots,\vec x_n) = \sum_{j=1}^N \frac{|p_j|^2}{2m_j} + \frac{1}{N} \sum_{1 \leq j < k \leq N} V_{jk}(q_j-q_k).

Now imagine a system of {N} indistinguishable particles. By this, we mean that all the state spaces {\Omega_1 = \ldots = \Omega_N} are identical, and all observables (including the Hamiltonian) are symmetric functions of the product space {\Omega = \Omega_1^N} (i.e. invariant under the action of the symmetric group {S_N}). In such a case, one may as well average over this group (since this does not affect any physical observable), and assume that all mixed states {\rho} are also symmetric. (One cost of doing this, though, is one has to largely give up pure states {(\vec x_1,\ldots,\vec x_N)}, since such states will not be symmetric except in the very exceptional case {\vec x_1=\ldots=\vec x_N}.)

A typical example of a symmetric Hamiltonian is

\displaystyle H(\vec x_1,\ldots,\vec x_n) = \sum_{j=1}^N \frac{|p_j|^2}{2m} + \frac{1}{N} \sum_{1 \leq j < k \leq N} V(q_j-q_k)

where {V} is even (thus all particles have the same individual Hamiltonian, and interact with the other particles using the same interaction potential). In many physical systems, it is natural to consider only short-range interaction potentials, in which the interaction between {q_j} and {q_k} is localised to the region {q_j-q_k=O(r)} for some small {r}. We model this by considering Hamiltonians of the form

\displaystyle H(\vec x_1,\ldots,\vec x_n) = \sum_{j=1}^N H(\vec x_j) + \frac{1}{N} \sum_{1 \leq j < k \leq N} \frac{1}{r^d} V(\frac{\vec x_j-\vec x_k}{r})

where {d} is the ambient dimension of each particle (thus in physical models, {d} would usually be {3}); the factor of {\frac{1}{r^d}} is a normalisation factor designed to keep the {L^1} norm of the interaction potential of size {O(1)}. It turns out that an interesting limit occurs when {r} goes to zero as {N} goes to infinity by some power law {r = N^{-\beta}}; imagine for instance {N} particles of “radius” {r} bouncing around in a box, which is a basic model for classical gases.

An important example of a symmetric mixed state is a factored state

\displaystyle \rho(\vec x_1,\ldots,\vec x_N) = \rho_1(\vec x_1) \ldots \rho_1(\vec x_N)

where {\rho_1} is a single-particle probability density function; thus {\rho} is the tensor product of {N} copies of {\rho_1}. If there are no interaction terms in the Hamiltonian, then Hamiltonian’s equation of motion will preserve the property of being a factored state (with {\rho_1} evolving according to the one-particle equation); but with interactions, the factored nature may be lost over time.

— 2. A quick review of quantum mechanics —

Now we turn to quantum mechanics. This theory is fundamentally rather different in nature than classical mechanics (in the sense that the basic objects, such as states and observables, are a different type of mathematical object than in the classical case), but shares many features in common also, particularly those relating to the Hamiltonian and other observables. (This relationship is made more precise via the correspondence principle, and more precise still using semi-classical analysis.)

The formalism of quantum mechanics for a given physical system can be summarised briefly as follows:

  • The physical system has a phase space {{\Bbb H}} of states {|\psi\rangle} (which is often parameterised as a complex-valued function of the position space). Mathematically, it has the structure of a complex Hilbert space, which is traditionally manipulated using bra-ket notation.
  • The complete state of the system at any given time {t} is given (in the case of pure states) by a unit vector {|\psi(t)\rangle} in the phase space {{\Bbb H}}.
  • Every physical observable {A} is associated to a linear operator on {{\Bbb H}}; real-valued observables are associated to self-adjoint linear operators. If one measures the observable {A} at time {t}, one will obtain the random variable whose expectation {\langle A \rangle} is given by {\langle \psi(t) | A | \psi(t) \rangle}. (The full distribution of {A} is given by the spectral measure of {A} relative to {|\psi(t)\rangle}.)
  • There is a special observable, the Hamiltonian {H: {\Bbb H} \rightarrow {\Bbb H}}, which governs the evolution of the state {|\psi(t)\rangle} through time, via Schrödinger’s equations of motion

    \displaystyle i\hbar \partial_t |\psi(t)\rangle = H |\psi(t) \rangle. \ \ \ \ \ (5)

Schrödinger’s equation of motion can also be expressed in a dual form in terms of observables {A}, as Heisenberg’s equation of motion

\displaystyle \partial_t \langle \psi | A | \psi \rangle = \frac{i}{\hbar} \langle \psi | [ H, A ] | \psi \rangle

or more abstractly as

\displaystyle \partial_t A = \frac{i}{\hbar} [ H, A ] \ \ \ \ \ (6)

where {[,]} is the commutator or Lie bracket (compare with (3)).

The states {|\psi\rangle} are pure states, analogous to the pure states {x} in Hamiltonian mechanics. One also has mixed states {\rho} in quantum mechanics. Whereas in classical mechanics, a mixed state {\rho} is a probability distribution (a non-negative function of total mass {\int_\Omega \rho = 1}), in quantum mechanics a mixed state is a non-negative (i.e. positive semi-definite) operator {\rho} on {{\Bbb H}} of total trace {\hbox{tr} \rho = 1}. If one measures an observable {A} at a mixed state {\rho}, one obtains a random variable with expectation {\hbox{tr} A \rho}. From (6) and duality, one can infer that the correct equation of motion for mixed states must be given by

\displaystyle \partial_t \rho = \frac{i}{\hbar} [H,\rho]. \ \ \ \ \ (7)

One can view pure states as the special case of mixed states which are rank one projections,

\displaystyle \rho= |\psi\rangle \langle \psi|.

Morally speaking, the space of mixed states is the convex hull of the space of pure states (just as in the classical case), though things are a little trickier than this when the phase space {{\Bbb H}} is infinite dimensional, due to the presence of continuous spectrum in the spectral theorem.

Pure states suffer from a phase ambiguity: a phase rotation {e^{i\theta} |\psi\rangle} of a pure state {|\psi \rangle} leads to the same mixed state, and the two states cannot be distinguished by any physical observable.

In a single particle system, modeling a (scalar) quantum particle in a {d}-dimensional position space {{\mathbb R}^d}, one can identify the Hilbert space {{\Bbb H}} with {L^2({\mathbb R}^d \rightarrow {\mathbb C})}, and describe the pure state {|\psi\rangle} as a wave function {\psi: {\mathbb R}^d \rightarrow {\mathbb C}}, which is normalised as

\displaystyle \int_{{\mathbb R}^d} |\psi(x)|^2\ dx = 1

as {|\psi\rangle} has to be a unit vector. (If the quantum particle has additional features such as spin, then one needs a fancier wave function, but let’s ignore this for now.) A mixed state is then a function {\rho: {\mathbb R}^d \times {\mathbb R}^d \rightarrow {\mathbb C}} which is Hermitian (i.e. {\rho(x,x') = \overline{\rho(x',x)}}) and positive definite, with unit trace {\int_{{\mathbb R}^d} \rho(x,x)\ dx = 1}; a pure state {\psi} corresponds to the mixed state {\rho(x,x') = \psi(x) \overline{\psi(x')}}.

A typical Hamiltonian in this setting is given by the operator

\displaystyle H \psi(x) := \frac{|p|^2}{2m} \psi(x) + V(x) \psi(x)

where {m > 0} is a constant, {p} is the momentum operator {p := -i \hbar \nabla_x}, and {\nabla_x} is the gradient in the {x} variable (so {|p|^2 = -\hbar^2 \Delta_x}, where {\Delta_x} is the Laplacian; note that {\nabla_x} is skew-adjoint and should thus be thought of as being imaginary rather than real), and {V: {\mathbb R}^d \rightarrow {\mathbb R}} is some potential. Physically, this depicts a particle of mass {m} in a potential well given by the potential {V}.

Now suppose one has an {N}-particle system of scalar particles. A pure state of such a system can then be given by an {N}-particle wave function {\psi: ({\mathbb R}^d)^N \rightarrow {\mathbb C}}, normalised so that

\displaystyle \int_{({\mathbb R}^d)^N} |\psi(x_1,\ldots,x_N)|^2\ dx_1 \ldots dx_N = 1

and a mixed state is a Hermitian positive semi-definite function {\rho: ({\mathbb R}^d)^N \times ({\mathbb R}^d)^N \rightarrow {\mathbb C}} with trace

\displaystyle \int_{({\mathbb R}^d)^N} \rho(x_1,\ldots,x_N; x_1,\ldots,x_N)\ dx_1 \ldots dx_N = 1,

with a pure state {\psi} being identified with the mixed state

\displaystyle \rho(x_1,\ldots,x_N; x'_1,\ldots,x'_N) := \psi(x_1,\ldots,x_N) \overline{\psi(x'_1,\ldots,x'_N)}.

In classical mechanics, the state of a single particle was the marginal distribution of the joint state. In quantum mechanics, the state of a single particle is instead obtained as the partial trace of the joint state. For instance, the state of the first particle is given as

\displaystyle \rho_1(x_1; x'_1) := \int_{({\mathbb R}^d)^{N-1}} \rho(x_1,x_2,\ldots,x_N; x'_1,x_2,\ldots,x_N)\ dx_2 \ldots dx_N,

the state of the first two particles is given as

\displaystyle \rho_{12}(x_1,x_2; x'_1,x'_2) := \int_{({\mathbb R}^d)^{N-2}} \rho(x_1,x_2,x_3,\ldots,x_N;

\displaystyle x'_1,x'_2,x_3,\ldots,x_N)\ dx_3 \ldots dx_N,

and so forth. (These formulae can be justified by considering observables of the joint state that only affect, say, the first two position coordinates {x_1,x_2} and using duality.)

A typical Hamiltonian in this setting is given by the operator

\displaystyle H \psi(x_1,\ldots,x_N) = \sum_{j=1}^N \frac{|p_j|^2}{2m_j} \psi(x_1,\ldots,x_N)

\displaystyle + \frac{1}{N} \sum_{1 \leq j < k \leq N} V_{jk}(x_j-x_k) \psi(x_1,\ldots,x_N)

where we normalise just as in the classical case, and {p_j := -i \hbar\nabla_{x_j}}.

An interesting feature of quantum mechanics – not present in the classical world – is that even if the {N}-particle system is in a pure state, individual particles may be in a mixed state: the partial trace of a pure state need not remain pure. Because of this, when considering a subsystem of a larger system, one cannot always assume that the subsystem is in a pure state, but must work instead with mixed states throughout, unless there is some reason (e.g. a lack of coupling) to assume that pure states are somehow preserved.

Now consider a system of {N} indistinguishable quantum particles. As in the classical case, this means that all observables (including the Hamiltonian) for the joint system are invariant with respect to the action of the symmetric group {S_N}. Because of this, one may as well assume that the (mixed) state of the joint system is also symmetric with respect to this action. In the special case when the particles are bosons, one can also assume that pure states {|\psi\rangle} are also symmetric with respect to this action (in contrast to fermions, where the action on pure states is anti-symmetric). A typical Hamiltonian in this setting is given by the operator

\displaystyle H \psi(x_1,\ldots,x_N) = \sum_{j=1}^N \frac{|p_j|^2}{2m} \psi(x_1,\ldots,x_N)

\displaystyle + \frac{1}{N} \sum_{1 \leq j < k \leq N} V(x_j-x_k) \psi(x_1,\ldots,x_N)

for some even potential {V}; if one wants to model short-range interactions, one might instead pick the variant

\displaystyle H \psi(x_1,\ldots,x_N) = \sum_{j=1}^N \frac{|p_j|^2}{2m} \psi(x_1,\ldots,x_N) + \frac{1}{N} \sum_{1 \leq j < k \leq N} r^d V(\frac{x_j-x_k}{r}) \psi(x_1,\ldots,x_N) \ \ \ \ \ (8)

for some {r>0}. This is a typical model for an {N}-particle Bose-Einstein condensate. (Longer-range models can lead to more non-local variants of NLS for the limiting equation, such as the Hartree equation.)

— 3. NLS —

Suppose we have a Bose-Einstein condensate given by a (symmetric) mixed state

\displaystyle \rho(t, x_1,\ldots,x_N; x'_1,\ldots,x'_N )

evolving according to the equation of motion (7) using the Hamiltonian (8). One can take a partial trace of the equation of motion (7) to obtain an equation for the state {\rho_1(t, x_1; x'_1)} of the first particle (note from symmetry that all the other particles will have the same state function). If one does take this trace, one soon finds that the equation of motion becomes

\displaystyle \partial_t \rho_1(t,x_1;x'_1) = \frac{i}{\hbar} [ (\frac{|p_1|^2}{2m} - \frac{|p'_1|^2}{2m}) \rho_1(t,x_1;x'_1)

\displaystyle + \frac{1}{N} \sum_{j=2}^N \int_{{\mathbb R}^d} \frac{1}{r^d} [ V( \frac{x_1 - x_j}{r} ) - V( \frac{x'_1 -x_j}{r} ) ] \rho_{1j}(t,x_1,x_j; x'_1,x_j)\ dx_j

where {\rho_{1j}} is the partial trace to the {1, j} particles. Using symmetry, we see that all the summands in the {j} summation are identical, so we can simplify this as

\displaystyle \partial_t \rho_1(t,x_1;x'_1) = \frac{i}{\hbar} [ (\frac{|p_1|^2}{2m} - \frac{|p'_1|^2}{2m}) \rho_1(t,x_1;x'_1)

\displaystyle + \frac{N-1}{N} \int_{{\mathbb R}^d} \frac{1}{r^d} [ V( \frac{x_1 - x_2}{r} ) - V( \frac{x'_1 -x_2}{r} ) ] \rho_{12}(t,x_1,x_2; x'_1,x_2)\ dx_2.

This does not completely describe the dynamics of {\rho_1}, as one also needs an equation for {\rho_{12}}. But one can repeat the same argument to get an equation for {\rho_{12}} involving {\rho_{123}}, and so forth, leading to a system of equations known as the BBGKY hierarchy. But for simplicity we shall just look at the first equation in this hierarchy.

Let us now formally take two limits in the above equation, sending the number of particles {N} to infinity and the interaction scale {r} to zero. The effect of sending {N} to infinity should simply be to eliminate the {\frac{N-1}{N}} factor. The effect of sending {r} to zero should be to send {\frac{1}{r^d} V(\frac{x}{r})} to the Dirac mass {\lambda \delta(x)}, where {\lambda := \int_{{\mathbb R}^d} V} is the total mass of {V}. Formally performing these two limits, one is led to the equation

\displaystyle \partial_t \rho_1(t,x_1;x'_1) = \frac{i}{\hbar} [ (\frac{|p_1|^2}{2m} - \frac{|p'_1|^2}{2m}) \rho_1(t,x_1;x'_1)

\displaystyle + \lambda (\rho_{12}(t,x_1,x_1;x'_1,x_1) - \rho_{12}(t,x_1,x'_1;x'_1,x'_1)) ].

One can perform a similar formal limiting procedure for the other equations in the BBGKY hierarchy, obtaining a system of equations known as the Gross-Pitaevskii hierarchy.

We next make an important simplifying assumption, which is that in the limit {N \rightarrow \infty} any two particles in this system become decoupled, which means that the two-particle mixed state factors as the tensor product of two one-particle states:

\displaystyle \rho_{12}(t,x_1,x_2; x'_1,x_2) \approx \rho_1(t,x_1;x'_1) \rho_1(t,x_2;x'_2).

One can view this as a mean field approximation, modeling the interaction of one particle {x_1} with all the other particles by the mean field {\rho_1}.

Making this assumption, the previous equation simplifies to

\displaystyle \partial_t \rho_1(t,x_1;x'_1) = \frac{i}{\hbar} [ (\frac{|p_1|^2}{2m} - \frac{|p'_1|^2}{2m})

\displaystyle + \lambda (\rho_1(t,x_1;x_1) - \rho_1(t,x'_1;x'_1))] \rho_1(t,x_1;x'_1).

If we assume furthermore that {\rho_1} is a pure state, thus

\displaystyle \rho_1(t,x_1;x'_1) = \psi(t,x_1) \overline{\psi(t,x'_1)}

then (up to the phase ambiguity mentioned earlier), {\psi(t,x)} obeys the Gross-Pitaevskii equation

\displaystyle \partial_t \psi(t,x) = \frac{i}{\hbar} [ (\frac{|p|^2}{2m} + \lambda |\psi(t,x)|^2 ] \psi(t,x)

which (up to some factors of {\hbar} and {m}, which can be renormalised away) is essentially (1).

An alternate derivation of (1), using a slight variant of the above mean field approximation, comes from studying the Hamiltonian (8). Let us make the (very strong) assumption that at some fixed time {t}, one is in a completely factored pure state

\displaystyle \psi(x_1,\ldots,x_N) = \psi_1(x_1) \ldots \psi_1(x_N),

where {\psi_1} is a one-particle wave function, in particular obeying the normalisation

\displaystyle \int_{{\mathbb R}^d} |\psi_1(x)|^2\ dx = 1.

(This is an unrealistically strong version of the mean field approximation. In practice, one only needs the two-particle partial traces to be completely factored for the discussion below.) The expected value of the Hamiltonian,

\displaystyle \langle \psi|H|\psi\rangle = \int_{({\mathbb R}^d)^N} \psi(x_1,\ldots,x_N) \overline{H \psi(x_1,\ldots,x_N)}\ dx_1 \ldots dx_N,

can then be simplified as

\displaystyle N \int_{{\mathbb R}^d} \psi_1(x) \overline{\frac{|p_1|^2}{2m} \psi_1(x)}\ dx

\displaystyle + \frac{N-1}{2} \int_{{\mathbb R}^d \times {\mathbb R}^d} r^{-d} V(\frac{x_1-x_2}{r}) |\psi_1(x_1)|^2 |\psi_1(x_2)|\ dx_1 dx_2.

Again sending {r \rightarrow 0}, this formally becomes

\displaystyle N \int_{{\mathbb R}^d} \psi_1(x) \overline{\frac{|p_1|^2}{2m} \psi_1(x)}\ dx + \frac{N-1}{2} \lambda \int_{{\mathbb R}^d \times {\mathbb R}^d} |\psi_1(x_1)|^4\ dx_1

which in the limit {N \rightarrow \infty} is asymptotically

\displaystyle N \int_{{\mathbb R}^d} \psi_1(x) \overline{\frac{|p_1|^2}{2m} \psi_1(x)} + \frac{\lambda}{2} |\psi_1(x_1)|^4\ dx_1.

Up to some normalisations, this is the Hamiltonian for the NLS equation (1).

There has been much progress recently in making the above derivations precise, by Erdös-Schlein-Yau, Klainerman-Machedon, Kirkpatrick-Schlein-Staffilani, Chen-Pavlovic, and others. A key step is to show that the Gross-Pitaevskii hierarchy necessarily preserves the property of being a completely factored state. This requires a uniqueness theory for this hierarchy, which is surprisingly delicate, due to the fact that it is a system of infinitely many coupled equations over an unbounded number of variables.

[Update, Dec 8: Interestingly, the above heuristic derivation only works when the interaction scale {r} is much larger than {N^{-1}}. For {r \sim N^{-1}}, the coupling constant {\lambda} acquires a nonlinear correction, becoming essentially the scattering length of the potential rather than its mean. See comments below.]