  *

# Random walk and non-Gaussianity of the 3D second-quantized Schrödinger-
Newton nonlocal soliton

Claudio Conti  Department of Physics, University Sapienza, Piazzale Aldo Moro
5, 00185 Rome, Italy  Institute for Complex Systems, National Research Council
(ISC-CNR), Via dei Taurini 19, 00185 Rome, Italy  Research Center Enrico
Fermi, Via Panisperna 89a, 00184 Rome, Italy<EMAIL_ADDRESS>
###### Abstract

Nonlocal quantum fluids emerge as dark-matter models and tools for quantum
simulations and technologies. However, strongly nonlinear regimes, like those
involving multi-dimensional self-localized solitary waves, are marginally
explored for what concerns quantum features. We study the dynamics of 3D+1
solitons in the second-quantized nonlocal nonlinear Schrödinger-Newton
equation. We theoretically investigate the quantum diffusion of the soliton
center of mass and other parameters, varying the interaction length. 3D+1
simulations of the Ito partial differential equations arising from the
positive P-representation of the density matrix validate the theoretical
analysis. The numerical results unveil the onset of non-Gaussian statistics of
the soliton, which may signal quantum-gravitational effects and be a resource
for quantum computing. The non-Gaussianity arises from the interplay between
the soliton parameter quantum diffusion and the stable invariant propagation.
The fluctuations and the non-Gaussianity are universal effects expected for
any nonlocality and dimensionality.

Keywords: Nonlocal solitons, dark matter, positive P-representation, quantum
gravity, non-Gaussian statistics, quantum simulations.

##  1 Introduction

Three-dimensional (3D) self-localized nonlinear waves enter various fields of
research [1, 2], but their quantum properties are unexplored. Classical three-
dimensional solitary waves (in short, 3D solitons) need to be stabilized
against catastrophic collapse. Nonlocality is a well-known mechanism for the
stabilization [3, 4, 5] and nonlocal soliton are a fascinating research
direction involving long-range Bose-Einstein condensates (BECs) [6, 7, 8, 9],
boson stars [10] and dark-matter models [11, 12]. However, a mean-field
description that overlooks quantum effects provides limited information on the
dynamics of self-trapped multidimensional waves. This limitation is
specifically relevant as recent investigations suggest the solitons as non-
classical sources for quantum technologies and fundamental studies [13, 14,
15, 16, 17, 18]. Results in 1D [19, 20] suggest that nonlocality frustrates
fluctuations. However, despite ab-initio investigations on long-range
interactions [21, 22], the quantum statistics of self-trapped 3D nonlocal
solitons is an open issue.

In addition, recent work on gravitational interaction in BEC predicts non-
Gaussian statistics [23]. Non-Gaussianity is a resource for continuous-
variable quantum information science [24, 25] and its understanding in quantum
fluids may enable new universal quantum processors. Also, emerging of non-
Gaussian statistics in table-top experiments may open the way to study - or
simulate - quantum gravity in the laboratory. Ref. [23] predicts that a BEC in
a trap, once prepared in a squeezed state or Schrödinger-cat state, triggers
the non-Gaussian statistics measured by a signal-to-noise ratio (SNR)
parameter, which reveals quantized gravity. However - so far - no experiments
or numerical simulations validate these theoretical predictions. Also, quantum
fluctuations and non-Gaussianity in multidimensional self-trapped solitonic
nonlocal condensates have never been considered before.

Here, we study theoretically and numerically the quantum dynamics 3D nonlocal
solitons. We use a perturbative approach and we analytically predict the
quantum diffusion of the soliton position and other parameters. We validate
our analytical results by ab-initio numerical simulations based on the 3D+1
positive P-representation [26, 21, 22]. We compute the SNR parameter
introduced in [23], which shows that non-Gaussianity arises in the quantum
dynamics of 3D+1 nonlocal solitons, starting from a coherent state.

##  2 Model and scaling

We consider the many-body Hamiltonian

$\hat{H}=\frac{\hbar^{2}}{2m}\int\nabla\hat{\psi}^{\dagger}\cdot\nabla\hat{\psi}\,\mathrm{d}^{3}{\bf
x}+\int U({\bf x}-{\bf x}^{\prime})\hat{\psi}^{\dagger}({\bf
x}^{\prime})\hat{\psi}({\bf x})^{\dagger}\hat{\psi}({\bf
x}^{\prime})\hat{\psi}({\bf x})\,\mathrm{d}^{3}{\bf x}\,\mathrm{d}^{3}{\bf
x}^{\prime}\;,$ (1)

with $m$ is the boson mass, and $U$ is the interaction potential. We adopt the
phase-space representation methods [27, 21] for studying the nonlocal
interaction. The quantum field model is equivalent to a Fokker-Planck
equation, which is mapped to Ito nonlinear partial differential equations
coupling two fields $\psi$, and $\psi^{+}$

$\begin{array}[]{lll}\imath\hbar\partial_{t}\psi&=&-\frac{\hbar^{2}}{2m}\nabla^{2}\psi+\psi\,U*\rho+\sqrt{\imath\hbar}\psi\,\xi_{U}\\\\[4.0pt]
-\imath\hbar\partial_{t}\psi^{+}&=&-\frac{\hbar^{2}}{2m}\nabla^{2}\psi^{+}+\psi^{+}\,U*\rho+\sqrt{-\imath\hbar}\psi^{+}\,\xi_{U}^{+}\end{array}$
(2)

where the asterisk denotes a convolution integral. In (2) $\rho=\psi^{+}\psi$,
$\xi_{U}$ and $\xi_{U}^{+}$ are independent noises such that

$\langle\xi_{U}({\bf x},t)\xi_{U}({\bf
x}^{\prime},t^{\prime})\rangle=\langle\xi_{U}^{+}({\bf x},t)\xi_{U}^{+}({\bf
x}^{\prime},t^{\prime})\rangle=U({\bf x}-{\bf
x}^{\prime})\delta(t-t^{\prime})\;.$ (3)

The total number of particles is $\int\psi\psi^{+}\,\mathrm{d}V$, its mean
value is $N_{T}=\int\langle\psi\psi^{+}\rangle\,\mathrm{d}V$; the brackets
here denote the mean-field solution obtained with $\xi_{U}=\xi_{U}^{+}=0$.

In our numerical calculations below, we consider self-gravitating screened
potential $U=-Gm^{2}e^{-r/\Lambda}/r$, where $\Lambda$ is the interaction
length. $G$ measures the coupling corresponding to the gravitational constant,
but $U$ also models other long-range interactions as, e.g., thermal effects in
photonic BEC [8]. The mean-field theory is obtained by
$\xi_{U}=\xi^{+}_{U}=0$, and $\psi^{+}=\psi^{*}$, and corresponds to the
Schrödinger-Newton equation [28].

Figure 1: Self-gravitating solitonic core. (a) Comparison of time dynamics
with and without interaction ($G=0$); 3D isodensity surfaces at different
instants for freely evolving fields (top panel) and in the presence of the
nonlinearity (bottom panel) with the time-invariant self-trapped wave-packet.
(b) Two-dimensional projection (average in the $z-$direction) of the density
profile (blue) and resulting long-range potential (yellow)  Figure 2:
Classical and quantum evolution of the 3D+1 nonlocal soliton. (a) Isosurfaces
of the density $\rho$ of the solitonic core at different instants. We show a
single run (top panel), an average of $10$ runs, and the classical propagation
invariant solution (bottom panel). (b) Diffusion coefficient $D_{X}$ after
(56) for various degrees of nonlocality $\sigma$. (c,d,e) Computed
trajectories for the $X^{1}(t)$ displacement for $100$ runs for three values
of $\sigma$. The thick line is the standard deviation $\langle
X^{1}(t)^{2}\rangle^{1/2}$, the dashed line is (56) for comparison with theory
without fitting parameters.

We write the stochastic equations in dimensionless units by letting

$\displaystyle(x,y,z)$ $\displaystyle\rightarrow$ $\displaystyle(x,y,z)r_{0}$
$\displaystyle t$ $\displaystyle\rightarrow$ $\displaystyle tt_{0}$
$\displaystyle\left(\psi,\psi^{+}\right)$ $\displaystyle\rightarrow$
$\displaystyle\psi_{0}\left(\psi,\psi^{+}\right)$

and

$\displaystyle t_{0}$ $\displaystyle=$ $\displaystyle 2mr_{0}^{2}/\hbar\;,$
(4) $\displaystyle\psi_{0}^{2}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}}{2Gm^{3}r_{0}^{4}}\;,$ (5) $\displaystyle n_{0}$
$\displaystyle=$ $\displaystyle\frac{\hbar^{2}}{2Gm^{3}r_{0}}=\frac{N_{T}}{M}$
(6)

being

$M=\langle\int\psi\psi^{+}\,\mathrm{d}^{3}{\bf x}\rangle\;.$ (7)

$n_{0}$ measures the number of particles in the condensate in units of $M$,
the norm of the numerically obtained bound state profile [Eq. (13) below].

In the dimensionless units, Eqs. (2) read

$\begin{array}[]{lll}+\imath\partial_{t}\psi+\nabla^{2}\psi-\psi\,U\ast\rho&=s\\\
-\imath\partial_{t}\psi^{+}+\nabla^{2}\psi^{+}-\psi^{+}\,U\ast\rho&=s^{+}\end{array}$
(8)

with

$U=-\frac{e^{-r/\sigma}}{r}\;,$ (9)

$\sigma=\Lambda/r_{0}$, and

$\begin{array}[]{lll}s=\sqrt{\frac{\imath}{n_{0}}}\xi_{U}\psi\;,\\\\[4.0pt]
s^{+}=\sqrt{\frac{-\imath}{n_{0}}}\xi^{+}_{U}\psi^{+}\;.\end{array}$ (10)

$\xi_{U}(x,y,z,t)$ and $\xi_{U}^{+}(x,y,z,t)$ are uncorrelated noise terms
such that

$\begin{array}[]{lll}\langle\xi_{U}({\bf x},t)\xi_{U}({\bf
x}^{\prime},t^{\prime})\rangle=\langle\xi^{+}({\bf x},t)\xi^{+}({\bf
x}^{\prime},t^{\prime})\rangle=U({\bf x}-{\bf
x}^{\prime})\delta(t-t^{\prime}).\end{array}$ (11)

$\sigma$ and $n_{0}$ are the dimensionless interaction length and particle
number, respectively.

According to eq. (6), one can either fix $r_{0}$ or $n_{0}$ to set all the
other normalization constants. We choose to use $n_{0}$ because it appears
explicitly in the normalized equation in a way such that the limit
$n_{0}\rightarrow\infty$ corresponds to the mean-field regime. Indeed, the
total mean particle number is $N_{T}=n_{0}M$.

Once we have the numerical solution of the bound state [Eq. (13) below], which
is determined by the scale $\sigma$, we study its quantum fluctuations by
numerically solving Eqs. (8). In this paper, we fix a specific value for
$n_{0}$, which allows us to perform numerical simulations with unitary time-
scale in our normalized scale (i.e., $t\simeq 1$ as in figure 2 ), and we
study the effects of a varying interaction length $\sigma$.

##  3 Self-gravitating non-local soliton

In the mean-field theory, equations (8) admit a stable radially-symmetric
bound-state solution: a self-localized three-dimensional solitary wave. We
write the solution with a Galileian boost as

$\psi=u(x^{a}-X^{a})\exp\left[\imath\theta-\imath
Et+\frac{i}{2}V^{a}(x_{a}-X_{a})\right]\;,$ (12)

$\psi^{+}=\psi^{*}$, with $a=1,2,3$, $x^{1}=x$, $x^{2}=y$, $x^{3}=z$, omitting
the sum symbol over repeated Latin indices. $u(x^{a})$ is the real-valued
soliton profile, such that

$\Delta u-U\ast u^{2}u=Eu\;.$ (13)

The soliton energy $E$ is time-independent. For the position $X^{a}=X^{a}(t)$,
we have (dot is the time-derivative)

$\begin{array}[]{lll}\dot{X}^{a}&=&V^{a}\\\ \dot{V}^{a}&=&0\\\
\dot{\theta}&=&\frac{1}{4}V^{2},\end{array}$ (14)

with $V^{2}=\delta_{ab}V^{a}V^{b}$, and $\delta_{ab}$ the Kronecker symbol
($b=1,2,3$). Equations (14) imply

$\begin{array}[]{lll}X^{a}&=X^{a}(t)=X^{a}(0)+V^{a}t\;,\\\
\theta&=\theta(t)=\theta(0)+\frac{1}{4}V^{2}\;.\end{array}$ (15)

Figure 1 shows the evolution of the time-invariant soliton profile obtained
from Eq. (13), compared with the evolution in the absence of nonlinearity
($U=0$). We also show the soliton compared with corresponding potential
$U\ast\rho$. The field profile and the potential are computed numerically. We
use a pseudo-spectral parallel relaxation procedure in a 3D Cartesian domain.
Figure 1 shows the calculated classical bound state $u$. In the absence of
interaction, the mass spreads upon evolution. In the presence of self-
attraction, the solitonic wave packet is invariant upon propagation.

##  4 Quantum effects on the 3D nonlocal soliton

In the quantum regime, with $\xi_{U}\neq 0$ and $\xi_{U}^{+}\neq 0$, the
soliton, initially prepared in a coherent state, evolves with fluctuations
depending on the interaction length $\sigma$. The 3D+1 stochastic partial
differential equations in (8) are solved by following Drummond and coworkers
[21]. We adopt an iterative stochastic solver with pseudospectral
discretization and parallelized with the FFTW [29] and the Message Passing
Interface (MPI) protocol. Figure 2a shows the numerical solution of the
stochastic equations (8), which unveils that the soliton undergoes a random
walk (figure 2c-e).

To study the quantum regime, we derive equations for the soliton parameters by
(8) using soliton perturbation theory. In the presence of noise, Eqs. (14) are
replaced by stochastic differential equations, which we derive by introducing
the vectorial notation

${\bm{\psi}}=\left(\begin{array}[]{c}\psi\\\ \psi^{+}\end{array}\right)\;.$
(16)

Equations (8) are written as

$\imath\sigma_{3}\partial_{t}{\bm{\psi}}+\Delta{\bm{\psi}}-U\ast(\psi\psi^{+}){\bm{\psi}}={\bm{s}}\;,$
(17)

with the Pauli matrix

$\sigma_{3}=\left(\begin{array}[]{ll}1&0\\\ 0&-1\end{array}\right)\;,$ (18)

and

${\bf s}=\left(\begin{array}[]{c}s\\\ s^{+}\end{array}\right)\;.$ (19)

We introduce the following vector

${\bf e}=\left(\begin{array}[]{lll}e^{\imath\theta-\imath
Et+\frac{i}{2}V^{a}(x_{a}-X_{a})}\\\ e^{-\imath\theta+\imath
Et-\frac{i}{2}V^{a}(x_{a}-X_{a})}\end{array}\right)\;,$ (20)

being $u=u(x^{a}-X^{a},E)$ solution of Eq. (13). We also define

$\displaystyle{\bf f}_{\theta}$ $\displaystyle=$ $\displaystyle{\bf e}u$ (21)
$\displaystyle{\bf f}_{E}$ $\displaystyle=$ $\displaystyle\imath\sigma_{3}{\bf
e}\frac{\partial u}{\partial E}\;,$ (22) $\displaystyle{\bf f}_{X}^{a}$
$\displaystyle=$ $\displaystyle-\imath\sigma_{3}{\bf e}\frac{\partial
u}{\partial x^{a}}\;,$ (23) $\displaystyle{\bf f}_{V}^{a}$ $\displaystyle=$
$\displaystyle{\bf e}\frac{1}{2}(x^{a}-X^{a})u\;.$ (24)

We introduce a scalar product for two vectors f and g such that

$\left({\bf f},{\bf g}\right)=2\Re\int{\bf f}^{*}\cdot{\bf g}\,\mathrm{d}V\;.$
(25)

By using this scalar product, we build a bi-orthogonal system by introducing
the conjugate vectors to (24)

$\displaystyle\hat{\bf f}_{\theta}$ $\displaystyle=$
$\displaystyle\imath\sigma_{3}{\bf f}_{E}$ (26) $\displaystyle\hat{\bf f}_{E}$
$\displaystyle=$ $\displaystyle-\imath\sigma_{3}{\bf f}_{\theta}$ (27)
$\displaystyle\hat{\bf f}_{X}^{a}$ $\displaystyle=$
$\displaystyle\imath\sigma_{3}{\bf f}_{V}^{a}$ (28) $\displaystyle\hat{\bf
f}_{V}^{a}$ $\displaystyle=$ $\displaystyle-\imath\sigma_{3}{\bf
f}_{X}^{a}\;.$ (29)

We have

$\left(\hat{{\bf f}}_{X}^{a},{\bf f}_{X}^{b}\right)=\left(\hat{{\bf
f}}_{V}^{a},{\bf f}_{V}^{b}\right)=M\delta_{ab}$ (30)

and

$\begin{array}[]{lll}\left(\hat{{\bf f}}_{E},{\bf
f}_{E}\right)=\left(\hat{{\bf f}}_{\theta},{\bf
f}_{\theta}\right)=\frac{\,\mathrm{d}M}{\,\mathrm{d}E}\\\\[4.0pt]
\left(\hat{{\bf f}}_{\theta},{\bf f}_{E}\right)=\left(\hat{{\bf f}}_{E},{\bf
f}_{\theta}\right)=0\;,\end{array}$ (31)

with

$M=\int u^{2}\,\mathrm{d}^{3}{\bf x}\;$ (32)

and all the other scalar products are vanishing.

In the presence of the quantum noise ${\bf s}$, we assume that all the soliton
parameters are time-dependent, and using (12) and we have after (8)

$\begin{array}[]{lll}\imath\sigma_{3}\partial_{t}{\bm{\psi}}+\Delta{\bm{\psi}}-U\ast(\psi\psi^{+}){\bm{\psi}}={\bf
s}\\\ ={\bf
f}_{\theta}\left(-\dot{\theta}+t\dot{E}+\frac{1}{2}V^{a}\dot{X}^{a}-\frac{V^{2}}{4}\right)+{\bf
f}_{E}\dot{E}+{\bf f}_{X}^{a}\left(\dot{X}^{a}-V^{a}\right)+{\bf
f}_{V}^{a}\left(-\dot{V}^{a}\right)\end{array}$ (33)

where the dot indicates the time derivative. Equations (33) are valid at the
lowest order of perturbation, higher orders can be determined by radiative
corrections to the soliton profile. By scalar multiplying by $\hat{{\bf
f}}_{X^{a}}$ and $\hat{{\bf f}}_{V^{a}}$, we obtain the stochastic equations
for the position and the velocity of the soliton

$\displaystyle M\dot{X}^{a}=MV^{a}+\left(\hat{{\bf f}}_{X}^{a},{\bf s}\right)$
(34) $\displaystyle M\dot{V}^{a}=-\left(\hat{{\bf f}}_{V}^{a},{\bf
s}\right)\;.$ (35)

Equations (34) and (35) describe the dynamics of the soliton position and
velocity with quantum noise. Seemingly, we get equations for $\theta$ and $E$.
At the lowest order in $t$ and $V^{a}=X^{a}=0$ at $t=0$, we have

$\displaystyle M^{\prime}\dot{E}$ $\displaystyle=$
$\displaystyle\left(\hat{{\bf f}}_{E},{\bf s}\right)$ (36) $\displaystyle
M^{\prime}\dot{\theta}$ $\displaystyle=$ $\displaystyle-\left(\hat{{\bf
f}}_{\theta}^{a},{\bf s}\right)\;.$ (37)

being

$M^{\prime}=\frac{dM}{dE}$ (38)

##  5 Quantum-induced parameter diffusion and random walk

We have for the perturbation vector ${\bf s}$,

${\bf s}=\sqrt{\frac{\imath}{n_{0}}}\left(\begin{array}[]{lll}\xi_{U}\\\
\imath\xi_{U}^{+}\end{array}\right)u{\bf e}\;.$ (39)

As detailed in A, by using (39), (34), (35) and (29), we obtain

$\displaystyle M\dot{X}^{a}$ $\displaystyle=$ $\displaystyle
MV^{a}+F^{a}_{X}({\bf X},t)$ (40) $\displaystyle M\dot{V}^{a}$
$\displaystyle=$ $\displaystyle F^{a}_{V}({\bf X},t)\;$ (41)

where $F^{a}_{X}$ and $F^{a}_{V}$, with $a=1,2,3$, are stochastic terms acting
on the position $X^{a}$ and velocity $V^{a}$ of the soliton with mass $M$. We
have

$\displaystyle F_{a}^{X}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{n_{0}}}\int\rho({\bf x}-{\bf
X})(x^{a}-X^{a})\xi_{+}({\bf x},t)\,\mathrm{d}^{3}{\bf x}$ (42) $\displaystyle
F_{a}^{V}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{n_{0}}}\int\frac{\partial\rho}{\partial
x^{a}}\left({\bf x}-{\bf X}\right)\xi_{+}({\bf x},t)\,\mathrm{d}^{3}{\bf x}\;$
(43)

with ${\bf X}=(X^{1},X^{2},X^{3})$ the soliton position, $\rho=u^{2}$, and
$\xi_{+}({\bf x},t)$ a real noise such that

$\langle\xi_{+}({\bf x}^{\prime},t^{\prime})\xi_{+}({\bf x},t)\rangle=-U({\bf
x}-{\bf x}^{\prime})\delta(t-t^{\prime})\;.$ (44)

For the stochastic terms in Eqs. (41) we have

$\displaystyle\langle F_{X}^{a}({\bf X},t)F_{X}^{b}({\bf
X},t^{\prime})\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{n_{0}}Q_{X}^{a}\delta_{ab}\delta(t-t^{\prime})\;,$ (45)
$\displaystyle\langle F_{V}^{a}({\bf X},t)F_{V}^{b}({\bf
X},t^{\prime})\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{n_{0}}Q_{V}^{a}\delta_{ab}\delta(t-t^{\prime})\;,$ (46)
$\displaystyle\langle F_{X}^{a}({\bf X},t)F_{V}^{b}({\bf
X},t^{\prime})\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{n_{0}}Q_{XV}^{a}\delta_{ab}\delta(t-t^{\prime})\;,$
(47)

with the correlation coefficients

$\displaystyle Q_{X}^{a}$ $\displaystyle=$ $\displaystyle-\int
x_{1}^{a}x_{2}^{a}\rho({\bf x}_{1})\rho({\bf x}_{2})U({\bf x}_{1}-{\bf
x}_{2})\,\mathrm{d}^{3}{\bf x}_{1}\,\mathrm{d}^{3}{\bf x}_{2}$ (48)
$\displaystyle Q_{V}^{a}$ $\displaystyle=$
$\displaystyle-\int\frac{\partial\rho({\bf x}_{1})}{\partial
x_{1}^{a}}\frac{\partial\rho({\bf x}_{2})}{\partial x_{2}^{a}}U({\bf
x}_{1}-{\bf x}_{2})\,\mathrm{d}^{3}{\bf x}_{1}\,\mathrm{d}^{3}{\bf x}_{2}$
(49) $\displaystyle Q_{XV}^{a}$ $\displaystyle=$ $\displaystyle-\int
x_{1}^{a}\rho({\bf x}_{1})\frac{\partial\rho({\bf x}_{2})}{\partial
x_{2}^{a}}U({\bf x}_{1}-{\bf x}_{2})\,\mathrm{d}^{3}{\bf
x}_{1}\,\mathrm{d}^{3}{\bf x}_{2}\;.$ (50)

Equations (40) and (41) with $X^{a}(0)=V^{a}(0)=0$ give for the moments

$\langle{\left[V^{a}(t)\right]}^{2}\rangle=D_{V}^{a}t\;,$ (51)

and

$\langle{\left[X^{a}(t)\right]}^{2}\rangle=D_{X}^{a}t+D_{XV}^{a}t^{2}+D_{V}^{a}\frac{t^{3}}{3}\;.$
(52)

The velocity and the position undergo a diffusive random walk with

$\displaystyle D_{X}^{a}$ $\displaystyle=$
$\displaystyle\frac{Q_{X}^{a}}{n_{0}M^{2}}$ (53) $\displaystyle D_{V}^{a}$
$\displaystyle=$ $\displaystyle\frac{Q_{V}^{a}}{n_{0}M^{2}}$ (54)
$\displaystyle D_{XV}^{a}$ $\displaystyle=$
$\displaystyle\frac{Q_{XV}^{a}}{n_{0}M^{2}}\;.$ (55)

The diffusion in the position in (52) arise from both the quantum noise and
the diffusion of the velocity. At the lowest order in $t$, we have

$\langle{\left[X^{a}(t)\right]}^{2}\rangle=D_{X}^{a}t+D_{XV}^{a}t^{2}+D_{V}^{a}\frac{t^{3}}{3}\simeq
D_{X}^{a}t\;.$ (56)

For the random walk of $\theta$ and $E$, we obtain after Eqs. (37)

$\begin{array}[]{lll}\langle E(t)^{2}\rangle\simeq D_{E}t\\\
\langle\theta(t)^{2}\rangle\simeq D_{\theta}t\\\ \end{array}$ (57)

with

$\displaystyle D_{E}$ $\displaystyle=$
$\displaystyle-\frac{4}{(M^{\prime})^{2}n_{0}}\int\rho({\bf x}_{1})\rho({\bf
x}_{2})U({\bf x}_{1}-{\bf x}_{2})\,\mathrm{d}^{3}{\bf
x}_{1}\,\mathrm{d}^{3}{\bf x}_{2}\;,$ (58) $\displaystyle D_{\theta}$
$\displaystyle=$
$\displaystyle-\frac{1}{(M^{\prime})^{2}n_{0}}\int\rho^{\prime}({\bf
x}_{1})\rho^{\prime}({\bf x}_{2})U({\bf x}_{1}-{\bf
x}_{2})\,\mathrm{d}^{3}{\bf x}_{1}\,\mathrm{d}^{3}{\bf x}_{2}\;,$ (59)

being $\rho^{\prime}({\bf x})=\partial\rho({\bf x})/\partial E$.

For the radially symmetric soliton, we show in figure 2b the diffusion
coefficient $D_{X}^{1}=D_{x}^{2}=D_{x}^{3}=D_{X}$, as obtained by the
numerical profile $u$ computed with the screened gravitational potential $U$.
One finds that for a growing $\sigma$ the quantum diffusion is frustrated, as
it happens in 1D [19, 20]. This can be deduced from (48), indeed, as
$\sigma\rightarrow\infty$, one has $U({\bf x}_{1}-{\bf x}_{2})\simeq$
constant, and $Q_{X}^{a}\rightarrow 0$, as for the soliton profile $\rho({\bf
x})=\rho(-{\bf x})$.

We compare (56) with the full 3D+1 stochastic simulations and we find
excellent agreement, as shown in figure 2c-e where we report the dynamics of
solitary waves with $n_{0}M\simeq 10^{6}$ atoms.

The diffusion constant $D=\hbar D_{X}/2m$ in physical units reads

$D=\left(\frac{Q_{X}^{2}}{2n_{0}M^{2}}\right)\frac{\hbar}{m}=\left(\frac{Q_{X}^{2}}{2M}\right)\frac{\hbar}{N_{T}m}\;.$
(60)

In Eq.(60) equation $N_{T}=n_{0}M$ is the total number of particles, and
$N_{T}m=n_{0}Mm$ the total mass of the condensate. $Q_{X}^{2}/2M$ is a
numerical constant that depends on the profile of the soliton and $\Lambda$.
We find that quantum fluctuations vanish when $N_{T}\rightarrow\infty$ or
$\hbar\rightarrow 0$. In the original units of Eq.(2), the diffusion constant
can be also cast as

$D=\frac{Gm^{2}}{\hbar}\int\int\frac{x_{1}\rho({\bf
x}_{1})}{\int\rho\,\mathrm{d}^{3}{\bf x}}\frac{x_{2}\rho({\bf
x}_{2})}{\int\rho\,\mathrm{d}^{3}{\bf x}}\frac{e^{-|{\bf x}_{1}-{\bf
x}_{2}|/\Lambda}}{|{\bf x}_{1}-{\bf x}_{2}|}\,\mathrm{d}^{3}{\bf
x}_{1}\,\mathrm{d}^{3}{\bf x}_{2}\;.$ (61)

We remark that Eq. (61) is written in the original physical units of Eq.(2),
such that in Eq. (61) $x$ is a length and the dimensions of $D$ are $m^{2}/s$
in the MKS system. Equation (61) shows the interplay of quantum and
gravitational effects through the ratio $Gm^{2}/\hbar$ and returns $D$ in
terms of the measurable density profile $\rho({\bf x})$.

##  6 Non-Gaussian statistics

In our stochastic simulations, the initial state is a coherent state, whose
statistical properties change upon evolution. Here we follow [23] to determine
if deviations from Gaussianity arise. We report in figure 3a the evolution of
the statistical distribution of the density $\rho({\bf x}=0)$ as computed by
Eqs. (8) at the center of the classical solitonic core. The initial state is
coherent, and the histogram is localized in the initial value of the peak.
Upon evolution, the distribution spreads and manifestly displays a bell-shaped
non-Gaussian profile. Similar behavior is also obtained for the quadratures of
the field (not reported).

To quantify the deviation from Gaussianity, we consider the SNR introduced in
[23]

$\rm{SNR}=\frac{|\kappa_{4}|}{\sqrt{\rm{var}\,k_{4}}}$ (62)

here $\kappa_{4}$ is the fourth cumulant of the statistical distribution.
$\rm{var}\,k_{4}$ is its uncertainty (see B). For Gaussian statistics, all the
cumulants higher than second order vanish, hence SNR measures deviation from
non-Gaussianity including the uncertainty $\rm{var}\,k_{4}$ due to a finite
number of samples. We compute SNR for the density and the field quadratures
with similar results.

At variance with [23], we account for the heterogeneous features of SNR, i.e.,
we measure SNR in different spatial locations. Figure 3b shows the 3D
isosurface of the SNR at different instants. The statistical distributions at
different positions become non-Gaussian with time. Figure 3c shows the
spatially averaged value of the SNR, which demonstrates that a self-trapped
solitonic wave packet develops non-Gaussian statistics. Results in figure 3
refer to a representative case with $n_{0}M\simeq 10^{4}$ atoms; we found
these dynamics for different interaction lengths and particle numbers.

To understand the physical origin of the non-Gaussianity, we observe that - at
the lowest order in $t$ \- the soliton parameters $X^{a}$, $V^{a}$, $\theta$,
and $E$, are the time-integral of white noise terms (i.e., Wiener processes).
Thus they are the sum of many independent variables and hence obey Gaussian
statistics. Non-Gaussianity arises from the fact that the soliton profile is a
nonlinear function of these parameters, and any observable depends on the
soliton profile. In general terms, the statistical distribution of a nonlinear
function of a Gaussian variable is expected to be non-Gaussian. Thus, as far
as the soliton is stable with respect to fluctuations, non-Gaussianity arises.
Nonlocal solitons are stable self-trapped nonlinear waves, and their
robustness against quantum fluctuations induces non-Gaussianity.

Figure 3: Time evolution and non-Gaussianity of the statistical distribution
of 3D nonlocal solitons. (a) Histogram (after $1000$ runs) at different
instants of the density $\rho$ computed at the classical soliton peak ${\bf
x}=0$; the vertical axis is truncated at $300$, the number of counts at $t=0$
and $\rho(0)=1.0$ is $1000$ as indicated by the arrow. (b) Temporal evolution
of the spatial distribution of the non-Gaussianity parameter SNR of the
density as 3D isosurfaces. (c) Mean value of the SNR computed on the spatial
profile versus time (parameters $\sigma=1,M=138$, $n_{0}=10^{2}$).

##  7 Summary

In conclusion, we studied theoretically and by first-principle numerical
simulations the 3D+1 dynamics of non-local self-gravitating boson fluids. The
quantum noise induces diffusion in the self-localized wave-packet position
determined by the degree of nonlocality and the particle number. The
theoretical results agree with ab-initio 3D+1 simulations with no fitting
parameters.

The quantum diffusion is due to the interplay of the quantum fluctuations and
the long-range self-interaction. This interplay causes non-Gaussian statistics
that spread in the solitonic core upon evolution. We remark that this is a
universal phenomenon that is not dependent on the specific interaction
potential $U$ but arises from the general stability properties of solitons.

Experimental investigations may involve long-range Bose-Einstein condensates
(see, e.g.,[8] and references therein), and also nonlinear optical systems,
where low-dimensional reductions of the Schrödinger-Newton equation have been
considered [30].

The results open the way to using non-Gaussian multidimensional solitary waves
as non-classical reservoirs for continuous-variable quantum information and as
quantum simulators for quantum gravity models. Notably enough, the numerical
simulations suggest that signatures of a quantized gravity may arise even
without careful preparation of the initial state as Schrödinger cat (or
squeezed state), but starting from a coherent solitonic state. Also, the
results show the relevance of quantum fluctuations in cold dark-matter models,
which can potentially impact the investigation of self-gravitating BEC and
enable tests within astrophysical observations.

We acknowledge support from the H2020 PhoQus project (Grant no.820392).

## Data availability statement

All data that support the findings of this study are included within the
article (and any supplementary files).

##  Appendix A Stochastic equations for the soliton parameter

The perturbation vector is written as

${\bf
s}=\frac{1}{\sqrt{n_{0}}}\left(\begin{array}[]{ccc}\sqrt{\imath}\xi_{U}\\\
\sqrt{-\imath}\xi_{U}^{+}\end{array}\right)u{\bf e}\;,$

where the two independent complex noises $\xi_{U}$ and $\xi_{U}^{+}$ are such
that

$\langle\xi_{U}({\bf x},t)\xi_{U}({\bf x}^{\prime},t^{\prime})\rangle=U({\bf
x}-{\bf x}^{\prime})\delta(t-t^{\prime})$ (63)

and

$\langle\xi^{+}_{U}({\bf x},t)\xi^{+}_{U}({\bf
x}^{\prime},t^{\prime})\rangle=U({\bf x}-{\bf
x}^{\prime})\delta(t-t^{\prime})\;,$ (64)

being $U({\bf x})=-\exp(-r/\sigma)/r<0$. We let

$\begin{array}[]{lll}\xi_{U}({\bf x},t)&=&\imath C({\bf x})\xi({\bf x},t)\\\
\xi^{+}_{U}({\bf x},t)&=&\imath C({\bf x})\xi^{+}({\bf x},t)\\\ \end{array}$
(65)

which satisfy (63) and (64), with $C({\bf x})$ a real-valued function such
that

$C({\bf x})\ast C(-{\bf x})=\int\,\mathrm{d}^{3}{\bf x}^{\prime}C({\bf x}-{\bf
x}^{\prime})C(-{\bf x}^{\prime})=-U({\bf x})>0\;,$ (66)

or, equivalently,

$\int C({\bf x}-{\bf x}^{\prime})C({\bf y}-{\bf
x}^{\prime})\,\mathrm{d}^{3}{\bf x}^{\prime}=-U({\bf x}-{\bf y})>0\;.$ (67)

Following Eq.(34), we need the scalar product $\left(\hat{{\bf
f}}_{X}^{a},{\bf s}\right)$ at $X^{a}=0$, that is

$\displaystyle\sqrt{n_{0}}\left(\hat{{\bf f}}_{X}^{a},{\bf s}\right)$
$\displaystyle=2\Re\int\,\mathrm{d}^{3}{\bf
x}\left[-\imath\sigma_{3}\left(\hat{\bf f}_{V}^{a}\right)^{*}\right]\cdot{\bf
s}$ (68) $\displaystyle=\Re\int\,\mathrm{d}^{3}{\bf x}\,x^{a}u^{2}({\bf
x})\left(-\imath\sqrt{\imath}\xi_{U}+\imath\sqrt{-\imath}\xi_{U}^{+}\right)$
$\displaystyle=\Re\int\,\mathrm{d}^{3}{\bf
x}\,x^{a}u^{2}C\ast(\sqrt{\imath}\xi-\sqrt{-\imath}\xi^{+})$
$\displaystyle=\int\,\mathrm{d}^{3}{\bf
x}\,x^{a}u^{2}C\ast\frac{\xi-\xi^{+}}{\sqrt{2}}$
$\displaystyle=\int\,\mathrm{d}^{3}{\bf x}\,x^{a}u^{2}C\ast\xi_{-}\;,$ (69)

with $\xi_{-}\equiv(\xi-\xi^{+})/\sqrt{2}$ a real noise such that
$\langle\xi_{-}({\bf x},t)\xi_{-}({\bf
x}^{\prime},t^{\prime})\rangle=\delta({\bf x}-{\bf
x}^{\prime})\delta(t-t^{\prime})$. Seemingly, we have in Eq.(35)

$\sqrt{n_{0}}\left(\hat{{\bf f}}_{V}^{a},{\bf
s}\right)=\int\,\mathrm{d}V\left(-\frac{\partial u^{2}}{\partial
x^{a}}\right)C\ast\xi_{-}\;.$ (70)

To solve the resulting Ito stochastic equations we define
$\xi_{+}=C\ast\xi_{-}$ [see (43)], and

$F_{X}^{a}(t)=\int\,\mathrm{d}^{3}{\bf x}\,\left(\hat{{\bf f}}_{X}^{a},{\bf
s}\right)=\frac{1}{\sqrt{n_{0}}}\int\,\mathrm{d}^{{3}}{\bf
x}x^{a}u^{2}C\ast\xi_{-}=\frac{1}{\sqrt{n_{0}}}\int\,\mathrm{d}^{3}{\bf
x}\,x^{a}u^{2}\xi_{+}\;,$ (71)

and

$F_{V}^{a}(t)=-\int\,\mathrm{d}^{3}{\bf x}\,\left(\hat{{\bf f}}_{V}^{a},{\bf
s}\right)=\frac{1}{\sqrt{n_{0}}}\int\,\mathrm{d}^{3}{\bf x}\,\frac{\partial
u^{2}}{\partial
x^{a}}\,C\ast\xi_{-}=\frac{1}{\sqrt{n_{0}}}\int\,\mathrm{d}^{3}{\bf
x}\frac{\partial u^{2}}{\partial x^{a}}\xi_{+}\;.$ (72)

Eqs. (34) and (35) read

$\eqalign{M\dot{X}^{a}=MV^{a}+F_{X}^{a}(t)\\\ M\dot{V}^{a}=F_{V}^{a}(t)\;.}$
(73)

Eqs.(73) are solved by quadratures as follows

$\displaystyle MV^{a}(t)=\int_{0}^{t}F_{V}^{a}(s)\,\mathrm{d}s$ (74)
$\displaystyle
MX^{a}(t)=M\int_{0}^{t}V^{a}(s)\,\mathrm{d}s+\int_{0}^{t}F_{X}^{a}(s)\,\mathrm{d}s=\int_{0}^{t}\int_{0}^{s}F_{V}(u)\,\mathrm{d}u\,\mathrm{d}s+\int_{0}^{t}F_{X}^{a}(s)\,\mathrm{d}s\;.$
(75)

From (74) we have

$M^{2}\langle
V^{a}(t)V^{b}(t^{\prime})\rangle=\int_{0}^{t}\int_{0}^{t^{\prime}}F_{V}^{a}(s)F_{V}^{b}(s^{\prime})\,\mathrm{d}s\,\mathrm{d}s^{\prime}\;.$
(76)

From (75)

$\displaystyle M^{2}\langle X^{a}(t)X^{b}(t^{\prime})\rangle$
$\displaystyle=\int_{0}^{t}\int_{0}^{s}\int_{0}^{t^{\prime}}\int_{0}^{s^{\prime}}\langle
F_{V}^{a}(u)F_{V}^{b}(u^{\prime})\rangle\,\mathrm{d}u\,\mathrm{d}u^{\prime}\,\mathrm{d}s\,\mathrm{d}s^{\prime}+$
(77) $\displaystyle+\int_{0}^{t}\int_{0}^{t^{\prime}}\int_{0}^{s}\langle
F_{V}^{a}(u)F_{X}^{b}(s^{\prime})\rangle\,\mathrm{d}u\,\mathrm{d}s^{\prime}\,\mathrm{d}s+$
$\displaystyle+\int_{0}^{t}\int_{0}^{t^{\prime}}\int_{0}^{s^{\prime}}\langle
F_{V}^{a}(u^{\prime})F_{X}^{b}(s)\rangle\,\mathrm{d}u^{\prime}\,\mathrm{d}s^{\prime}\,\mathrm{d}s+$
$\displaystyle+\int_{0}^{t}\int_{0}^{t^{\prime}}\langle
F_{X}^{a}(s)F_{X}^{b}(s^{\prime})\rangle\,\mathrm{d}s\,\mathrm{d}s^{\prime}$

We also have from Eqs. (71) and (72) the following

$\displaystyle\langle
F_{V}^{a}(t)F_{V}^{b}(t^{\prime})\rangle=Q_{V}\delta_{ab}\delta(t-t)$ (78)
$\displaystyle\langle
F_{X}^{a}(t)F_{V}^{b}(t^{\prime})\rangle=Q_{XV}\delta_{ab}\delta(t-t)$ (79)
$\displaystyle\langle
F_{X}^{a}(t)F_{X}^{b}(t^{\prime})\rangle=Q_{X}\delta_{ab}\delta(t-t)$ (80)
(81)

where we accounted for the fact that $u^{2}({\bf x})=u^{2}(-{\bf x})$, and
(48),(50) and (50) hold. By using (81) in (76) and (77) and letting
$t=t^{\prime}$ we have Eq. (51) and Eq. (56). Similar arguments lead to (57).

##  Appendix B Non-Gaussianity parameter

The fourth order cumulant $\kappa_{4}$ is computed by using the value of the
density $\rho({\bf x},t)$, or of the field quadratures. Denoting as $q$ a
value of a single run, we first determine the non-central moments
($m=0,1,2,\ldots$)

$\mu^{\prime}_{m}=\langle q^{m}\rangle\,.$ (82)

Then we compute the first $8$ cumulants $\kappa_{n}$ with
$\kappa_{1}=\mu^{\prime}_{1}$, and ($n>1$)

$\kappa_{n}=\mu_{n}^{\prime}-\sum_{m=1}^{n-1}\kappa_{n-m}\mu_{m}^{\prime}\;.$
(83)

For the $k-$statistics, we have $\langle k_{4}\rangle=\kappa_{4}$, and,
letting $\mathcal{M}$ the number of runs,

$\begin{array}[]{lll}\rm{var}(k_{4})&=\langle(k_{4}-\kappa_{4})^{2}\rangle=+\frac{\kappa_{8}}{\mathcal{M}}+16\frac{\kappa_{2}\kappa_{6}}{\mathcal{M}-1}+48\frac{\kappa_{3}\kappa_{5}}{\mathcal{M}-1}+34\frac{\kappa_{4}^{2}}{\mathcal{M}-1}+\\\
&72\frac{\mathcal{M}\kappa_{2}^{2}\kappa_{4}}{(\mathcal{M}-1)(\mathcal{M}-2)}+144\frac{\mathcal{M}\kappa_{2}\kappa_{3}^{2}}{(\mathcal{M}-1)(\mathcal{M}-2)}+24\frac{\mathcal{M}(\mathcal{M}+1)\kappa_{2}^{2}}{(\mathcal{M}-1)(\mathcal{M}-2)(\mathcal{M}-3)}\;.\end{array}$
(84)

## References

## References

  * [1] Kivshar Y and Agrawal G P 2003 Optical solitons (New York: Academic Press)
  * [2] Malomed B 2021 arXiv:2111.00547
  * [3] Turitsyn S K 1985 Teor. Mat. Fiz. 64 226
  * [4] Pérez-García V M, Konotop V V and García-Ripoll J J 2000 Phys. Rev. E 62 4300–4308
  * [5] Bang O, Krolikowski W, Wyller J and Rasmussen J J 2002 Phys Rev E 66 046619
  * [6] Klaers J, Schmitt J, Vewinger F and Weitz M 2010 Nature 468 545–548
  * [7] Carusotto I and Ciuti C 2013 Rev. Mod. Phys. 85 299–366
  * [8] Calvanese Strinati M and Conti C 2014 Phys. Rev. A 90 043853
  * [9] Defenu N, Donner T, Macrì T, Pagano G, Ruffo S and Trombettoni A 2021 arXiv:2109.01063
  * [10] O’Dell D, Giovanazzi S, Kurizki G and Akulin V M 2000 Phys. Rev. Lett. 84 5687–5690
  * [11] Paredes A and Michinel H 2016 Phys. Dark Universe 12 50 – 55
  * [12] Garnier J, Baudin K, Fusaro A and Picozzi A 2021 arXiv:2108.13250
  * [13] Conti C 2014 Phys. Rev. A 89 061801
  * [14] Liang Q Y, Venkatramani A V, Cantu S H, Nicholson T L, Gullans M J, Gorshkov A V, Thompson J D, Chin C, Lukin M D and Vuletić V 2018 Science 359 783–786
  * [15] Villari L D M, Faccio D, Biancalana F and Conti C 2018 Phys. Rev. A 98 043859
  * [16] Marchukov O V, Malomed B A, Dunjko V, Ruhl J, Olshanii M, Hulet R G and Yurovsky V A 2020 Phys Rev Lett 125 050405
  * [17] Conti C 2022 Phys. Rev. A 106 013518
  * [18] Alodjants A, Tsarev D, Ngo T V and Lee R K 2022 Physical Review A 105 012606
  * [19] Folli V and Conti C 2010 Phys. Rev. Lett. 104 193901
  * [20] Batz S and Peschel U 2011 Phys. Rev. A 83 033826
  * [21] D Drummond P and Chaturvedi S 2016 Phys Scripta 91 073007
  * [22] Wüster S, Corney J F, Rost J M and Deuar P 2017 Phys Rev E 96 013309
  * [23] Howl R, Vedral V, Naik D, Christodoulou M, Rovelli C and Iyer A 2021 PRX Quantum 2 010325
  * [24] Hughes C, Genoni M G, Tufarelli T, Paris M G A and Kim M S 2014 Phys. Rev. A 90 013810
  * [25] Zhuang Q, Shor P W and Shapiro J H 2018 Phys. Rev. A 97 052317
  * [26] Drummond P D, Deuar P and Kheruntsyan K V 2004 Phys. Rev. Lett. 92 040405
  * [27] Gardiner C W and Zoller P 2004 Quantum Noise 3rd ed (Berlin: Springer-Verlag)
  * [28] Ruffini R and Bonazzola S 1969 Phys. Rev. 187 1767–1783
  * [29] Frigo M and Johnson S G 2005 Proceedings of the IEEE 93 216–231 special issue on “Program Generation, Optimization, and Platform Adaptation”
  * [30] Rogert T, Maitland C, Wilson K, Westerberg N, Vocke D, Wright E W and Faccio D 2016 Nat. Commun. 7 13492