Centre-of-mass motion in multi-particle Schrödinger–Newton dynamics

How does a quantum system in a non-classical state gravitate? There is no unanimously accepted answer to this seemingly obvious question. If we assume that gravity is fundamentally quantum, as most physicists assume, the fairest answer is simply that we do not know. If gravity stays fundamentally classical, a perhaps less likely but not altogether outrageous possibility [1, 2], we also do not know; but we can guess. One such guess is that semi-classical gravity stays valid, beyond the realm it would be meant for if gravity were quantum [1, 2]. Semi-classical gravity in that extended sense is the theory that we wish to pursue in this paper. Since eventually we are aiming for the characterization of experimentally testable consequences of such gravitational self-interaction through matter-wave interferometry, we focus attention on the centre-of-mass motion.

Note that by 'quantum system' we refer to the possibility for the system to assume states that have no classical counterpart, like superpositions of spatially localized states. We are not primarily interested in matter under extreme conditions (energy, pressure, etc). Rather we are interested in ordinary laboratory matter described by non-relativistic quantum mechanics, whose states will source a classical gravitational field according to semi-classical equations. Eventually we are interested in the question concerning the range of validity of such equations. Since we do not exclude the possibility that gravity might stay classical at the most fundamental level, we explicitly leave open the possibility that these equations stay valid even for strongly fluctuating states of matter.

Now, if we assume that a one-particle state $\psi$ gravitates like a classical mass density $\tilde{\rho }\left( {\bf x} \right)=m{{\left| \psi \left( {\bf x} \right) \right|}^{2}}$, we immediately get the coupled equations (neglecting other external potentials for simplicity)

$$\begin{eqnarray}&&{\rm i}\hbar {{\partial }_{t}}\psi (t;{\bf x})=\left( -\frac{{{\hbar }^{2}}}{2\;m}\Delta +{{V}_{{\rm g}}}(t;{\bf x}) \right)\psi (t;{\bf x}),\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&\Delta {{V}_{{\rm g}}}(t;{\bf x})=4\pi \;G\;{{m}^{2}}\;{{\left| \psi (t;{\bf x}) \right|}^{2}}.\end{eqnarray} \tag{ 1b }$$

These equations are known as the (one-particle) Schrödinger–Newton system. This system can be transformed into a single, nonlinear and non-local equation for $\psi$ by first solving (1b ) with the boundary condition that $\phi$ be zero at spatial infinity, which leads to

$$\begin{eqnarray}&&{{V}_{{\rm g}}}\left( t;{\bf x} \right)=-G{{m}^{2}}\int \frac{{{\left| \psi \left( t;{\bf x}^{\prime} \right) \right|}^{2}}}{\parallel {\bf x}-{\bf x}^{\prime} \parallel }\;{{{\rm d}}^{3}}{\bf x}^{\prime} \;.\end{eqnarray} \tag{ 2 }$$

Inserting (2) into (1a ) results in the one-particle Schrödinger–Newton equation:

$$\begin{eqnarray}&&{\rm i}\hbar {{\partial }_{t}}\psi \left( t;{\bf x} \right)=\left( -\frac{{{\hbar }^{2}}}{2\;m}\Delta -G{{m}^{2}}\int \frac{{{\left| \psi \left( t;{\bf x}^{\prime} \right) \right|}^{2}}}{\parallel {\bf x}-{\bf x}^{\prime} \parallel }\;{{{\rm d}}^{3}}{\bf x}^{\prime} \right)\psi \left( t;{\bf x} \right).\end{eqnarray} \tag{ 3 }$$

Concerning the theoretical foundation of (3), the nonlinear self-interaction should essentially be seen as a falsifiable hypothesis on the gravitational interaction of matter fields, where the reach of this hypothesis delicately depends on the kind of 'fields' it is supposed to cover. For example, (3) has been shown to follow in a suitable non-relativistic limit from the Einstein–Klein–Gordon or Einstein–Dirac systems [3], i.e., systems where the energy–momentum tensor ${{T}_{\mu \nu }}$ on the right-hand side of Einsteinʼs equations,

$$\begin{eqnarray}&&{{R}_{\mu \nu }}-\frac{1}{2}{{g}_{\mu \nu }}R=\frac{8\pi G}{{{c}^{4}}}\;{{T}_{\mu \nu }}\;,\end{eqnarray} \tag{ 4 }$$

is built from classical Klein–Gordon or classical Dirac fields. Such an expression for ${{T}_{\mu \nu }}$ results from the expectation value $\langle \psi |{{\hat{T}}_{\mu \nu }}|\psi \rangle$ in quantum-field theory, where $\psi$ labels the amplitude (wave function) of a one-particle state, ${{\hat{T}}_{\mu \nu }}$ is the operator-valued energy–momentum tensor, which has been suitably regularized ¹ . The non-relativistic limit is then simply the (regularized) mass density operator whose expectation value in a one-particle state is $m{{\left| \psi \right|}^{2}}$; see e.g. [5].

Now, if we believe that there exists an underlying quantum theory of gravity of which the semi-classical Einstein equation (4) with ${{T}_{\mu \nu }}$ replaced by $\langle \psi |{{\hat{T}}_{\mu \nu }}|\psi \rangle$ is only an approximation, then this will clearly only make sense in situations where the source-field for gravity, which is an operator, may be replaced by its mean-field approximation. This is the case in many-particle situations, i.e., where $\psi$ is a many-particle amplitude, and then only in the limit as the particle number tends to infinity. From that perspective it would make little sense to use one-particle expectation values on the right-hand side of Einsteinʼs equation, for their associated classical gravitational field according to (4) will not be any reasonable approximation of the (strongly fluctuating) fundamentally quantum gravitational field. This has been rightfully stressed recently [5, 6].

On the other hand, if we consider the possibility that gravity stays fundamentally classical, as we wish to do so here, then we are led to contemplate the strict (and not just approximate) sourcing of gravitational fields by expectation values rather than operators. In this case we do get nonlinear self-interactions due to gravity in the equations, even for the one-particle amplitudes. Note that it would clearly not be proper to regard these amplitudes as classical fields and once more (second) quantize them. This is an important conceptual point that seems to have caused some confusion recently. We will therefore briefly return to this issue at the end of section 2. Also recall that the often alleged existing evidence, experimental [7] or conceptual [8], is generally found inconclusive, e.g., [9–11].

Taken as a new hypothesis for the gravitational interaction of matter, the Schrödinger–Newton equation has attracted much attention in recent years. First of all, it raises the challenge to experimentally probe the consequences of the nonlinear gravitational self-interaction term [12]. More fundamentally, the verification of the existence of this semi-classical self-interaction could shed new light on the holy grail of theoretical physics: $quantum\:gravity$ and its alleged necessity; compare [2]. And even though the original numerical estimates made in [12] were too optimistic by many orders of magnitude, there is now consensus as to the prediction of (3) concerning gravity-induced inhibition of quantum-mechanical dispersion [13].

However, concerning the current and planned interference experiments, it must be stressed that they are made with extended objects, like large molecules or tiny 'nanospheres' [14], and that the so-called 'large superpositions' concern only the centre-of-mass part of the overall multi-particle wavefunction. But even if we assume the elementary constituents in isolation to obey (3), there is still no obvious reason why the centre of mass of a compound object would obey a similar equation. These equations are nonlinear and 'separating off' degrees of freedom is not as obvious a procedure as in the linear case. The study of this issue is the central concern of this paper. For this we start afresh from a multi-particle version of the Schrödinger–Newton equation.

In this paper we consider the $\left( N+1 \right)$-particle Schrödinger–Newton equation for a function $\Psi :{{\mathbb{R}}^{1+3\left( N+1 \right)}}\to \mathbb{C}$, where $3\left( N+1 \right)$ arguments correspond to the three coordinates each of $\left( N+1 \right)$ particles of masses ${{m}_{0}},{{m}_{1}},\cdots ,{{m}_{N}}$, and one argument is given by the (Newtonian absolute) time t. In the presence of non-gravitational two-body interactions represented by potentials ${{V}_{ab}}(\parallel {{{\bf x}}_{a}}-{{{\bf x}}_{b}}\parallel )$, where ${{V}_{ab}}={{V}_{ba}}$, for the pair labelled by $\left( ab \right)$, the $\left( N+1 \right)$-particle Schrödinger–Newton equation reads in full glory

$$\begin{eqnarray}&&\begin{array}{l} {\rm i}\hbar {{\partial }_{t}}\Psi \left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right)=\left( -\mathop{\sum }\limits_{a=0}^{N}\frac{{{\hbar }^{2}}}{2{{m}_{a}}}{{\Delta }_{a}}+\mathop{\sum }\limits_{a=0}^{N}\mathop{\sum }\limits_{b>a}^{N}{{V}_{ab}}\left( \parallel {{{\bf x}}_{a}}-{{{\bf x}}_{b}}\parallel \right) \right. \\ \quad \left. -G\mathop{\sum }\limits_{a=0}^{N}\mathop{\sum }\limits_{b=0}^{N}{{m}_{a}}{{m}_{b}}\left\{ \int \mathop{\prod }\limits_{c=0}^{N}{{{\rm d}}^{3}}{\bf x}_{c}^{^{\prime} } \right\}\frac{{{\left| \Psi \left( t;{\bf x}_{0}^{^{\prime} },\cdots ,{\bf x}_{N}^{^{\prime} } \right) \right|}^{2}}}{\parallel {{{\bf x}}_{a}}-{\bf x}_{b}^{^{\prime} }\parallel } \right)\Psi \left( t;{{{\bf x}}_{0}},\cdots {{{\bf x}}_{N}} \right).\qquad \; \\ \end{array}\end{eqnarray} \tag{ 5 }$$

Here and in the following, we write

$$\begin{eqnarray}&&{{{\rm d}}^{3}}{{{\bf x}}_{c}}:={\rm d}x_{c}^{1}\wedge {\rm d}x_{c}^{2}\wedge {\rm d}x_{c}^{3}\quad {\rm and}\quad \mathop{\prod }\limits_{c=0}^{N}{{{\rm d}}^{3}}{{{\bf x}}_{c}}:={{{\rm d}}^{3}}{{{\bf x}}_{0}}\wedge \cdots \wedge {{{\rm d}}^{3}}{{{\bf x}}_{N}}.\end{eqnarray} \tag{ 6 }$$

The second, nonlinear and non-local potential term is meant to represent the gravitational interaction according to a suggestion first made in [15]. The structure of this term seems rather complicated, but the intuition behind it is fairly simple:

Assumption 1. Each particle represents a mass distribution in physical space that is proportional to its marginal distribution derived from $\Psi \left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right)$. More precisely, the mass distribution represented by the bth particle is

$$\begin{eqnarray}\begin{array}{rcl} {{\tilde{\rho }}_{b}}\left( t;{\bf x} \right) & = & {{m}_{b}}\left\{ \int \mathop{\prod }\limits_{\begin{array}{ccccccccccccccc} c=0 \\ c\ne b \\ \end{array}}^{N}{{{\rm d}}^{3}}{{{\bf x}}_{c}} \right\}\ {{\left| \Psi \left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{b-1}},{\bf x},{{{\bf x}}_{b+1}},\cdots ,{{{\bf x}}_{N}} \right) \right|}^{2}} \\ {} & = & {{m}_{b}}\left\{ \int \mathop{\prod }\limits_{c=0}^{N}{{{\rm d}}^{3}}{{{\bf x}}_{c}} \right\}\ {{\delta }^{\left( 3 \right)}}\left( {\bf x}-{{{\bf x}}_{b}} \right){{\left| \Psi \left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right) \right|}^{2}}. \\ \end{array}\end{eqnarray} \tag{ 7 }$$

Assumption 2. The total gravitational potential $\Phi$ at ${\bf x}$ in physical space is that generated by the sum of the mass distributions (7) according to Newtonian gravity. More precisely, the Newtonian gravitational potential is given by

$$\begin{eqnarray}&&\Phi \left( t;{\bf x} \right)=-G\int {{{\rm d}}^{3}}{\bf x}^{\prime} \frac{\mathop{\sum }\limits_{b=0}^{N}\;{{\tilde{\rho }}_{b}}\left( t;{\bf x}^{\prime} \right)}{\parallel {\bf x}-{\bf x}^{\prime} \parallel }.\end{eqnarray} \tag{ 8 }$$

Assumption 3. The gravitational contribution ${{V}_{{\rm g}}}\left( {{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right)$ that enters the Hamiltonian in the multi-particle Schrödinger equation

$$\begin{eqnarray}\begin{array}{rcl} {\rm i}\hbar {{\partial }_{t}}\psi \left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right) & = & \left( -\mathop{\sum }\limits_{a=0}^{N}\frac{{{\hbar }^{2}}}{2{{m}_{a}}}{{\Delta }_{a}}+{{V}_{{\rm other}}}\left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right)+{{V}_{{\rm g}}}\left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right) \right) \\ {} & {} & \times \Psi \left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right).\qquad \qquad \qquad \qquad \; \\ \end{array}\end{eqnarray} \tag{ 9 }$$

is the sum of the gravitational potential energies of $\left( N+1 \right)$ point-particles (sic!) of masses ${{m}_{a}}$ situated at positions ${{{\bf x}}_{a}}$. More precisely, the total gravitational contribution to the Hamiltonian is

$$\begin{eqnarray}&&{{V}_{{\rm g}}}\left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right)=\mathop{\sum }\limits_{a=0}^{N}{{m}_{a}}\Phi \left( t;{{{\bf x}}_{a}} \right),\end{eqnarray} \tag{ 10 }$$

where $\Phi$ is given by (8).

Taken together, all three assumptions result in a gravitational contribution to the Hamiltonian of

$$\begin{eqnarray}&&{{V}_{{\rm g}}}\left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right)=-G\mathop{\sum }\limits_{a=0}^{N}\mathop{\sum }\limits_{b=0}^{N}{{m}_{a}}{{m}_{b}}\left\{ \int \mathop{\prod }\limits_{c=0}^{N}{{{\rm d}}^{3}}{\bf x}_{c}^{^{\prime} } \right\}\frac{{{\left| \Psi \left( t;{\bf x}_{0}^{^{\prime} },\cdots ,{\bf x}_{N}^{^{\prime} } \right) \right|}^{2}}}{\parallel {{{\bf x}}_{a}}-{\bf x}_{b}^{^{\prime} }\parallel }\end{eqnarray} \tag{ 11 }$$

just as in (5). We note that (5) can be derived from a Lagrangian

$$\begin{eqnarray}&&L=T-U,\end{eqnarray} \tag{ 12 }$$

where the kinetic part ² , T, is

$$\begin{eqnarray}\begin{array}{rcl} T & = & \frac{{\rm i}\hbar }{2}\left\{ \int \mathop{\prod }\limits_{a=0}^{N}{{{\rm d}}^{3}}{{{\bf x}}_{a}} \right\}\left( \bar{\Psi }{{\partial }_{t}}\Psi -\Psi {{\partial }_{t}}\bar{\Psi } \right) \\ {} & {} & +{{\hbar }^{2}}\left\{ \int \mathop{\prod }\limits_{a=0}^{N}{{{\rm d}}^{3}}{{{\bf x}}_{a}} \right\}\mathop{\sum }\limits_{b=0}^{N}\frac{1}{{{m}_{b}}}{{\nabla }_{b}}\bar{\Psi }\cdot {{\nabla }_{b}}\Psi .\qquad \qquad \qquad \qquad \qquad \qquad \; \\ \end{array}\end{eqnarray} \tag{ 13 }$$

Here all functions are taken at the same argument $\left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right)$, which we suppressed. The potential part, U, consists of a sum of two terms. The first term represents possibly existent two-body interactions, such as, e.g., electrostatic energy:

$$\begin{eqnarray}&&{{U}^{{\rm local}2-{\rm body}}}=\left\{ \int \mathop{\prod }\limits_{c=0}^{N}{{{\rm d}}^{3}}{{{\bf x}}_{c}} \right\}\mathop{\sum }\limits_{a=0}^{N}\mathop{\sum }\limits_{b>a}^{N}{{V}_{ab}}\left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right){{\left| \Psi \left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right) \right|}^{2}}.\end{eqnarray} \tag{ 14 }$$

The second contribution is that of gravity:

$$\begin{eqnarray}\begin{array}{rcl} {{U}^{{\rm grav}}} & = & -\frac{G}{2}\left\{ \int \mathop{\prod }\limits_{c=0}^{N}{{{\rm d}}^{3}}{{{\bf x}}_{c}} \right\}\left\{ \int \mathop{\prod }\limits_{d=0}^{N}{{{\rm d}}^{3}}{\bf x}_{d}^{^{\prime} } \right\}\mathop{\sum }\limits_{a=0}^{N}\mathop{\sum }\limits_{b=0}^{N}{{m}_{a}}{{m}_{b}} \\ {} & {} & \times \frac{{{\left| \Psi \left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right) \right|}^{2}}{{\left| \Psi \left( t;{\bf x}_{0}^{^{\prime} },\cdots ,{\bf x}_{N}^{^{\prime} } \right) \right|}^{2}}}{\parallel {{{\bf x}}_{a}}-{\bf x}_{b}^{^{\prime} }\parallel } \\ {} & = & -\frac{G}{2}\mathop{\sum }\limits_{a=0}^{N}\mathop{\sum }\limits_{b=0}^{N}\int {{{\rm d}}^{3}}{\bf x}\int {{{\rm d}}^{3}}{\bf x}^{\prime} \quad \frac{{{\tilde{\rho }}_{a}}\left( {\bf x} \right){{\tilde{\rho }}_{b}}\left( {\bf x}^{\prime} \right)}{\parallel {\bf x}-{\bf x}^{\prime} \parallel }.\qquad \qquad \qquad \qquad \qquad \qquad \; \\ \end{array}\end{eqnarray} \tag{ 15 }$$

The last line shows that the gravitational energy is just the usual binding energy of $\left( N+1 \right)$ lumps of matter distributed in physical space according to (7). Note that the sum not only contains the energies for the mutual interactions between the lumps, but also the self-energy of each lump. The latter are represented by the diagonal terms in the double sum, i.e. the terms where a = b. These self-energy contributions would diverge for pointlike mass distributions, i.e. if ${{\tilde{\rho }}_{a}}\left( {\bf x} \right)={{m}_{a}}{{\delta }^{\left( 3 \right)}}\left( {\bf x}-{{{\bf x}}_{a}} \right)$, as in the case of electrostatic interaction (see below). Here, however, the hypotheses underlying the three assumptions above imply that gravitationally the particles interact differently, resulting in finite self-energies. Because of these self-energies we already obtain a modification of the ordinary Schrödinger equation in the one-particle case, which is just given by (3). Explicit expressions for the double integrals over ${{\tilde{\rho }}_{a}}\left( {\bf x} \right){{\tilde{\rho }}_{b}}\left( {\bf x}^{\prime} \right)/\parallel {\bf x}-{\bf x}^{\prime} \parallel$ can, e.g., be found in [16] for some special cases where ${{\tilde{\rho }}_{a}}$ and ${{\tilde{\rho }}_{b}}$ are spherically symmetric.

Finally we wish to come back to the fundamental issue already touched upon in the introduction, namely how to relate the interaction term (15) to known physics as currently understood. As already emphasized in the context of (3), i.e. for just one particle, the gravitational interaction contains self-energy contributions. In the multi-particle scheme they just correspond to the diagonal terms a = b in (15). These terms are certainly finite for locally bounded ${{\tilde{\rho }}_{a}}$.

This would clearly not be the case in a standard quantum field-theoretic treatment, like QED, outside the mean-field limit. In non-relativistic quantum field theory the interaction Hamiltonian would be a double integral over ${{\Psi }^{\dagger }}\left( {\bf x} \right)\Psi \left( {\bf x} \right){{\Psi }^{\dagger }}\left( {\bf x}^{\prime} \right)\Psi \left( {\bf x}^{\prime} \right)/\parallel {\bf x}-{\bf x}^{\prime} \parallel$, where $\Psi$ is the (non-relativistic) field operator. (See, e.g., chapter 11 of [17] for a text-book account of non-relativistic QFT.) This term will lead to divergent self-energies, which one renormalizes through normal ordering, and pointwise Coulomb interactions of pairs. This is just the known and accepted strategy followed in deriving the multi-particle Schrödinger equation for charged point-particles from QED. This procedure has a long history. In fact, it can already be found in the appendix of Heisenbergʼs 1929 Chicago lectures [18] on quantum mechanics.

It has therefore been frequently complained that the Schrödinger–Newton equation does not follow from 'known physics' [5, 6, 19, 20]. This is true, of course. But note that this does not imply the sharper argument according to which the Schrödinger–Newton equation even contradicts known physics. Such sharper arguments usually beg the question by assuming some form of quantum gravity to exist. But this hypothetical theory is not yet part of 'known physics' either, and may never be! Similarly, by rough analogy of the classical fields in gravity and electromagnetism, the Schrödinger–Newton equation is sometimes argued to contradict known physics because the analogous nonlinear 'Schrödinger–Coulomb' equation yields obvious nonsense, like a grossly distorted energy spectrum for hydrogen. In fact, this was already observed in 1927 by Schrödinger, who wondered about this factual contradiction with what he described as a natural demand (self-coupling) from a classical field-theoretic point of view [21]. Heisenberg in his 1929 lectures also makes this observation, which he takes as irrefutable evidence for the need to (second) quantize the Schrödinger field, thereby turning a nonlinear 'classical' field theory into a linear quantum version of it.

To say it once more: all this is only an argument against the Schrödinger–Newton equation provided we assume an underlying theory of quantum gravity to exist and whose effective low energy approximation can be dealt with in full analogy to, say, QED. But our attitude here is different! What we have is a hypothesis that is essentially based on the assumption that gravity behaves differently as regards its coupling to matter and, in particular, its need for quantization. The interesting aspect of this is that it gives rise to potentially observable consequences that render this hypothesis falsifiable.

Instead of the $\left( N+1 \right)$ positions ${{{\bf x}}_{a}}$, $a=0,\cdots ,N$, in absolute space, we introduce the centre of mass and N positions relative to it. We write

$$\begin{eqnarray}&&M:=\mathop{\sum }\limits_{a=0}^{N}{{m}_{a}}\end{eqnarray} \tag{ 16 }$$

for the total mass and adopt the convention that greek indices $\alpha ,\beta ,\cdots$ take values in $\left\{ 1,\cdots ,N \right\}$, in contrast to latin indices $a,b,\cdots$, which we already agreed to take values in $\left\{ 0,1,\cdots ,N \right\}$. The centre-of-mass and the relative coordinates of the N particles labelled by $1,\cdots ,N$ are given by (thereby distinguishing the particle labelled by 0)

$$\begin{eqnarray}&&{\bf c}:=\frac{1}{M}\mathop{\sum }\limits_{a=0}^{N}{{m}_{a}}\;{{{\bf x}}_{a}}=\frac{{{m}_{0}}}{M}{{{\bf x}}_{0}}+\mathop{\sum }\limits_{\beta =1}^{N}\frac{{{m}_{\beta }}}{M}{{{\bf x}}_{\beta }}\;,\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&{{{\bf r}}_{\alpha }}\;:={{{\bf x}}_{\alpha }}-{\bf c}=-\frac{{{m}_{0}}}{M}\;{{{\bf x}}_{0}}+\mathop{\sum }\limits_{\beta =1}^{N}\left( {{\delta }_{\alpha \beta }}-\frac{{{m}_{\beta }}}{M} \right){{{\bf x}}_{\beta }}.\end{eqnarray} \tag{ 17b }$$

The inverse transformation is obtained by simply solving (17) for ${{{\bf x}}_{0}}$ and ${{{\bf x}}_{\alpha }}$:

$$\begin{eqnarray}&&{{{\bf x}}_{0}}={\bf c}-\mathop{\sum }\limits_{\beta =1}^{N}\frac{{{m}_{\beta }}}{{{m}_{0}}}\;{{{\bf r}}_{\beta }}\;,\end{eqnarray} \tag{ 18a }$$

$$\begin{eqnarray}&&{{{\bf x}}_{\alpha }}={\bf c}+{{{\bf r}}_{\alpha }}\;.\end{eqnarray} \tag{ 18b }$$

All this may be written in a self-explanatory $\left( 1+N \right)$ split matrix form

$$\begin{eqnarray}&&(\begin{array}{ccccccccccccccc} {\bf c} \\ {{{\bf r}}_{\alpha }} \\ \end{array})=(\begin{array}{ccccccccccccccc} \frac{{{m}_{0}}}{M} & \frac{{{m}_{\beta }}}{M} \\ -\frac{{{m}_{0}}}{M} & {{\delta }_{\alpha \beta }}-\frac{{{m}_{\beta }}}{M} \\ \end{array})(\begin{array}{ccccccccccccccc} {{{\bf x}}_{0}} \\ {{{\bf x}}_{\beta }} \\ \end{array})\;,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&(\begin{array}{ccccccccccccccc} {{{\bf x}}_{0}} \\ {{{\bf x}}_{\alpha }} \\ \end{array})=(\begin{array}{ccccccccccccccc} 1 & -\frac{{{m}_{\beta }}}{{{m}_{0}}} \\ 1 & {{\delta }_{\alpha \beta }} \\ \end{array})(\begin{array}{ccccccccccccccc} {\bf c} \\ {{{\bf r}}_{\beta }} \\ \end{array})\;.\end{eqnarray} \tag{ 20 }$$

For the wedge product of the $\left( N+1 \right)$ 1-forms ${\rm d}x_{a}^{1}$ for $a=0,1,\cdots ,N$ we easily get from (18)

$$\begin{eqnarray}\begin{array}{rcl} {\rm d}x_{0}^{1}\wedge \cdots \wedge {\rm d}x_{N}^{1} & = & \left( {\rm d}{{c}^{1}}-\mathop{\sum }\limits_{\beta =1}^{N}\frac{{{m}_{\beta }}}{{{m}_{0}}}{\rm d}r_{\beta }^{1} \right)\wedge \left( {\rm d}{{c}^{1}}+{\rm d}r_{1}^{1} \right)\wedge \cdots \wedge \left( {\rm d}{{c}^{1}}+{\rm d}r_{N}^{1} \right) \\ {} & = & \frac{M}{{{m}_{0}}}\left( {\rm d}{{c}^{1}}\wedge {\rm d}r_{1}^{1}\wedge \cdots \wedge {\rm d}r_{N}^{1} \right). \\ \end{array}\end{eqnarray} \tag{ 21 }$$

Hence, writing ${{{\rm d}}^{3}}{{{\bf r}}_{\alpha }}:={\rm d}r_{\alpha }^{1}\wedge {\rm d}r_{\alpha }^{2}\wedge {\rm d}r_{\alpha }^{3}$ and $\prod _{\alpha =1}^{N}$ for the N-fold wedge product, we have

$$\begin{eqnarray}&&\mathop{\prod }\limits_{a=0}^{N}{{{\rm d}}^{3}}{{{\bf x}}_{a}}={{\left( \frac{M}{{{m}_{0}}} \right)}^{3}}\left( {{{\rm d}}^{3}}{\bf c}\wedge \mathop{\prod }\limits_{\alpha =1}^{N}{{{\rm d}}^{3}}{{{\bf r}}_{\alpha }} \right).\end{eqnarray} \tag{ 22 }$$

Note that the sign changes that may appear in rearranging the wedge products on both sides coincide and hence cancel. From (22) we can just read off the determinant of the Jacobian matrix for the transformation (18):

$$\begin{eqnarray}&&\left| \frac{\partial \left( {{{\bf x}}_{0}},{{{\bf x}}_{\alpha }} \right)}{\partial \left( {\bf c},{{{\bf r}}_{\beta }} \right)} \right|:={\rm det} \left\{ \frac{\partial \left( {{{\bf x}}_{0}},{{{\bf x}}_{\alpha }} \right)}{\partial \left( {\bf c},{{{\bf r}}_{\beta }} \right)} \right\}={{\left( \frac{M}{{{m}_{0}}} \right)}^{3}}\;.\end{eqnarray} \tag{ 23 }$$

Equation (18) also allows to simply rewrite the kinetic-energy metric

$$\begin{eqnarray}&&G=\mathop{\sum }\limits_{a=0}^{N}\mathop{\sum }\limits_{b=0}^{N}{{G}_{ab}}\;{\rm d}{{{\bf x}}_{a}}\otimes {\rm d}{{{\bf x}}_{b}}:=\mathop{\sum }\limits_{a=0}^{N}{{m}_{a}}\;{\rm d}{{{\bf x}}_{a}}\otimes {\rm d}{{{\bf x}}_{a}}\end{eqnarray} \tag{ 24 }$$

in terms of the new coordinates: it is given by

$$\begin{eqnarray}\begin{array}{rcl} G & = & {{m}_{0}}\left( {\rm d}{\bf c}-\mathop{\sum }\limits_{\alpha =1}^{N}\frac{{{m}_{\alpha }}}{{{m}_{0}}}{\rm d}{{{\bf r}}_{\alpha }} \right)\otimes \left( {\rm d}{\bf c}-\mathop{\sum }\limits_{\beta =1}^{N}\frac{{{m}_{\beta }}}{{{m}_{0}}}{\rm d}{{{\bf r}}_{\beta }} \right)+\mathop{\sum }\limits_{\alpha =1}^{N}{{m}_{\alpha }}\left( {\rm d}{\bf c}+{\rm d}{{{\bf r}}_{\alpha }} \right)\otimes \left( {\rm d}{\bf c}+{\rm d}{{{\bf r}}_{\alpha }} \right) \\ {} & = & M\;{\rm d}{\bf c}\otimes {\rm d}{\bf c}+\mathop{\sum }\limits_{\alpha =1}^{N}\mathop{\sum }\limits_{\beta =1}^{N}{{H}_{\alpha \beta }}\;{\rm d}{{{\bf r}}_{\alpha }}\otimes {\rm d}{{{\bf r}}_{\beta }}. \\ \end{array}\end{eqnarray} \tag{ 25 }$$

The first thing to note is that there are no off-diagonal terms, i.e. terms involving tensor products between ${\rm d}{\bf c}$ and ${\rm d}{{{\bf r}}_{\alpha }}$. This means that the degrees of freedom labelled by our ${{{\bf r}}_{a}}$ coordinates are perpendicular (with respect to the kinetic-energy metric) to the centre-of-mass motion. The restriction of the kinetic-energy metric to the relative coordinates has the components

$$\begin{eqnarray}&&{{H}_{\alpha \beta }}=\left( \frac{{{m}_{\alpha }}{{m}_{\beta }}}{{{m}_{0}}}+{{m}_{\alpha }}{{\delta }_{\alpha \beta }} \right).\end{eqnarray} \tag{ 26 }$$

The determinant of $\{{{H}_{\alpha \beta }}\}$ follows from taking the determinant of the transformation formula for the kinetic-energy metric (taking due account of the three-fold multiplicities hidden in the inner products in ${{\mathbb{R}}^{3}}$)

$$\begin{eqnarray}&&{{\left( {\rm det} \left\{ {{G}_{ab}} \right\} \right)}^{3}}\times {{\left| \frac{\partial \left( {{{\bf x}}_{0}},{{{\bf x}}_{\alpha }} \right)}{\partial \left( {\bf c},{{{\bf r}}_{\beta }} \right)} \right|}^{2}}={{M}^{3}}\times {{\left( {\rm det} \left\{ {{H}_{\alpha \beta }} \right\} \right)}^{3}}\end{eqnarray} \tag{ 27 }$$

which, using (23) and ${\rm det} \{{{G}_{ab}}\}=\prod _{a=0}^{N}({{m}_{a}}/2)$, results in

$$\begin{eqnarray}&&{\rm det} \left\{ {{H}_{\alpha \beta }} \right\}=\frac{M}{m_{0}^{2}}\;\mathop{\prod }\limits_{a=0}^{N}{{m}_{a}}.\end{eqnarray} \tag{ 28 }$$

Finally we consider the inverse of the kinetic-energy metric:

$$\begin{eqnarray}&&{{G}^{-1}}=\mathop{\sum }\limits_{a=0}^{N}\mathop{\sum }\limits_{b=0}^{N}{{G}^{ab}}\frac{\partial }{\partial {{{\bf x}}_{a}}}\otimes \frac{\partial }{\partial {{{\bf x}}_{b}}}=\mathop{\sum }\limits_{a=0}^{N}\frac{1}{{{m}_{a}}}\frac{\partial }{\partial {{{\bf x}}_{a}}}\otimes \frac{\partial }{\partial {{{\bf x}}_{a}}}\end{eqnarray} \tag{ 29 }$$

Using (17) we have

$$\begin{eqnarray}&&\frac{\partial }{\partial {{{\bf x}}_{0}}}=\frac{{{m}_{0}}}{M}\left( \frac{\partial }{\partial {\bf c}}-\mathop{\sum }\limits_{\alpha =1}^{N}\frac{\partial }{\partial {{{\bf r}}_{\alpha }}} \right),\end{eqnarray} \tag{ 30a }$$

$$\begin{eqnarray}&&\frac{\partial }{\partial {{{\bf x}}_{\alpha }}}=\frac{\partial }{\partial {{{\bf r}}_{\alpha }}}+\frac{{{m}_{\alpha }}}{M}\left( \frac{\partial }{\partial {\bf c}}-\mathop{\sum }\limits_{\beta =1}^{N}\frac{\partial }{\partial {{{\bf r}}_{\beta }}} \right).\end{eqnarray} \tag{ 30b }$$

Inserting this into (29) we obtain the form

$$\begin{eqnarray}&&{{G}^{-1}}=\frac{1}{M}\frac{\partial }{\partial {\bf c}}\otimes \frac{\partial }{\partial {\bf c}}+\mathop{\sum }\limits_{\alpha =1}^{N}\mathop{\sum }\limits_{\beta =1}^{N}{{H}^{\alpha \beta }}\frac{\partial }{\partial {{{\bf r}}_{\alpha }}}\otimes \frac{\partial }{\partial {{{\bf r}}_{\beta }}},\end{eqnarray} \tag{ 31 }$$

where $\{{{H}^{\alpha \beta }}\}$ is the inverse matrix to $\{{{H}_{\alpha \beta }}\}$, which turns out to be surprisingly simple:

$$\begin{eqnarray}&&{{H}^{\alpha \beta }}=\left( m_{\alpha }^{-1}\;{{\delta }_{\alpha \beta }}-{{M}^{-1}} \right).\end{eqnarray} \tag{ 32 }$$

In fact, the relation $\sum _{\beta =1}^{N}{{H}_{\alpha \beta }}{{H}^{\beta \gamma }}=\delta _{\alpha }^{\gamma }$ is easily checked from the given expressions.

Note that the kinetic part in (5) is just $(-{{\hbar }^{2}}/2)$ times the Laplacian on ${{\mathbb{R}}^{3\left( N+1 \right)}}$ with respect to the kinetic-energy metric. Since ${\rm det} \left( G \right)$ and ${\rm det} \left( H \right)$ are constant, this Laplacian is just:

$$\begin{eqnarray}\begin{array}{rcl} {{\Delta }_{G}} & = & \mathop{\sum }\limits_{a=0}^{N}\mathop{\sum }\limits_{b=0}^{N}{{G}^{ab}}\frac{\partial }{\partial {{{\bf x}}_{a}}}\cdot \frac{\partial }{\partial {{{\bf x}}_{b}}}=\mathop{\sum }\limits_{a=0}^{N}\frac{1}{{{m}_{a}}}\frac{\partial }{\partial {{{\bf x}}_{a}}}\cdot \frac{\partial }{\partial {{{\bf x}}_{a}}} \\ {} & = & \frac{1}{M}\frac{\partial }{\partial {\bf c}}\cdot \frac{\partial }{\partial {\bf c}}+\mathop{\sum }\limits_{\alpha =1}^{N}\mathop{\sum }\limits_{\beta =1}^{N}{{H}^{\alpha \beta }}\frac{\partial }{\partial {{{\bf r}}_{\alpha }}}\cdot \frac{\partial }{\partial {{{\bf r}}_{\beta }}} \\ {} & =: & {{\Delta }_{c}}+{{\Delta }_{r}}. \\ \end{array}\end{eqnarray} \tag{ 33 }$$

Here ${{\Delta }_{c}}$ is the part just involving the three centre-of-mass coordinates ${\bf c}$ and ${{\Delta }_{r}}$ the part involving the derivatives with respect to the $3N$ relative coordinates ${{{\bf r}}_{\alpha }}$. Note that there are no terms that mix the derivatives with repect to ${\bf c}$ and ${{{\bf r}}_{\alpha }}$, but that ${{\Delta }_{r}}$ mixes any two derivatives with respect to ${{{\bf r}}_{\alpha }}$ due to the second term on the right-hand side of (32). Clearly, a further linear redefinition of the relative coordinates ${{{\bf r}}_{\alpha }}$ could be employed to diagonalize ${{H}_{\alpha \beta }}$ and ${{H}^{\alpha \beta }}$, but that we will not need here.

Having introduced the centre-of-mass coordinates, one can consider the possibility that the wavefunction separates into a centre-of-mass and a relative part ³ ,

$$\begin{eqnarray}&&\Psi \left( t;{{{\bf x}}_{0}},\cdots ,{{{\bf x}}_{N}} \right)={{\left( \frac{{{m}_{0}}}{M} \right)}^{3/2}}\psi \left( t;{\bf c} \right)\;\chi \left( t;{{{\bf r}}_{1}},\cdots ,{{{\bf r}}_{N}} \right).\end{eqnarray} \tag{ 34 }$$

In order to obtain an independent equation for just the centre-of-mass dynamics one is, however, left with the necessity to show that equation (5) also separates for this ansatz. This is true for the kinetic term, as shown in (33), and it is also obvious for the non-gravitational contribution ${{V}_{ab}}$, which depends on the relative distances, and therefore the relative coordinates, only.

As long as non-gravitational interactions are present these are presumably much stronger than any gravitational effects. Hence, the latter can be ignored for the relative motion, which leads to a usually complicated but well-known equation: the ordinary, linear Schrödinger equation whose solution becomes manifest in the inner structure of the present lump of matter.

However, while separating the linear multi-particle Schrödinger equation in the absence of external forces (i.e. equation (5) with the gravitational constant G set to zero) yields a free Schrödinger equation for the evolution of the centre of mass, the $\left( N+1 \right)$-particle Schrödinger–Newton equation (5) will comprise contributions of the gravitational potential to the centre-of-mass motion. The reason for these to appear is the non-locality of the integral term in the equation (and not the mere existence of the diagonal term a = b as one could naively assume).

Let us take a closer look at the gravitational potential (11). Using the results from the previous section, in centre-of-mass coordinates it reads:

$$\begin{eqnarray}\begin{array}{rcl} {{V}_{{\rm g}}}\left( t;{\bf c},{{{\bf r}}_{1}},\cdots ,{{{\bf r}}_{N}} \right) & = & -G\;\int {{{\rm d}}^{3}}{\bf c}^{\prime} {{\left| \psi \left( t;{\bf c}^{\prime} \right) \right|}^{2}}\left\{ \int \mathop{\prod }\limits_{\gamma =1}^{N}{{{\rm d}}^{3}}{\bf r}_{\gamma }^{^{\prime} } \right\} \\ {} & {} & \times [m_{0}^{2}\frac{{{\left| \chi \left( t;{\bf r}_{1}^{^{\prime} },\cdots ,{\bf r}_{N}^{^{\prime} } \right) \right|}^{2}}}{\parallel {\bf c}-{\bf c}^{\prime} -\mathop{\sum }\limits_{\delta =1}^{N}\frac{{{m}_{\delta }}}{{{m}_{0}}}\left( {{{\bf r}}_{\delta }}-{\bf r}_{\delta }^{^{\prime} } \right)\parallel } \\ {} & {} & +{{m}_{0}}\mathop{\sum }\limits_{\alpha =1}^{N}{{m}_{\alpha }}\frac{{{\left| \chi \left( t;{\bf r}_{1}^{^{\prime} },\cdots ,{\bf r}_{N}^{^{\prime} } \right) \right|}^{2}}}{\parallel {\bf c}-{\bf c}^{\prime} -\mathop{\sum }\limits_{\delta =1}^{N}\frac{{{m}_{\delta }}}{{{m}_{0}}}{{{\bf r}}_{\delta }}-{\bf r}_{\alpha }^{^{\prime} }\parallel } \\ {} & {} & +{{m}_{0}}\mathop{\sum }\limits_{\alpha =1}^{N}{{m}_{\alpha }}\frac{{{\left| \chi \left( t;{\bf r}_{1}^{^{\prime} },\cdots ,{\bf r}_{N}^{^{\prime} } \right) \right|}^{2}}}{\parallel {\bf c}-{\bf c}^{\prime} +{{{\bf r}}_{\alpha }}+\mathop{\sum }\limits_{\delta =1}^{N}\frac{{{m}_{\delta }}}{{{m}_{0}}}{\bf r}_{\delta }^{^{\prime} }\parallel } \\ {} & {} & +\mathop{\sum }\limits_{\alpha =1}^{N}\mathop{\sum }\limits_{\beta =1}^{N}{{m}_{\alpha }}{{m}_{\beta }}\frac{{{\left| \chi \left( t;{\bf r}_{1}^{^{\prime} },\cdots ,{\bf r}_{N}^{^{\prime} } \right) \right|}^{2}}}{\parallel {\bf c}-{\bf c}^{\prime} +{{{\bf r}}_{\alpha }}-{\bf r}_{\beta }^{^{\prime} }\mathop{\sum }\limits_{\delta =1}^{N}\parallel }].\qquad \qquad \qquad \qquad \; \\ \end{array}\end{eqnarray} \tag{ 35 }$$

The ${{m}_{0}}$ dependent terms in the second, third, and fourth lines are more intricate than those in the last line; but they are only $\left( 2N+1 \right)$ out of ${{\left( N+1 \right)}^{2}}$ terms and therefore can be neglected for large N ⁴ . In this 'large N'-approximation only the last double-sum in (35) survives. All ${\bf r}_{\gamma }^{^{\prime} }$ integrations except that where $\gamma =\beta$ can be carried out (obtaining the $\beta$th marginal distributions for $|\chi (t;{\bf r}_{1}^{^{\prime} },\cdots ,{\bf r}_{N}^{^{\prime} }){{|}^{2}}$). Because of the remaining integration over ${\bf r}_{\beta }^{^{\prime} }$ we may rename the integration variable ${\bf r}_{\beta }^{^{\prime} }\to {\bf r}^{\prime}$, thereby removing its fictitious dependence on $\beta$. All this leads to the expression

$$\begin{eqnarray}&&{{V}_{{\rm g}}}\left( t;{\bf c},{{{\bf r}}_{1}},\cdots ,{{{\bf r}}_{N}} \right)=-G\mathop{\sum }\limits_{\alpha =1}^{N}{{m}_{\alpha }}\int {{{\rm d}}^{3}}{\bf c}^{\prime} \int {{{\rm d}}^{3}}{\bf r}^{\prime} \frac{{{\left| \psi \left( t;{\bf c}^{\prime} \right) \right|}^{2}}{{\rho }_{c}}\left( {\bf r}^{\prime} \right)}{\parallel {\bf c}-{\bf c}^{\prime} +{{{\bf r}}_{\alpha }}-{\bf r}^{\prime} \parallel }\;,\end{eqnarray} \tag{ 36 }$$

where we defined

$$\begin{eqnarray}&&{{\rho }_{c}}\left( t;{\bf r} \right):=\mathop{\sum }\limits_{\beta =1}^{N}{{m}_{\beta }}\{\int \prod \gamma \ne \beta N\gamma =1{{{\rm d}}^{3}}{{{\bf r}}_{\gamma }}\mathop{\sum }\limits_{\beta =1}^{N}\}\ {{\left| \chi \left( t;{{{\bf r}}_{1}},\cdots ,{{{\bf r}}_{\beta -1}},{\bf r},{{{\bf r}}_{\beta +1}},\cdots ,{{{\bf r}}_{N}} \right) \right|}^{2}}\;.\end{eqnarray} \tag{ 37 }$$

This 'relative' mass distribution is built analogously to (7) from the marginal distributions, here involving only the relative coordinates of all but the zeroth particle. In the large N approximation this omission of ${{m}_{0}}$ should be neglected and ${{\rho }_{c}}\left( t;{\bf r} \right)$ should be identified as the mass distribution relative to the centre of mass. Given a (stationary) solution $\chi$ of the Schrödinger equation for the relative motion, ${{\rho }_{c}}$ is then simply the mass density of the present lump of matter (e.g. a molecule) relative to the centre of mass. Although for the following discussion the time dependence of ${{\rho }_{c}}$ makes no difference, we will omit it. This may be justified by an adiabatic approximation, since the typical frequencies involved in the relative motions are much higher than the frequencies involved in the centre-of-mass motion.

Note that the only approximation that entered the derivation of (36) so far is that of large N. For the typical situations we want to consider, where N is large indeed, this will be harmless. However, the analytic form taken by the gravitational potential in (36) is not yet sufficiently simple to allow for a separation into centre-of-mass and relative motion. In order to perform such a separation we have to get rid of the ${{{\bf r}}_{\alpha }}$-dependence. This can be achieved if further approximations are made, as we shall explain now.

5.1. Wide wavefunctions

As long as the centre-of-mass wavefunction is much wider than the extent of the considered object one can assume that it does not change much over the distance ${{{\bf r}}_{\alpha }}$, i.e. $|\psi (t;{\bf c}^{\prime} +{{{\bf r}}_{\alpha }})|\approx \left| \psi \left( t;{\bf c}^{\prime} \right) \right|$. Substituting ${\bf c}^{\prime}$ by ${\bf c}^{\prime} +{{{\bf r}}_{\alpha }}$ in (36) then yields the following potential, depending only on the centre-of-mass coordinate:

$$\begin{eqnarray}&&V_{{\rm g}}^{\left( A \right)}\left( t;{\bf c} \right)\approx -GM\int {{{\rm d}}^{3}}{\bf c}^{\prime} \int {{{\rm d}}^{3}}{\bf r}^{\prime} \ \frac{{{\left| \psi \left( t;{\bf c}^{\prime} \right) \right|}^{2}}{{\rho }_{c}}\left( {\bf r}^{\prime} \right)}{\parallel {\bf c}-{\bf c}^{\prime} -{\bf r}^{\prime} \parallel }\;.\end{eqnarray} \tag{ 38 }$$

As a result, the equation for the centre of mass is now indeed of type (1) with ${{V}_{{\rm g}}}=V_{{\rm g}}^{\left( A \right)}$ being given by M times the convolution of ${{\left| \psi \right|}^{2}}$ with the Newtonian gravitational potential for the mass density ${{\rho }_{c}}$. Case A has been further analysed in [22].

5.2. Born–Oppenheimer-type approximation

An alternative way to get rid of the dependence of (36) on the relative coordinates, i.e., the ${{{\bf r}}_{\alpha }}$-dependence on the right-hand side, is to just replace ${{V}_{{\rm g}}}$ with its expectation value in the state $\chi$ of the relative motion ⁵ . This procedure corresponds to the Born–Oppenheimer approximation in molecular physics where the electronic degrees of freedom are averaged over in order to solve the dynamics of the nuclei. The justification for this procedure in molecular physics derives from the much smaller timescales for the motion of the fast and lighter electrons as compared to the slow and heavier nuclei. Hence the latter essentially move only according to the averaged potential sourced by the electrons. In case of the Schrödinger–Newton equation, the justification is formally similar, even though it is clear that there is no real material object attached to the centre of mass. What matters is that the relative interactions (based on electrodynamic forces) are much stronger than the gravitational ones, so that the characteristic frequencies of the former greatly exceed those of the latter; compare, e.g., the discussion in [23].

Now, the expectation value is easily calculated:

$$\begin{eqnarray}\begin{array}{rcl} V_{{\rm g}}^{(B)}(t;{\bf c}) & = & \left\{ \int \mathop{\prod }\limits_{\beta =1}^{N}{{{\rm d}}^{3}}{\bf r}_{\beta }^{^{\prime\prime} } \right\}{{\left| \chi ({\bf r}_{1}^{^{\prime\prime} },\cdots ,{\bf r}_{N}^{^{\prime\prime} }) \right|}^{2}}\ {{V}_{{\rm g}}}(t;{\bf c},{\bf r}_{1}^{^{\prime\prime} },\cdots ,{\bf r}_{N}^{^{\prime\prime} }) \\ {} & = & -G\mathop{\sum }\limits_{\alpha =1}^{N}{{m}_{\alpha }}\int {{{\rm d}}^{3}}{\bf c}^{\prime} \int {{{\rm d}}^{3}}{\bf r}^{\prime} \left\{ \int \mathop{\prod }\limits_{\beta =1}^{N}{{{\rm d}}^{3}}{\bf r}_{\beta }^{^{\prime\prime} } \right\}\frac{{{\left| \psi (t;{\bf c}^{\prime} ) \right|}^{2}}{{\rho }_{c}}({\bf r}^{\prime} ){{\left| \chi ({\bf r}_{1}^{^{\prime\prime} },\cdots ,{\bf r}_{N}^{^{\prime\prime} }) \right|}^{2}}}{\parallel {\bf c}-{\bf c}^{\prime} -{\bf r}^{\prime} +{\bf r}_{\alpha }^{^{\prime\prime} }\parallel } \\ {} & = & -G\int {{{\rm d}}^{3}}{\bf c}^{\prime} \int {{{\rm d}}^{3}}{\bf r}^{\prime} \int {{{\rm d}}^{3}}{\bf r}^{\prime\prime} \ \frac{{{\left| \psi (t;{\bf c}^{\prime} ) \right|}^{2}}{{\rho }_{c}}({\bf r}^{\prime} ){{\rho }_{c}}({\bf r}^{\prime\prime} )}{\parallel {\bf c}-{\bf c}^{\prime} -{\bf r}^{\prime} +{\bf r}^{\prime\prime} \parallel }. \\ \end{array}\end{eqnarray} \tag{ 39 }$$

Note that this expression involves one more ${{\mathbb{R}}^{3}}$ integration than (38).

In [22] we studied two simple models for the matter density ${{\rho }_{c}}$: a solid and a hollow sphere. The solid sphere suffers from some peculiar divergence issues, which we explain in appendix B, and is also mathematically slightly more difficult to handle than the hollow sphere whose radial mass distribution is just a $\delta$-function. We therefore use the hollow sphere as a model to compare the two approximation ansätze given above.

While in [2, 12, 13, 22] the expression 'collapse mass' was used in a rather loosely defined manner, here we define as critical mass the mass value for which at t = 0 the second-order time derivative of the second moment $Q\left( t \right)=\int {{{\rm d}}^{3}}{\bf c}\;{{\left| {\bf c} \right|}^{2}}{{\left| \psi \left( t;{\bf c} \right) \right|}^{2}}$ vanishes, i.e. $\ddot{Q}\left( t=0 \right)=0$. (Note that for a real-valued initial wave-packet the first-order time derivative always vanishes.) For the one-particle Schrödinger–Newton equation and a Gaussian wave packet of $0.5\;\mu {\rm m}$ width this yields a critical mass of $6.5\times {{10}^{9}}$ u, which fits very well with the numerical results obtained in [13].

For the hollow sphere we then obtain the analytic expression

$$\begin{eqnarray}&&{{m}_{{\rm crit}}}={{\left( \sqrt{\frac{\pi }{2}}\;\frac{3{{\hbar }^{2}}}{G\;a\;f\left( R/a \right)} \right)}^{1/3}}\approx 5.153\times {{10}^{9}}\;{\rm u}{{\left( \left( a/\mu {\rm m} \right)\times f\left( \frac{R}{a} \right) \right)}^{-1/3}}\end{eqnarray} \tag{ 40 }$$

for the critical mass. This expression is derived in appendix A. The function f is constantly 1 in the case of the one-particle Schrödinger–Newton equation and shows exponential dependence on R in the case of the wide wavefunction approximation. In the case of the Born–Oppenheimer approximation, f is a rather complicated function that can be found in the appendix.

The resulting critical mass for a width of the centre-of-mass wavefunction of $0.5\;\mu {\rm m}$ is plotted as a function of the hollow-sphere radius in figure 1. The curve that the figure shows for the wide wavefunction approximation coincides well with the results we obtained in the purely numerical analysis in [22]. For the Born–Oppenheimer-type approximation the plot shows a radius dependence of the collapse mass that is almost linear. This is in agreement with the result by Diósi [15], who estimates the width of the ground state for a solid sphere to be proportional to ${{\left( R/M \right)}^{3/4}}$.

Zoom In Zoom Out Reset image size

5.3. Narrow wavefunctions in the Born–Oppenheimer scheme

With the Born–Oppenheimer-type approximation scheme just derived we now possess a tool with which we can consider the opposite geometric situation to that in Case A, namely for widths of the centre-of-mass wavefunction $\psi$ that are much smaller than the extensions (diameters of the support) of the matter distribution ${{\rho }_{c}}$, i.e., for well localized mass centres inside the bulk of matter.

Let us recall that in Newtonian gravitational physics the overall gravitational self-energy of a mass distribution $\tilde{\rho }$ is given by

$$\begin{eqnarray}&&{{U}_{{\rm g}}}\left( \tilde{\rho } \right):=-\frac{G}{2}\int \;{{{\rm d}}^{3}}{\bf x}\int \;{{{\rm d}}^{3}}{\bf x}^{\prime} \ \frac{\tilde{\rho }\left( {\bf x} \right)\tilde{\rho }\left( {\bf x}^{\prime} \right)}{\parallel {\bf x}-{\bf x}^{\prime} \parallel }\;.\end{eqnarray} \tag{ 41 }$$

If $\tilde{\rho }=\rho +\rho ^{\prime}$, we have by the simple quadratic dependence on $\tilde{\rho }$

$$\begin{eqnarray}&&{{U}_{{\rm g}}}\left( \rho +\rho ^{\prime} \right):={{U}_{{\rm g}}}\left( \rho \right)+{{U}_{{\rm g}}}\left( \rho ^{\prime} \right)+{{I}_{{\rm g}}}\left( \rho ,\rho ^{\prime} \right),\end{eqnarray} \tag{ 42 }$$

where

$$\begin{eqnarray}&&{{I}_{{\rm g}}}\left( \rho ,\rho ^{\prime} \right):=-G\int \;{{{\rm d}}^{3}}{\bf x}\int \;{{{\rm d}}^{3}}{\bf x}^{\prime} \ \frac{\rho \left( {\bf x} \right)\rho ^{\prime} \left( {\bf x}^{\prime} \right)}{\parallel {\bf x}-{\bf x}^{\prime} \parallel }\;.\end{eqnarray} \tag{ 43 }$$

represents the mutual gravitational interaction of the matter represented by $\rho$ with that represented by $\rho ^{\prime}$. In the special case $\rho ^{\prime} ={{T}_{{\bf d}}}\rho$, where ${{T}_{{\bf d}}}$ denotes the operation of translation by the vector ${\bf d}$,

$$\begin{eqnarray}&&\left( {{T}_{{\bf d}}}\rho \right)\left( {\bf x} \right):=\rho \left( {\bf x}-{\bf d} \right),\end{eqnarray} \tag{ 44 }$$

we set

$$\begin{eqnarray}&&{{I}_{\rho }}\left( {\bf d} \right):={{I}_{{\rm g}}}\left( \rho ,{{T}_{{\bf d}}}\rho \right).\end{eqnarray} \tag{ 45 }$$

It is immediate from (43) that ${{I}_{\rho }}:{{\mathbb{R}}^{3}}\to \mathbb{R}$ has a zero derivative at the origin ${\bf 0}\in {{\mathbb{R}}^{3}}$,

$$\begin{eqnarray}&&I_{\rho }^{^{\prime} }\left( {\bf 0} \right)=0,\end{eqnarray} \tag{ 46 }$$

and that it satisfies the following equivariance

$$\begin{eqnarray}&&{{I}_{\rho }}\left( R{\bf d} \right)={{I}_{\rho {}^\circ R}}\left( {\bf d} \right)\end{eqnarray} \tag{ 47 }$$

for any orthogonal 3 × 3 matrix R. The latter implies the rather obvious result that the function ${\bf d}\mapsto {{I}_{\rho }}\left( {\bf d} \right)$ is rotationally invariant if $\rho$ is a rotationally invariant distribution, i.e., the interaction energy depends only on the modulus of the shift, not its direction.

For example, given that $\rho$ is the matter density of a homogeneous sphere of radius R and mass M,

$$\begin{eqnarray}\rho \left( {\bf x} \right)=\left\{ \begin{array}{ccccccccccccccc} \frac{3\;M}{4\pi {{R}^{3}}} & {\rm for} & \parallel {\bf x}\parallel \leqslant R \\ 0 & {\rm for} & \parallel {\bf x}\parallel >R, \\ \end{array} \right.\end{eqnarray} \tag{ 48 }$$

the gravitational interaction energy between two such identical distributions a distance $d:=\parallel {\bf d}\parallel$ apart is

$$\begin{eqnarray}{{I}_{\rho }}\left( d \right)=-\frac{G{{M}^{2}}}{R}\times \left\{ \begin{array}{ccccccccccccccc} \frac{6}{5}-2{{\left( \frac{d}{2R} \right)}^{2}}+\frac{3}{2}{{\left( \frac{d}{2R} \right)}^{3}}-\frac{1}{5}{{\left( \frac{d}{2R} \right)}^{5}} & {\rm for} & d\leqslant 2R, \\ \frac{R}{d} & {\rm for} & d\geqslant 2R. \\ \end{array} \right.\end{eqnarray} \tag{ 49 }$$

The second line is obvious, whereas the first line follows, e.g., from specializing the more general formula (42) of [16] to equal radii (${{R}_{p}}={{R}_{t}}$) and making the appropriate redefinitions in order to translate their electrostatic to our gravitational case. This formula also appears in [24].Using the definitions (43) and (45), we can rewrite the right-hand side of (39) as convolution of ${{\left| \psi \right|}^{2}}$ with ${{I}_{{{\rho }_{c}}}}$:

$$\begin{eqnarray}&&V_{{\rm g}}^{\left( B \right)}\left( t;{\bf c} \right)=\int \;{{{\rm d}}^{3}}{\bf c}^{\prime} \;{{I}_{{{\rho }_{c}}}}\left( {\bf c}-{\bf c}^{\prime} \right){{\left| \psi \left( t;{\bf c}^{\prime} \right) \right|}^{2}}\;.\end{eqnarray} \tag{ 50 }$$

Since in equation (3) this potential is multiplied with $\psi \left( t,{\bf c} \right)$, we see that only those values of ${{I}_{{{\rho }_{c}}}}\left( {\bf c}-{\bf c}^{\prime} \right)$ will contribute where ${{\left| \psi \left( t;{\bf c}^{\prime} \right) \right|}^{2}}\psi \left( t;{\bf c} \right)$ appreciably differs from zero. Hence if $\psi$ is concentrated in a region of diameter D then we need to know ${{I}_{{{\rho }_{c}}}}\left( {\bf c}-{\bf c}^{\prime} \right)$ only for $\parallel {\bf c}-{\bf c}^{\prime} \parallel <D$. Assuming D to be small we expand ${{I}_{{{\rho }_{c}}}}$ in a Taylor series. Because of (46) there is no linear term, so that up to and including the quadratic terms we have (using that ${{\left| \psi \left( t;{\bf c}^{\prime} \right) \right|}^{2}}$ is normalized with respect to the measure ${{{\rm d}}^{3}}{\bf c}^{\prime}$)

$$\begin{eqnarray}&&V_{{\rm g}}^{\left( B \right)}\left( t;{\bf c} \right)\approx {{I}_{{{\rho }_{c}}}}\left( {\bf 0} \right)+\frac{1}{2}I_{{{\rho }_{c}}}^{^{\prime\prime} }\left( {\bf 0} \right)\cdot \left( {\bf c}\otimes {\bf c}-2\;{\bf c}\otimes \left\langle {\bf c} \right\rangle +\left\langle {\bf c}\otimes {\bf c} \right\rangle \right).\end{eqnarray} \tag{ 51 }$$

Here $I_{{{\rho }_{c}}}^{^{\prime\prime} }\left( {\bf 0} \right)$ denotes the second derivative of the function ${{I}_{{{\rho }_{c}}}}:{{\mathbb{R}}^{3}}\to \mathbb{R}$ at ${\bf 0}\in \mathbb{R}$ (which is a symmetric bilinear form on ${{\mathbb{R}}^{3}}$) and $\left\langle \cdot \right\rangle$ denotes the expectation value with respect to $\psi$. We stress that the nonlinearity in $\psi$ is now entirely encoded into this state dependence of the expectation values which appear in the potential. If, for simplicity, we only consider centre-of-mass motions in one dimension, the latter being coordinatized by $c\in \mathbb{R}$, then (51) simplifies to

$$\begin{eqnarray}&&V_{{\rm g}}^{(B)}(t;c)\;\approx \;{{I}_{{{\rho }_{c}}}}(0)+\frac{1}{2}I_{{{\rho }_{c}}}^{^{\prime\prime} }(0)({{c}^{2}}-2c\;\langle c\rangle +\langle {{c}^{2}}\rangle )\end{eqnarray} \tag{ 52a }$$

$$\begin{eqnarray}&&\qquad \qquad \ =\;{{I}_{{{\rho }_{c}}}}(0)+\frac{1}{2}I_{{{\rho }_{c}}}^{^{\prime\prime} }(0){{(c-\langle c\rangle )}^{2}}+\frac{1}{2}I_{{{\rho }_{c}}}^{^{\prime\prime} }(0)(\langle {{c}^{2}}\rangle -{{\left\langle c \right\rangle }^{2}}).\qquad \qquad \qquad \qquad \quad \end{eqnarray} \tag{ 52b }$$

The first term, ${{I}_{{{\rho }_{c}}}}\left( 0 \right)$, just adds a constant to the potential, which can be absorbed by adding $-\left( {\rm i}/\hbar \right){{I}_{{{\rho }_{c}}}}\left( 0 \right)t$ to the phase $\psi$. The second term is the crucial one and has been shown in [23] to give rise to interesting and potentially observable effects for Gaussian states.

More precisely, consider a one-dimensional nonlinear Schrödinger evolution of the form (1a ) with ${{V}_{{\rm g}}}$ given by the second term in (52) and an additional external harmonic potential for the centre of mass, then we get the following nonlinear Schrödinger–Newton equation for the centre-of-mass wavefunction,

$$\begin{eqnarray}&&{\rm i}\hbar {{\partial }_{t}}\psi \left( t;c \right)=\left( -\frac{{{\hbar }^{2}}}{2\;M}\frac{{{\partial }^{2}}}{\partial {{c}^{2}}}+\frac{1}{2}M\omega _{c}^{2}{{c}^{2}}+\frac{1}{2}M\omega _{{\rm SN}}^{2}{{\left( c-\langle c\rangle \right)}^{2}} \right)\psi \left( t;c \right),\end{eqnarray} \tag{ 53 }$$

where ${{\omega }_{{\rm SN}}}:=\sqrt{I_{{{\rho }_{c}}}^{^{\prime\prime} }\left( 0 \right)/M}$ is called the Schrödinger–Newton frequency. This equation has been considered in [23], where the last term on the right-hand side of (52b) has been neglected for a priori no good reason. Note that $\langle {{c}^{2}}\rangle$ and ${{\langle c\rangle }^{2}}$ contain the wavefunction and hence are therefore not constant (in time). Now, in the context of [23] the consequences of interest were the evolution equations for the first and second moments in the canonical phase-space variables, and it shows that for them only spatial derivatives of the potential contribute. As a consequence, the term in question makes no difference. The relevant steps in the computation are displayed in appendix C.

Based on the observation that equation (53) evolves Gaussian states into Gaussian states, it has then been shown that the covariance ellipse of the Gaussian state rotates at frequency ${{\omega }_{q}}:=\sqrt{\omega _{c}^{2}+\omega _{{\rm SN}}^{2}}$, whereas the centre of the ellipse orbits the origin in phase with frequency ${{\omega }_{c}}$. This asynchrony results from a difference between first- and second-moment evolution and is a genuine effect of self-gravity. It has been suggested that it may be observable via the output spectra of optomechanical systems [23].

Although the many-particle Schrödinger–Newton equation (5) does not exactly separate into centre-of-mass and relative motion, we could show that for some well-motivated approximations such a separation is possible. As long as the extent of an object is negligible in comparison to the uncertainty in localization of its centre of mass the one-particle equation (3) is a good model in both approximation schemes considered.

In the opposite case of a well localized object, i.e. one that has a narrow wavefunction compared to its extent, the gravitational potential takes the form (52), which yields a closed system of equations for the first and second moments and therefore the effects described in [23]. The nonlinear Schrödinger equation resulting from the potential (52) is also considered in [24], where it is used for comparison of Schrödinger–Newton dynamics with models of quantum state reduction and decoherence.

The modification (38) provides a valid correction of the one-particle Schrödinger–Newton equation for objects of finite but small radii. This equation was considered in [25] and studied numerically in [22]. It remains unclear for which ratio of the objectʼs extent to the width of the wavefunction the Born–Oppenheimer-type approximation (39) starts to be superior to the wide wavefunction approximation. It may even be the better approximation throughout the whole range of possible object sizes and wavefunctions, since a Born–Oppenheimer-like approximation is implicitly assumed also for the wide wavefunction when the mass density is taken to be that of a solid object.

In passing we make the final technical remark that the analysis of the critical mass for the hollow sphere shows that this mass increases linearly with the radius R of the sphere. Given a fixed mass, this implies that the width of the stationary solution increases like ${{R}^{3/4}}$, a relation already found by Diósi [15].

The interface between quantum mechanics and gravity theory remains one of the most interesting and profound challenges with hopefully revealing experimental consequences, which we are only beginning to explore. In this context one should also mention that nonlinear one-particle Schrödinger equations are of course also considered for Einstein–Bose condensates, in which case inclusion of self-gravity adds a Schrödinger–Newton term in addition to that nonlinear term obtained from the effective potential within the Hartree–Fock approximation (Gross–Pitaevskii–Newton equation). Such equations are derivable for particle numbers $N\to \infty$ without further hypotheses and may open up the possibility to test self-gravity effects on large quantum systems. Recent experiments have demonstrated the high potential of atom interferometry on freely falling Einstein–Bose condensates [26] and it seems an interesting question whether this may be used to see self-gravity effects on such systems.

We gratefully acknowledge funding and support through the Center for Quantum Engineering and Space-Time Research (QUEST) at the Leibniz University Hannover and the Center of Applied Space Technology and Microgravity (ZARM) at the University of Bremen. AG is supported by the John Templeton foundation (grant 39530).

For both the wide wavefunction approximation (38) and the Born–Oppenheimer-type approximation (39) one must solve integrals of the type

$$\begin{eqnarray}&&I\left( {\bf a} \right)=\int {{{\rm d}}^{3}}{\bf r}\ \frac{{{\rho }_{c}}\left( {\bf r} \right)}{\parallel {\bf r}-{\bf a}\parallel }\;.\end{eqnarray} \tag{ A.1 }$$

In a spherically symmetric situation these take the form

$$\begin{eqnarray}\begin{array}{rcl} I(a) & = & \int _{0}^{\infty }{{r}^{2}}\;{\rm d}r\int _{-1}^{1}{\rm d}{\rm cos} \theta \int _{0}^{2\pi }{\rm d}\varphi \;\frac{{{\rho }_{c}}(r)}{\sqrt{{{r}^{2}}+{{a}^{2}}-2ra{\rm cos} \theta }} \\ {} & = & \frac{4\pi }{a}\int _{0}^{a}{\rm d}r\;{{r}^{2}}\;{{\rho }_{c}}(r)+4\pi \int _{a}^{\infty }{\rm d}r\;r\;{{\rho }_{c}}(r), \\ \end{array}\end{eqnarray} \tag{ A.2 }$$

where we write a for the absolute value $\left| {\bf a} \right|$, etc. If now we assume that ${{\rho }_{c}}$ is the mass density of a hollow sphere of radius R, i.e.

$$\begin{eqnarray}&&{{\rho }_{c}}\left( r \right)=\frac{M}{4\pi \;{{r}^{2}}}\;\delta \left( r-R \right)\;,\end{eqnarray} \tag{ A.3 }$$

these integrals simplify to

$$\begin{eqnarray}{{I}_{R}}\left( a \right)=\left\{ \begin{array}{ccccccccccccccc} \frac{M}{R} & {\rm if}\;a<R \\ \frac{M}{a} & {\rm if}\;a\geqslant R. \\ \end{array} \right.\;\end{eqnarray} \tag{ A.4 }$$

With this the wide wavefunction approximation (38) results in

$$\begin{eqnarray}&&V_{{\rm g}}^{\left( A \right)}\left( t;c;R \right)=-GM\int {{{\rm d}}^{3}}{\bf c}^{\prime} \;{{\left| \psi \left( t;{\bf c}^{\prime} \right) \right|}^{2}}\;{{I}_{R}}\left( \parallel {\bf c}-{\bf c}^{\prime} \parallel \right).\end{eqnarray} \tag{ A.5 }$$

On the other hand, the Born–Oppenheimer approximation (39) leads to

$$\begin{eqnarray}&&V_{{\rm g}}^{\left( B \right)}\left( t;c;R \right)=-G\int {{{\rm d}}^{3}}{\bf c}^{\prime} \;{{\left| \psi \left( t;{\bf c}^{\prime} \right) \right|}^{2}}\int {{{\rm d}}^{3}}{\bf r}^{\prime} \;{{\rho }_{c}}\left( {\bf r}^{\prime} \right)\;{{I}_{R}}\left( \parallel {\bf c}-{\bf c}^{\prime} -{\bf r}^{\prime} \parallel \right)\;.\end{eqnarray} \tag{ A.6 }$$

In order to be able to obtain an analytical result we consider the initial Gaussian wave packet

$$\begin{eqnarray}&&\psi \left( t=0;{\bf c} \right)={{\left( \pi {{a}^{2}} \right)}^{-3/4}}\;{\rm exp} \left( -\frac{{{c}^{2}}}{2{{a}^{2}}} \right)\;,\end{eqnarray} \tag{ A.7 }$$

for which these potentials take the form

$$\begin{eqnarray}\begin{array}{rcl} V_{A}^{0}(c;R) & = & -\frac{G{{M}^{2}}}{2}\left\{ \frac{a}{\sqrt{\pi }\;c\;R}\left[ {\rm exp} \left( -\frac{{{\left( c+R \right)}^{2}}}{{{a}^{2}}} \right)-{\rm exp} \left( -\frac{{{\left( c-R \right)}^{2}}}{{{a}^{2}}} \right) \right] \right. \\ {} & {} & +\frac{1}{c}\left[ {\rm erf}\left( \frac{c+R}{a} \right)+{\rm erf}\left( \frac{c-R}{a} \right) \right] \\ {} & {} & +\left. \frac{1}{R}\left[ {\rm erf}\left( \frac{c+R}{a} \right)-{\rm erf}\left( \frac{c-R}{a} \right) \right] \right\} \\ \end{array}\end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}\begin{array}{rcl} V_{B}^{0}(c;R) & = & -\frac{G{{M}^{2}}}{2}\left\{ \frac{a}{\sqrt{\pi }\;c}\left( \frac{1}{R}+\frac{8R}{3{{a}^{2}}} \right){\rm exp} \left( -\frac{{{\left( c+2R \right)}^{2}}}{{{a}^{2}}} \right) \right. \\ {} & {} & +\frac{a}{\sqrt{\pi }\;R}\left( \frac{1}{2c}-\frac{1}{4R} \right)\left[ {\rm exp} \left( -\frac{{{\left( c+2R \right)}^{2}}}{{{a}^{2}}} \right)-{\rm exp} \left( -\frac{{{\left( c-2R \right)}^{2}}}{{{a}^{2}}} \right) \right] \\ {} & {} & +\frac{1}{c}\left[ {\rm erf}\left( \frac{c+2R}{a} \right)+{\rm erf}\left( \frac{c-2R}{a} \right) \right] \\ {} & {} & +\left. \frac{1}{R}\left( 1-\frac{c}{4R}-\frac{{{a}^{2}}}{8cR} \right)\left[ {\rm erf}\left( \frac{c+2R}{a} \right)-{\rm erf}\left( \frac{c-2R}{a} \right) \right] \right\}\;. \\ \end{array}\end{eqnarray} \tag{ A.9 }$$

Note that both potentials agree in the limits

$$\begin{eqnarray}&&\mathop{{\rm lim} }\limits_{R\to 0}V_{A,B}^{0}\left( c;R \right)=-\frac{G{{M}^{2}}}{c}\;{\rm erf}\left( \frac{c}{a} \right)\quad {\rm and}\quad \mathop{{\rm lim} }\limits_{R\to \infty }V_{A,B}^{0}\left( c;R \right)=0.\end{eqnarray} \tag{ A.10 }$$

As a measure to compare these potentials with each other and the one-particle Schrödinger–Newton equation we use the second moment $Q\left( t \right)=\int {{{\rm d}}^{3}}{\bf c}\;{{\left| {\bf c} \right|}^{2}}{{\left| \psi \left( t;{\bf c} \right) \right|}^{2}}$. For a real wave packet its first-order time derivative can be shown to vanish. Therefore the sign of the second-order time derivative $\ddot{Q}$ at t = 0 determines if a wave packet initially shrinks or increases in width. In general, the second-order time derivative is

$$\begin{eqnarray}&&\ddot{Q}\left( t \right)=\int {{{\rm d}}^{3}}{\bf c}\left( \frac{2{{\hbar }^{2}}}{{{M}^{2}}}{{\left| \nabla \psi \left( t;{\bf c} \right) \right|}^{2}}+\frac{2}{M}{{V}_{{\rm g}}}\left( t;{\bf c} \right)\left( 3\;{{\left| \psi \left( t;{\bf c} \right) \right|}^{2}}+{\bf c}\cdot \nabla {{\left| \psi \left( t;{\bf c} \right) \right|}^{2}} \right) \right),\end{eqnarray} \tag{ A.11 }$$

which for the spherically symmetric Gaussian state (A.7) takes the form

$$\begin{eqnarray}&&\ddot{Q}(t=0)=\frac{3{{\hbar }^{2}}}{{{M}^{2}}{{a}^{2}}}+\frac{8}{\sqrt{\pi }M{{a}^{5}}}\int _{0}^{\infty }{\rm d}c{\rm exp} \left( -\frac{{{c}^{2}}}{{{a}^{2}}} \right){{V}^{0}}(c)\left( 3{{a}^{2}}{{c}^{2}}-2{{c}^{4}} \right)\end{eqnarray} \tag{ A.12 }$$

$$\begin{eqnarray}&&\qquad \quad \ \ \;=\;\frac{3{{\hbar }^{2}}}{{{M}^{2}}{{a}^{2}}}-\sqrt{\frac{2}{\pi }}\frac{G\;M}{a}\;f\left( \frac{R}{a} \right)\;.\end{eqnarray} \tag{ A.13 }$$

The critical mass defined by $\ddot{Q}\left( t=0 \right)=0$ is then given by

$$\begin{eqnarray}&&{{m}_{{\rm crit}}}={{\left( \sqrt{\frac{\pi }{2}}\;\frac{3{{\hbar }^{2}}}{G\;a\;f\left( R/a \right)} \right)}^{1/3}}\approx 5.153\times {{10}^{9}}\;{\rm u}{{\left( \left( a/\mu {\rm m} \right)\times f\left( \frac{R}{a} \right) \right)}^{-1/3}}\;.\end{eqnarray} \tag{ A.14 }$$

The function f is $f\equiv 1$ for the one-particle Schrödinger–Newton equation. For the hollow sphere potential in the wide wavefunction and Born–Oppenheimer-type approximations, (A.8) and (A.9), respectively, this function can be calculated as

$$\begin{eqnarray}&&{{f}_{A}}\left( \frac{R}{a} \right)={\rm exp} \left( -\frac{{{R}^{2}}}{2{{a}^{2}}} \right)\end{eqnarray} \tag{ A.15 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{f}_{B}}\left( \frac{R}{a} \right) & = & \frac{2}{3}\sqrt{\frac{2}{\pi }}{\rm exp} \left( -\frac{4{{R}^{2}}}{{{a}^{2}}} \right)\frac{R}{a}\left( 1-{{\left( \frac{2R}{a} \right)}^{2}} \right) \\ {} & {} & +{\rm exp} \left( -\frac{2{{R}^{2}}}{{{a}^{2}}} \right)\left( 1-{\rm erf}\left( \sqrt{2}\frac{R}{a} \right) \right)\left( 1+\frac{1}{3}{{\left( \frac{2R}{a} \right)}^{4}} \right) \\ {} & {} & +\frac{{{a}^{2}}}{2{{R}^{2}}}\left( \frac{1}{\sqrt{2}}{\rm erf}\left( 2\frac{R}{a} \right)-{\rm exp} \left( -\frac{2{{R}^{2}}}{{{a}^{2}}} \right){\rm erf}\left( \sqrt{2}\frac{R}{a} \right) \right)\;. \\ \end{array}\end{eqnarray} \tag{ A.16 }$$

Given a spherically symmetric situation the wide wavefunction approximation (38) takes the form

$$\begin{eqnarray}&&V\left( c \right)=\left( {{\left| \psi \right|}^{2}}\;*\;\Phi \right)\left( c \right)=4\pi \ \int _{0}^{\infty }{\rm d}c^{\prime} \ c{{^{\prime} }^{2}}\;{{\left| \psi \left( c^{\prime} \right) \right|}^{2}}\;\Phi \left( \left| c-c^{\prime} {\mkern 1mu} \right|\; \right)\;,\end{eqnarray} \tag{ B.1 }$$

where for the potential $\Phi$ we want to consider the following three cases:

• Coulomb potential (i.e. the case of the Schrödinger–Newton equation (3)):
  $$\begin{eqnarray}&&{{\Phi }_{c}}\left( c \right)=-\frac{1}{c}\;,\end{eqnarray} \tag{ B.2 }$$
• hollow sphere of radius R:
  $$\begin{eqnarray}{{\phi }_{h}}\left( c \right)=\left\{ \begin{array}{ccccccccccccccc} -\frac{1}{R} & {\rm if}\;c<R \\ {{\Phi }_{c}}\left( c \right) & {\rm if}\;c\geqslant R, \\ \end{array} \right.\end{eqnarray} \tag{ B.3 }$$
• solid sphere of radius R:
  $$\begin{eqnarray}{{\Phi }_{s}}\left( c \right)=\left\{ \begin{array}{ccccccccccccccc} -\frac{3}{2R}+\frac{{{c}^{2}}}{2{{R}^{3}}} & {\rm if}\;c<R \\ {{\Phi }_{c}}\left( c \right) & {\rm if}\;c\geqslant R. \\ \end{array} \right.\;\end{eqnarray} \tag{ B.4 }$$

First we want to study the behaviour of ${{V}_{0}}=V\left( c=0 \right)$ for a Gaussian wave packet ${{\left| \psi \left( c \right) \right|}^{2}}=\frac{1}{4\pi }{\rm exp} (-{{c}^{2}})$. For convenience we omit all pre-factors. Equation (B.1) then reads:

$$\begin{eqnarray}&&{{V}_{0}}=\int _{0}^{\infty }{\rm d}c\ {{c}^{2}}\;{\rm exp} (-{{c}^{2}})\;\Phi (c)\end{eqnarray} \tag{ B.5 }$$

$$\begin{eqnarray}&&\quad =\int _{0}^{R}{\rm d}c\ {{c}^{2}}\;{\rm exp} (-{{c}^{2}})\;\Phi (c)+\int _{R}^{\infty }{\rm d}c\ {{c}^{2}}\;{\rm exp} (-{{c}^{2}})\;{{\Phi }_{c}}(c)\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&\quad =\int _{0}^{R}{\rm d}c\ {{c}^{2}}\;{\rm exp} (-{{c}^{2}})\;\Phi (c)-\frac{{\rm exp} (-{{R}^{2}})}{2}\;.\end{eqnarray} \tag{ B.7 }$$

For the three different potentials one obtains

$$\begin{eqnarray}\begin{array}{rcl} {{V}_{0,c}} & = & \int _{0}^{R}{\rm d}c\ {{c}^{2}}\;{\rm exp} (-{{c}^{2}})\;{{\Phi }_{c}}(c)-\frac{{\rm exp} (-{{R}^{2}})}{2} \\ {} & = & -\frac{1}{2}+\frac{{\rm exp} (-{{R}^{2}})}{2}-\frac{{\rm exp} (-{{R}^{2}})}{2} \\ {} & = & -\frac{1}{2}, \\ \end{array}\end{eqnarray} \tag{ B.8 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{V}_{0,h}} & = & \int _{0}^{R}{\rm d}c\ {{c}^{2}}\;{\rm exp} (-{{c}^{2}})\;{{\Phi }_{h}}(c)-\frac{{\rm exp} (-{{R}^{2}})}{2} \\ {} & = & \frac{{\rm exp} (-{{R}^{2}})}{2}-\frac{\sqrt{\pi }}{4R}{\rm erf}(R)-\frac{{\rm exp} (-{{R}^{2}})}{2} \\ {} & = & -\frac{\sqrt{\pi }}{4R}{\rm erf}(R), \\ \end{array}\end{eqnarray} \tag{ B.9 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{V}_{0,s}} & = & \frac{3}{2}\int _{0}^{R}{\rm d}c\ {{c}^{2}}\;{\rm exp} (-{{c}^{2}})\;{{\Phi }_{h}}(c)+\frac{1}{2{{R}^{3}}}\int _{0}^{R}{\rm d}c\ {{c}^{4}}{\rm exp} (-{{c}^{2}})-\frac{{\rm exp} (-{{R}^{2}})}{2} \\ {} & = & -\frac{3\sqrt{\pi }}{8R}{\rm erf}(R)-\frac{3}{8{{R}^{2}}}{\rm exp} (-{{R}^{2}})-\frac{3\sqrt{\pi }}{16{{R}^{2}}}{\rm erf}(R). \\ \end{array}\end{eqnarray} \tag{ B.10 }$$

In the limit $R\to 0$ the function ${\rm erf}\left( R \right)/R$ converges to $2/\sqrt{\pi }$. Thus, (B.9) converges to $-1/2$ and yields the same value as one gets for ${{\Phi }_{c}}$. For (B.10) both the second and third terms diverge but the sum of both terms converges and altogether ${{V}_{0,s}}$ also converges to the value of $-1/2$. So everything seems fine.

But now consider the behaviour of V(c) in a small neighbourhood of c = 0, i.e. ${{V}_{\varepsilon }}=V\left( c=\varepsilon \right)$. For the hollow sphere this changes nothing of course, since the potential is constant within radius R. The potentials ${{\Phi }_{c}}$ and ${{\Phi }_{s}}$ can be expanded around $\varepsilon =0$ and yield

$$\begin{eqnarray}&&{{\Phi }_{c}}(c+\varepsilon )={{\Phi }_{c}}(c)+\varepsilon \;\frac{1}{{{c}^{2}}}+\mathcal{O}({{\varepsilon }^{2}}),\end{eqnarray} \tag{ B.11 }$$

$$\begin{eqnarray}&&{{\Phi }_{s}}(c+\varepsilon )={{\Phi }_{s}}(c)+\varepsilon \;\frac{c}{{{R}^{3}}}+\mathcal{O}({{\varepsilon }^{2}})\;.\end{eqnarray} \tag{ B.12 }$$

This gives the additional contributions

$$\begin{eqnarray}&&{{V}_{\varepsilon ,c}}={{V}_{0,c}}+\varepsilon \;\int _{0}^{R}{\rm d}c\ {\rm exp} (-{{c}^{2}})\end{eqnarray} \tag{ B.13 }$$

$$\begin{eqnarray}&&\quad \ {\mkern 1mu} =\;{{V}_{0,c}}+\varepsilon \;\frac{\sqrt{\pi }}{2}\;{\rm erf}(R),\end{eqnarray} \tag{ B.14 }$$

$$\begin{eqnarray}&&{{V}_{\varepsilon ,s}}={{V}_{0,s}}+\frac{\varepsilon }{{{R}^{3}}}\;\int _{0}^{R}{\rm d}c\ {{c}^{3}}\;{\rm exp} (-{{c}^{2}})\end{eqnarray} \tag{ B.15 }$$

$$\begin{eqnarray}&&\quad \ {\mkern 1mu} =\;{{V}_{0,s}}-\varepsilon \;\frac{1+{{R}^{2}}}{2}\;{\rm exp} (-{{R}^{2}})+\frac{\varepsilon }{2{{R}^{3}}}\end{eqnarray} \tag{ B.16 }$$

to the potentials. For the Coulomb potential everything is fine since ${\rm erf}\left( R \right)\to 0$ for $R\to 0$. Hence, both the Coulomb and the hollow sphere potential obtain no further contributions at this order and it can be easily checked that this also holds for all higher orders in $\varepsilon$.

For the solid sphere potential, however, things are not fine at all. Not only does the term proportional to ${\rm exp} (-{{R}^{2}})$ in the limit $R\to 0$ yield a contribution $-\varepsilon /2$, which already makes it differ from the Coulomb potential; the last term is even worse because it diverges in this limit. Therefore, we cannot take this model seriously for small radii of the solid sphere and we are better off taking the hollow sphere potential as a toy model for the density of a molecule.

Here we will explicitly derive the self-contained system of evolution equations for the first and second moments given in [23]. It has been noted there that since this system is closed, Gaussian states will remain Gaussian under evolution. We will show that the difference of our equation (52) from equation (53) given in [23] has no influence on this set of equations.

For this we consider the Schrödinger equation

$$\begin{eqnarray}&&{\rm i}\dot{\psi }=\frac{{{p}^{2}}}{2\;M}\psi +{{H}_{1}}\psi ,\end{eqnarray} \tag{ C.1 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl} {{H}_{1}} & = & \frac{k}{2}{{{\bf x}}^{2}}-{{k}_{{\rm SN}}}{\bf x}\cdot \left\langle {\bf x} \right\rangle +\alpha {{\left\langle {\bf x} \right\rangle }^{2}}+\beta \left\langle {{{\bf x}}^{2}} \right\rangle \\ {{p}_{i}} & = & -{\rm i}{{\partial }_{i}} \\ k & = & {{k}_{{\rm CM}}}+{{k}_{{\rm SN}}}=M\omega _{{\rm CM}}^{2}+M\omega _{{\rm SN}}^{2}. \\ \end{array}\end{eqnarray*}$$

In principle, in the case of equation (52) we have $\alpha =0$, $\beta ={{k}_{{\rm SN}}}/2$, while in the case of equation (53) $\alpha ={{k}_{{\rm SN}}}/2$, $\beta =0$. But note that

$$\begin{eqnarray}&&{{\partial }_{i}}{{\left\langle {\bf x} \right\rangle }_{j}}=0,\end{eqnarray} \tag{ C.2a }$$

$$\begin{eqnarray}&&{{\partial }_{i}}\langle {{{\bf x}}^{2}}\rangle =0.\end{eqnarray} \tag{ C.2b }$$

Therefore, independent of the choice of $\alpha$ and $\beta$ the derivatives of ${{H}_{1}}$ are

$$\begin{eqnarray}&&{{\partial }_{i}}{{H}_{1}}=k{{x}_{i}}-{{k}_{{\rm SN}}}{{\left\langle {\bf x} \right\rangle }_{i}},\end{eqnarray} \tag{ C.3a }$$

$$\begin{eqnarray}&&\Delta {{H}_{1}}=3k.\end{eqnarray} \tag{ C.3b }$$

We will see that ${{H}_{1}}$ will enter into the evolution equations for the first and second moments only through these derivatives. Thus, for the different equations (52) and (53) we obtain the same evolution equations for the first and second moments, which are:

$$\begin{eqnarray}&&{{\partial }_{t}}{{\left\langle {\bf x} \right\rangle }_{i}}=\int {{{\rm d}}^{3}}x\;{{x}_{i}}\left( {{\psi }^{*}}\dot{\psi }+\psi {{{\dot{\psi }}}^{*}} \right)=\frac{{\rm i}}{2\;M}\int {{{\rm d}}^{3}}x\;{{x}_{i}}\left( {{\psi }^{*}}\Delta \psi -\psi \Delta {{\psi }^{*}} \right)\end{eqnarray} \tag{ C.4a }$$

$$\begin{eqnarray}&&\qquad \ {\mkern 1mu} =\frac{{\rm i}}{2\;M}\int {{{\rm d}}^{3}}x\left( -({{\partial }_{j}}\psi ){{\partial }_{j}}({{x}_{i}}{{\psi }^{*}})+({{\partial }_{j}}{{\psi }^{*}}){{\partial }_{j}}({{x}_{i}}\psi ) \right)\end{eqnarray} \tag{ C.4b }$$

$$\begin{eqnarray}&&\qquad \ {\mkern 1mu} =\frac{{\rm i}}{2\;M}\int {{{\rm d}}^{3}}x\left( -{{\psi }^{*}}{{\partial }_{i}}\psi +\psi {{\partial }_{i}}{{\psi }^{*}} \right)=\frac{1}{M}\int {{{\rm d}}^{3}}x\;{{\psi }^{*}}(-{\rm i}{{\partial }_{i}})\psi \end{eqnarray} \tag{ C.4c }$$

$$\begin{eqnarray}&&\qquad \ {\mkern 1mu} =\frac{{{\left\langle {\bf p} \right\rangle }_{i}}}{M}\end{eqnarray} \tag{ C.4d }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{\left\langle {\bf p} \right\rangle }_{i}}=-{\rm i}\int {{{\rm d}}^{3}}x\left( {{{\dot{\psi }}}^{*}}{{\partial }_{i}}\psi +{{\psi }^{*}}{{\partial }_{i}}\dot{\psi } \right)\end{eqnarray} \tag{ C.5a }$$

$$\begin{eqnarray}&&\qquad \ \ =\;\frac{1}{2\;M}\int {{{\rm d}}^{3}}x\left( -(\Delta {{\psi }^{*}}){{\partial }_{i}}\psi +{{\psi }^{*}}{{\partial }_{i}}\Delta \psi \right)-\left\langle {{\partial }_{i}}{{H}_{1}} \right\rangle \end{eqnarray} \tag{ C.5b }$$

$$\begin{eqnarray}&&\qquad \ \ =\;\frac{1}{2\;M}\int {{{\rm d}}^{3}}x\left( ({{\partial }_{j}}{{\psi }^{*}}){{\partial }_{i}}{{\partial }_{j}}\psi +{{\psi }^{*}}{{\partial }_{i}}\Delta \psi \right)-\left\langle {{\partial }_{i}}{{H}_{1}} \right\rangle \end{eqnarray} \tag{ C.5c }$$

$$\begin{eqnarray}&&\qquad \ \ =\;\frac{1}{2\;M}\int {{{\rm d}}^{3}}x\left( -{{\psi }^{*}}{{\partial }_{i}}\Delta \psi +{{\psi }^{*}}{{\partial }_{i}}\Delta \psi \right)-\left\langle {{\partial }_{i}}{{H}_{1}} \right\rangle \end{eqnarray} \tag{ C.5d }$$

$$\begin{eqnarray}&&\qquad \ \ =-k{{\left\langle {\bf x} \right\rangle }_{i}}+{{k}_{{\rm SN}}}{{\left\langle {\bf x} \right\rangle }_{i}}\end{eqnarray} \tag{ C.5e }$$

$$\begin{eqnarray}&&\qquad \ \ =-{{k}_{{\rm CM}}}{{\left\langle {\bf x} \right\rangle }_{i}}\end{eqnarray} \tag{ C.5f }$$

$$\begin{eqnarray}&&{{\partial }_{t}}\left\langle {{{\bf x}}^{2}} \right\rangle =\int {{{\rm d}}^{3}}x\;{{{\bf x}}^{2}}\left( {{\psi }^{*}}\dot{\psi }+\psi {{{\dot{\psi }}}^{*}} \right)=\frac{{\rm i}}{2\;M}\int {{{\rm d}}^{3}}x\;{{{\bf x}}^{2}}\left( {{\psi }^{*}}\Delta \psi -\psi \Delta {{\psi }^{*}} \right)\end{eqnarray} \tag{ C.6a }$$

$$\begin{eqnarray}&&\qquad \ \ {\mkern 1mu} =\frac{{\rm i}}{2\;M}\int {{{\rm d}}^{3}}x\left( -({{\partial }_{j}}\psi ){{\partial }_{j}}({{x}_{i}}{{x}_{i}}{{\psi }^{*}})+({{\partial }_{j}}{{\psi }^{*}}){{\partial }_{j}}({{x}_{i}}{{x}_{i}}\psi ) \right)\end{eqnarray} \tag{ C.6b }$$

$$\begin{eqnarray}&&\qquad \ \ {\mkern 1mu} =\frac{{\rm i}}{M}\int {{{\rm d}}^{3}}x\left( -{{x}_{i}}{{\psi }^{*}}{{\partial }_{i}}\psi +{{x}_{i}}\psi {{\partial }_{i}}{{\psi }^{*}} \right)\end{eqnarray} \tag{ C.6c }$$

$$\begin{eqnarray}&&\qquad \ \ {\mkern 1mu} =\frac{1}{M}\int {{{\rm d}}^{3}}x\left( {{\psi }^{*}}{{x}_{i}}(-{\rm i}{{\partial }_{i}})\psi +{{\psi }^{*}}(-{\rm i}{{\partial }_{i}})({{x}_{i}}\psi ) \right)\end{eqnarray} \tag{ C.6d }$$

$$\begin{eqnarray}&&\qquad \ \ {\mkern 1mu} =\frac{1}{M}\left( \left\langle {\bf x}\cdot {\bf p} \right\rangle +\left\langle {\bf p}\cdot {\bf x} \right\rangle \right)\end{eqnarray} \tag{ C.6e }$$

$$\begin{eqnarray}\begin{array}{rcl} {{\partial }_{t}}\left\langle {{{\bf p}}^{2}} \right\rangle & = & -\int {{{\rm d}}^{3}}x\left( {{{\dot{\psi }}}^{*}}\Delta \psi +{{\psi }^{*}}\Delta \dot{\psi } \right) \\ {} & = & \frac{{\rm i}}{2\;M}\int {{{\rm d}}^{3}}x\left( (\Delta {{\psi }^{*}})\Delta \psi -{{\psi }^{*}}\Delta \Delta \psi \right)+2{\rm i}\int {{{\rm d}}^{3}}x\;{{\psi }^{*}}({{\partial }_{i}}{{H}_{1}}){{\partial }_{i}}\psi +{\rm i}\left\langle \Delta {{H}_{1}} \right\rangle \\ \end{array}\end{eqnarray} \tag{ C.7a }$$

$$\begin{eqnarray}&&\qquad \ \ \;=\;2{\rm i}\int {{{\rm d}}^{3}}x\;{{\psi }^{*}}\left( k{{x}_{i}}-{{k}_{{\rm SN}}}{{\left\langle {\bf x} \right\rangle }_{i}} \right){{\partial }_{i}}\psi +3{\rm i}k\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \end{eqnarray} \tag{ C.7b }$$

$$\begin{eqnarray}&&\qquad \ \ \;=-2k\int {{{\rm d}}^{3}}x\;{{\psi }^{*}}{{x}_{i}}(-{\rm i}{{\partial }_{i}})\psi +2{{k}_{{\rm SN}}}{{\left\langle {\bf x} \right\rangle }_{i}}\int {{{\rm d}}^{3}}x\;{{\psi }^{*}}(-{\rm i}{{\partial }_{i}})\psi +3{\rm i}k\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{eqnarray} \tag{ C.7c }$$

$$\begin{eqnarray}&&\qquad \ \ \;=-2k\left\langle {\bf x}\cdot {\bf p} \right\rangle +2{{k}_{{\rm SN}}}\left\langle {\bf x} \right\rangle \cdot \left\langle {\bf p} \right\rangle +k\left\langle {\bf x}\cdot {\bf p} \right\rangle -k\left\langle {\bf p}\cdot {\bf x} \right\rangle \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{eqnarray} \tag{ C.7d }$$

$$\begin{eqnarray}&&\qquad \ \ \;=-k\left( \left\langle {\bf x}\cdot {\bf p} \right\rangle +\left\langle {\bf p}\cdot {\bf x} \right\rangle \right)+2{{k}_{{\rm SN}}}\left\langle {\bf x} \right\rangle \cdot \left\langle {\bf p} \right\rangle \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{eqnarray} \tag{ C.7e }$$

$$\begin{eqnarray}&&{{\partial }_{t}}\left\langle {\bf x}\cdot {\bf p} \right\rangle ={{\partial }_{t}}\left\langle {\bf p}\cdot {\bf x} \right\rangle =\int {{{\rm d}}^{3}}x\left( {{{\dot{\psi }}}^{*}}{{x}_{i}}(-{\rm i}{{\partial }_{i}})\psi +{{\psi }^{*}}{{x}_{i}}(-{\rm i}{{\partial }_{i}})\dot{\psi } \right)\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{eqnarray} \tag{ C.8a }$$

$$\begin{eqnarray}&&\qquad \quad \ \ \;=\;-\frac{1}{2\;M}\int {{{\rm d}}^{3}}x\;{{x}_{i}}\left( (\Delta {{\psi }^{*}}){{\partial }_{i}}\psi -{{\psi }^{*}}{{\partial }_{i}}\Delta \psi \right)-\left\langle {{x}_{i}}{{\partial }_{i}}{{H}_{1}} \right\rangle \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{eqnarray} \tag{ C.8b }$$

$$\begin{eqnarray}&&\qquad \quad \ \ \;=\;\frac{1}{2\;M}\int {{{\rm d}}^{3}}x({{\partial }_{j}}{{\psi }^{*}}){{\partial }_{j}}({{x}_{i}}{{\partial }_{i}}\psi )+\frac{1}{2\;M}\int {{{\rm d}}^{3}}x\;{{x}_{i}}{{\psi }^{*}}{{\partial }_{i}}\Delta \psi -\left\langle {{x}_{i}}{{\partial }_{i}}{{H}_{1}} \right\rangle \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{eqnarray} \tag{ C.8c }$$

$$\begin{eqnarray}&&\qquad \quad \ \ \;=\;\frac{1}{2\;M}\int {{{\rm d}}^{3}}x\left( ({{\partial }_{i}}{{\psi }^{*}}){{\partial }_{i}}\psi +{{x}_{i}}({{\partial }_{j}}{{\psi }^{*}}){{\partial }_{i}}{{\partial }_{j}}\psi +{{x}_{i}}{{\psi }^{*}}{{\partial }_{i}}\Delta \psi \right)-\left\langle {{x}_{i}}{{\partial }_{i}}{{H}_{1}} \right\rangle \end{eqnarray} \tag{ C.8d }$$

$$\begin{eqnarray}&&\qquad \quad \ \ \;=\;\frac{1}{2\;M}\int {{{\rm d}}^{3}}x\left( -{{\psi }^{*}}\Delta \psi -{{\psi }^{*}}{{\partial }_{j}}({{x}_{i}}{{\partial }_{i}}{{\partial }_{j}}\psi )+{{x}_{i}}{{\psi }^{*}}{{\partial }_{i}}\Delta \psi \right)-\left\langle {{x}_{i}}{{\partial }_{i}}{{H}_{1}} \right\rangle \end{eqnarray} \tag{ C.8e }$$

$$\begin{eqnarray}&&\qquad \quad \ \ \;=\;-\frac{1}{M}\int {{{\rm d}}^{3}}x\;{{\psi }^{*}}\Delta \psi -\left\langle {{x}_{i}}{{\partial }_{i}}{{H}_{1}} \right\rangle \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{eqnarray} \tag{ C.8f }$$

$$\begin{eqnarray}&&\qquad \quad \ \ \;=\;\frac{\left\langle {{{\bf p}}^{2}} \right\rangle }{M}-k\left\langle {{{\bf x}}^{2}} \right\rangle +{{k}_{{\rm SN}}}{{\left\langle {\bf x} \right\rangle }^{2}}\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{eqnarray} \tag{ C.8g }$$

The same evolution equations are obtained by operators ${\bf x}$ and ${\bf p}$ that in the Heisenberg picture fulfil

$$\begin{eqnarray}&&{{\partial }_{t}}{\bf x}=\frac{{\bf p}}{M}\end{eqnarray} \tag{ C.9a }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{\bf p}=-{{k}_{{\rm CM}}}{\bf x}-{{k}_{{\rm SN}}}({\bf x}-\left\langle {\bf x} \right\rangle ).\end{eqnarray} \tag{ C.9b }$$

This was used in [23] to describe the effect of the Schrödinger–Newton equation on Gaussian states.