# Dark matter condensates as highly nonlocal solitons: instability in the
Schwarzschild metric and laboratory analog

Ludovica Dieli<EMAIL_ADDRESS>Claudio Conti
<EMAIL_ADDRESS>Dipartimento di Fisica, Università di Roma ”La
Sapienza”, 00185 Rome, Italy

###### Abstract

Theories on the bosonic nature of dark matter are a promising alternative to
the cold dark matter model. Here we consider a dark matter halo in the state
of a Bose-Einstein condensate, subject to the gravitation of a black hole. In
the low energy limit, we bring together the general relativity in the
Schwarzschild metric and the quantum description of the Bose-Einstein
condensate. The model is solvable in the Fermi normal coordinates with the so
called highly nonlocal approximation and describes tidal deformations in the
condensate wave function. The black hole deforms the localized condensate
until the attraction of the compact object overcomes the self-gravitation and
destabilizes the solitonic dark matter. Moreover, the model can be implemented
as a gravitational analog in the laboratory; the time-dependent potential
generated by the galactic black hole can be mimicked by an optical trap acting
on a conventional condensate. The results open the way to new laboratory
simulators for quantum gravitational effects.

††preprint: APS/123-QED

##  I Introduction

One of the main unsolved questions about the Universe is the nature of dark
matter (DM). Many astrophysical observations have revealed the presence of a
new kind of matter, which gravitationally interacts with the visible matter
but does not emit light [1]. Among the proposed theories, one of the most
promising is the cold dark matter model (CDM), based on collisionless
particles. Despite CDM explains anomalous galaxy rotation curves, new
discrepancies, such as the cuspy problem regarding the flat density
distribution in the galactic core, remain not yet understood [1]. Authors
proposed the introduction of particle self-interaction to solve the
inconsistencies in the CDM [2]. Inspired by early Universe theories, bosonic
DM with ultrasmall masses $\sim 10^{-24}-10^{-20}\,\text{eV}$ [3] allows the
formation of galactic Bose-Einstein condensates (BECs) [4, 5, 6, 7, 9]. One
key aspect of BEC DM halos is the role of spacetime curvature first considered
in [4]. In [10], the curved spacetime appears as a perturbative potential
inducing an instability dependent on the condensate mass. However, the effect
of an external gravitational field on the dynamics of the galactic condensate
is unexplored.

In this manuscript, we study a DM condensate subject to a rotationally
symmetric gravitational field described by the Schwarzschild metric. The Bose-
condensed DM experiences the gravitational field of a black hole (BH) in the
form of an effective potential in the non-relativistic limit. By a simple
closed-form theory, we show how the interaction induces the excitation of
high-energy states. These states are required to fit the galaxy rotation
curves in [3]. Also, we predict the existence of a breaking time for the
destabilization of the BEC, which may explain the absence of DM in specific
galaxies.
Moreover, we consider the feasibility of BEC DM simulations in a laboratory
[11]. We detail an experimental setup to address the gravitational instability
in an optically trapped Bose gas. Indeed, in the context of analog gravity
[12], physical systems such as in polaritonics, atomic gases, and nonlinear
optics are used to mimic the cosmological models.
The manuscript is organized as follows. In Sec. II, we introduce the model:
$(i)$ the Fermi normal coordinates used to make the expansion of the
d’Alembert operator and $(ii)$ the multiple scale expansion of the solution of
the Klein-Gordon equation in the non-relativistic limit.
Section III is about the effective potential in the Schwarzschild metric,
simplified by the highly-nonlocal approximation in Sec. IV.
In Sec. V, we analyze the asymptotic model and predict the existence of a
critical time for the instability dependent on the masses of the DM and of the
BH.
In Sec. VI, we specifically address the significant timescales and other
parameters in the interaction with BHs; in Sec. VII, we described an
experiment for the laboratory analog. Conclusions are drawn in Sec. VIII.

##  II Model and methods

The hypothesis of BEC dark matter halos opens the question about the effects
of the spacetime curvature. We analyze the scenario in which a BEC DM is
freely falling in a curved metric. We take into account two aspects: $(i)$ the
BEC with its quantum description and $(ii)$ the relativistic spacetime.
Concerning point $(i)$, we describe the evolution of the self-gravitating BEC
by the Gross-Pitaevskii equation (GPE) [13] with a gravitational potential

$\displaystyle
i\hbar\partial_{t}\psi(t,\vec{r})=-\frac{\hbar^{2}}{2m}\nabla^{2}\psi(t,\vec{r})+m\Phi(t,\vec{r})\psi(t,\vec{r});$
(1) $\displaystyle\nabla^{2}\Phi(t,\vec{r})-K^{2}\Phi(t,\vec{r})=4\pi
Gm|\psi(t,\vec{r})|^{2};$ (2)

where $K$ is inversely proportional to the interaction length of the potential
$\Phi$ and $G$ is the universal gravitational constant. $\Phi(t,\vec{r})$ is a
nonlinear long-range gravitational potential solution of the Poisson equation
with the mass density distribution $m|\psi(t,\vec{r})|^{2}$. The BEC
description is clearly non-relativistic since the Poisson equation concerns
Newtonian gravitational field. However, eq.(2) is naturally generalised to a
relativistic scenario in the form of Einstein’s equations [14].
Concerning point $(ii)$, the BEC experiences an external gravitational field
described by a metric $g_{\mu\nu}$. However, the GPE is valid in a non-
relativistic regime, while the curved metric requires a four-dimensional
formalism.
To unify the two descriptions in a non-relativistic regime, we expand the
d’Alembert operator by means of the Fermi normal coordinates suitable to
describe a freely falling object. The expansion is applied to the Klein-Gordon
equation recovering a Schroedinger-like equation in the non-relativistic limit
[15]. We also use a multiple scale expansion and treat the nonlinearity as a
mass-dependent potential within the highly nonlocal approximation.

###  II.1 Curved metric in Fermi normal coordinates

According to the Equivalence Principle, a freely falling object is locally
described by the flat metric. The Fermi normal coordinates enlarge this
description to a world tube around the timelike trajectory [16]. The Fermi
normal coordinates satisfy the condition $\Gamma^{\alpha}_{\mu\nu}=0$ along
the object trajectory. The basis vectors are defined in the following way: the
time coordinate $\tau$ is such that $\bm{e_{0}}=\frac{\partial}{\partial\tau}$
is tangent to the trajectory. From every point of the trajectory, an
orthogonal hypersurface $\Sigma$ is defined, and $3$ independent vectors are
chosen on $\Sigma$ as the spatial vectors of the basis
$\bm{e_{j}}=\frac{\partial}{\partial x_{j}}$ with $j=1,2,3$. A $(2+1)D$
simplified representation is reported in Fig. 1. $(c\tau,x_{1},x_{2},x_{3})$
are the Fermi normal coordinates. They describe the spacetime around the
falling object as an expansion of the flat metric

$\begin{split}g_{\mu\nu}(c\tau,\vec{x})&=\eta_{\mu\nu}+g_{\mu\nu,\alpha}(c\tau,0)x^{\alpha}\\\
&+\frac{1}{2}g_{\mu\nu,\alpha\beta}(c\tau,0)x^{\alpha}x^{\beta}+o(x^{2})\end{split}$

where the comma stands for the partial derivative of the metric tensor with
respect to the Fermi normal coordinates.
The linear term vanishes due to $\Gamma^{\alpha}_{\mu\nu}=0$ on the timelike
geodesic, and the first corrections are quadratic. According to [16], we found
the following expressions for the metric elements

$\begin{split}&g_{00}=-1-R_{0l0m}x^{l}x^{m}+o(x^{2});\\\
&g_{0k}=-\frac{2}{3}R_{0jki}x^{i}x^{j}+o(x^{2});\\\
&g_{lm}=\delta_{lm}-\frac{1}{3}R_{iljm}x^{i}x^{j}+o(x^{2});\end{split}$ (3)

where $R_{\alpha\beta\mu\nu}$ is the Riemann tensor. The quadratic terms are
small in the neighborhood of the geodesic, defining $\varepsilon\ll 1$:
$g_{\mu\nu}=\eta_{\mu\nu}+\varepsilon^{2}g^{(1)}_{\mu\nu}$. Thus, given the
spacetime metric $g_{\mu\nu}$, its expansion is defined by the Riemann tensor
in Fermi normal coordinates.

Figure 1: The timelike geodesic $G$ (black line) is surrounded by a world tube
in which Fermi normal coordinates are defined, and special relativity laws are
valid. The green line represents a spacelike geodesic $S$ on an orthogonal
hypersurface $\Sigma$. The time component is tangent to the geodesic $G$:
$\bm{e}_{0}=\bm{u}$, while the space components lie on the surface $\Sigma$:
$\bm{t}=t^{j}\bm{e}_{j}$, $j=1,2,3$. In the world tube, the metric is flat up
to quadratic corrections.

We use this metric expansion to calculate the d’Alembert operator in the
curved spacetime as a series expansion defined in the neighbor of the
geodesic.
The d’Alembertian in general coordinates is

$\square=\frac{1}{\sqrt{g}}\partial_{\mu}(\sqrt{g}g^{\mu\nu}\partial_{\nu})$

where $g=-\text{det}(g_{\mu\nu})$, after (3)

$\begin{split}g&=-\det(g_{\mu\nu})=1+\varepsilon^{2}\frac{2}{3}R_{0l0m}\,x^{l}x^{m}+o(\varepsilon^{2})\\\
&=1+\varepsilon^{2}g^{(1)}+o(\varepsilon^{2}).\end{split}$ (4)

The inverse matrix for the metric can be found from
$g^{\mu\nu}g_{\nu\xi}=\delta^{\mu}_{\xi}$

$\begin{split}&g^{\mu\nu}=\eta^{\mu\nu}+\varepsilon^{2}g^{(1)\mu\nu};\\\
&g_{\nu\xi}=\eta_{\nu\xi}+\varepsilon^{2}g^{(1)}_{\nu\xi}.\end{split}$

We have

$\begin{split}&g^{(1)00}=-g^{(1)}_{00}=R_{0l0m}\,x^{l}x^{m};\\\
&g^{(1)0i}=-g^{(1)}_{0i}=\frac{2}{3}R_{0lim}\,x^{l}x^{m};\\\
&g^{(1)ij}=-g^{(1)}_{ij}=\frac{1}{3}R_{iljm}\,x^{l}x^{m}.\end{split}$ (5)

and

$\begin{split}\square&=\left(1-\frac{\varepsilon^{2}}{2}g^{(1)}\right)\,\partial_{\mu}\left[\left(1+\frac{\varepsilon^{2}}{2}g^{(1)}\right)\left(\eta^{\mu\nu}+\varepsilon^{2}g^{(1)\mu\nu}\right)\partial_{\nu}\right]=\\\
&=-\partial^{2}_{0}+\nabla^{2}+\varepsilon^{2}\partial_{\mu}(g^{(1)\mu\nu}\partial_{\nu})+\frac{\varepsilon^{2}}{2}\eta^{\mu\nu}\partial_{\mu}(g^{(1)})\partial_{\nu}+o(\varepsilon^{2}).\end{split}$
(6)

This expansion in the Fermi normal coordinates needs a specific choice of the
metric describing the falling object.

###  II.2 Multiple scale expansion

We start with the Klein-Gordon (KG) equation

$(\square-m^{2})\phi=0,$ (7)

which reduces to the Schroedinger equation in the non-relativistic limit, and
use Eq. (6)

$\begin{split}\square&=\square_{\text{flat}}+\varepsilon^{2}\partial_{\mu}(g^{(1)\mu\nu}\partial_{\nu})\\\
&+\frac{\varepsilon^{2}}{2}\eta^{\mu\nu}(\partial_{\mu}g^{(1)})\partial_{\nu}+o(\varepsilon^{2}).\end{split}$

Where we adopted units with $c=1$ and $\hbar=1$.
Since the effect of the curvature appears as a second-order perturbation of
the flat d’Alembertian, we add the nonlinear potential of the DM condensate as
a term of the same order

$(\square-m^{2})\phi+\varepsilon^{2}\Phi[\phi]\phi=0.$ (8)

We define a multiple scale expansion by the slow variables
$x_{1}^{\mu}=\varepsilon x^{\mu}$ and $x_{2}^{\mu}=\varepsilon^{2}x^{\mu}$

$\begin{split}\phi(x^{\mu},x_{1}^{\mu},x_{2}^{\mu})&=\phi_{0}(x^{\mu},x_{1}^{\mu},x_{2}^{\mu})+\varepsilon\phi_{1}(x^{\mu},x_{1}^{\mu},x_{2}^{\mu})\\\
&+\varepsilon^{2}\phi_{2}(x^{\mu},x_{1}^{\mu},x_{2}^{\mu})+o(\varepsilon^{2}).\end{split}$
(9)

Due to the presence of the slow variables, the flat d’Alembertian is

$\begin{split}\square_{\text{flat}}&=\square_{0}+2\varepsilon\square_{01}+\varepsilon^{2}\square_{1}+2\varepsilon^{2}\square_{02}+o(\varepsilon^{2})\\\
&=\eta^{\mu\nu}\frac{\partial^{2}}{\partial x^{\mu}\partial
x^{\nu}}+2\varepsilon\eta^{\mu\nu}\frac{\partial^{2}}{\partial x^{\mu}\partial
x_{1}^{\nu}}\\\ &+\varepsilon^{2}\eta^{\mu\nu}\frac{\partial^{2}}{\partial
x_{1}^{\mu}\partial
x_{1}^{\nu}}+2\varepsilon^{2}\eta^{\mu\nu}\frac{\partial^{2}}{\partial
x^{\mu}\partial x_{2}^{\nu}}+o(\varepsilon^{2}).\end{split}$

For the self-gravitating potential we have

$\begin{split}\Phi[\phi]&=\int
U(\vec{r}-\vec{r}\,^{\prime})|\phi(\vec{r}\,^{\prime})|^{2}d^{3}r\,^{\prime}\\\
&=\int
U(\vec{r}-\vec{r}\,^{\prime})[|\phi_{0}(\vec{r}\,^{\prime})|^{2}+\varepsilon
2\phi_{0}\phi_{1}\\\
&+\varepsilon^{2}|\phi_{1}(\vec{r}\,^{\prime})|^{2}+\varepsilon^{2}2\phi_{0}\phi_{2}\\\
&+\varepsilon^{3}2\phi_{1}\phi_{2}+\varepsilon^{4}|\phi_{2}(\vec{r}\,^{\prime})|^{2}]d^{3}r\,^{\prime}.\end{split}$
(10)

Equation (8) is

$\begin{split}&\left[\square_{0}+2\varepsilon\square_{01}+\varepsilon^{2}\square_{1}+2\varepsilon^{2}\square_{02}+\varepsilon^{2}\partial_{\mu}(g^{(1)\mu\nu}\partial_{\nu})+\varepsilon^{2}\frac{\eta^{\mu\nu}}{2}(\partial_{\mu}g^{(1)})\partial_{\nu}\right](\phi_{0}+\varepsilon\phi_{1}+\varepsilon^{2}\phi_{2})\\\
&-m^{2}(\phi_{0}+\varepsilon\phi_{1}+\varepsilon^{2}\phi_{2})+\varepsilon^{2}\Phi[|\phi_{0}|^{2}](\phi_{0}+\varepsilon\phi_{1}+\varepsilon^{2}\phi_{2})=0.\end{split}$
(11)

At order $O(1)$ we have

$(\square_{0}-m^{2})\phi_{0}=0$ (12)

which is the KG equation in the flat spacetime. The solution is a free
particle wave function modulated by an amplitude depending on the slow
variables

$\begin{split}\phi_{0}(x^{\mu},x_{1}^{\mu},x_{2}^{\mu})&=e^{ik_{\mu}x^{\mu}}\psi(x_{1}^{\mu},x_{2}^{\mu})\\\
&+e^{-ik_{\mu}x^{\mu}}\psi^{*}(x_{1}^{\mu},x_{2}^{\mu}).\end{split}$

At the order $O(\varepsilon)$ in Eq. (11) we have

$(\square_{0}-m^{2})\phi_{1}+2\square_{01}\phi_{0}=0$ (13)

that is

$(\square_{0}-m^{2})\phi_{1}=-2\square_{01}\phi_{0}\\\
=-2\eta^{\mu\nu}\frac{\partial^{2}\phi_{0}}{\partial x^{\mu}\partial
x_{1}^{\nu}}\\\ =-2i\omega
e^{-i\omega\tau+i\vec{k}\cdot\vec{x}}\left(-\frac{\partial\psi}{\partial\tau_{1}}+\frac{k^{j}}{\omega}\frac{\partial\psi}{\partial
x_{1}^{j}}\right)+c.c.$ (14)

We notice that the forcing term is resonant with the solution of the
homogeneous equation: it has to be null to avoid secular terms

$\left(-\frac{\partial\psi}{\partial\tau_{1}}+\frac{k^{j}}{\omega}\frac{\partial\psi}{\partial
x_{1}^{j}}\right)=0.$ (15)

As a consequence, Eq. (14) reduces to the KG equation, and we choose
$\phi_{1}=0$.
After Eq. (15) we have

$\psi(x_{1}^{\mu},x_{2}^{\mu})=\psi(\xi_{1}^{j},x_{2}^{\mu});\\\
\xi_{1}^{j}=k^{j}\left(x_{1}^{j}+\frac{k^{j}}{\omega}\tau_{1}\right)\quad
j=1,2,3;$ (16)

with $k^{j}$ arbitrary functions. The upper index $j$ identifies the component
of a $3$-dimensional vector.
At the order $O(\varepsilon^{2})$ in Eq. (11) we find

$\begin{split}&\left[\square_{1}+2\square_{02}+\partial_{\mu}(g^{(1)\mu\nu}\partial_{\nu})+\frac{\eta^{\mu\nu}}{2}(\partial_{\mu}g^{(1)})\partial_{\nu}\right]\phi_{0}\\\
&+m\Phi[\phi_{0}]\phi_{0}+(\square_{0}-m^{2})\phi_{2}=0.\end{split}$ (17)

The gravitational potential is a function of $\phi_{0}$, thus

$\begin{split}\Phi[\phi_{0}]&=\int
G(\vec{r}-\vec{r}\,^{\prime})\left[2|\psi(\vec{r}\,^{\prime})|^{2}+\right.\\\
&\left.+2\real\left(e^{2ik_{\mu}x^{\mu}}\psi^{2}(\vec{r}\,^{\prime})\right)\right]d^{3}r^{\prime}\end{split}$
(18)

this means that, solving the equation (17) with respect to $\phi_{2}$, it
needs the resonant term to vanish, obtaining the following

$\begin{split}&\left[\square_{1}+2\square_{02}+\partial_{\mu}(g^{(1)\mu\nu}\partial_{\nu})+\frac{\eta^{\mu\nu}}{2}(\partial_{\mu}g^{(1)})\partial_{\nu}\right.\\\
&+2m\left.\int
d^{3}r^{\prime}\,G(\vec{r}-\vec{r}\,^{\prime})|\psi(\vec{r}\,^{\prime})|^{2}\right]\phi_{0}(x^{\mu},x_{1}^{\mu},x_{2}^{\mu});\end{split}$
(19)

and

$(\square_{0}-m^{2})\phi_{2}\\\ +2m\phi_{0}\int
d^{3}r^{\prime}G(\vec{r}-\vec{r}\,^{\prime})\,\real\left(e^{2ik_{\mu}x^{\mu}}\psi^{2}\right)=0.$
(20)

We consider the non-relativistic limit of the leading Eq. (19) by assuming a
low energy regime with $\omega^{2}=m^{2}+k^{2}\simeq m^{2}$

$\begin{split}&\phi_{0}(x^{\mu},x_{1}^{\mu},x_{2}^{\mu})=e^{-im\tau}\psi(\xi_{1}^{j},x_{2}^{\mu})+c.c.;\\\
&|\vec{k}\cdot\vec{x}|\ll 1,\quad\omega^{2}\sim m^{2}.\end{split}$ (21)

As a consequence, the third and fourth terms of equation (19) contain only
time derivatives

$\begin{split}&\partial_{\mu}(g^{(1)\mu\nu}\partial_{\nu}\phi_{0})=(\partial_{0}g^{(1)00})\partial_{0}\phi_{0}+g^{(1)00}\partial_{0}^{2}\phi_{0}\\\
&=[-im(\partial_{0}g^{(1)00})-m^{2}g^{(1)00}]e^{-im\tau}\psi;\end{split}$

and

$\begin{split}&\frac{\eta^{\mu\nu}}{2}(\partial_{\mu}g^{(1)})\partial_{\nu}\phi_{0}=-\frac{1}{2}(\partial_{0}g^{(1)})\partial_{0}\phi_{0}\\\
&+\frac{1}{2}\delta^{kj}(\partial_{k}g^{(1)})\partial_{j}\phi_{0}\\\
&=\frac{im}{2}(\partial_{0}g^{(1)})e^{-im\tau}\psi.\end{split}$

From Eq. (5), $g^{(1)}$ and $g^{(1)00}$ depend on time $\tau$ only through the
Riemann tensor, which is a function of the coordinate system centered on the
BH. Noting that the time variation of the distance between the BEC and the BH
is much smaller than the distance itself, namely $\partial_{0}g^{(1)00}\ll
g^{(1)00}$ and $\partial_{0}g^{(1)}\ll g^{(1)00}$ we neglect the corresponding
terms (as detailed for Schwarzschild metric in Sec. III).
The first two terms of the equation (19) are

$\begin{split}&(\square_{1}+2\square_{02})\phi_{0}\\\
&=\left(\eta^{\mu\nu}\frac{\partial^{2}}{\partial x_{1}^{\mu}\partial
x_{1}^{\nu}}+2\eta^{\mu\nu}\frac{\partial^{2}}{\partial x^{\mu}\partial
x_{2}^{\nu}}\right)e^{-im\tau}\psi\\\
&=-e^{-im\tau}\frac{\partial^{2}\psi}{\partial\tau_{1}^{2}}+e^{-im\tau}\nabla_{1}^{2}\psi+2ime^{-im\tau}\frac{\partial\psi}{\partial\tau_{2}}.\end{split}$

Due to the constraint imposed on the slow time derivative, Eq. (15), we have

$\begin{split}\frac{\partial^{2}\psi}{\partial\tau_{1}^{2}}&=\frac{\partial}{\partial\tau_{1}}\left(\frac{k^{j}}{\omega}\frac{\partial\psi}{\partial
x_{1}^{j}}\right)\\\
&=\frac{\partial}{\partial\tau_{1}}\frac{k^{2}}{\omega}\frac{\partial\psi}{\partial\xi_{1}}\end{split}$

and we can neglect this term due to the approximation $\omega^{2}\simeq m^{2}$
and $|k|^{2}\sim 0$.
Concluding the expansion by posing $\varepsilon=1$ and reinserting all the
natural constants, the final equation results to be

$i\hbar\frac{\partial\psi}{\partial\tau}=-\frac{\hbar^{2}}{2m}\nabla^{2}\psi+\frac{mc^{2}}{2}g^{(1)00}\psi+m\Phi[|\psi|^{2}]\psi.$
(22)

We highlight that, so far, the only assumption on the nonlinear potential is
its quadratic form in the field, thus the result is completely general with
respect to any kind of quadratic nonlinearity.

##  III Schwarzschild metric

In the Schwarzschild metric generated by a non-rotating massive object like a
static BH

$ds^{2}=-X\,dT^{2}+X^{-1}\,dR^{2}+R^{2}\,d\Theta^{2}\\\
+R^{2}\,\sin^{2}\Theta\,d\varphi^{2}$ (23)

where $T$, $R$, $\Theta$ and $\varphi$ are Schwarzschild coordinates and
$X=\left(1-\frac{2GM}{c^{2}R}\right)=\left(1-\frac{R_{s}}{R}\right)$, with $G$
universal gravitational constant, $c$ light speed, $M$ the BH mass and
$R_{s}=\frac{2GM}{c^{2}}$ Schwarzschild radius.
The scenario is sketched in Fig. 2: the self-gravitating BEC is radially
falling toward the BH. We are considering a small DM halo with respect to the
BH, thus the spacetime curvature will be the one generated by the BH. The
hypothesis of small DM halo is supported by recent studies about space
extension of DM halos [8].
Eq. (22) describes the BEC time evolution, where the effective potential is
evaluated using the Schwarzschild metric and its expansion in Fermi normal
coordinates. By the metric element expansion in (5) and the Riemann tensor in
the Schwarzschild metric, we find the effective potential

$\begin{split}V_{\text{eff}}&=\frac{mc^{2}}{2}g^{(1)00}=\frac{mc^{2}}{2}R_{0l0m}\,x^{l}x^{m}\\\
&=\frac{mc^{2}}{2}(R_{0101}\,x^{2}+R_{0202}\,y^{2}+R_{0303}\,z^{2})\\\
&=\frac{mc^{2}R_{s}}{4R^{3}(\tau)}(-2x^{2}+y^{2}+z^{2}).\end{split}$ (24)

We notice that $V_{\text{eff}}$ is time-dependent through $R(\tau)$, which is
the radial distance between the BEC and the BH, expressed as a function of the
BEC proper time $\tau$. $R(\tau)$ is the solution of the geodesic equation

$\begin{split}\dot{R}(\tau)=-c\sqrt{E^{2}-V(R)}\end{split}$ (25)

where $E$ is a constant and
$V(R)=\left(1+\frac{l^{2}}{R^{2}}\right)\left(\frac{R_{s}}{R(\tau)}\right)$.
We verify for this particular choice the validity of the approximation
$\dot{R}(\tau)\ll R(\tau)$.
Moreover the potential (24) is quadratic, being attractive in the directions
$y$ and $z$ and repulsive along $x$. This difference in sign expresses the
presence of tidal forces exerted by the BH on the DM halo.

Having found the explicit form of the effective potential generated by the
Schwarzschild metric, in order to make the problem analytically tractable, we
consider the self-gravitating potential.

Figure 2:  Schematic illustration of the studied phenomenon. A dark matter
halo in the form of a self-gravitating Bose-Einstein condensate (BEC) is
radially falling toward a Schwarzschild black hole. a) At the beginning of the
evolution, the BEC is localized due to self-gravitational interaction. b) At
the end of the evolution, the BEC is deformed and, at a critical proper time,
destabilizes.

##  IV The highly nonlocal approximation

The nonlocality of the self-gravitating potential $\Phi[|\psi|^{2}]$ is the
key property we exploit. We assume the self-gravitating BEC self-localized;
the nonlocal effect of the self-gravitating potential is much wider than the
BEC mass distribution (Fig. 3). We hence replace the nonlocal nonlinear
potential with a parabolic one. This approach is the so called highly nonlocal
approximation (HNA) [17]. The potential $\Phi[|\psi|^{2}]$ is given by the
convolution of the mass density with the Green function

$\displaystyle\Phi(\vec{r})=\int
G(\vec{r}-\vec{r}\,^{\prime})|\psi(t,\vec{r}\,^{\prime})|^{2}d^{3}r^{\prime};$
(26) $\displaystyle G(\vec{r})=-\frac{Gm}{|\vec{r}|}e^{-K|\vec{r}|}.$ (27)

The HNA reduces the nonlinear long-range potential to a parabolic one

$\Phi\simeq\Phi_{M}+\tilde{\Phi}r^{2}=\Phi_{M}+\frac{\Omega^{2}}{2}r^{2}.$
(28)

Following the variational method presented in [17], we find the form for the
frequency of the harmonic approximation

$\Omega^{2}=\frac{(4G)^{4}}{(3\pi^{1/2})^{4}\hbar^{6}}m^{10}N^{4}.$ (29)

We underline that the potential $\Phi$ depends on the BEC total number of
particles $N$.
Within an additive constant potential term, equation (22) becomes

$\begin{split}i\hbar\partial_{\tau}\psi(\tau,\vec{r})&=-\frac{\hbar^{2}}{2m}\nabla^{2}\psi(\tau,\vec{r})+\frac{m\Omega^{2}}{2}r^{2}\,\psi(\tau,\vec{r})\\\
&+\frac{mc^{2}R_{s}}{4R^{3}(\tau)}(-2x^{2}+y^{2}+z^{2})\psi(\tau,\vec{r}).\end{split}$
(30)

Equation (30) is characterized by two parabolic terms of potential energy
causing the deformation of the BEC DM halo.

Figure 3:  Comparison between the spatial density distribution of the BEC
$|\psi(r)|^{2}$ and the intensity of the self-gravitational potential
$|\Phi(r)|$. We highlight the nonlocality of the interaction, which spreads
beyond the BEC mass distribution.

##  V Time evolution of the DM condensate

A critical element of the effective potential is its time dependence. In our
simplified model, it is due to the variation of the radial distance $R(\tau)$,
given by the solution of the geodesic equation

$\begin{split}\dot{R}(\tau)&=-c\sqrt{E^{2}-V(R)}\\\
&=-c\sqrt{\gamma^{2}-1+\frac{R_{s}}{R(\tau)}},\end{split}$ (31)

where $E$ is a constant equal to $\gamma=1/\sqrt{1-(v/c)^{2}}$, and assuming
the BEC has constant initial velocity $v$, with $\dot{R}(0)=-\gamma v$ for
$R(0)\rightarrow\infty$. We also have
$V(R)=\left(1+\frac{l^{2}}{R^{2}}\right)\left(\frac{R_{s}}{R(\tau)}\right)=1-\frac{R_{s}}{R(\tau)}$
for a non-rotating geometry.
Moreover, if we take zero initial velocity $v=0$, the solution of the geodesic
equation is ($\gamma=1$)

$\begin{split}R(\tau)=\left(-\frac{3}{2}\sqrt{R_{s}}c\tau+R^{3/2}(0)\right)^{2/3}.\end{split}$
(32)

From Eq. (32), we derive $3$ different time scales: (I) the time when the BEC
reaches the position $R(\tau_{0})=0$

$\tau_{0}=\frac{2R_{0}}{3c}\sqrt{\frac{R_{0}}{R_{s}}};$ (33)

(II) the time when the BEC crosses the event horizon $R(\tau_{s})=R_{s}$, i.e.
the Schwarzschild time

$\tau_{s}=\tau_{0}\left[1-\left(\frac{R_{s}}{R_{0}}\right)^{3/2}\right];$ (34)

and (III) the time at which the effective potential dominates over the self-
gravitation, i.e. the breaking time

$\tau_{B}=\tau_{0}-\frac{2}{3\Omega}.$ (35)

When $\tau=\tau_{B}$ the effective potential along $x$ becomes equal in
strength to the self-gravitating one

$\begin{split}\frac{1}{2}m\Omega^{2}x^{2}-\frac{mR_{s}c^{2}}{2R^{3}(\tau_{B})}x^{2}=0.\end{split}$
(36)

As a consequence, the sign of the total potential along $x$ changes, and the
BEC experiences a repulsion in the $x$-direction.
To analytically study the evolution, we consider the dimensionless version of
our model equation (30)

$\begin{split}i\,\partial_{\zeta}u(\zeta;\xi,\eta,\theta)&=-\nabla^{2}u(\zeta;\xi,\eta,\theta)\\\
&+\frac{1}{2}[\nu^{2}-V_{\xi}(\zeta)\,\xi^{2}+V_{\eta}(\zeta)\,\eta^{2}\\\
&+V_{\theta}(\zeta)\,\theta^{2}]\,u(\zeta;\xi,\eta,\theta)\;,\end{split}$ (37)

where we introduce the dimensionless wave function $u=\psi\sqrt{x_{p}^{3}}$
and variables

$\begin{split}\zeta=\frac{\tau}{\tau_{p}};\,\,\,\,\xi=\frac{x}{x_{p}};\,\,\,\,\eta=\frac{y}{x_{p}};\,\,\,\,\theta=\frac{z}{x_{p}};\end{split}$
(38)

with $\nu^{2}=\xi^{2}+\eta^{2}+\theta^{2}$, defined by the constants

$x_{p}^{2}=\frac{\hbar}{\sqrt{2}m\Omega};\qquad\tau_{p}=\frac{\sqrt{2}}{\Omega}.$
(39)

The dimensionless effective potential terms are

$\begin{split}&V_{\xi}=-\frac{2}{9}\frac{1}{(\zeta-\zeta_{0})^{2}};\\\
&V_{\eta,\theta}=\frac{1}{9}\frac{1}{(\zeta-\zeta_{0})^{2}}.\end{split}$

As the BEC Hamiltonian is the sum of spatially independent terms, a generic
solution is written as the product of three wave functions

$\begin{split}u(\zeta,\xi,\eta,\theta)=u_{\xi}(\zeta,\xi)\,u_{\eta}(\zeta,\eta)\,u_{\theta}(\zeta,\theta).\end{split}$

The motion equation splits into three different equations

$\begin{split}&i\partial_{\zeta}u_{\xi}=-\partial^{2}_{\xi}u_{\xi}+\frac{1}{2}\xi^{2}\left[1-\frac{2}{9(\zeta-\zeta_{0})^{2}}\right]u_{\xi};\\\
&i\partial_{\zeta}u_{\eta}=-\partial^{2}_{\eta}u_{\eta}+\frac{1}{2}\eta^{2}\left[1+\frac{1}{9(\zeta-\zeta_{0})^{2}}\right]u_{\eta};\\\
&i\partial_{\zeta}u_{\theta}=-\partial^{2}_{\theta}u_{\theta}+\frac{1}{2}\theta^{2}\left[1+\frac{1}{9(\zeta-\zeta_{0})^{2}}\right]u_{\theta}.\end{split}$
(40)

Eqs. (40) have a Gaussian solution

$\begin{split}&u_{j}(\zeta,j)=\frac{1}{\pi^{1/4}\sqrt{q_{j}(\zeta)}}\,\exp\left[\frac{i}{4}\frac{\dot{q}_{j}(\zeta)}{q_{j}(\zeta)}j^{2}+iP_{j}(\zeta)\right]\\\
&j=\xi,\eta,\theta.\end{split}$

where $q_{j}(\zeta)$ and $P_{j}(\zeta)$ are complex functions:
$\frac{\dot{q}(\zeta)}{q(\zeta)}=\frac{1}{R(\zeta)}+\frac{2i}{W^{2}(\zeta)}$,
where $R(\zeta)$ is the curvature radius and $W(\zeta)$ is the width of the
wave function. The functions $q_{j}(\zeta)$ are given by the equations

$\begin{split}\ddot{q}_{j}(\zeta)+\omega^{2}_{j}(\zeta)q_{j}(\zeta);\end{split}$
(41)

with different frequencies in the spatial directions

$\begin{split}&\omega^{2}_{\xi}(\zeta)=2\left[1-\frac{2}{9(\zeta-\zeta_{0})^{2}}\right];\\\
&\omega^{2}_{\theta}(\zeta)=2\left[1+\frac{1}{9(\zeta-\zeta_{0})^{2}}\right];\\\
&\omega^{2}_{\eta}(\zeta)=2\left[1+\frac{1}{9(\zeta-\zeta_{0})^{2}}\right].\end{split}$
(42)

Thus, the width in each direction evolves oscillating with a $\zeta$-dependent
frequency. In the unperturbed case, the frequency is constant and the
equations have oscillating solutions and a stationary one for the initial
waist $W_{0}^{2}=\sqrt{2}$.
Considering specifically the $\xi$-direction, the frequency becomes imaginary,
as $\omega_{\xi}^{2}<0$ when $\zeta>\zeta_{B}$, with

$\begin{split}\zeta_{B}=\frac{\tau_{B}}{\tau_{p}}=\zeta_{0}-\frac{\sqrt{2}}{3}.\end{split}$
(43)

Correspondingly, the solution of the width equation $q_{\xi}(\zeta)$ from
oscillating becomes exponential. As a consequence, the wave function $u_{\xi}$
ceases to be a Gaussian, the BEC is no more localized and a phase transition
occurs.
We show in Figure 4 the width of the $\xi$-dependent solution during the
evolution, including the critical case when the initial waist corresponds to
the stationary solution of the unperturbed Eq. (41). $W$ diverges for times
greater than the breaking time $\tau_{B}$.
We consider specifically the dependence of the breaking time $\tau_{B}$ on the
BEC mass (Fig. 5). Given the BH mass, which fixes the Schwarzschild time
$\tau_{s}$, the breaking time increases with the BEC mass. Below a BEC
critical mass, the breaking time is smaller than the Schwarzschild time, thus
the DM halo becomes unstable before crossing the event horizon giving rise to
an observable instability. we remark that at a time comparable with the
breaking time, the perturbative approach is expected not to provide an
accurate description of the condensate, and also the mean field approximation
is expected to fail, while opening the road to novel and unexplored quantum-
gravitational effects.

Figure 4:  a) BEC wave function width along the $\xi$ direction during the
evolution. Before the breaking time, the mass density distribution has a
constant width or an oscillating one according to its initial value. The total
potential is attractive and the BEC is in a bound state. After the breaking
time, the width follows an exponential trend independently from the previous
history. Since the potential is repulsive, the BEC delocalizes. b) Detail of
the unstable evolution. Figure 5:  Breaking time versus the BEC mass and
Schwarzschild time (green line).

##  VI The astrophysical scenario

In the astrophysical context, the relevant parameters are the masses: the
single particle $m$, the DM halo $M$, and the BH $M_{BH}$, chosen as follows

$\begin{split}&\text{BEC:}\quad m\sim 10^{-22}\text{eV};\qquad M\sim
10^{6}M_{\odot};\\\ &\text{BH:}\quad M_{BH}\sim 10^{9}M_{\odot};\end{split}$
(44)

where $M_{\odot}$ is the solar mass.
The ultralight single particle provides the formation of galactic BECs, in
accordance with previous work [3, 6]. The supermassive BH suits the
description of galactic core phenomena.
We have the characteristic frequency of the self-gravitational potential,
$\Omega$

$\Omega=\frac{16G^{2}m^{5}N^{2}}{9\pi\hbar^{3}}\sim 10^{-11}\,\text{rad/s},$

and the width of the BEC wave function

$w_{0}=\sqrt{\frac{\hbar}{m\Omega}}\sim 1.2\,\times 10^{7}\text{m},$

extending over galactic distances.
From the BH mass, we find the Schwarzschild radius

$R_{s}=\frac{2GM_{BH}}{c^{2}}\sim 10^{12}\,\text{m}.$

Assuming a starting distance between the center of the BEC and the BH of order
$R_{0}\sim 10^{5}\,R_{s}\sim 10^{17}\text{m}$, we determine the relevant time
scales

$\begin{split}&\tau_{0}=\frac{2R_{0}}{3c}\sqrt{\frac{R_{0}}{R_{s}}}\sim
2\times 10^{11}\,\text{s}\sim 6600\text{y};\\\
&\tau_{s}=\tau_{0}\left(1-\frac{R_{s}^{3/2}}{R_{0}^{3/2}}\right)\sim 2\times
10^{11}\,\text{s}\sim 6600\text{y};\\\
&\tau_{B}=\tau_{0}-\frac{2}{3\Omega}\sim 1.7\times 10^{11}\,\text{s}\sim
5000\,\text{y}.\end{split}$

DM instability occurs when the breaking time is shorter than the Schwarzschild
time, in a way dependent on the BEC mass and the interacting massive objects.
The onset of instability may explain some astrophysical observations,
including the DM absence in galaxies of different sizes.

##  VII The laboratory analog

The laboratory realization of an analog system for our model equation (22) is
based on the production of a BEC from a gas of bosons like Sodium atoms. The
condensation is achieved by cooling and confining a Bose gas with tools like a
magneto-optic trap [18, 19].
By using the explicit form for the radius $R(\tau)$ in (32), we find that the
total potential Eq. (24) is

$\frac{m\Omega^{2}}{2}\vec{r}^{\,2}+\frac{m}{9(\tau-\tau_{0})^{2}}(-2x^{2}+y^{2}+z^{2}).$
(45)

This potential can be realized by an optical trap with a monochromatic beam of
frequency $\omega$ and intensity $I(\vec{r})$ that generates a spatial
dependent potential of the form

$V(\vec{r})=\frac{3\pi c^{2}\Gamma}{2\omega_{0}^{3}\Delta}I(\vec{r})$ (46)

where $\omega_{0}$ is the characteristic resonant frequency of the trapped
gas, $\Delta=\omega-\omega_{0}$ is the laser detuning, and $\Gamma$ is the
damping rate due to the radiative energy loss.
The constant term of the potential (45) is obtained by a $3$-dimensional
optical trap, where $3$ orthogonal laser beams have constant intensity and
negative detuning $(\Delta<0)$, namely red-detuned with respect to the
condensate characteristic frequency. For the time-dependent term of the
potential, we need time-varying intensity beams and a repulsive trap along the
$x$-direction. With optical traps, a repulsive potential is obtained by a
blue-detuned $(\Delta>0)$ laser beam [18]. The whole potential is a composite
optical trap (Figure 7). The red-detuned laser beams at
$\lambda=700\,\text{nm}$ comprise the constant parabolic term analog to the
self-gravitation, while the red laser at $\lambda=760\,\text{nm}$ and the
green one at $\lambda=540\,\text{nm}$ are the analogs of the time-dependent
potential generated by the BH and felt by the BEC during its radial fall. The
laser beam intensities vary in time as reported in Figure 6.
For a BEC of Sodium atoms

$\begin{split}&\text{Atomic mass:}\quad m_{A}\sim 23\,\text{u.m.a.}\sim
10^{-26}\,\text{kg};\\\ &\text{Resonant frequency:}\quad\omega_{0}\sim 3\times
10^{15}\,\text{rad/s}.\\\ \end{split}$

Red-detuned laser beams with constant power generate the potential well
mimicking the BEC DM self-gravitational potential

$\begin{split}\text{Laser wavelength:}\quad&\lambda\sim 700\,\text{nm};\\\
\text{Laser beam waist:}\quad&w_{L}\sim 1\,\text{mm};\\\ \text{Laser
power:}\quad&P\sim 20\,\text{mW};\\\ \text{Detuning:}\quad&\Delta\sim-3\times
10^{14}\text{rad/s}.\end{split}$

In the laboratory simulation, the final value $\tau_{0}$ corresponds to the
duration of the experiment. A reasonable value for according to BEC lifetime
is $\tau_{0}\sim 1\,\text{s}$. To be consistent with our model we assume the
DM particle to have mass $m\sim 10^{-22}\,\text{eV}$. With these parameters we
can mimic the following astrophysical scenario:

$\begin{split}\text{BEC DM waist:}\quad&w_{0}\sim 10^{9}\text{m};\\\ \text{BEC
total mass:}\quad&M_{\text{BEC}}\sim 10^{4}\,M_{\odot};\\\ \text{BH
mass:}\quad&M_{\text{BH}}\sim 10^{7}M_{\odot};\\\ \text{Schwarzschild
radius:}\quad&R_{s}\sim 10^{10}\,\text{m};\\\ \text{Breaking
time:}\quad&\tau_{B}\sim 10^{9}\,\text{s}.\end{split}$

The corresponding breaking time in the laboratory is $\tau_{L}\sim
0.2\,\text{s}$. At this temporal scale, the BEC becomes unstable, and one
observes the wave breaking caused by the interaction with the BH.

Figure 6:  Laser power of the optical trap of Sodium atoms to mimic the total
potential exerted by a Schwarzschild black hole on a galactic BEC. The red-
detuned laser beams at $\lambda=700\,\text{nm}$ generate the constant term of
the potential, while the time-varying beams mimic the potential felt by the
halo during the radial fall. The beam at $\lambda=760\,\text{nm}$ produces an
attractive potential while the green laser at $\lambda=540\,\text{nm}$ induces
a repulsive potential.  Figure 7:  Sketch of the optical trap to mimic the
gravitational interactions in the astrophysical DM BEC. a) An optical trap
with lasers of constant intensity realizes the self-gravitational potential;
b) two additional time-dependent laser beams generate the time-varying
potential of the BH on the approaching BEC. The detuning between the lasers
and the BEC atoms induces deforming tidal forces. The red laser realizes an
attractive potential, while the green one is repulsive.

##  VIII Conclusions

Following the hypothesis of the solitonic bosonic nature of interacting dark
matter, we investigated the evolution of a Bose-Einstein condensate in curved
spacetime. We considered a self-gravitating dark matter halo condensate in the
Fermi normal coordinates, which express the d’Alembert operator as the flat
Laplacian with parabolic terms generated by the curvature. By a multiple scale
expansion, we studied the gravitational interaction with a black hole within
highly nonlocal approximation. The resulting effective potential in the GPE is
parabolic with a time-dependent frequency.
As observed in the literature [3], the collective bosonic wave function
explains many of the DM anomalous observations only if the BEC is in its
excited states. In our model, the tidal forces in the effective potential
induce exponential deformations and instability in the dark matter halo, which
excite higher-order internal modes. The BEC, initially in the ground state of
the attractive self-gravitating potential, undergoes a perturbation due to the
curved spacetime and evolves into a superposition of excited states.
In the direction of the repulsive effective potential, the BEC undergoes a
transition from the stable harmonic oscillator to an unstable inverted
oscillator [20]. The instability is reached at a breaking time, which grows
with the DM and decreases with the BH mass. This may be used to explain the
presence, or the absence, of DM halos in galaxies of different sizes.
Our model is not only a simple closed-form analysis of the interaction of DM
with BH but suggests simulating the physics of DM interaction with a curved
spacetime in the laboratory by time-varying traps. This may open the way to
study first and second quantized dark matter dynamics [21, 22] in simulated
gravity with atoms, polaritons, and optical systems [23, 24, 25].

###### Acknowledgements.

We thank prof. Paolo Maria Santini for the fruitful discussions.

## References

  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Particle dark matter: Evidence, candidates and constraints, Physics reports 405, 279 (2005).
  * Suárez  _et al._ [2014] A. Suárez, V. H. Robles, and T. Matos, A review on the scalar field/bose-einstein condensate dark matter model, in  _Accelerated Cosmic Expansion: Proceedings of the Fourth International Meeting on Gravitation and Cosmology_ (Springer, 2014) pp. 107–142.
  * Sin [1994] S.-J. Sin, Late-time phase transition and the galactic halo as a bose liquid, Physical Review D 50, 3650 (1994).
  * Boehmer and Harko [2007] C. Boehmer and T. Harko, Can dark matter be a bose–einstein condensate?, Journal of Cosmology and Astroparticle Physics 2007 (06), 025.
  * Crăciun and Harko [2020] M. Crăciun and T. Harko, Testing bose–einstein condensate dark matter models with the sparc galactic rotation curves data, The European Physical Journal C 80, 735 (2020).
  * Schive  _et al._ [2014] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Cosmic structure as the quantum interference of a coherent dark wave, Nature Physics 10, 496 (2014).
  * Paredes and Michinel [2016] A. Paredes and H. Michinel, Interference of dark matter solitons and galactic offsets, Physics of the Dark Universe 12, 50 (2016).
  * Wang and  _et al._ [2020] J. Wang and S. Bose, C. S. Frenk, L. Gao, A. Jenkins, V. Springel, W. Simon, Universal structure of dark matter haloes over a mass range of 20 orders of magnitude, Nature 585, 39–42 (2020).
  * Garnier  _et al._ [2021] J. Garnier, K. Baudin, A. Fusaro, and A. Picozzi, Incoherent localized structures and hidden coherent solitons from the gravitational instability of the schrödinger-poisson equation, Physical Review E 104, 054205 (2021).
  * Barranco  _et al._ [2012] J. Barranco, A. Bernal, J. C. Degollado, A. Diez-Tejedor, M. Megevand, M. Alcubierre, D. Núñez, and O. Sarbach, Schwarzschild black holes can wear scalar wigs, Physical review letters 109, 081102 (2012).
  * O’dell  _et al._ [2000] D. O’dell, S. Giovanazzi, G. Kurizki, and V. M. Akulin, Bose-einstein condensates with 1/r interatomic attraction: Electromagnetically induced “gravity”, Physical Review Letters 84, 5687 (2000).
  * Barceló  _et al._ [2011] C. Barceló, S. Liberati, and M. Visser, Analogue gravity, Living reviews in relativity 14, 1 (2011).
  * Dalfovo  _et al._ [1999] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of bose-einstein condensation in trapped gases, Reviews of modern physics 71, 463 (1999).
  * Thorne  _et al._ [2000] K. S. Thorne, J. A. Wheeler, C. W. Misner, Gravitation, Freeman San Francisco, CA (2000).
  * Exirifard and Karimi [2022] Q. Exirifard and E. Karimi, Schrödinger equation in a general curved spacetime geometry, International Journal of Modern Physics D 31, 2250018 (2022).
  * Manasse and Misner [1963] F. Manasse and C. W. Misner, Fermi normal coordinates and some basic concepts in differential geometry, Journal of mathematical physics 4, 735 (1963).
  * Aleksić  _et al._ [2014] B. N. Aleksić, N. B. Aleksić, M. S. Petrović, A. I. Strinić, and M. R. Belić, Variational and accessible soliton approximations to multidimensional solitons in highly nonlocal nonlinear media, Optics Express 22, 31842 (2014).
  * Grimm  _et al._ [2000] R. Grimm, M. Weidemüller, and Y. B. Ovchinnikov, Optical dipole traps for neutral atoms, in  _Advances in atomic, molecular, and optical physics_ , Vol. 42 (Elsevier, 2000) pp. 95–170.
  * Metcalf and van der Straten [2003] H. J. Metcalf and P. van der Straten, Laser cooling and trapping of atoms, JOSA B 20, 887 (2003).
  * Marcucci and Conti [2016] G. Marcucci and C. Conti, Irreversible evolution of a wave packet in the rigged-hilbert-space quantum mechanics, Physical Review A 94, 052136 (2016).
  * Howl  _et al._ [2021] R. Howl, V. Vedral, D. Naik, M. Christodoulou, C. Rovelli, and A. Iyer, Non-gaussianity as a signature of a quantum theory of gravity, PRX Quantum 2, 010325 (2021).
  * Conti [2023] C. Conti, Random walk and non-gaussianity of the 3d second-quantized schrödinger–newton nonlocal soliton, New Journal of Physics 25, 023026 (2023).
  * Carusotto and Ciuti [2013] I. Carusotto and C. Ciuti, Quantum fluids of light, Reviews of Modern Physics 85, 299 (2013).
  * Calvanese Strinati and Conti [2014] M. Calvanese Strinati and C. Conti, Bose-einstein condensation of photons with nonlocal nonlinearity in a dye-doped graded-index microcavity, Phys. Rev. A 90, 043853 (2014).
  * Klaers  _et al._ [2010] J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz, Bose–einstein condensation of photons in an optical microcavity, Nature 468, 545 (2010).