# Coherent and incoherent structures in fuzzy dark matter halos

I-Kang Liu, Nick P. Proukakis, and Gerasimos Rigopoulos
School of Mathematics, Statistics and Physics, Newcastle University, Newcastle
upon Tyne, NE1 7RU, United Kingdom
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ. Comments:)

###### Abstract

We show that fuzzy dark matter halos exhibit spatial differentiation in the
coherence of the field configuration, ranging from completely coherent in the
central solitonic core to incoherent outside it, with a significant drop of
the dimensionless phase-space density with increasing radius. The core is a
pure condensate, overlapping perfectly with the Penrose-Onsager mode
corresponding to the largest eigenvalue of the one-particle density matrix.
The virialized outer halo exhibits no clear coherence as a whole upon radial
and temporal averaging but can be described as a collection of local, short-
lived quasi-condensate lumps, exhibiting suppressed fluctuations, which can be
identified with the structures commonly referred to as granules. These
localized regions are separated by vortices that form a dynamical web
inhibiting phase coherence across the entire halo. We further examine the core
oscillations, finding that they are accurately described by two time-dependent
parameters characterizing the size of the core, $r_{c}(t)$, and a crossover
region, $r_{t}(t)$. For the halos in our merger simulations this feature is
reflected in the (anti-)correlated oscillation of the peak value of the
density power-spectrum. The halo’s turbulent vortex tangle appears to reach a
quasi-equilibrium state over probed timescales, with the incompressible
component of the kinetic energy exhibiting a characteristic $k^{-3}$ tail in
its spectrum, indicative of a $\rho\sim r^{2}$ density profile around the
quantum vortex cores. Comparison of the peak wavenumbers in the corresponding
power-spectra shows the inter-vortex spacing and the granule length scale in
the outer halo to be very similar and slightly above the core size.

###### keywords:

methods: numerical – galaxies: halos – dark matter

††pubyear: 2022††pagerange: Coherent and incoherent structures in fuzzy dark
matter halos–LABEL:lastpage

##  1 Introduction

The idea that dark matter may be made up of bosonic particles, light enough to
exhibit coherent dynamics and wave-like behaviour on galactic scales, was
introduced some time ago [Hu et al. (2000)] as a solution to a number of
perceived problems exhibited by CDM on sub-galactic length scales [Weinberg et
al. (2015); Bullock & Boylan-Kolchin (2017); Del Popolo & Le Delliou (2017)]
and has attracted significant attention recently, especially after the
pioneering simulations of [Schive et al. (2014a)] - see [Marsh (2016); Hui
(2021); Ferreira (2021)] for recent reviews. Sometimes referred to as Fuzzy
Dark Matter (FDM), Wave Dark Matter or $\psi$-Dark Matter ($\psi$DM), in
contrast to the prevailing Cold Dark Matter model (CDM), it presents a
different phenomenology to that of CDM on short scales but is
indistinguishable from it on longer scales [Hui et al. (2017)] and can thus
reproduce CDM’s remarkable ability to explain the whole history of structure
formation in the universe [Frenk & White (2012)]. The main parameter
delineating these two regimes is the boson’s mass with a usually quoted value
within an order of magnitude of $m\sim 10^{-22}$ eV/$c^{2}$ in FDM models;
particles of such mass, and with typical velocities found in halos hosting
Milky Way-sized galaxies, acquire a de Broglie wavelength of $\mathcal{O}(1)$
kpc and can therefore exhibit wave-like features at galactic scales. Various
large scale structure observations have begun placing constraints on the
boson’s mass [Schive et al. (2016); Rogers & Peiris (2021); Dentler et al.
(2022)] but a lot of work remains to be done, especially to properly include
effects on nonlinear scales which can only be approached via numerical
simulations [Schive et al. (2014a); Schwabe et al. (2016); Mocz et al. (2017);
May & Springel (2021)].

There is ongoing debate on whether the aforementioned short scale features of
CDM cosmological structures pose a fundamental problem for the model,
requiring an alteration of dark matter’s fundamental properties, or if it is a
matter of proper modeling of complex sub-galactic baryonic physics (gas, SN
explosions etc) that feed back on the dark matter short scale distribution
[Bullock & Boylan-Kolchin (2017)]. Even if the latter view ends up explaining
all discrepancies, it is nevertheless interesting to explore the implications
of the FDM idea with its different phenomenology acting as an extension of CDM
that can be constrained, allowing for example to place limits on the mass (or
self-coupling) of scalar Dark Matter particles.111Here we set the self-
coupling to zero and focus on the Schrödinger-Poisson system of equations. We
will raise this restriction by introducing a quartic self-interaction and
study the Gross-Pitaevskii-Poisson system in upcoming work.

Beyond observationally motivated considerations, examining the Schrödinger
equation with the inclusion of the gravitational force is an interesting
physics exercise in itself. The phenomenon of Bose-Einstein condensation is
universal [Proukakis et al. (2017)] and has been studied for a long time in
various weakly-interacting laboratory-controlled settings, such as atomic
(ultracold quantum matter [Anderson et al. (1995); Davis et al. (1995);
Pitaevskii & Stringari (2003); Pethick & Smith (2008)]), solid-state/optical
(exciton-polariton [Kasprzak et al. (2006)] and photon [Klaers et al. (2010)])
systems, alongside studies in superfluid liquid Helium [Kapitza (1938); Allen
& Misener (1938); London (1938)], in which the strength of interactions leads
to a significantly reduced condensate fraction, even in a pure superfluid
[Penrose & Onsager (1956)]. Unlike the more strongly interacting helium
superfluids, atomic/optical condensates lie in parameter regimes where nearly
pure condensates (with a condensate particle fraction approaching $100\%$) can
be achieved, and facilitate control of many parameters, including trapping
potential and interaction strength, with various ultracold atomic species
featuring a strong dipolar component and facilitating a controllable interplay
between contact and long-range interactions [Stuhler et al. (2005); Ferrier-
Barbut et al. (2016)]. Therefore, examining the features introduced by the
long range gravitational interaction may lead to behaviours that could bear
similarities to atomic condensates forming in external potential traps [O’Dell
et al. (2000); Papadopoulos et al. (2007)]. Accordingly, insights from the
field of cold atom physics may be used to further understand the structure and
dynamics of FDM halos.

In this paper we examine the structure and coherence of a single virialised
FDM halo: in our simulations, such a halo is formed from the merger of a
number of initially coherent soliton spheres placed at random initial
positions. After the resulting halo has virialised, the well known central
soliton has formed and the Navarro–Frenk–White (NFW) density profile has been
established outside it, we examine the coherence of the field configuration in
the spatial domain by measuring (through temporal and generally also angular
averaging) the first and second order spatial correlation functions
characterizing phase coherence and density fluctuations respectively. We also
directly measure the coherent, Penrose-Onsager mode, corresponding to the
single-particle density matrix’s eigenmode with the largest eigenvalue, which
ends up overlapping with the empirical soliton density profile. Such analysis
allows us to clearly distinguish a spatial structure consisting of the inner
coherent (solitonic) core, a crossover region and a largely incoherent outer
halo, which can be parameterized by two characteristic radii $r_{c}$ and
$r_{t}$ delineating these regimes. Our simulations also reveal the previously
reported core oscillations [Li et al. (2021); Chiang et al. (2021); Dutta
Chowdhury et al. (2021)] which we find here to be very well captured by
assigning time dependence to $r_{c}$.

Bose-Einstein condensation corresponds to the existence of macroscopic
occupation in a particular state, exhibiting off-diagonal-long-range-order
(ODLRO) and implying a two-point correlation function which does not tend to
zero as $|\mathbf{r}-\mathbf{r}^{\prime}|\rightarrow\infty$, a feature absent
in the normal phase. In textbook statistical mechanics, three-dimensional (3D)
homogeneous, non-interacting systems exhibit Bose-Einstein condensation when
the so-called (dimensionless) phase space density ${\cal
D}=(\rho/m)\lambda_{\rm dB}^{3}>2.612$ [Pethick & Smith (2008)] (or more
generally ${\cal D}\gtrsim O(1)$). In other words, the emergence of Bose-
Einstein condensation arises when the characteristic de Broglie wavelength
$\lambda_{\rm dB}$ exceeds the average interparticle distance. In this work,
we thus also investigate the value of ${\cal D}$, as a relevant probe of the
wave-mechanical nature of the system.

In a typical 3D harmonically-confined atomic system, the central high-density
region is dominated by an inhomogeneous coherent condensate component, with
lower-density regions around such a core consisting of incoherent (thermal)
particles. Moreover, the existence of strong phase fluctuations can also lead
to the emergence of a so-called quasi-condensate state [Popov (1983); Kagan &
Svistunov (1997)], which features (relatively) suppressed density
fluctuations, but sustained phase fluctuations. Such a state, associated with
long-wavelength fluctuations destroying the phase coherence, is particularly
pronounced in equilibrium lower-dimensional (quasi-2D, quasi-1D) laboratory
settings [Prokof’ev et al. (2001); Petrov et al. (2000, 2001); Dettmer et al.
(2001); Al Khawaja et al. (2002); Proukakis (2006); Cockburn et al. (2011);
Hadzibabic & Dalibard (2011)] (which facilitate a decoupling between
characteristic temperatures for density and phase fluctuations) or close to
the critical region, and emerges as a transient state during the formation of
Bose-Einstein condensation [Berloff & Svistunov (2002); Proukakis et al.
(2017)].

Fuzzy dark matter halos show a combination of these features, with all four
coherence measures we employ resulting in quantitatively compatible radial
profiles (as shown in subsequent figure 3). Firstly, such systems typically
satisfy the phase-space density condition (see also subsequent parameter
regime discussion), thus justifying a classical field description throughout,
as implemented in this work. Moreover, such systems exhibit behaviour rather
analogous to the harmonically-confined atomic systems: in particular, they
consist of a coherent solitonic core surrounded by fluctuating regions of
significantly lower density and coherence, mediated by an intermediate region
over which the system density fluctuations become gradually enhanced. The
solitonic core is established by the attractive gravitational potential, which
plays the role of the trap concentrating the particles. Density and phase
fluctuations increase gradually, but significantly, over the immediate
surrounding region, with such region further embedded within a fluctuating
halo state which, on average, features no phase, or density, coherence.
Nonetheless, closer inspection of such outer halo region reveals a much richer
structure featuring the co-existence of patches of slowly-varying density and
variable (but often significant) phase coherence, separated by vortices, which
are themselves dynamically propagating (but not dissipating) through the halo,
ensuring any local coherence does not extend into the surrounding halo
regions.

The coherence of both phase and density in the soliton core region remains
intact throughout the post-virialization evolution of the halo. In contrast,
the outer halo exhibits both phase and density fluctuations222Interestingly,
our findings appear compatible with the qualitative suggestion of Guth et al.
(2015) for the nature of the equilibrium state of axionic particles with
attractive self-interactions, including gravity.; we explicitly show that they
peak at the same characteristic ‘granular’ length scale and are suppressed
below it, thus giving rise to a rapidly-fluctuating tangle of quantum vortices
intertwined with the granular density distribution, as first suggested in [Hui
et al. (2021)]. We further consolidate this picture through studying the
relevant power spectra, the second object of this study. It should be noted
that, interestingly, for the halos in our merger numerical experiments, the
power spectrum of the density field also carries the signature of the
oscillating core

Beyond the formation of a stable, density and phase coherent solitonic core in
virialized halos, the other main difference of FDM compared to CDM is the way
vorticity manifests itself and is what we also focus on in this work. Any
initial vorticity present in the early universe is expected to decay like
$1/a$ in the matter-dominated era, as long as density perturbations remain
small and linear perturbation theory is sufficient for describing the
evolution of the density field [Durrer (2020)]. Hence, a vorticity-free
initial condition is a standard consideration in late universe cosmology
[Jelic-Cizmek et al. (2018)]. As long as CDM can be modeled by a pressure-less
perfect fluid, vorticity will not reappear. As perturbations grow however, CDM
enters the so called multi-stream regime, the perfect fluid approximation
breaks down and vorticity is generated [Pueblas & Scoccimarro (2009)]. In the
FDM model vorticity also arises during halo formation but, unlike corpuscular
CDM, a rotational velocity field can only be accommodated through the
appearance of quantized vortices, a well-studied subject in the field of
laboratory condensates. The presence of vortices in the central region of a
rotating FDM halo could deform it from spherical symmetry to a disk-spiral
shape [Zinner (2011)]. However, it has been theoretically found (for Riemann-S
ellipsoid solutions) that a vortex could not be stable without a strong enough
repulsive self-coupling between FDM particles and a heavy enough bosonic mass
such that the core density is/becomes low enough to energetically support the
presence of a vortex [Rindler-Daller & Shapiro (2012, 2014); Dmitriev et al.
(2021); Schobesberger et al. (2021); Nikolaieva et al. (2021)]. The dynamical
instability of a vortex sitting at the centre of a halo leads to its migration
out of the core which would temporarily split [Dmitriev et al. (2021)] into a
number of fragments dependent on the charge of the vortex [Nikolaieva et al.
(2021)]. These instabilities lead to the expectation that in a dynamical halo
formation no vortices will reside in the core. Indeed, we find that the core
region remains shielded from vortices, retaining its density and phase
coherence throughout our simulations.

As mentioned above, the density beyond the confines of the core is
characterised by dynamical granular structures and vortex filaments that form
topologically closed loops which deform, stretch and shrink during the time
evolution [Hui et al. (2021)] - see Fig. 1 and section 5. The typical inter-
vortex separation scale was suggested in [Mocz et al. (2017); Hui et al.
(2021)] to be comparable to the soliton core size. We find it to be a bit
larger than that (corresponding to the extended core including the crossover
region discussed below) and roughly equal to the granule size.

Furthermore, [Mocz et al. (2017)] has numerically demonstrated that vortices
can exhibit turbulent features driven by their reconnections within the vortex
tangle that permeates the halo. The existence of a tangle of dynamically-
evolving vortices separating regions of near-constant phase – and thus
facilitating short-range superfluid order – has been previously discussed in
the context of condensate formation from an initial highly non-equilibrium
state, proceeding through regimes of superfluid turbulence with quasi-
condensate local correlations [Berloff & Svistunov (2002)]. In such works,
which dealt with homogeneous systems, the emerging vortex tangle was
numerically observed to relax over long times, with such decaying behaviour
also qualitatively persisting in the context of equilibration in an elongated
3D harmonic trap [Liu et al. (2018)], where coherence becomes established from
the trap centre outwards [Liu et al. (2020)]. In our present setting, we do
find the emergence of a phase coherent central solitonic core, surrounded by a
halo featuring short-range order across various locations within our
simulation box, but with long-range order inhibited by the random motion of
the quantum vortices. Unlike such earlier simulations, however, the vortex
tangle in the gravitational case does not seem to be noticeably decaying, at
least for the rather long timescales probed: such difference could be due to
the absence of efficient mode-mixing between different modes of the system in
the absence of a self-interaction [Chiueh et al. (2011)]. Although the
stability of such vortex tangle remains an open question, the present case
could constitute a cosmological-scale example of a self-sustained state of
strong turbulence which has no current analogue in laboratory condensates.

The incompressible (rotational) kinetic energy spectrum carries a clear
signature of the vortices and, through the position of its peak, provides the
inter-vortex distance length scale. As mentioned above, we find the latter to
be equal to the size of the granular structures that give a “wavy” appearance
to fuzzy dark matter halos. Beyond the peak, we find the incompressible
kinetic energy spectrum in the halo to exhibit a $k^{-3}$ tail, also found in
[Nore et al. (1997a)]. This spectrum demonstrates that vortices are not ideal
line-like defects but have internal structure and a characteristic size in the
superfluid [Dalfovo et al. (1999)]. Indeed, the $k^{-3}$ tail of the
incompressible kinetic energy spectrum implies that $\rho v^{2}\approx{\rm
const}$ in the vincinity of the vortex centre. Since the velocity profile of a
quantized vortex is $v\sim r^{-1}$[Nore et al. (1997a); Dalfovo et al. (1999)]
it directly follows that $\rho\approx r^{2}$ with the centrifugal contribution
in the quantum kinetic energy depleting the density around the vortex centre.
Using the spectrum of the total velocity field, Mocz et al. (2017) has
verified the velocity profile of a quantized vortex to be close to $v\sim
r^{-1}$. On the other hand, Hui et al. (2021) has demonstrated that the
density profile exhibits a $\rho\sim r^{2}$ behaviour in the vicinity of an
isolated vortex core in the Schrödinger-Poisson case. In our simulations we
verify the density profile around FDM vortices in the fully dynamical vortex
tangle using the spectrum of the incompressible kinetic energy, a measure more
commonly employed in quantum turbulence studies.

Figure 1: Rendered illustration of the core (inner yellow region) and outer
halo (blue region) density profiles, on which the co-existing topologically-
closed tangled vortical structures (purple lines) are superimposed, based on
our Schrödinger-Poisson simulations. The vortical structures, whose
characteristic separation scale is a bit larger than that of the soliton core,
only appear outside the outer red dashed circle that demarcates a transition
zone from coherent to incoherent regions, and do not enter the soliton core
region (within the dashed blue circle) once the virialized halo has formed.

A schematic of the core-outer halo features that are at the heart of the study
presented in this work is depicted in Fig. 1. It is a rendered image for the
mass density and spatial distribution of vortices and vortex rings in the
final virialized halo. The blue-yellow hues in the figure correspond to the
density amplitude while the purple-red is the volume-rendering profile of the
modulus square of the velocity field’s rotational part in the outer region of
the halo. The inner high-density region supports the purely coherent core of
condensed axions for $r\leq r_{c}$ (within the dashed blue circle); this is
surrounded by a broader radial ring delineating the crossover region (dotted
red circle), where phase and density fluctuations start appearing; at larger
radii from the solitonic core the phase coherence is reduced and the system
becomes chaotic, where a large number of tangled quantum vortices lie between
regions of otherwise nearly-constant density which also exhibit variable local
coherence. We will see that the inter-vortex distance and ring size are
slightly larger than the extended core diameter $\sim 3r_{t}$. We discuss the
details of this image below, the aim being to characterize these distinct
features using a combination of characterization techniques from statistical
physics, quantum optics, cold quantum matter, condensed matter and turbulence
theory.

The paper is structured as follows: In section 2, we discuss the model
formulation of FDM and its hydrodynamic mapping and outline our numerical
method, modelling scheme and procedure of the performed soliton merger
simulation which allows us to characterize the static and dynamical core-halo
properties. In section 3 we discuss the interplay between coherent and
incoherent features through various appropriate characterizations, and present
the corresponding energy distributions. In section 4 we discuss the breathing
oscillation exhibited in the soliton core and its correlation to the peak
wavenumber of the power spectrum for a full halo density profile in the cases
we study. In section 5, we demonstrate the existence of a large number of
vortices in the outer halo and use the granule power and vortex energy spectra
to obtain the relation of the granule-size and the inter-vortex distance
length scales. Our conclusions and final remarks are given in section 6.
Numerical and other technical details, along with some statistical
measurements supporting the universality of the presented picture are given in
the Appendices.

##  2 Theoretical background and Numerical Implementation

###  2.1 Schrödinger-Poisson equations and corresponding hydrodynamic
description

The ultra-light scalar field making up fuzzy dark matter can be described in
the non-relativistic limit by the Schrödinger-Poisson system of equations
(SPE) [Ruffini & Bonazzola (1969); Widrow & Kaiser (1993); Sin (1994); Lee &
Koh (1996); Coles & Spencer (2003); Uhlemann et al. (2014)]

$i\hbar\frac{\partial}{\partial
t}\Psi(\mathbf{r},t)=\left[-\frac{\hbar^{2}\nabla^{2}}{2m}+m\Phi(\mathbf{r},t)\right]\Psi(\mathbf{r},t)$
(1)

where $\Psi=\Psi(\mathbf{r},t)$ is the wavefunction of fuzzy dark matter in
the physical coordinate $\mathbf{r}$, and the gravitational potential
$\Phi(\mathbf{r},t)$ is determined via the Poisson equation

$\nabla^{2}\Phi(\mathbf{r},t)=4\pi
G\left[\rho(\mathbf{r},t)-\bar{\rho}\right]$ (2)

with $\rho(\mathbf{r},t)=|\Psi(\mathbf{r},t)|^{2}$ the mass density, and
$\bar{\rho}$ its spatially averaged, constant value. As we are focusing on the
dynamics of a virialized halo, we would not expect the cosmic expansion to
have a strong influence on the halo structures we wish to study here, so we do
not account for it in the present work. The subtraction of $\bar{\rho}$ from
the r.h.s. of (2) is reminiscent of the Jeans swindle for a fluid of mass
density coupled to Newtonian gravity [Binney & Tremaine (2011)] but actually
has a rigorous mathematical basis for Newtonian gravity [Kiessling (1999);
Chavanis (2011)]. It is also associated with the periodic boundary condition
for the numerical solution found by the implementation of the pseudo-Fourier
spectrum method [Dabo et al. (2008)], which is commonly used in most studies
and also in this work.

The system involves three main conserved quantities [Chavanis (2011); Chavanis
& Delfini (2011); Nikolaieva et al. (2021); Schobesberger et al. (2021)]: the
total mass

$M=\int d^{3}\mathbf{r}\,\rho(\mathbf{r},t)$ (3)

the energy

$E=\int
d^{3}\mathbf{r}\,\,\Psi^{\ast}(\mathbf{r},t)\left[-\frac{\hbar^{2}\nabla^{2}}{2m^{2}}+\frac{1}{2}\Phi(\mathbf{r},t)\right]\Psi(\mathbf{r},t)$
(4)

obtained from the usual expression after integration by parts and ignoring the
boundary terms, and the angular momentum

$\mathbf{L}=\int
d^{3}\mathbf{r}\,\,\Psi^{\ast}(\mathbf{r},t)\left[\mathbf{r}\times\left(\frac{\hbar}{im}\nabla\right)\right]\Psi(\mathbf{r},t)$
(5)

which is zero in the case we examine here (as also confirmed numerically).
Expression (4) for the energy functional can also be used as an energy
functional of a Bose-Einstein condensate system with the Hartree variational
principle to derive the SPE [Dalfovo et al. (1999); Chavanis (2011); Chavanis
& Delfini (2011)]. The source of gravitational potential in the SPE is the
collective potential of the mass distribution $\rho(\mathbf{r},t)$ which
arises from the long-range, two-body interaction
$U_{2B}(\mathbf{r},\mathbf{r}^{\prime})=Gm/|\mathbf{r}-\mathbf{r}^{\prime}|$
in the context of mean-field description for a Bose-Einstein condensate
[O’Dell et al. (2000); Papadopoulos et al. (2007); Chavanis (2011); Chavanis &
Delfini (2011); Veltmaat et al. (2018); Soto et al. (2022)]. The SPE satisfies
well-known scaling properties for the spatial scale factor $\Lambda$
($x\rightarrow\Lambda\hat{x}$) and boson mass $m\rightarrow\alpha\hat{m}$, see
e.g. [Mocz et al. (2017)], which can be found in Appendix A.

The SPE can be mapped into hydrodynamic equations by means of the Madelung
transformation,
$\Psi(\mathbf{r},t)=\sqrt{\rho(\mathbf{r},t)}e^{i\varphi(\mathbf{r},t)}$ and
$\mathbf{v}(\mathbf{r},t)=(\hbar/m)\nabla\varphi(\mathbf{r},t)$. One thus
obtains

$\frac{\partial}{\partial t}\rho+\nabla\cdot(\rho\mathbf{v})=0$ (6)

and

$\frac{\partial}{\partial
t}\mathbf{v}+\frac{\nabla}{m}\left[\frac{m|\mathbf{v}|^{2}}{2}+m\Phi-\frac{\hbar^{2}}{2m}\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}\right]=0.$
(7)

The equation of motion of $\rho$ is the continuity equation while that of
$\mathbf{v}$ contains the gravitational potential and quantum pressure terms.
The second hydrodynamic equation for the velocity field $\mathbf{v}$ is in
fact identical to the Euler equation for a classical irrotational inviscid
fluid in the limit of $\hbar\rightarrow 0$. This indicates that the quantum
pressure term can be a source of crucial differences between the classical and
quantum fluid but also suggests that classical turbulence can also be present
in quantum fluids. In spite of the SPE and the hydrodynamic equations being
mathematically equivalent, the latter cannot fully resolve numerically regions
of vanishing density that occur due to destructive interference without
further appropriate numerical implementation [Jelic-Cizmek et al. (2018)].
Since sites of vanishing density are places where vortices can form [Li et al.
(2021)], the SPE system is the appropriate framework for studying them.333It
is worth noting that Chiueh et al. (2011) has suggested supplementing the
hydrodynamic equations by equations of motion for the vortices to handle low
density regions of destructive interference. With the Madelung transformation,
the energy of the system, Eq. (4), can be rewritten as a sum of three
different parts,

$E=\displaystyle E_{\mathrm{ke}}+E_{\mathrm{qp}}+E_{\Phi}=\int
d^{3}\textsterling\mathbf{r}\left[\varepsilon_{\mathrm{ke}}(\mathbf{r})+\varepsilon_{\mathrm{qp}}(\mathbf{r})+\varepsilon_{\Phi}(\mathbf{r})\right]$
(8)

where $E_{\mathrm{ke}}$, $E_{\mathrm{qp}}$ and $E_{\Phi}$ are classical
kinetic, quantum pressure and gravitational potential energies respectively,
with the energy densities given by

$\varepsilon_{\mathrm{ke}}(\mathbf{r})=\frac{1}{2}\rho(\mathbf{r})|\mathbf{v}(\mathbf{r})|^{2},\quad\varepsilon_{\mathrm{qp}}(\mathbf{r})=\frac{\hbar^{2}}{2m^{2}}\left|\nabla\sqrt{\rho(\mathbf{r})}\right|^{2},$
(9)

contributed by the Laplacian term in the SPE, and

$\varepsilon_{\Phi}(\mathbf{r})=\frac{1}{2}\Phi(\mathbf{r})\rho(\mathbf{r}).$
(10)

Since the velocity field is defined through the gradient of a scalar field,
one would expect the system to be irrotational, namely
$\nabla\times\mathbf{v}=0$, which would leave no room for vorticity. This
conclusion only holds however if $\varphi$ is continuously differentiable to
second order. If there are points of vanishing mass density, the phase becomes
ill-defined there but the wavefunction can remain single-valued by considering
a phase winding around the singularity being $\Delta\varphi=2\pi\kappa$ for
integer $\kappa$. This suggests that a vortical structure in $\Psi$ appears as
a vortex line with a divergence in $\mathbf{v}$ along the vortex core, leading
to the tangential velocity being

$|\mathbf{v}(r)|=\frac{\hbar}{m}\frac{\kappa}{r}\,.$ (11)

The circulation of the velocity field, specifically for the rotational
(incompressible) part over a closed contour on the plane normal to the vortex
is $\Gamma=\oint_{C}d\mathbf{l}\cdot\mathbf{v}^{i}=(2\pi\hbar/m)\kappa$ ,
yielding $v\propto 1/r$, where $r$ is the distance from the vortex core, while
the vorticity right on a vortex core at $\mathbf{r}_{q}$ is
$\mathbf{\sigma}(\mathbf{r})=\nabla\times\mathbf{v}(\mathbf{r})=\kappa\delta(\mathbf{r}-\mathbf{r}_{q}).$
The divergence in $|\mathbf{v}^{i}(r\rightarrow 0)|$ requires a density
depletion at the center of the vortical structure, forming a topological
defect known as a quantized vortex. Such a defect localized around the vortex
core is characterized by a length scale $\xi$ which, in the absence of self-
interactions, is expected to be of the order of the local de Broglie
wavelength [Hui et al. (2021)].

The existence of vortices in FDM halos was numerically demonstrated in [Mocz
et al. (2017)] and was further examined in [Hui et al. (2021)], where it was
estimated that there should be approximately one vortex line in a unit volume
with length scale $\lambda_{\mathrm{dB}}=2\pi\hbar/m|\mathbf{v}|$, which
should also be the granule length scale. Beyond the velocity field
$\mathbf{v}$ which diverges at a vortex core, a key role in the study of
superfluid turbulence is played by the density-associated current

$\mathbf{F}=\sqrt{\rho}\mathbf{v}$ (12)

which tends to a constant at the vortex core. This is the quantity we will use
to study turbulence in the following. A vector field can always be decomposed
into compressible (irrotational) and incompressible (rotational) parts via the
Helmholtz decomposition, namely,

$\mathbf{v}=\mathbf{v}^{c}+\mathbf{v}^{i}\textrm{ and
}\mathbf{F}=\mathbf{F}^{c}+\mathbf{F}^{i}$ (13)

satisfying

$\nabla\times\mathbf{v}^{c}=\nabla\times\mathbf{F}^{c}=0\textrm{ and
}\nabla\cdot\mathbf{v}^{i}=\nabla\cdot\mathbf{F}^{i}=0$ (14)

and can be computed via Fourier transformation [Kobayashi & Tsubota (2005)].
The rotational component of velocity field is linked to the velocity profile
Eq. (11), and the extremely large $|\mathbf{v}|^{i}$ localized around the
vortex cores can also be used as a probe of vortices. As a result, we can
further decompose the classical kinetic energy into compressible and
incompressible parts,

$E_{\mathrm{ke}}=E_{\mathrm{ke}}^{c}+E_{\mathrm{ke}}^{i}$ (15)

with the corresponding energy densities
$\varepsilon_{\mathrm{ke}}(\mathbf{r})=\varepsilon_{\mathrm{ke}}^{c}(\mathbf{r})+\varepsilon_{\mathrm{ke}}^{i}(\mathbf{r})=(1/2)|\mathbf{F}^{c}(\mathbf{r})|^{2}+(1/2)|\mathbf{F}^{i}(\mathbf{r})|^{2}$.
Such decomposition, which is standard in quantum fluid studies [Nore et al.
(1997b, a); Kobayashi & Tsubota (2005); Numasato et al. (2010); Stagg et al.
(2016)], is rather crucial in the context of superfluid turbulence and the
study of its arising spectra. Note that
$\mathbf{F}^{i}\neq\sqrt{\rho}\mathbf{v}^{i}$ and
$\mathbf{F}^{c}\neq\sqrt{\rho}\mathbf{v}^{c}$ essentially contain similar
information as the full wavefunction $\Psi$ since the Helmholtz decomposition
for $\mathbf{F}$ also involves the spatial variation of the density.

###  2.2 Dimensionless Equations

To solve the SPE numerically we adopt characteristic reference scales for
energy, time and length, respectively

$E_{\mathrm{ref}}=\hbar\sqrt{G\rho_{\mathrm{ref}}},\hskip
2.84544pt\tau_{\mathrm{ref}}=\displaystyle\frac{1}{\sqrt{G\rho_{\mathrm{ref}}}},\hskip
2.84544pt\textrm{and}\hskip
2.84544ptl_{\mathrm{ref}}=\displaystyle\left(\frac{\hbar^{2}}{m^{2}G\rho_{\mathrm{ref}}}\right)^{1/4}.$
(16)

The dimensionless form of the SPE then reads (we suppress labelling the
dimensionless $t$ and $\mathbf{r}$ with different symbols, see Appendix B)

$\begin{array}[]{rl}\displaystyle i\frac{\partial}{\partial
t^{\prime}}\Psi^{\prime}(\mathbf{r},t)=&\displaystyle\left[-\frac{\nabla^{\prime
2}}{2}+\Phi^{\prime}(\mathbf{r},t)\right]\Psi^{\prime}(\mathbf{r}^{\prime},t^{\prime})\\\
\\\ \nabla^{\prime
2}\Phi^{\prime}(\mathbf{r},t)=&4\pi\varrho^{\prime}\left(|\Psi^{\prime}(\mathbf{r}^{\prime},t^{\prime})|^{2}-1\right)\end{array}$
(17)

where $\varrho^{\prime}=\rho_{0}/\rho_{\rm ref}$ is introduced to tune the
total mass of simulation, $M_{\rm tot}=\rho_{0}V=(\rho_{\rm ref}l^{3}_{\rm
ref})\varrho^{\prime}V^{\prime}$ with the dimensionless volume of the
computational box $V^{\prime}$. We have additionally scaled the density to
$\rho_{\rm 0}$ and wavefunction to $1/\sqrt{\rho_{0}}$ giving the
dimensionless wavefunction normalized to $V^{\prime}$. Note that $E_{\rm ref}$
is a very small unit of energy, more appropriate for a single particle. A more
useful reference energy unit for the configuration is
$\mathcal{E}_{\mathrm{ref}}\equiv\rho_{0}E_{\mathrm{ref}}l_{\mathrm{ref}}^{3}/m=(NE_{\mathrm{ref}}/V)l_{\mathrm{ref}}^{3}$
which refers to the energy contained in a reference volume. These units can be
scaled to different physical configurations according to the scaling
invariance properties of SPE, Eq. (48-51). Without any loss of generality we
thus present our results in the units of $\mathcal{E}_{\mathrm{ref}}$,
$\tau_{\mathrm{ref}}$ and $l_{\mathrm{ref}}$. Further details of our
simulation scheme and parameters are given in Appendix B.

Figure 2: (a) An example of a halo forming in a soliton merger experiment in
one of the $M=100M_{\mathrm{ref}}$ realizations, initiated by 10 randomly
distributed solitons of comparable masses. The projection images on each side
of the box plot the volume rendering of the density and the yellow and cyan
isosurfaces respectively depict the regions bounded within one half and one
tenth of the peak density. As is expected in an FDM simulation of such
parameters, a soliton core forms at the halo’s centre. (b) The energy
evolution for the soliton merger simulation shown in (a). The various curves
correspond to, in order of appearance in the legend, the total energy ($E_{\rm
tot}$, black) which remains constant to better than 1.6% accuracy, the total
kinetic ($E_{\rm ke}$, blue) and gravitational potential ($E_{\rm\Phi}$,
purple) energies, the energy associated with the quantum pressure term
($E_{\rm qp}$, green), the incompressible ($E_{\rm ke}^{i}$, red) and
compressible ($E_{\rm ke}^{c}$, yellow) components of the kinetic energy. Note
that $E^{\rm i}_{\rm ke}$, associated with vorticity, is always slightly above
$E^{\rm c}_{\rm ke}$, indicating the importance of quantum vortex contribution
to the energy budget of the halo. (c) Corresponding virialized radial density
profiles, $\rho_{\mathrm{avg}}(r)$, for 30 different realizations with 3
different total masses and averaged over 390 snapshots, presented in scaled
axes (in terms of $\rho_{c}$ and $r_{c}$ respectively). The selected data in
(a) and (b) are presented in black, and its core-halo fit, Eq. (20), is shown
in the dashed grey line. The inset presents the relative error, $[\rho_{\rm
avg}-\rho_{\rm cNFW}]/\rho_{\rm avg}$, showing that the core-halo fit can
sufficiently capture the core profile with error not exceeding 20% at any
point. (d) The core density - radius relation in log-log scale for 30
individual simulations. The error bars here represent the standard deviations
for the measured quantities after $t=\tau_{\mathrm{virial}}$ (see Sec. 2.5 for
more details). The black dashed line plots the core density-radius relation,
Eq. (24), without any free parameters. The highlighted data point refers to
the data presented in (a) and (b) with $E/M_{\rm tot}\approx-70.1$ and
corresponds to the primary dataset analyzed in subplots (a)-(c) and the
remainder of this work.

###  2.3 Simulation protocol

We start with a number of randomly distributed idealized solitons with the
profile given by Eq. (18) of comparable masses (determined by a set of uniform
random numbers), and allow them to coalesce through their gravitational
attraction [Schive et al. (2014b); Schwabe et al. (2016); Mocz et al. (2017);
Chan et al. (2022)] in order to obtain a single virialized halo. This allows
us to investigate both (i) the time-averaged spatial and coherence properties
of the virialized core-halo system, and (ii) the fluctuating post-
virialization dynamics of the system, giving us access to the combined core-
halo oscillations and the associated power spectrum. We simulate the
dimensionless SPE in a $V^{\prime}=L^{\prime 3}=10^{3}$ box discretized by
$288^{3}$ grid points imposing a periodic boundary condition from the
implementation of Fourier transformation. Such simulations are repeated 10
times each for 3 different total masses, $50M_{\mathrm{ref}}$,
$100M_{\mathrm{ref}}$ and $150M_{\mathrm{ref}}$ ($\varrho^{\prime}=0.05,\;0.1$
and $0.15$ respectively) with
$M_{\mathrm{ref}}=\rho_{\mathrm{ref}}(l_{\mathrm{ref}})^{3}$, which is about
$1.26\times 10^{6}M_{\odot}$ for our reference parameters (see section 6). The
$50M_{\mathrm{ref}}$ simulations begin with 5 initial solitons, while 10
initial solitons are chosen for each of the $100M_{\mathrm{ref}}$ and
$150M_{\mathrm{ref}}$ simulations. The initial solitons are randomly
distributed within the simulation box such that for each soliton centre
$|x_{\mathrm{soliton}}|$, $|y_{\mathrm{soliton}}|$,
$|z_{\mathrm{soliton}}|\leq 0.85L^{\prime}$ and under the constraint that the
configuration’s centre of mass is close to the centre of our numerical grid.
The solitons have very similar masses with a small variance
$\sigma_{M}=0.01/\ln(N_{\mathrm{soliton}}$). Due to the slightly different
initial conditions, this gives rise to a range of system energies for each
total mass configuration, and so we have access to a number of closely-
related, yet different, configurations. Specifically, these sets of initial
conditions correspond to the system having total energies that lie in the
region of
$E_{\mathrm{tot}}/M_{\mathrm{tot}}\in(-217,-57)\,\mathcal{E}_{\mathrm{ref}}/M_{\mathrm{ref}}$.
All simulations proceed to a merger of the initially placed masses, eventually
resulting in a virialized halo and exhibiting only a very small energy loss at
the level of 1.6%. This way we can build a complete picture of the physical
state and properties of the core-halo system in a general manner through
confirmation of our findings over a number of similar mass-energy
configurations.

Having analyzed all individual simulations in detail – and ensured that we
observe consistent findings – the main part of this paper focuses primarily on
a single numerical realization shown in Fig. 2 based on
$M=100M_{\mathrm{ref}}$ and $E\approx-7006.4\,\mathcal{E}_{\mathrm{ref}}$,
marked by the arrow in Fig. 2 (d), with our findings applicable to all
conducted simulations. The evolution of 10 initially randomly located solitons
is shown in 3 characteristic images, which also depict the volume rendering
density profiles, projected to different planes of the simulation box along
each perpendicular axis, to give a comprehensive view of the wave-like
evolution of density in time. The depicted 3D isosurfaces correspond to a
density value equal to one half (yellow), and one tenth (cyan) of the peak
value density at a given time, with the former providing an indication of the
soliton core size according to the core width definition, Eq. (18). A movie of
this simulation can be accessed via the supplementary material [SM Movie 1].

To better appreciate the dynamical virialization process, and following [Mocz
et al. (2017)], we first consider the evolution of the different energy
contributions for the presented example in Fig. 2 (b). While the total system
energy remains constant (horizontal straight black line), we see significant
transfer of energy from gravitational (bottom purple curve) to kinetic (top,
blue line) over a typical timescale of $\sim(1-2)\tau_{\rm ref}$, as the
initially isolated idealized self-bound solitons approach each other and
coalesce due to gravity. The generated core is heavily excited after the
merger, exhibiting a range of dynamical density oscillations, during which
heavily reduced, but nonetheless non-negligible, amounts of energy are
periodically exchanged between gravitational and kinetic energy, with most
exchange concerning the quantum pressure term acting to oppose gravitational
attraction. Gradually, such modulations die out, with the different energy
components effectively executing only minor residual noisy oscillations around
a converged final value. We interpret the time of such occurrence as the time
of dynamical equilibration, which we refer to as the virialization time
$\tau_{\mathrm{vir}}$. For the probed configurations this tends to happen over
a time scale $\tau_{\mathrm{vir}}\sim(4-5)\tau_{\rm ref}$, with the particular
depicted example revealing a $\tau_{\mathrm{vir}}\sim(4.5)\tau_{\rm ref}$
(dashed vertical line in Fig. 2(b)). In order to accurately perform the
temporal averaging discussed in subsequent sections, after
$\tau_{\mathrm{vir}}$ we record successive field configurations every
$0.0125\tau_{\mathrm{ref}}$.

Examining in more detail the evolution of the different kinetic energy
components $E_{\mathrm{ke}}^{\mathrm{qp}}$, $E_{\mathrm{ke}}^{i}$ and
$E_{\mathrm{ke}}^{c}$, we note that the quantum pressure energy mostly
dominates over the classical kinetic one, and provides support against
gravitational collapse on small scales. In addition, the incompressible part
of the classical kinetic energy is always higher than the compressible one, by
at least 10% (with such difference being much more striking at the point of
first merging of all input solitons, where it approaches 85%). The final
configuration consists of a single halo containing a quasi-spherical, high-
density soliton core at its centre, surrounded by the granular density
distribution whose average results in the NFW profile evident at large $r$
[Dutta Chowdhury et al. (2021)]. We analyze this further in the sections
below.

###  2.4 Bimodal Characterization of core-halo profile

We now discuss the arising density profiles in more detail, using previously
documented ‘bimodal’ core-halo profiles. This picture is in a manner
reminiscent of condensate-thermal cloud density analysis in confined
inhomogeneous ultracold quantum matter.

A dynamically-equilibrated FDM halo contains at its centre what is often
referred to as a soliton, characterised by a flat density at small $r$, which
can be described by an empirical formula [Schive et al. (2014a); Mocz et al.
(2017); Chan et al. (2018); Chiang et al. (2021)],

$\rho_{\mathrm{soliton}}(r)=\rho_{c}\left[1+\lambda\left(\frac{r}{r_{c}}\right)^{2}\right]^{-8},$
(18)

giving the total mass of the soliton as

$M_{\mathrm{soliton}}=\frac{33\pi^{2}}{1024\lambda^{3/2}}\rho_{c}r_{c}^{3}\approx
11.6811\rho_{c}r_{c}^{3}\,,$ (19)

where $\rho_{c}$ and $r_{c}$ are mutually dependent parameters - see the
following discussion in Sec. 2.5. The central peak density is $\rho_{c}$ and
the length scale $r_{c}$ is defined as the radius at which the density drops
to half its central value (i.e. $\rho_{c}/2)$), giving $\lambda\equiv
2^{1/8}-1\approx 0.091$ [Schive et al. (2014a)].

The soliton profile, Eq. (18), has been implemented to probe the solitonic
cores of FDM halos, showing good agreement and deviating only by a few percent
for $r\lesssim 3r_{c}$ [Mocz et al. (2017); Schive et al. (2014a); Chiang et
al. (2021)], which may be taken as approximately the crossover radius $r_{t}$.
The recent theoretical studies [Chan et al. (2022); Yavetz et al. (2022)] and
observational surveys [Salucci et al. (2007); Di Paolo et al. (2019)] unveil
the possibility that the core profile may only extend to a distance less than
$3r_{c}$. Besides, the mass relation for the core and NFW halo of Schrödinger-
Poison system is not yet fully determined [Chan et al. (2022)].

Even though a lot of attention in the literature has focused on the pure
solitonic core, which is indeed a novel feature of FDM compared to the cuspy
central region of CDM halos, the whole halo is not just composed of the core
but also exhibits the well-known NFW density profile at larger $r$, predicted
in CDM and already verified in many FDM simulations. Thus, one can approach
such a core-outer halo density profile as consisting of the soliton core and
NFW profile with a crossover radius $r_{t}$ [Lin et al. (2018); Bernal et al.
(2018); Zagorac et al. (2022); Chan et al. (2022); Chiang et al. (2021)],
defining an inner and outer halo via

$\rho_{\mathrm{cNFW}}(r)=\Theta(r-r_{t})\rho_{\mathrm{soliton}}(r)+\Theta(r_{t}-r)\rho_{\mathrm{NFW}}(r)$
(20)

where $\Theta(r)$ is the Heaviside step function, and the NFW profile is given
by

$\rho_{\mathrm{NFW}}(r)=\rho_{h}\left(\frac{r}{r_{h}}\right)^{-1}\left[1+\left(\frac{r}{r_{h}}\right)\right]^{-2}$
(21)

featuring the halo length scale $r_{h}$ and an $r^{-3}$ trend at large $r$.
The continuity of $\rho$ at $r_{t}$ limits the number of free parameters,
determining the value of $\rho_{h}$ as

$\rho_{h}=\rho_{c}\left[1+\lambda\left(\frac{r_{t}}{r_{c}}\right)^{2}\right]^{-8}\left(\frac{r_{t}}{r_{c}}\right)\left[1+\left(\frac{r_{t}}{r_{h}}\right)\right]^{2}\,.$
(22)

It is worth noting that the continuity of $\rho$ does not necessarily give a
smooth differentiable density profile; a smooth transition has been tested in
[Chan et al. (2022)], giving a rather small validity region of the core
region, namely, a rather small $r_{t}$ which could be smaller than $r_{c}$.
Besides, the core density can be typically one or even above two orders of
magnitude higher than the density at the crossover radius [Schive et al.
(2014a); Schwabe et al. (2016); Mocz et al. (2017); Lin et al. (2018); Schive
et al. (2020); Li et al. (2021); Dutta Chowdhury et al. (2021); Schwabe &
Niemeyer (2022); Chiang et al. (2021); Zagorac et al. (2022)].

We now investigate how close the core-halo profile given by Eq. (20) matches
the actual radial density distribution in our simulations. To identify the
halo’s centre we note that the positions of the peak density and potential
minimum are nearly identical or only marginally different [Schive et al.
(2020)], and thus we associate the halo center to the position of peak density
without loss of generality. We also observe random-like motions in both of
these points, with the mean displacement being less than a fraction of the
core radius, as already reported in [Li et al. (2021); Schive et al. (2020);
Dutta Chowdhury et al. (2021)]. The actual position of the core’s centre is
actively traced in the simulation, and, after compensating for this random
motion by re-centering the halo, we average the density profiles over 390
snapshots for the presented primary data after $\tau_{\mathrm{vir}}$ with a
time spacing of $0.0125\tau_{\mathrm{ref}}$ to obtain the time averaged
density $\rho_{\mathrm{avg}}(\mathbf{r})$. Such averaged profiles are shown in
Fig. 2 (c) for all of our simulations in scaled axes, $\rho/\rho_{c}-r/r_{c}$.
The radial profile is then obtained by a further angular average
$f(r)=\sum_{r\leq|\mathbf{x}_{i}|<r+\Delta x}f(\mathbf{x}_{i})$ with
$f(\mathbf{r})=\rho(\mathbf{r},t)=|\Psi(\mathbf{r},t)|^{2}$. The best fit to
Eq. (20) is obtained in log-log scale by examining the least value of the
$\chi^{2}$ error among all possible choices of $r_{t}\in[1,N_{x}/2]\Delta x$
with the space resolution $\Delta x$. Such a process uniquely identifies the
spatial soliton extent ($r_{c}$), the outer crossover region ($r_{t}$) and the
NFW profile length scale ($r_{h}$).

Later, in Sec. 4, we will extend such fitting to the time-domain, making both
$r_{c}$ and $r_{t}$ (and thus also $\rho_{c}$) time-dependent. To our best
knowledge the fitting has only been applied to the time-averaged halo profile
so far. The fitting for the dynamical halos gives the time-dependent length
scales and, as we will discuss below, we also find that the oscillation in the
core radius drives an oscillatory power spectrum in our halos.

###  2.5 Verification of the $\rho_{c}-r_{c}$ relation

Inspired by [Chavanis (2011)] and following [Indjin et al. (2022)], the static
soliton energy (ignoring the classical kinetic energy part) can be
analytically computed as

$\begin{array}[]{rl}E_{\mathrm{soliton}}=&E_{\mathrm{qp}}^{\mathrm{soliton}}+E_{\Phi}^{\mathrm{soliton}}\\\
\\\
=&\displaystyle\frac{143}{1024}\frac{\hbar^{2}\pi^{2}\rho_{c}r_{c}}{m^{2}\lambda^{1/2}}-\frac{152}{96115}\frac{\pi^{3}G\rho_{c}^{2}r_{c}^{5}}{\lambda^{5/2}}\\\
\\\
=&\displaystyle\frac{13}{3}\frac{\hbar^{2}\lambda}{m^{2}}\frac{M_{\mathrm{soliton}}}{r_{c}^{2}}-\frac{8506}{2793}\frac{\lambda^{1/2}}{\pi}\frac{GM_{\mathrm{soliton}}^{2}}{r_{c}}~{}\;.\end{array}$
(23)

Considering the energetically stable state for an isolated soliton, free from
fluctuations, and taking $r_{c}$ as the variational parameter, we set
$dE_{\mathrm{tot}}(r_{c})/dr_{c}=0$ to obtain analytical equilibrium relations
connecting the spatial extent of the solitonic core $r_{c}$ to the soliton
mass $M_{\mathrm{soliton}}$ and peak density $\rho_{c}$.

$\begin{array}[]{rl}r_{c}=&\displaystyle\frac{1107}{389}\frac{\hbar^{2}\pi\lambda^{1/2}}{m^{2}GM_{\mathrm{soliton}}}=\left(\frac{11303\lambda^{2}}{128}\frac{\hbar^{2}}{m^{2}G\pi}\right)^{1/4}\rho_{c}^{-1/4}\\\
\\\ \approx&\displaystyle 23.686\left(\frac{2.5\times
10^{-23}\textrm{eV}/\textrm{c}^{2}}{m}\right)^{1/2}\left(\frac{10^{3}M_{\odot}\textrm{kpc}^{-3}}{\rho_{c}}\right)^{1/4}\textrm{
kpc}.\end{array}$ (24)

Such relations are in agreement with many independent numerical simulations
[Schive et al. (2014a); Mocz et al. (2017)], and also satisfy the SPE’s
scaling invariance [Seidel & Suen (1990); Guzmán & Ureña-López (2004)] (see
also Appendix A). Recalling Eq. (19), the mass of a soliton with radius
$r_{c}$ can be written as

$\begin{array}[]{rl}M_{\mathrm{soliton}}=&\displaystyle\frac{1107}{389}\frac{\sqrt{\lambda}\pi\hbar^{2}}{m^{2}G}r_{c}^{-1}\\\
\\\ \approx&\displaystyle 3.68\times 10^{7}\times\left(\frac{2.5\times
10^{-23}\textrm{eV}/\textrm{c}^{2}}{m}\right)^{2}\left(\frac{r_{c}}{\textrm{kpc}}\right)^{-1}M_{\odot}\end{array}$
(25)

which is compatible with previous numerical studies.

Having obtained $r_{c}$ from the fit, we examine the core-radius relation for
the soliton, Eq. (24). The core density $\rho_{c}$ after $\tau_{\mathrm{vir}}$
oscillates by about 50% around the time-averaged value which we take as the
value of $\rho_{c}$. The core radius as function of $\rho_{c}$ obtained from
our simulations is plotted in Fig. 2 (d) and the trend matches the derived
relation Eq. (24) without any free parameters.

The nice agreement with the expected $\rho_{c}-r_{c}$ trend indicates that the
empirical core profile is a suitable one for the core region of the FDM halo.
There is a small deviation for $r_{c}<0.08l_{\mathrm{ref}}\lessapprox 3\Delta
x$ in $M=150M_{\mathrm{ref}}$ simulations, suggesting a soliton needs to be at
least 3 grid points across to be well resolved. Nevertheless, the overall
trend is robust. It is worth emphasising that in our simulations multiple
excitations are generated during the merger history of the final halo, and it
is consequently hard to define a precise soliton mass. However, the core
density and radius relations above can still be taken as useful and relatively
accurate guides.

##  3 Coherent and Incoherent Features of the core-halo system

Having discussed our core-halo system generation protocol and fitting schemes,
we next proceed to discuss key features of such virialized systems in terms of
spatial coherence, power spectra and dynamics.

###  3.1 Condensation and spatial correlation functions

The phenomenon of Bose-Einstein condensation has been extensively studied both
theoretically and experimentally in laboratories using ultracold atomic gases
(among other systems). In most configurations, the system evolves towards
(potentially through highly non-equilibrium states, e.g. as those occurring
during the condensation process [Berloff & Svistunov (2002); Proukakis et al.
(2017); Beugnon & Navon (2017); Liu et al. (2018)], and eventually relaxes to,
a state of (at least local) thermal equilibrium, described by the associated
Bose-Einstein distribution with a fixed temperature and chemical potential, a
situation applicable in many experimental setups. However, this assumption
cannot be made in the case of a collection of self-gravitating bosons such as
those that make up cosmological halos in FDM, because such systems are neither
thermalized nor have a well-defined chemical potential in the expanding
Universe. Nonetheless, condensation will still occur in this circumstance and
it is expected that, like with any collection of bosons, ultralight dark
matter will Bose-Einstein condense when the particle number in a volume
characterized by the de Broglie wavelength scale is much larger than order
unity [Sikivie & Yang (2009); Banik & Sikivie (2017)]. In cosmological
settings, this criterion has already been associated with FDM exhibiting
distinctive quantum, or “wavy”, behaviour [Schive et al. (2014a); Berezhiani &
Khoury (2015); Hui (2021)].

A standard criterion delineating the emergence of Bose-Einstein condensation
in an equilibrium setting is the value of the dimensionless phase-space
density ${\cal D}$, with the textbook critical value for a non-interacting,
homogeneous, 3D condensate being ${\cal D}>\zeta(3/2)\approx 2.612$ [Huang
(2000); Pitaevskii & Stringari (2003); Pethick & Smith (2008)]. Since in the
FDM cosmological context the stationary-states are not thermal (unlike many
laboratory condensate systems), one is restricted to using the more general
relation $\lambda_{\rm dB}=h/p$ as opposed to the thermal de Broglie
wavelength $\lambda_{\rm dB}=\lambda_{\rm T}=\sqrt{2\pi\hbar^{2}/mk_{B}T}$
used in ultracold superfluid studies. Obviously, both density and velocity
are, in general, functions of ${\bf r}$ in the halo, so here $p({\bf
r})=mv({\bf r})$. After averaging over both azimuthal directions (to remove
angular dependence) and time (during the post-virialized evolution), and under
the ergodicity assumption within which temporal averages in steady-state
approach the ensemble average [Blakie et al. (2008); Liu et al. (2018, 2020)],
denoted by $\langle\cdots\rangle$ – and taking for the sake of our present
arguments the simplest known limit of a 3D non-interacting homogeneous
equilibrium condensate – we thus obtain the radially-dependent condensation
criterion

${\cal
D}(r)=\frac{\rho(r)}{m}\lambda_{\mathrm{dB}}^{3}=\frac{h^{3}}{m^{4}}\frac{\langle\rho(r)\rangle}{\langle|{v}(r)|\rangle^{3}}\gtrsim\zeta(3/2)\approx
2.612\,.$ (26)

Furthermore, in terms of the dimensionless density $\rho^{\prime}(r)$ and
velocity ${v}^{\prime}(r)$ provided by the numerical simulation we can write

${\cal
D}(r)=\frac{h^{3}\rho_{\mathrm{ref}}}{m^{4}v_{\mathrm{ref}}^{3}}\frac{\langle\rho^{\prime}(r)\rangle}{\langle|{v}^{\prime}(r)|\rangle^{3}}={\cal
D}_{\mathrm{ref}}\frac{\langle\rho^{\prime}(r)\rangle}{\langle|{v}^{\prime}(r)|\rangle^{3}}~{},$
(27)

where

${\cal
D}_{\mathrm{ref}}=\displaystyle\frac{h^{3}\rho_{\mathrm{ref}}}{m^{4}v_{\mathrm{ref}}^{3}}~{}\;.$
(28)

As we discuss later, it is useful to consider the quantity ${\cal D}(r)/{\cal
D}_{\rm ref}$ which has a universal character that can be applied for
different boson masses. This exhibits a significant drop of many orders of
magnitude with increasing radius, with its value dropping to $\sim O(1)$ at
the crossover radius $r_{t}$, beyond which point its rate of decrease becomes
noticeably reduced. Regardless of the high occupation numbers for small boson
masses across our entire simulation grid (see section 6), there are in fact
qualitative differences in the state of the field between the core and the
outer halo. Furthermore, the term "Bose-Einstein condensate" typically refers
to a state which not only exhibits a high number of particles per de Broglie
volume but also a suppression of phase and density fluctuations [Huang (2000);
Pitaevskii & Stringari (2003); Pethick & Smith (2008)]. We therefore introduce
below commonly used, and appropriately normalized, spatial correlation
functions which characterize the state of the system with respect to such
variables. Although, as shown later, the core-halo system oscillates
dynamically, once it has virialized such dynamics are perturbations around a
steady-state. This makes our time-averaging protocol, used for extracting
smoother spatial correlation functions, appropriate for the inner regions.
Nonetheless, the surrounding region features intense fluctuations in both the
spatial and temporal domain, and we comment further on that below.

The textbook definition of Bose-Einstein condensation relies on the emergence
of Off-Diagonal Long-Range Order (ODLRO) [Pitaevskii & Stringari (2003)], i.e.
the fact that the equal-time one-body density matrix does not vanish at large
distances (as $(|{\bf r}-{\bf r^{\prime}}|)\rightarrow\infty$) but remains
finite [Penrose & Onsager (1956)]. This is best quantified in the normalized
first order spatial correlation function444Note that the general
characterization in terms of appropriate field operators, is here replaced by
corresponding expressions in terms of classical fields.

$g_{1}(\mathbf{r})=g_{1}(\mathbf{0},\mathbf{r})=\frac{\langle\Psi^{\ast}(\mathbf{0})\Psi(\mathbf{r})\rangle}{\sqrt{\langle|\Psi(\mathbf{0})|^{2}\rangle\langle|\Psi(\mathbf{r})|^{2}\rangle}}$
(29)

quantifying the change of _phase coherence_ from the (numerically identified)
centre of the solitonic core (${\bf r}={\bf 0}$) to a general radial position
${\bf r}$ away from it. For a Gaussian fluctuating field, full coherence and
complete incoherence are characterized by $g_{1}\approx 1$ and $0$
respectively (in the limit ${\bf r}\rightarrow\infty$). We also consider the
equal-position second-order correlation function

$g_{2}(\mathbf{r})=g_{2}(\mathbf{r},\mathbf{r},\mathbf{r},\mathbf{r})=\frac{\langle|\Psi(\mathbf{r})|^{4}\rangle}{\langle|\Psi(\mathbf{r})|^{2}\rangle^{2}}=\frac{\langle[\rho(\mathbf{r})]^{2}\rangle}{\langle\rho(\mathbf{r})\rangle^{2}},$
(30)

which characterizes _local density fluctuations_ in the system, charting the
spatial distribution of the degree of local coherence. Completely suppressed
local density fluctuations, consistent with a pure condensate, imply
$g_{2}=1$, while $g_{2}=2$ for an incoherent (chaotic) Gaussian field. In what
follows we present the second-order correlation in the form of $2-g_{2}({\bf
r})$ in order to provide a consistent visualization for the two coherence
measures.

Figure 3:  Characterization of core-halo coherence. (a)-(b) Different
condensation and coherence measures for the simulation data of Fig 2(a),
highlighting the crossover from coherent to incoherent field configurations as
a function of radius, scaled to numerically-obtained radius $r_{c}$ of the
solitonic core. All radial cases show angular- and time-averaged results,
taken over $390$ snapshots during the
$t\geq\tau_{\mathrm{vir}}=4.5\tau_{\mathrm{ref}}$ evolution period in our
reference simulation. For comparison, both subplots also show the averaged
density profile ($\rho_{\rm avg}$, black line), and locations of fitting
parameters $r_{c}$ (dashed blue vertical line) and $r_{t}$ (dotted red
vertical line) of Fig 2(c). (a) Penrose-Onsager (PO), or condensate, density
($\rho_{\rm PO}=|\Psi_{\rm PO}|^{2}$ ) (green) and numerical core-halo fit of
$\rho_{\rm avg}$ (dashed grey line) as a function of scaled radius $r/r_{c}$.
Inset plots scaled densities $(\rho/\rho_{0})$ as a function of $x$ for
$y=z=0$, depicting the core-halo fit (dashed grey line), and the
decompositions of its exact, numerically-obtained, density, $\rho_{\rm avg}$,
into the condensate PO mode (central peak, solid green) and the residual lower
density $(\rho_{\rm avg}-\rho_{\rm PO})$ (broad red line), clearly
highlighting our decomposition into coherent (green) and incoherent (red)
contributions to the total density, in qualitative analogy to condensed and
thermal features in harmonically-trapped ultracold atomic superfluids. The
grey band shows the standard deviation of $\rho(r,t)$ over times beyond
$\tau_{\mathrm{vir}}$, and the green band depicts the standard deviation of
the PO density at a given $r$. (b) The radial profiles for the coherence
measures: dimensionless reduced phase space density ${\cal D}(r)$, scaled to
the reference value ${\cal D}_{\rm ref}$ (light blue, left axis); $g_{1}(r)$
(green), $2-g_{2}(r)$ (purple) and condensate fraction
$\rho_{{}_{\mathrm{PO}}}(r)/\rho_{\mathrm{avg}}(r)$. A horizontal reference
line is drawn at ${\cal D}(r)=2.612$ which is the critical phase space density
for condensation for a 3D homogeneous non-interacting Bose gas in the
thermodynamic limit. Purple band labels the variance of the time-averaged
$(2-g_{2}(r))$, with such band becoming significant at distances $r>r_{t}$
when $g_{2}$ approaches the value of 2, valid for a Gaussian fluctuating
field. (c) Planar profiles for $z=0$ of (i) the time-averaged density
$\rho_{\rm avg}$, (ii) the PO density corresponding to the largest eigenvalue,
(iii)-(iv) the next two dominant eigenmodes of the single-particle density
matrix (with eignevalues $N_{2}\approx N_{3}\approx 0.1N_{\rm PO}$, obtained
by temporal averaging over entire post-virialized evolution. (d) Corresponding
$z=0$ planar profiles of $g_{2}$ averaged over two different timescales and
plotted on two different colourbars to highlight [top] range of values $1\leq
g_{2}\lesssim 2$, and [bottom] regions of significantly suppressed density
fluctuations, values $g_{2}\lesssim 1.2$.: Left plots [(i), (iii)] correspond
to an averaging timescale $t\in[8,8.5]\tau_{\mathrm{ref}}$, broadly consistent
with timescale for noticeable vortex evolution, highlighting local regions of
significantly suppressed density fluctuations; Right plots [(ii), (iv)]
indicate corresponding results when averaged over the entire post-
virialization evolution (i.e. same averaging timescale as used in obtaining
plots (a)-(c)).

Having specified our distinct measures of coherence, we now determine them
numerically to characterize the state of our representative virialized core-
halo system shown in Fig. 3(a). Such analysis, and its confirmation across
different $(E,M)$ combinations [see Appendix C] provide clear evidence of the
properties of the obtained spatially-distinct regions. The averaged radial
dependence of the aforementioned quantities is shown in Fig. 3 (b) as a
function of $r/r_{c}$, where $r_{c}$ is extracted from the fit to
$\rho_{\mathrm{avg}}(\mathbf{r})\equiv\langle\rho(\mathbf{r})\rangle$, as
previously discussed.

Firstly, we note that although the phase space density ${\cal D}(r)$ can be
very high for typical FDM boson masses, the reduced phase space density ${\cal
D}(r)/{\cal D}_{\rm ref}\sim O(10^{5})$ in the solitonic core centre and
decreases rather rapidly when moving radially outwards, reaching the critical
value $\sim 2.612$ and exhibiting a clear change in its slope around the edge
of the crossover region, $r_{t}$. It then drops below unity in the outer parts
of the halo, not far beyond $r_{t}$.

Looking at the spatial coherence measures, $g_{1}$ and $(2-g_{2})$, we note
that they are almost perfectly flat and $\approx 1$ within the entire inner
solitonic core region ($r<r_{c}$), signifying the near-complete suppression of
both density and phase fluctuations, consistent with the solitonic core being
a pure condensate. Within the region $r_{c}\lesssim r\lesssim r_{t}$ both
coherence measures decrease noticeably, and, for $r\gtrsim r_{t}$, the field
becomes incoherent ($g_{1}\approx 2-g_{2}\approx 0$). Our numerical results
point towards the potential existence of regions locally exhibiting values of
$g_{2}>2$ (marginally), with large error bars in its determination, shifting
concrete conclusions beyond our present numerical investigations: while such
values would arise in the case of non-Gaussian fluctuations, this could also
here be evidence that the averaging time is not enough to completely wash out
such overdensities in the slower-evolving outer parts of the halo (see also
Appendix C). Either way, the combination of angular and extensive temporal
averaging obscure information about the existence of smaller-scale, localized,
partly coherent density patches of variable (non-radially-symmetric) shapes in
the outer halo, which are further discussed below – see also Fig. 3(d).

An alternative way to extract information about the presence of the condensate
and relate to the mass density is by directly identifying those parts of the
density contributions which exhibit both full phase, and density coherence, in
direct analogy to the previously identified correlation functions $g_{1}$ and
$g_{2}$. In particular, one can directly extract from the numerics of the
entire field the corresponding condensate mode. This can be done according to
the Penrose-Onsager criterion [Penrose & Onsager (1956); Blakie & Davis
(2005); Blakie et al. (2008)], which identifies the condensate mode
$\Psi_{\mathrm{PO}}(\mathbf{r})$ as that mode of the single-particle density
matrix555We remind the reader that our wavefunction $\Psi$ normalizes to the
system mass, rather than particle number.

$\varrho(\mathbf{r},\mathbf{r}^{\prime})\equiv\frac{1}{m}\langle\Psi^{\ast}(\mathbf{r})\Psi(\mathbf{r}^{\prime})\rangle$
(31)

which has the largest eigenvalue.

Mathematically, the dominant eigenvalues and corresponding eigenfunctions are
identified from the solution of

$\int{\rm d}{\bf r^{\prime}}\varrho({\bf r},{\bf
r^{\prime}})\psi_{n}(\mathbf{r^{\prime}})=N_{n}\psi_{n}(\mathbf{r})$ (32)

where the index $n$ labels the $n$-th eigenfunction $\psi_{n}$ of $\varrho$.
Such identification only makes interpretative sense if one, or few,
eigenvalues are dominant: in the former case the system is said to exhibit a
well-defined Bose-Einstein condensate, with $\psi_{0}$ the condensate mode,
and $N_{\mathrm{PO}}=N_{0}$ the condensate particle number. A system with no
(significantly) prevailing eigenvalue, but exhibiting co-existence of few
macroscopically occupied eigenmodes is characteristic of a quasi-condensate
state, in which there is no overall phase coherence (except locally) due to
the competition between dominant eigenmodes (exhibiting different random
phases).

As the exact numerical diagonalization of such a large density matrix is
numerically impractical, we instead evaluate the largest eigenvalue of the
dynamical steady-state system (post-virialization) by replacing the ensemble
average in Eq. (31) by a time average: Strictly speaking this is true for an
ergodic system, and this is a standard technique in the field of ultracold
quantum matter666In ultracold atomic systems, the condensate ultimately
reaches an equilibrium state over relevant probing timescales, i.e. before any
losses, not typically included, become important and lead to the ultimate
solidification of the quantum gas in such systems. [Blakie et al. (2008); Liu
et al. (2018, 2020); Proukakis et al. (2013)]. The temporal duration of such
averaging should be chosen to be long enough to wash out random fluctuations,
but short enough not to obscure any key underlying physical dynamics being
probed (see, e.g. [Liu et al. (2018)]). Post virialization, the solitonic core
acquires a near-equilibrium profile (modulo small oscillations - see below),
and so the exact averaging duration beyond $\tau_{\rm vir}$ has practically no
impact on its properties, provided enough samples are averaged over. However,
more attention is needed when addressing the halo regions, where the quantum
vortices and granular structures propagate dynamically; in this case our
findings become sensitive on whether the averaging timescale is shorter than,
or longer than the typical evolution timescales of such structures, as
discussed below. In order to obtain most averaged results presented in this
work, we typically probe the entire post-virialization timescale of
$4.875\tau_{\rm ref}$. However, when probing short-time local density
correlations [left panels in Fig. 3(d)], we also consider a shorter averaging
time $t_{\mathrm{aver}}=0.5\tau_{\mathrm{ref}}$. In practice we average over a
certain number of snapshots taken $0.0125\tau_{\mathrm{ref}}$ apart, with
respective numbers being 390 (long-time) and 40 (short-time) snaphsots \- see
Appendix B for further details.

Our simulations clearly reveal the existence of one macroscopically-occupied
eigenmode – found to have an eigenvalue which is $\approx 10$ times larger
than that of any other modes – followed by a range of few approximately
equally important eigenmodes (e.g. for $n=1$ to $5$ we found that
$N_{n}/N_{0}\approx 0.1$ – see Appendix B.2 for details). The three first
eigenmodes are shown in Fig. 3(c)(ii)-(iv), alongside the corresponding time-
averaged density [Fig. 3(c)(i)]. The largest eigenvalue mode, or Penrose-
Onsager (PO) mode [Fig. 3(c)(ii)] clearly coincides with the central soliton
(including the crossover region); contrary to this, as expected, the next two
eigenmodes have $N_{2}\approx N_{3}\approx 0.1N_{0}$ and a much broader radial
extent. The radially-averaged PO mode is shown by the green line in Fig. 3(a),
and corresponds to that part of the density with suppressed density and phase
fluctuations.

At this point, it is perhaps instructive to draw a closer analogy between our
system and ultracold atomic gases. In the latter system, under a harmonic
confining potential, one finds a dominant central coherent profile, i.e. a
condensate (in the shape of an inverted parabola) embedded within an
incoherent Gaussian thermal cloud distribution (which is used to identify the
system temperature), thus giving rise to a bimodal configuration777Note that
despite the presence of (typically) repulsive interactions, careful modelling
in terms of a condensate and a thermal cloud (see, e.g. Proukakis & Jackson
(2008)) reveals that a small (approximately constant density) incoherent
component persists even within the condensate profile.. This is highly
analogous to the spatial distribution obtained in our present cosmological
simulations. To make this ‘two subsystems’ analogy clearer, we impose here a
separation of the numerically-obtained total system density into a condensed
(PO mode) and an uncondensed component, obtained by subtracting the PO mode
from the average density; these are shown respectively by the green and red
lines in the inset to Fig. 3(a). We stress that such determination is based
only on the numerically-obtained density and does not in any way invoke the
soliton and NFW profiles. To highlight the excellent agreement of our PO-based
analysis for the density components to the bimodal fit based on Eq. (20),
which uses a solitonic core profile in the inner region ($r<r_{t}$) and an NFW
profile for $r>r_{t}$, we overlay in the inset the total density obtained via
Eq. (20) by a grey dashed line. Such near-perfect agreement implies that the
well-studied solitonic core of FDM is a _true inhomogeneous condensate_
embedded within a dynamical, and on average incoherent, NFW halo. The
condensate fraction is defined by
$f_{{}_{\mathrm{PO}}}(r)\equiv\rho_{{}_{\mathrm{PO}}}(r)/\rho(r)$, showing
within $r<r_{c}$ the halo density is nearly a pure condensate with
$f_{{}_{\mathrm{PO}}}(r)\sim 1$ in agreement with the above coherence
measurements.

Having established that the core is a proper BEC with suppressed density and
phase fluctuations, the question naturally arises to what extent any such
features can be discerned in the outer halo. A useful quantity characterising
an ‘intermediate’ state exhibiting suppressed density fluctuations but no true
ODLRO, which is actively discussed in the ultracold atomic community, is the
quasi-condensate density. This is closely related to $g_{2}$ and is defined
according to [Prokof’ev et al. (2001); Berloff & Svistunov (2002); Proukakis
(2006); Blakie et al. (2008); Cockburn et al. (2011); Liu et al. (2020)]

$\rho_{\mathrm{qc}}(\mathbf{r})=\sqrt{2\langle|\Psi(\mathbf{r})|^{2}\rangle^{2}-\langle|\Psi(\mathbf{r})|^{4}\rangle}=\rho_{\mathrm{avg}}(\mathbf{r})\sqrt{2-g_{2}(\mathbf{r})}\;.$
(33)

As we have already seen in Fig. 3(b), the value of $g_{2}$ extracted as an
average over the entire post-virialized evolution can still (marginally)
exceed the value of 2, while exhibiting large fluctuations – a finding
consistent across all our simulations (details in Appendix C). As a result, a
radial quasi-condensate density can only be defined within our numerics up to
$r\lesssim r_{t}$, where it practically coincides with the PO mode. Features
beyond such point are obscured by the temporal and angular averaging used to
obtain the value of $g_{2}$ in Fig. 3 (b). Therefore, it is very instructive
to consider the full spatial dependence of $g_{2}({\bf r})$, which will allow
us to identify spatially- and temporally-localized features relating to
structures in the halo often referred to in the literature as ‘granules’ (e.g.
[Schive et al. (2014a)]).

To visualize this, we look at slices of $g_{2}$ in the $x$-$y$ plane, taken at
$z=0$ [Fig. 3(d)] and consider two different time averages (left/right
columns) to examine any dynamically evolving features which are washed away by
long-time averaging. Moreover, in Fig. 3(d) we plot the same information on
two different colourbars [top and bottom] to more clearly highlight the key
arising features. Looking first at the same (long) time averaging used to
extract $\rho_{\rm avg}(r)$, $\rho_{\rm PO}(r)$ and $g_{2}(r)$ in Fig.
3(a)-(c), we see from Fig. 3(d)(ii) the existence of the prominent solitonic
core with completely suppressed density fluctuations $g_{2}\approx 1$ (yellow
core). This is surrounded by regions of intense density fluctuations
($g_{2}\approx 2$, dark red regions) which dominate the picture, with only a
very small number of regions exhibiting suppressed density fluctuations
($g_{2}\lesssim 1.5$, light blue regions). Such rich structure explains why
the quasi-condensate cannot be numerically identified (with our current
averaging scheme at least) in the radial direction beyond $r_{t}$. Over the
extended time-averaging, there are in fact no regions, other than the extended
solitonic core, with $g_{2}<1.2$ [Fig. 3(d)(iv)].

However, such picture changes drastically when considering shorter time
averages of 0.5$\tau_{\rm ref}$ [Fig. 3(d), left column]: In this case, we see
clear evidence of widespread regions of locally (partially) suppressed density
fluctuations in the halo [yellow/green regions in Fig. 3(d)(i)]. Although
density fluctuations in such regions (where the radially-integrated density is
well-described by the NFW profile) are evidently locally suppressed
($g_{2}<1.2$), they are nonetheless not as suppressed as in the core itself
($g_{2}\approx 1$) [Fig. 3(d)(iii)]. We have further found that, at a given
time within the probed interval, the phase over a region of suppressed density
fluctuations remains largely constant. This is a strong indication of the
local existence of quasi-condensate lumps in the halo, with such lumps
destroyed on longer timescales through the chaotic passage of quantum
vortices. We stress that consolidation of the above quasi-condensate
discussion would require extensive further analysis, which we defer to future
work.

Figure 4:  (a) The energy distribution as a function of radius for the quantum
pressure energy, $\varepsilon_{\mathrm{qp}}$ (green), and classical kinetic
energy $\varepsilon_{\mathrm{ke}}$ (dashed purple), decomposed into
incompressible $\varepsilon_{\mathrm{ke}}^{i}$ (red) and compressible
$\varepsilon_{\mathrm{ke}}^{c}$ (orange) components with $\rho_{\mathrm{avg}}$
as a reference for the data presented in Fig. 3. (b) The radial energy density
profile for the incompressible and compressible parts of the kinetic energy
per density, a measure of the velocity. The bands represent the uncertainty
due to the variance in time of the radial profiles for
$t\geq\tau_{\mathrm{vir}}$. Dotted green lines show the corresponding
prediction for $\varepsilon_{\mathrm{qp}}$ based on (34)

###  3.2 Energy density distribution

Next, we consider the spatially-dependent contributions of the different
(decomposed) energy densities. According to Eq. (18), the quantum energy
density of an isolated soliton core can be calculated to be

$\varepsilon_{\mathrm{qp}}^{\mathrm{soliton}}(r)=32\frac{\hbar^{2}\lambda^{2}r^{2}}{m^{2}r_{c}^{4}}\left[1+\lambda\left(\frac{r}{r_{c}}\right)^{2}\right]^{-2}\rho_{\mathrm{soliton}}(r)$
(34)

suggesting a quadratic increase of the energy associated to the quantum
pressure term in the flat core region.

In Fig. 4 (a), we present the distribution of the classical kinetic energy,
broken down into incompressible and compressible parts, as well as the energy
associated with the quantum pressure, as a function of radius for the example
halo described in this work. One can readily see the clear dominance of
quantum pressure in the core region while the classical kinetic components are
at least an order of magnitude smaller (but not zero). The trend matches well
with the estimation from the empirical soliton form up to $r\approx r_{t}$
where the density transits to the NFW profile for the outer halo region.
However, the empirical formula tends to underestimate the quantum pressure
contribution, indicating the slightly more dynamical nature of the soliton
embedded in a halo as oppossed to an isolated one. In the outer halo, we find
that both quantum pressure and classical kinetic energies are comparable with
${\varepsilon_{\mathrm{qp}}/{\varepsilon}_{\mathrm{ke}}}\approx 1$ [Mocz et
al. (2017)]. As can be clearly seen in Fig. 4 (a) and (b), in the outer halo
region the incompressible energy contributes slightly more than the classical
kinetic energy.

By investigating the energy distribution divided by the mass density,
$\rho^{-1}\varepsilon$ in Fig. 4 (b), which is essentially the modulus square
of effective velocity amplitude per particle, we also find that, as expected,
the quantum pressure remains dominant in the core region while the classical
kinetic velocities are much smaller. The incompressible part of the classical
kinetic energy component is smaller than the compressible part for $r<r_{c}$;
notably $\rho^{-1}\varepsilon^{i}$ is even less than order unity. As there are
no vortices within the core region (see below), its small but non-zero value
suggests that the whole core participates in a slow rotational motion relative
to points outside it. The compressible part can be larger and is due to the
non-vortical motion of the core itself or excitations propagating through it.
Visual inspection of the core region in our simulations supports this
interpretation: as mentioned above, the core exhibits a random-like walk
confined in a small region around the centre and is itself pulsating to a
small degree on top of the underlying solitonic density profile. For $r>r_{t}$
the classical kinetic part becomes comparable to the quantum pressure one,
apparently satisfying equipartition.

Figure 5: (a) Two dynamical radial density profiles with corresponding core-
halo fits for $t=5.8875\tau_{\mathrm{ref}}$ and $7.1375\tau_{\mathrm{ref}}$,
selected from the temporal evolution of $r_{c}(t)$ and $r_{t}(t)$ shown in
(c), such that they correspond to the extrema of the $r_{t}$ oscillation. The
blue and red vertical dashed and dotted lines mark the length scales of
$r_{c}(t)$ and $r_{t}(t)$ with the circle and triangle marking those points on
the density profile. The insets show the relative error,
$[\rho(r,t)-\rho_{\mathrm{fit}}(r)]/\rho(r,t)$, with the core-halo fit
sufficiently capturing the overall profile particularly for $r<r_{c}(t)$ with
less than 20% error. (b) The corresponding power spectra for the
configurations in (a) show a displacement of the peak location of the power
spectrum, marked by the vertical dashed lines with hollow circles for the two
times $t\geq\tau_{\mathrm{vir}}$. (c) The temporal evolution of $r_{c}(t)$ and
$r_{t}(t)$ extracted from the dynamical density profiles - the times shown in
(a) are marked by the black dashed-dotted lines and correspond to two extreme
variations. The shaded background shows the subtle density deformation in
time. (d) The evolution of $r_{c}(t)$ and the peak momentum of the power
spectrum $k_{\mathrm{max}}^{\Delta^{2}}(t)$ for the same time windows in (c).
A clear anti-correlation for these two quantities can be observed.

##  4 Core and power spectrum Oscillations

As we have already seen in Sec. 2.5, the (static) central core density and the
(static) soliton core width $r_{c}$ are related. The oscillation of the core
density and its frequency have been reported in previous studies [Guzmán &
Ureña-López (2004); Veltmaat et al. (2018); Li et al. (2021); Schive et al.
(2020); Chiang et al. (2021); Schwabe & Niemeyer (2022)]. We find the dominant
oscillation frequency, inferred by the Fourier spectra of both $\rho_{c}(t)$
and $r_{c}(t)$, to be

$\omega\propto\rho_{c}^{1/2}$ (35)

agreeing with the literature and an analytical prediction based on the
Gaussian density packet approach [Chavanis (2011)]. These density oscillations
may be driven by interactions between the soliton core and the outer halo [Li
et al. (2021); Zagorac et al. (2022)], resulting in fluctuating anisotropy of
the central density [Dutta Chowdhury et al. (2021)]. This relation can also be
seen from the scaling invariant relation, Eq. (48).

In addition, the Fourier spectra in both $r_{c}(t)$ and $\rho_{c}(t)$ (and
also the peak wavenumber of the overdensity power spectrum
$k_{\mathrm{max}}(t)$ \- see below) exhibit multiple peak-frequency features
[Li et al. (2021); Dutta Chowdhury et al. (2021)]. This suggests that there
are multiple excitations created in the merger simulation with the
nonlinearity in the equations of motion resulting in mode mixing [Morgan et
al. (1998); Li et al. (2021); Yavetz et al. (2022); Zagorac et al. (2022)].
All these lead to fluctuations in the precise core oscillation frequencies, a
situation that could also be generic for actual halos, coupled to the baryonic
matter. We will report more details about the topic of oscillations elsewhere
[Indjin et al. (2022)].

Going beyond the already studied oscillation of the central density, we now
discuss how the core-halo fitting formula (20) does not just capture the time-
averaged density profile but can also model the dynamical one fairly well. In
Fig 5 (a) (i) and (ii) we plot two radial density profiles with their core-
halo fits at two different times, chosen at the extrema of the post-virialized
oscillations to demonstrate the broad applicability of our fits. The core-halo
fit successfully captures the overall profile, particularly around the "neck"
at $r\approx r_{t}$, as illustrated in Fig. 5 (a) (i). Very occasionally in
our simulation the extracted value of $r_{t}$ approaches that of $r_{c}$, with
such an instance shown in Fig. (a) (ii) in which the typical head-neck-
shoulder profile is not clear. In Fig. 5 (c), we show an example of the
evolution of $r_{c}$ and $r_{t}$ in time with the two instances depicted in
Fig. (a) (i) and (ii) marked by the two vertical dashed-dotted lines. Clearly,
$r_{c}$ exhibits far less fluctuations than $r_{t}$.

We also examine the dimensionless power spectrum of the overdensity in the
field configuration, written as

$\begin{array}[]{rl}\Delta^{2}(k_{r},t)=&\displaystyle\frac{k_{r}^{3}}{(2\pi^{2})4\pi
k_{r}^{2}}\int d\mathbf{\Omega}_{k}P(\mathbf{k},t)\\\ \\\
\approx&\displaystyle\frac{k_{r}^{3}}{2\pi^{2}\mathcal{N}_{k_{r}}}\sum_{k_{r}\leq|\mathbf{k}|<k_{r}+\Delta
k}P(\mathbf{k},t)\,.\end{array}$ (36)

Here $P(\mathbf{k},t)=\left|\tilde{\eta}(\mathbf{k},t)\right|^{2}$ is the
power with the Fourier transformation of the scaled density
$\eta(\mathbf{r},t)=[\rho(\mathbf{r},t)-\rho_{0}]/\rho_{0}$ with respect to
the spatially averaged density $\rho_{0}$ which is a constant in time;
$\mathbf{\Omega}_{k}$ the solid angle in the momentum space and
$\mathcal{N}_{k_{r}}$ is the number of momentum grids within
$k_{r}\leq|\mathbf{k}|<k_{r}+\Delta k$ and $\Delta k_{r}=2\pi/L\approx
0.63l_{\mathrm{ref}}^{-1}$.

The power spectrum is time dependent and always exhibits a clear peak, the
height and position of which oscillates. Two instances of the power spectrum
are shown in Fig. 5 (b), corresponding to the profiles shown in (a-i) and
(a-ii) of the same figure and also indicate the range over which the power
spectrum varies in the virialized halo of our chosen example - at any given
instant the power spectrum lies between these two curves. Note that when the
peak is located at the maximum wavenumber the core and the crossover region
are at their smallest extent, while the minimum peak wavenumber corresponds to
the most extended phase of the core during its oscillation. This
correspondence between the peak wavenumber of the power spectrum and the
oscillation phase of the core is evident in Fig. 5 (d) where we plot both
$r_{c}(t)$ and the peak wavenumber of the power spectrum. It is clear from
that figure that the two quantities are almost perfectly anti-correlated. The
natural question to ask is whether this power spectrum oscillation is driven
by the core alone or whether the granules oscillate too. The dynamical state
of the core with its oscillatory and random walk motion has prevented us from
subtracting it cleanly from the halo to study the granule spectrum alone, but
preliminary results, to be reported elsewhere [Indjin et al. (2022)], indicate
that at least most of the oscillatory contribution to the power spectrum comes
from the core and the granule spectrum is much less time dependent. Thus, the
core’s time dependent size shows up in the power spectrum, at least for the
halos in the mass and energy regime formed in our simulations.

##  5 Vortices and Granules

Figure 6:  (a) A composite visualization of the halo at
$t=8\tau_{\mathrm{ref}}$, with peak density at the centre of the computational
box. The red filaments are the isosurfaces of
$|\mathbf{v}_{ke}^{i}|^{2}=800l_{\mathrm{ref}}^{2}/\tau_{\mathrm{ref}}^{2}\approx(\hbar^{2}/2m^{2}\Delta
x^{2})$, representing vortex structures in closed loop topologies. The slices
on the right and bottom of the sliced sphere show the density profile
$\rho(\mathbf{r})/\rho_{0}$ on them with colours ranging from blue (low
density) to yellow (high density). The phase profile of the wavefunction is
presented on the left slice with a pale pink-purple-blue-green-yellow cyclic
colour scheme for phases $-\pi\rightarrow\pi$. The plot also highlights the
characteristic radii $r_{c}$ (dashed blue) and $r_{t}$ (dotted red),
demonstrating the absence of vorticity within such radial regions – a feature
which persists during the entire probed evolution (see supplementary movie [SM
Movie 2] (b) The evolution of the local integral of the modulus square of the
incompressible velocity field, $\zeta(x,t)$, for $x=0$ (black), $x=r_{c}$
(blue), $x=r_{t}$ (red), and $x=10r_{c}$ (dashed green) for $r_{c}$ computed
from the time averaged density profile shown in Fig. 3.

###  5.1 Vortical structures outside the core

The energy distributions illustrated in Fig. 4 show that, opposite to the core
region where the quantum pressure dominates, in the outer halo region, or more
specifically for $r>r_{t}$, the classical kinetic energy is as important as
the quantum pressure one, namely
$\varepsilon_{\mathrm{ke}}(r)\approx\varepsilon_{\mathrm{qp}}(r)$. In
addition, the contribution of the kinetic energy’s incompressible component is
slightly higher than that of the compressible component, as discussed in Sec.
3.2. The rotational/incompressible part of the velocity $\mathbf{v}^{i}$ and
density current $\mathbf{F}^{i}$ in a scalar field can only receive
contributions from quantized vortices. In the vortex core the velocity
exhibits a divergence, see Eq. (11), while the density current converges to a
constant. Here we use the high incompressible velocity amplitude as a probe in
order to visualize vortex features in the outer halo [Hui et al. (2021)]. The
vortices identified are plotted in Fig. 6 (a) where the isosurface of
$|\mathbf{v}^{i}(\mathbf{r})|^{2}=800l_{\mathrm{ref}}^{2}/\tau_{\mathrm{ref}}^{2}\approx(\hbar^{2}/2m^{2})(\Delta
x)^{-2}$ is illustrated in red, vividly showcasing the tangled filamentary
structures that feature prominently outside the core.

To verify that these filaments are indeed quantized vortices, we also plot
slices of the phase and density profiles. The vortical filaments pierce
through the points at the centre of a $2\pi$ phase gradient and extremely low
density regions, providing support to the interpretation of these structures
as quantized vortices. This vortex tangle is made up from deformed vortex
rings or, more accurately, topologically closed loops. In our simulations,
vortex rings are formed which can expand, shrink and reconnect throughout the
dynamical evolution [Chiueh et al. (2011); Hui et al. (2021)]. In a
supplementary movie [SM Movie 2], we present the full dynamical evolution for
the halo displayed in Fig 6 (a). One can also readily observe there that the
vortices in the inner halo move faster than the outer ones, which reconnect to
each other rapidly.

In [Mocz et al. (2017); Hui et al. (2021)] it was mentioned that vortices
appear outside the core. To quantify this statement further, in Fig 6 (a) we
mark the region bounded by the scales of $r_{c}(t)$ and $r_{t}(t)$ and find
that throughout the evolution the phase within $r<r_{c}$ is almost fully
coherent, with the phase fluctuation increasing as we move outward to $r_{t}$
and beyond. From this observation, we indeed expect that there will be no
vortices passing within the core region since a vortex would be a source of
phase incoherence and density fluctuations. To corroborate such result,
supported by visual observation of our simulations, we introduce the local
integral of the modulus square of the incompressible velocity field,

$\zeta(x,t)=\int_{V_{c}(x)}d\mathbf{r}|\mathbf{v}^{i}(\mathbf{r},t)|^{2}$ (37)

where $V_{c}(x)$ stands for the volume of a spherical region of radius
$r_{c}$, placed at different positions along the $x$ axis. This quantity can
qualitatively capture the incident rate of vortices crossing $V_{c}(x)$, a
volume equal to that of the core but placed at various distances from the
centre. Fig. 6 (b) shows $\zeta(x,t)$ for 4 different positions. For $r=0$,
$\zeta(0,t)$ fluctuates around the values of $0.0344\pm 0.0185$, while
$\zeta(r_{c},t)$ exhibits a similar pattern ($0.0498\pm 0.0312$) (occasionally
rising briefly to the order of $0.1$). Moving towards the crossover region,
the fluctuation rises significantly with values $0.5446\pm 0.6960$, where the
large variance indicates a considerable contribution to $\zeta$ when vortices
pass through. A similar pattern appears for $\zeta(x=10r_{c},t)=0.4352\pm
0.5811$ where we observe a peak value of comparable magnitude but a longer
period of low values compared to $\zeta(x=r_{t},t)$. These features show that
there are no vortices entering the soliton core region, or more precisely, the
region $|\mathbf{r}|<r_{c}$ through time.

Around the crossover region there are several incidents of vortices crossing
our sample volume, indicating a considerable population there. The vortex-
crossing incident rate drops compared to the crossover region as we move
towards the outer halo, in agreement with the incompressible energy density
distribution in Fig. 4 (a).

###  5.2 Superfluid turbulence in the outer halo

Figure 7:  (a) 3D volume rendering for the density and the high vortical
velocity field, displaying densities in the range $[2.5\times
10^{-5},0.5]\rho_{c}(t)$; this demonstrates that the typical inter-vortex
distance coincides with the typical granule size with vortices winding in
between locally coherentdensity lumps. (b) The time-averaged incompressible
kinetic energy spectra for halos with 3 different total masses. The spectrum
peak indicates the typical inter-vortex distance. The black dashed line
indicates the $k^{-3}$ scaling behaviour at high momenta which probes the
vortex core structure - see text. (c) The overdensity power spectrum for the
same configurations as in (b). In (b) and (c) we show the momenta at which the
total density field $\rho(\mathbf{r},t)$ (hollow squares), overdensity
$\delta\rho(\mathbf{r},t)$ (filled diamonds) and incompressible power spectra
(hollow circles) respectively peak. Note that the power spectrum for the total
field $\rho(\mathbf{r},t)$ peaks at higher $k$ values (hollow squares)
compared to the overdensity (filled diamonds). This is due to the somewhat
smaller size of the core which contributes to the $\rho$ but not $\delta\rho$.
The error bands in (b) illustrate the standard deviation from the time
average, and the small horizontal error bars in both (b) and (c) show the
standard deviation of the peak location over the averaged times. (d) The
linear relation between $k_{\mathrm{max}}^{\Delta^{2}_{\delta}}$ and
$k_{\mathrm{max}}^{\mathrm{ke},i}$ is shown with the black dashed line being
the best fit from the scattered data giving
$k_{\mathrm{max}}^{\mathrm{ke},i}=(0.9300\pm
0.0274)k_{\mathrm{max}}^{\Delta^{2}_{\delta}}$. The outlier
$M=50M_{\mathrm{ref}}$ simulation (bottom right) is the most bounded among our
simulations with $E/M\approx-200$ and with a very dense core. The highlghted
data points correspond to the simulations represented in (b) and (c) with our
primary simulation represented by the black coloured data point.

Turbulence in a quantum/superfluid, carried by quantized vortices, is a
subject that has been extensively studied in non-self-gravitating atomic Bose-
Einstein condensate systems [Nore et al. (1997a); Tsubota (2014); Barenghi et
al. (2014); Tsatsos et al. (2016)], and in superfluid helium. Turbulence is
often regarded as a chaotic state but it can contain coherent structures in
the tangle of vortex lines [Tsubota (2014)]. In fact, during the equilibration
of an ultracold quantum fluid from a highly non-equilibrium initial state,
Berloff & Svistunov (2002) demonstrated the co-existence of dynamically-
evolving vortices and coherent structures, with the long-range coherence of
the system suppressed by the chaotic motions of the vortices (features also
supported during the controlled cooling of a trapped ultracold superfluid [Liu
et al. (2018)]). This view agrees qualitatively with the results depicted in
Fig. 3 and discussed in section 3.1, suggesting that the outer halo is indeed
made up of locally coherent (density) structures.

Vorticity in FDM halos has been studied in [Mocz et al. (2017)] and [Hui et
al. (2021)] with the main conclusion of the earlier work [Chiueh et al.
(2011)], which did not include gravity, surviving the inclusion of the
gravitational long-range force: vortices in such a system occur at places of
destructive interference where the density can become zero, giving rise to the
possibility of the phase winding by $2\pi$ around such regions. Interestingly,
we note that, unlike cases of dynamically equilibrating superfluids in which
efficient mode-mixing through the nonlinearity induced by the self-interaction
ensures vortices can transfer energy into sound waves, the vortex structures
in the non-self-interacting FDM case discussed here do not dissipate. We
further comment on this below.

One of the key quantities characterising quantum turbulence is the
incompressible energy spectrum which coincides with the vortex energy
spectrum, [Nore et al. (1997a, b); Tsubota & Kasamatsu (2012); Baggaley et al.
(2012a, b); Tsubota (2014); Tsubota et al. (2017)]. Here we compute the
incompressible energy spectrum in the form of [Nore et al. (1997a)]

$\tilde{\epsilon}_{\mathrm{ke}}^{i}(k_{r})\equiv\int
d\mathbf{\Omega}\,k^{2}\tilde{\epsilon}_{\mathrm{ke}}^{i}(\mathbf{k})=\sum_{k_{r}\leq|\mathbf{k}|\leq
k_{r}+\Delta k}\tilde{\epsilon}_{\mathrm{ke}}^{i}(\mathbf{k})$ (38)

which is a cumulative measure of the amplitude of
$\varepsilon_{\mathrm{ke}}^{i}(\mathbf{r})$’s (incompressible energy density)
Fourier components for $k_{r}\leq|\mathbf{k}|<k_{r}+\Delta k$.

In Fig. 7 (b) we present the time-averaged
$\tilde{\varepsilon}^{i}_{\mathrm{ke}}(k_{r})$, together with the variance
over time in shaded bands, for the three different total halo masses that we
have simulated. We can see that after $t=\tau_{\mathrm{vir}}$, the variance of
$\tilde{\varepsilon}_{\mathrm{ke}}^{i}$ is marginal in time, suggesting a good
convergence and that the system is in a steady state of turbulence, as first
examined for the Schrödinger equation in [Chiueh et al. (2011)]. The location
of the peak apparent in the incompressible kinetic energy spectrum, marked by
hollow circles in Fig. 7 (b), signifies the characteristic intervortex length
scale [Stagg et al. (2016); Baggaley et al. (2012a)]. As we discuss in the
following subsection, this scale is also correlated to the typical granule
size. In the high-$k$ region we also find a clear $k^{-3}$ behaviour,
exhibiting a slight tilting around the grid resolution length scale. This last
feature is an artefact of being too close to the smallest characteristic
length scales of the system. The
$\tilde{\varepsilon}^{i}_{\mathrm{k}}(k_{r})\propto k^{-3}$ behaviour carries
physical information and directly probes the vortex core structure [Nore et
al. (1997a); Stagg et al. (2016)], signifying that $\rho v^{2}\sim{\rm
constant}$ on short scales888A power law behaviour $r^{-s}$ leads to a $\sim
k^{s-3}$ Fourier transform and, correspondingly, to a $\sim k^{2s-3}$ power
spectrum.. Given that quantized vortices must exhibit a
$|\mathbf{v}^{i}|\propto r^{-1}$ velocity profile around the vortex cores (see
Eqn. (11)) it follows that our measured spectrum signifies a $\rho\sim r^{2}$
profile for the local density around the core region as theoretically
expected, see e.g. [Nore et al. (1997a); Chiueh et al. (2011)]. We note here
that FDM turbulence was first numerically studied by Mocz et al. (2017)999The
turbulent features in FDM were also discussed in Woo (2009). through the
spectrum of the full velocity for which they quoted a $k^{-1.1}$ behaviour.
Although the authors of that work did not directly comment on this, this
appears consistent with the expected velocity profile around the vortex cores.

A conspicuous feature of the incompressible kinetic energy spectrum is the
absence of any segment exhibiting the Kolmogorov $k^{-5/3}$ scaling, usually
associated with turbulent states. Furthermore, unlike most turbulence studies
in laboratory Bose-Einstein condensates which show a decay of the turbulent
state101010The decay can be associated with an energy cascade [Madeira et al.
(2020)] resulting in energy transfer from the system’s incompressible
component to the compressible one. together with the decay of the number and
line length of vortices via annihilations between vortex-antivortex pair and
the emission of sound waves/phonons through Kelvin waves [Leadbeater et al.
(2003); Kobayashi & Tsubota (2005); Stagg et al. (2016)], FDM turbulence does
not appear to decay in time as seen for example from the evolution of
incompressible energy, see Fig. 2 (b). The emergence of a steady state that
does not exhibit the Kolmogorov spectrum can be attributed to the absence of
both a self interaction coupling all length scales via a nonlinear term, and
of any source of dissipation, and was already recognised in [Chiueh et al.
(2011)] which studied the turbulence emerging in numerical solutions to the
Schrödinger equation. It appears that the addition of gravity only provides
spatial differentiation of the density and does not alter the basic picture of
this type of turbulence discussed in that work.

###  5.3 Granule spectrum

Lastly, we look at the interplay of turbulence and the structural opposite of
vortices, those regions of constructive interference that, as we discussed
above, exhibit local density coherence and are commonly referred to as
granules in the FDM literature. To characterize the density
fluctuations/granules, we reformulate the density as [Chan et al. (2018); Lin
et al. (2018); Dutta Chowdhury et al. (2021)]

$\rho(\mathbf{r},t)=\bar{\rho}(r)\left[1+\delta\rho(\mathbf{r},t)\right]$ (39)

where now $\bar{\rho}(r)$ is the spherically symmetric core-halo density
profile, obtained by the azimuthal average of the time averaged density,
$\langle\rho(\mathbf{r})\rangle_{t}$.111111$\bar{\rho}$ in this section should
not be confused with the use of the same symbol in earlier sections which
denotes a different average. Here the overdensity is defined by

$\delta\rho(\mathbf{r},t)=\frac{\rho(\mathbf{r},t)-\bar{\rho}(r)}{\bar{\rho}(r)}$
(40)

where $\rho(\mathbf{r},t)$ is re-centred by setting $\rho(\mathbf{r}=0)$ as
the peak density and gives the density fluctuations on top of the spherical
$\bar{\rho}(r)$. The azimuthal avaraging aims to reduce artefacts on the
identification of the granule size from very slow moving local structures that
could remain in $\langle\rho(\mathbf{r})\rangle_{t}$ after temporal averaging
as well as from the random motions of the soliton core.

In order to to probe the characteristc granule length scales [Chan et al.
(2018); Dutta Chowdhury et al. (2021); Lin et al. (2018)] and similar to the
discussion on the core, we consider the power spectrum of the overdensity,

$\begin{array}[]{rl}\Delta^{2}_{\delta}(k_{r},t)=&\displaystyle\frac{k_{r}^{3}}{(2\pi^{2})4\pi
k_{r}^{2}}\int d\mathbf{\Omega}_{k}P_{\delta}(\mathbf{k},t)\\\ \\\
\approx&\displaystyle\frac{k_{r}^{3}}{2\pi^{2}\mathcal{N}_{k_{r}}}\sum_{k_{r}\leq|\mathbf{k}|<k_{r}+\Delta
k}|\mathbf{k}|P_{\delta}(\mathbf{k},t)\end{array}$ (41)

where $P_{\delta}(\mathbf{k},t)=|\tilde{\eta}_{\delta}(\mathbf{k},t)|^{2}$ and
$\tilde{\eta}_{\delta}(\mathbf{k},t)$ is the Fourier transformation of the
overdensity $\delta(\mathbf{r},t)$. In Fig. 7 (c) we show the time averaged
power spectrum of the overdensity for the same set of data as in (b). The peak
locations of the overdensity power spectrum are marked as filled diamonds in
(c) and are also marked in (b). In Fig. 7 (d) we plot the relation of
$k_{\mathrm{max}}^{\Delta^{2}_{\delta}}$ and
$k_{\mathrm{max}}^{\mathrm{ke},i}$ for our 30 simualtions exhibiting a clear
linear relation between these two quantities. The best fit for the scattered
data is

$k_{\mathrm{max}}^{\mathrm{ke},i}=(0.9300\pm
0.0274)k_{\mathrm{max}}^{\Delta^{2}_{\delta}}$ (42)

showing the typical granule size to be essentially equal or marginally smaller
than the typical intervortex distance. This is consistent with the picture of
granules being regions of constructive interference with vortices surrounding
them, residing in low density regions of destructive interference. Fig. 7 (a)
shows a 3D volume rendering for the density and the regions identified as
vortices via their high velocity values, visually illustrating that idea that
vortices surround the granules in a FDM halo, as first pointed out in [Hui et
al. (2021)].

##  6 Discussion and Conclusions

In this work, we have carried out 30 simulations of virialized FDM halos,
generated by the merging of randomly located coherent solitons. The final
products of our merger numerical experiments contain the characteristic FDM
solitonic cores at the centre and exhibit the striking FDM granular structure
outside them, which, when averaged over time and angles, follows an $r^{-3}$
trend consistent with the NFW profile. We then examined these halos through
the prism of the field’s coherence, employing measures widely used in the
study of laboratory Bose-Einstein condensates: the lowest eigenmodes of the
single-particle density matrix $\rho(\mathbf{r},\mathbf{r^{\prime}})$, the
correlation functions $g_{1}$ and $g_{2}$ as well as the reduced phase space
density.

All coherence measures are compatible with a bimodal spatial configuration
that involves a central, fully coherent Bose-Einstein condensate that aligns
perfectly with the soliton, embedded in a turbulent halo that is on average
incoherent but appears to exhibit local (spatially and temporally) coherent
quasi-condensate patches, with an identifiable crossover region between the
solitonic core and (outer) halo. The configuration is in many respects
analogous to finite temperature laboratory BECs confined in a (harmonic)
potential where the condensate sits at the bottom of the potential and is
surrounded by a cloud of incoherent (thermal) particles. Here, the role of the
trap is played by the gravitational potential of the whole halo. It is also
worth stressing that, although the particle number is very large and fully
justifies a classical field approach, the _reduced_ phase space density beyond
the crossover region drops below unity, indicating the possibility that a
description via distinguishable, individual particles may capture some
features of the halo dynamics outside the core.

The identification of the central soliton with a perfect BEC state was
obtained via the numerical determination of the Penrose-Onsager mode of the
field configuration (mode with largest eigenvalue of the single-particle
density matrix) which coincides with the empirically derived soliton profile.
The outer halo was deemed incoherent on average. However, a more detailed
inspection of the $g_{2}$ correlation function, when averaged over differing
timescales, shows that some density coherence remains over short temporal
intervals and localised in space. These features may be tentatively identified
with the "granular texture" often employed to describe the appearance of the
outer parts of FDM halos. Precise characterization of these features though
would require more extensive numerics and should be left for future work.

We found that the "core-crossover region-outer halo" spatial structure implied
by the measures of field coherence doesn’t only reflect static or average
properties but also describes the dynamical oscillations of the core.121212We
have also confirmed the random motion of the central core reported in [Li et
al. (2021); Schive et al. (2020); Dutta Chowdhury et al. (2021)].
Interestingly, these oscillations were found to be anti-correlated to the peak
location of the power spectrum of the whole halo. This would suggest that
either the core leaves its imprint on the field fluctuation power spectrum, at
least for the halos in our numerical experiments, or that the typical size of
the density fluctuations, the granules, varies with time. Clarification of
this would be another topic for future work but preliminary results indicate
that the granule size appears to remain fixed with time, suggesting the core
is likely driving the power spectrum oscillation.

The phase incoherence of the halo outside the core is expressed via the
formation of a tangle of quantum vortices, surrounding the phase and density
semi-coherent granules and exhibiting sustained quantum turbulence with
persistent, non-decaying, incompressible kinetic energy. We identified
filamentary structures as regions of high velocity and confirmed that they are
indeed quantum vortices by cross-checking with the local phase profiles around
them. These vortices are confined outside the radius of the solitonic core. We
then examined the power spectrum of the incompressible kinetic energy which we
found to exhibit a characteristic peak that marks the intervortex distance.
Its location was found to essentially coincide with the peak of the granule
power spectrum, which in turn corresponds to a characteristic granule size.
Thus, we have explicitly found that vortices do indeed form in between the
granules. Furthermore, vortices are not ideal line-like defects but have a
specific width over which density and velocity develop characteristic
profiles. We found the incompressible kinetic energy spectrum to exhibit a
$k^{-3}$ scaling in the ultraviolet momentum region which, given that the
quantum vortex velocity profile has to be $|\mathbf{v}^{i}|\propto r^{-1}$,
implies a density $\rho\sim r^{2}$ around the vortex cores. It might be worth
noting that the vortex energy spectrum of the halos in our simulations do not
exhibit the well known Kolmogorov $k^{-5/3}$ scaling, which is associated with
the bundles of vortex lines on length scales larger than the vortex size
[Baggaley et al. (2012b)].

Before closing let us make a few comments regarding the physical scales for
which our simulations can be applied. Considering a fiducial boson mass of
$m=2.5\times 10^{-22}$ eV/$c^{2}$ and a fiducial mass density of
$10^{3}M_{\odot}\textrm{kpc}^{-3}$, comparable to the mean cosmic value (and
noting that the estimated average dark matter density in a halo hosting a
milky way-type galaxy is approximately $\rho_{\mathrm{Halo}}\approx
10^{7}M_{\odot}\textrm{kpc}^{-3}$) our reference units can be written as

$E_{\mathrm{ref}}\approx 5.38\times
10^{-27}\left(\frac{\rho_{\mathrm{ref}}}{10^{3}M_{\odot}\textrm{kpc}^{-3}}\right)^{1/2}\mathrm{eV}\;,$
(43) $l_{\mathrm{ref}}\approx 10.81\left(\frac{2.5\times
10^{-22}\mathrm{eV}}{mc^{2}}\right)^{1/2}\left(\frac{10^{3}M_{\odot}\textrm{kpc}^{-3}}{\rho_{\mathrm{ref}}}\right)^{1/4}\mathrm{kpc}\;,$
(44)

and

$\tau_{\mathrm{ref}}\approx
14.9\left(\frac{10^{3}M_{\odot}\textrm{kpc}^{-3}}{\rho_{\mathrm{ref}}}\right)^{1/2}\mathrm{Gyr}\;,$
(45)

with

$\begin{array}[]{rl}M_{\mathrm{ref}}=&\rho_{\mathrm{ref}}l_{\mathrm{ref}}^{3}\\\
\\\ \approx&\displaystyle 1.26\times 10^{6}M_{\odot}\left(\frac{2.5\times
10^{-22}\mathrm{eV}}{mc^{2}}\right)^{3/2}\left(\frac{\rho_{\mathrm{ref}}}{10^{3}M_{\odot}\textrm{kpc}^{-3}}\right)^{1/4}\;.\end{array}$
(46)

The reference velocity in our simulations,
$v_{\mathrm{ref}}=l_{\mathrm{ref}}/\tau_{\mathrm{ref}}$, is $0.71$ km/s for
$\rho_{\mathrm{ref}}=10^{3}M_{\odot}\mathrm{kpc}^{-3}$, or e.g. $3.99$ km/s
for $\rho_{\mathrm{ref}}=10^{6}M_{\odot}\mathrm{kpc}^{-3}$, and the velocities
of the FDM particles in the outer parts of the halo are measured to be $5$ to
$10$ times that reference velocity in the total mass range of
$50M_{\mathrm{ref}}$ to $150M_{\mathrm{ref}}$. We thus see that the halos in
our merger simulations are "colder" than what would be expected for typical
dark matter halos hosting a Milky Way-type galaxy and therefore expect the
features discussed in this paper to be present at relatively smaller scales
within the halo in more realistic conditions. The reference phase space
density is

$\begin{array}[]{rl}{\cal
D}_{\mathrm{ref}}=&\displaystyle\frac{h^{3}\rho_{\mathrm{ref}}}{m^{4}v_{\mathrm{ref}}^{3}}\\\
\\\ \approx&\displaystyle 1.25\times
10^{90}\left(\frac{10^{-22}\textrm{eV}/c^{2}}{m}\right)^{4}\left(\frac{\rho_{\mathrm{ref}}}{10^{3}M_{\odot}\mathrm{kpc}^{-3}}\right)\left(\frac{250\mathrm{km/s}}{v_{\mathrm{ref}}}\right)^{3}\end{array}$
(47)

where, the value of 250 km/s was chosen as the typical velocity of a particle
in a Milky-way-like environment. Because the boson mass relevant to FDM is so
small, the resulting number of particles is enormous for our parameters, and
the condensation criterion (26) is amply satisfied throughout our halos. This
of course amply justifies the use of a wave equation.

Although we performed 30 randomly generated simulations, the configurations we
examined gave very similar sizes of core radii while the initial energy
configurations involved only about a factor of 5 difference. Furthermore, as
discussed above all our configurations are less energetic than typical dark
matter halos of comparable masses. In addition, the periodic boundary
condition may play a role in the overall halo while the propagating waves
could reenter the system. The natural next step would therefore be to apply
the coherence and quantum turbulence analysis of this work to more realistic
halos obtained from a bigger simulation with cosmological initial conditions.
Moreover, it would be interesting to include self-coupling and also
investigate possible observational signatures of these features which would be
specific to FDM halos. We will focus on these topics in future works.

## Acknowledgements

I.K.L. acknowledges funding from European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement No.
897324 (upgradeFDM), while G.R. and N.P.P. acknowledge funding from the
Leverhulme Trust (Grant no. RPG-2021-010). I.K.L. gratefully thanks Professor
Hsi-Yu Schive and Professor Tzihong Chiueh for their hospitality and
enlightening discussions during his extended visit to their group which
contributed towards the precise identification of some of the questions
addressed here and which allowed us to crystallize our thinking on them. We
also acknowledge indirect contributions by Milos Indjin, and discussions with
Carlo Barenghi, Piotr Deuar, Luca Galantucci, Cora Uhlemann and Alex Gough.
The simulations were conducted on the HPC Cluster Rocket of Newcastle
University.

## Data Availability

Data supporting this publication can be openly accessed under an ‘Open Data
Commons Open Database License’ on the site
[https://doi.org/10.25405/data.ncl.22004978].

## References

  * Ajaib (2013) Ajaib M. A., 2013, arXiv preprint arXiv:1302.5601
  * Al Khawaja et al. (2002) Al Khawaja U., Andersen J. O., Proukakis N. P., Stoof H. T. C., 2002, Phys. Rev. A, 66, 013615
  * Allen & Misener (1938) Allen J. F., Misener A., 1938, Nature, 141, 75
  * Anderson et al. (1995) Anderson M. H., Ensher J. R., Matthews M. R., Wieman C. E., Cornell E. A., 1995, Science, 269, 198
  * Baggaley et al. (2012a) Baggaley A. W., Sherwin L. K., Barenghi C. F., Sergeev Y. A., 2012a, Phys. Rev. B, 86, 104501
  * Baggaley et al. (2012b) Baggaley A. W., Laurie J., Barenghi C. F., 2012b, Phys. Rev. Lett., 109, 205304
  * Banik & Sikivie (2017) Banik N., Sikivie P., 2017, in Proukakis N., Snoke D. W., Littlewood P. B., eds, , Universal Themes of Bose-Einstein Condensation. Cambridge University Press, Cambridge, Chapt. 31, pp 608–619
  * Barenghi et al. (2014) Barenghi C., Skrbek L., Sreenivasan K. R., 2014, Proceedings of the National Academy of Science, 111 (Supplement 1), 4647
  * Berezhiani & Khoury (2015) Berezhiani L., Khoury J., 2015, Phys. Rev. D, 92, 103510
  * Berloff & Svistunov (2002) Berloff N. G., Svistunov B. V., 2002, Phys. Rev. A, 66, 013603
  * Bernal et al. (2018) Bernal T., Fernández-Hernández L. M., Matos T., Rodríguez-Meza M. A., 2018, MNRAS, 475, 1447
  * Beugnon & Navon (2017) Beugnon J., Navon N., 2017, Journal of Physics B: Atomic, Molecular and Optical Physics, 50, 022002
  * Binney & Tremaine (2011) Binney J., Tremaine S., 2011, Galactic dynamics. Vol. 13, Princeton university press
  * Blakie & Davis (2005) Blakie P. B., Davis M. J., 2005, Phys. Rev. A, 72, 063608
  * Blakie et al. (2008) Blakie P. B., Bradley A. S., Davis M. J., Ballagh R. J., Gardiner C. W., 2008, Advances in Physics, 57, 363
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Chan et al. (2018) Chan J. H. H., Schive H.-Y., Woo T.-P., Chiueh T., 2018, MNRAS, 478, 2686
  * Chan et al. (2022) Chan H. Y. J., Ferreira E. G. M., May S., Hayashi K., Chiba M., 2022, MNRAS, 511, 943
  * Chavanis (2011) Chavanis P.-H., 2011, Phys. Rev. D, 84, 043531
  * Chavanis & Delfini (2011) Chavanis P.-H., Delfini L., 2011, Phys. Rev. D, 84, 043532
  * Chen et al. (2021) Chen J., Du X., Lentz E. W., Marsh D. J. E., Niemeyer J. C., 2021, Phys. Rev. D, 104, 083022
  * Chiang et al. (2021) Chiang B. T., Schive H.-Y., Chiueh T., 2021, Phys. Rev. D, 103, 103019
  * Chiueh et al. (2011) Chiueh T., Woo T.-P., Jian H.-Y., Schive H.-Y., 2011, Journal of Physics B Atomic Molecular Physics, 44, 115101
  * Cockburn et al. (2011) Cockburn S. P., Negretti A., Proukakis N., Henkel C., 2011, Physical Review A, 83, 043619
  * Coles & Spencer (2003) Coles P., Spencer K., 2003, MNRAS, 342, 176
  * Dabo et al. (2008) Dabo I., Kozinsky B., Singh-Miller N. E., Marzari N., 2008, Phys. Rev. B, 77, 115139
  * Dalfovo et al. (1999) Dalfovo F., Giorgini S., Pitaevskii L. P., Stringari S., 1999, Reviews of Modern Physics, 71, 463
  * Davis et al. (1995) Davis K. B., Mewes M. O., Andrews M. R., van Druten N. J., Durfee D. S., Kurn D. M., Ketterle W., 1995, Phys. Rev. Lett., 75, 3969
  * Del Popolo & Le Delliou (2017) Del Popolo A., Le Delliou M., 2017, Galaxies, 5, 17
  * Dentler et al. (2022) Dentler M., Marsh D. J. E., Hložek R., Laguë A., Rogers K. K., Grin D., 2022, MNRAS, 515, 5646
  * Dettmer et al. (2001) Dettmer S., et al., 2001, Phys. Rev. Lett., 87, 160406
  * Di Paolo et al. (2019) Di Paolo C., Salucci P., Erkurt A., 2019, MNRAS, 490, 5451
  * Dmitriev et al. (2021) Dmitriev A. S., Levkov D. G., Panin A. G., Pushnaya E. K., Tkachev I. I., 2021, Phys. Rev. D, 104, 023504
  * Durrer (2020) Durrer R., 2020, The cosmic microwave background. Cambridge University Press
  * Dutta Chowdhury et al. (2021) Dutta Chowdhury D., van den Bosch F. C., Robles V. H., van Dokkum P., Schive H.-Y., Chiueh T., Broadhurst T., 2021, ApJ, 916, 27
  * Ferreira (2021) Ferreira E. G. M., 2021, A&ARv, 29, 7
  * Ferrier-Barbut et al. (2016) Ferrier-Barbut I., Kadau H., Schmitt M., Wenzel M., Pfau T., 2016, Phys. Rev. Lett., 116, 215301
  * Frenk & White (2012) Frenk C. S., White S. D. M., 2012, Annalen der Physik, 524, 507
  * Guth et al. (2015) Guth A. H., Hertzberg M. P., Prescod-Weinstein C., 2015, Phys. Rev. D, 92, 103513
  * Guzmán & Ureña-López (2004) Guzmán F. S., Ureña-López L. A., 2004, Phys. Rev. D, 69, 124033
  * Hadzibabic & Dalibard (2011) Hadzibabic Z., Dalibard J., 2011, Riv. Nuovo Cimento, 34, 389
  * Hu et al. (2000) Hu W., Barkana R., Gruzinov A., 2000, Phys. Rev. Lett., 85, 1158
  * Huang (2000) Huang K., 2000, Statistical Mechanics. John Wiley and Sons
  * Hui (2021) Hui L., 2021, ARA&A, 59, 247
  * Hui et al. (2017) Hui L., Ostriker J. P., Tremaine S., Witten E., 2017, Phys. Rev. D, 95, 043541
  * Hui et al. (2021) Hui L., Joyce A., Landry M. J., Li X., 2021, J. Cosmology Astropart. Phys., 2021, 011
  * Indjin et al. (2022) Indjin M., Liu I.-K. G., Proukakis N. P., Rigopoulos G., 2022, To Appear
  * Jelic-Cizmek et al. (2018) Jelic-Cizmek G., Lepori F., Adamek J., Durrer R., 2018, J. Cosmology Astropart. Phys., 2018, 006
  * Kagan & Svistunov (1997) Kagan Y., Svistunov B. V., 1997, Phys. Rev. Lett., 79, 3331
  * Kapitza (1938) Kapitza P., 1938, Nature, 141, 74
  * Kasprzak et al. (2006) Kasprzak J., et al., 2006, Nature, 443, 409
  * Kiessling (1999) Kiessling M. K.-H., 1999, arXiv, p. 9910247
  * Klaers et al. (2010) Klaers J., Schmitt J., Vewinger F., Weitz M., 2010, Nature, 468, 545
  * Kobayashi & Tsubota (2005) Kobayashi M., Tsubota M., 2005, Journal of the Physical Society of Japan, 74, 3248
  * Leadbeater et al. (2003) Leadbeater M., Samuels D. C., Barenghi C. F., Adams C. S., 2003, Phys. Rev. A, 67, 015601
  * Lee & Koh (1996) Lee J.-W., Koh I.-G., 1996, Phys. Rev. D, 53, 2236
  * Li et al. (2021) Li X., Hui L., Yavetz T. D., 2021, Phys. Rev. D, 103, 023508
  * Lin et al. (2018) Lin S.-C., Schive H.-Y., Wong S.-K., Chiueh T., 2018, Phys. Rev. D, 97, 103523
  * Liu et al. (2018) Liu I. K., Donadello S., Lamporesi G., Ferrari G., Gou S. C., Dalfovo F., Proukakis N. P., 2018, Communications Physics, 1, 24
  * Liu et al. (2020) Liu I. K., Dziarmaga J., Gou S.-C., Dalfovo F., Proukakis N. P., 2020, Physical Review Research, 2, 033183
  * London (1938) London F., 1938, Nature, 141, 643
  * Madeira et al. (2020) Madeira L., Cidrim A., Hemmerling M., Caracanhas M. A., dos Santos F. E. A., Bagnato V. S., 2020, AVS Quantum Science, 2, 035901
  * Marsh (2016) Marsh D. J. E., 2016, Phys. Rep., 643, 1
  * May & Springel (2021) May S., Springel V., 2021, MNRAS, 506, 2603
  * Mocz et al. (2017) Mocz P., Vogelsberger M., Robles V. H., Zavala J., Boylan-Kolchin M., Fialkov A., Hernquist L., 2017, MNRAS, 471, 4559
  * Morgan et al. (1998) Morgan S. A., Choi S., Burnett K., Edwards M., 1998, Phys. Rev. A, 57, 3818
  * Nikolaieva et al. (2021) Nikolaieva Y. O., Olashyn A. O., Kuriatnikov Y. I., Vilchynskii S. I., Yakimenko A. I., 2021, Low Temperature Physics, 47, 684
  * Nore et al. (1997a) Nore C., Abid M., Brachet M. E., 1997a, Physics of Fluids, 9, 2644
  * Nore et al. (1997b) Nore C., Abid M., Brachet M. E., 1997b, Phys. Rev. Lett., 78, 3896
  * Numasato et al. (2010) Numasato R., Tsubota M., L’Vov V. S., 2010, Phys. Rev. A, 81, 063630
  * O’Dell et al. (2000) O’Dell D., Giovanazzi S., Kurizki G., Akulin V. M., 2000, Phys. Rev. Lett., 84, 5687
  * Papadopoulos et al. (2007) Papadopoulos I., Wagner P., Wunner G., Main J., 2007, Phys. Rev. A, 76, 053604
  * Penrose & Onsager (1956) Penrose O., Onsager L., 1956, Physical Review, 104, 576
  * Pethick & Smith (2008) Pethick C. J., Smith H., 2008, Bose-Einstein Condensation in Dilute Gases
  * Petrov et al. (2000) Petrov D. S., Holzmann M., Shlyapnikov G. V., 2000, Phys. Rev. Lett., 84, 2551
  * Petrov et al. (2001) Petrov D. S., Shlyapnikov G. V., Walraven J. T. M., 2001, Phys. Rev. Lett., 87, 050404
  * Pitaevskii & Stringari (2003) Pitaevskii L. P., Stringari S., 2003, Bose-Einstein Condensation. Clarendon Press
  * Popov (1983) Popov V. N., 1983, Functional integrals in quantum field theory and statistical physics. Mathematical Physics and Applied Mathematics, D. Reidel Publishing Company, Dordrecht, Holland
  * Prokof’ev et al. (2001) Prokof’ev N., Ruebenacker O., Svistunov B., 2001, Phys. Rev. Lett., 87, 270402
  * Proukakis (2006) Proukakis N. P., 2006, Phys. Rev. A, 74, 053617
  * Proukakis & Jackson (2008) Proukakis N. P., Jackson B., 2008, J. Phys. B: At. Mol. Opt. Phys, 41, 203002
  * Proukakis et al. (2013) Proukakis N., Gardiner S., Davis M., Szymańska M., eds, 2013, Quantum Gases: Finite Temperature and Non-Equilibrium Dynamics: 1 (Cold Atoms). ICP
  * Proukakis et al. (2017) Proukakis N., Snoke D. W., Littlewood P. B., 2017, Universal Themes of Bose-Einstein Condensation. Cambridge University Press, doi:10.1017/9781316084366
  * Pueblas & Scoccimarro (2009) Pueblas S., Scoccimarro R., 2009, Phys. Rev. D, 80, 043504
  * Rindler-Daller & Shapiro (2012) Rindler-Daller T., Shapiro P. R., 2012, MNRAS, 422, 135
  * Rindler-Daller & Shapiro (2014) Rindler-Daller T., Shapiro P. R., 2014, Modern Physics Letters A, 29, 1430002
  * Rogers & Peiris (2021) Rogers K. K., Peiris H. V., 2021, Phys. Rev. Lett., 126, 071302
  * Ruffini & Bonazzola (1969) Ruffini R., Bonazzola S., 1969, Physical Review, 187, 1767
  * Salucci et al. (2007) Salucci P., Lapi A., Tonini C., Gentile G., Yegorova I., Klein U., 2007, MNRAS, 378, 41
  * Schive et al. (2014a) Schive H.-Y., Chiueh T., Broadhurst T., 2014a, Nature Physics, 10, 496
  * Schive et al. (2014b) Schive H.-Y., Liao M.-H., Woo T.-P., Wong S.-K., Chiueh T., Broadhurst T., Hwang W. Y. P., 2014b, Phys. Rev. Lett., 113, 261302
  * Schive et al. (2016) Schive H.-Y., Chiueh T., Broadhurst T., Huang K.-W., 2016, ApJ, 818, 89
  * Schive et al. (2020) Schive H.-Y., Chiueh T., Broadhurst T., 2020, Phys. Rev. Lett., 124, 201301
  * Schobesberger et al. (2021) Schobesberger S. O., Rindler-Daller T., Shapiro P. R., 2021, MNRAS, 505, 802
  * Schwabe & Niemeyer (2022) Schwabe B., Niemeyer J. C., 2022, Phys. Rev. Lett., 128, 181301
  * Schwabe et al. (2016) Schwabe B., Niemeyer J. C., Engels J. F., 2016, Phys. Rev. D, 94, 043513
  * Schwabe et al. (2020) Schwabe B., Gosenca M., Behrens C., Niemeyer J. C., Easther R., 2020, Phys. Rev. D, 102, 083518
  * Seidel & Suen (1990) Seidel E., Suen W.-M., 1990, Phys. Rev. D, 42, 384
  * Sikivie & Yang (2009) Sikivie P., Yang Q., 2009, Phys. Rev. Lett., 103, 111301
  * Sin (1994) Sin S.-J., 1994, Phys. Rev. D, 50, 3650
  * Soto et al. (2022) Soto A., Proukakis N. P., Rigopoulos G., 2022, To Appear
  * Stagg et al. (2016) Stagg G. W., Parker N. G., Barenghi C. F., 2016, Phys. Rev. A, 94, 053632
  * Stuhler et al. (2005) Stuhler J., Griesmaier A., Koch T., Fattori M., Pfau T., Giovanazzi S., Pedri P., Santos L., 2005, Phys. Rev. Lett., 95, 150406
  * Tsatsos et al. (2016) Tsatsos M. C., Tavares P. E., Cidrim A., Fritsch A. R., Caracanhas M. A., dos Santos F. E. A., Barenghi C. F., Bagnato V. S., 2016, Physics Reports, 622, 1
  * Tsubota (2014) Tsubota M., 2014, Journal of Statistical Mechanics: Theory and Experiment, 2014, 02013
  * Tsubota & Kasamatsu (2012) Tsubota M., Kasamatsu K., 2012, arXiv e-prints, p. arXiv:1202.1863
  * Tsubota et al. (2017) Tsubota M., Fujimoto K., Yui S., 2017, Journal of Low Temperature Physics, 188, 119
  * Uhlemann et al. (2014) Uhlemann C., Kopp M., Haugg T., 2014, Phys. Rev. D, 90, 023517
  * Veltmaat et al. (2018) Veltmaat J., Niemeyer J. C., Schwabe B., 2018, Phys. Rev. D, 98, 043509
  * Weinberg et al. (2015) Weinberg D. H., Bullock J. S., Governato F., Kuzio de Naray R., Peter A. H. G., 2015, Proceedings of the National Academy of Science, 112, 12249
  * Widrow & Kaiser (1993) Widrow L. M., Kaiser N., 1993, ApJ, 416, L71
  * Woo (2009) Woo T.-P., 2009, Astrophys. J., 697, 850
  * Yavetz et al. (2022) Yavetz T. D., Li X., Hui L., 2022, Phys. Rev. D, 105, 023512
  * Zagorac et al. (2022) Zagorac J. L., Sands I., Padmanabhan N., Easther R., 2022, Phys. Rev. D, 105, 103506
  * Zinner (2011) Zinner N. T., 2011, arXiv e-prints, p. arXiv:1108.4290

##  Appendix A Scaling Transformations of Schrödinger–Poisson Equations

The SPE are invariant under the Lifshtz-type transformation Sin (1994); Guzmán
& Ureña-López (2004),

$\left\\{t,x,\Phi,\psi,\rho\right\\}\rightarrow\left\\{\Lambda^{-2}\hat{t},\Lambda^{-1}\hat{x},\Lambda^{2}\hat{\Phi},\Lambda^{2}\hat{\Psi},\Lambda^{4}\hat{\rho}\right\\}\,,$
(48)

With the total mass, energy and angular momentum of the system transforming as

$\left\\{M,E,L\right\\}\rightarrow\left\\{\Lambda\hat{M},\Lambda^{3}\hat{E},\Lambda\hat{L}\right\\}\,.$
(49)

Furthermore, the SPE can also be scaled by the boson mass, $m\rightarrow\alpha
m^{\prime}$ as [Mocz et al. (2017); Li et al. (2021)],

$\left\\{t,x,\Phi,\psi,\rho\right\\}\rightarrow\left\\{\alpha\hat{t},\hat{x},\alpha^{-2}\hat{\Phi},\alpha^{-1}\hat{\Psi},\alpha^{-2}\hat{\rho}\right\\}$
(50)

and, accordingly,

$\left\\{M,E,L\right\\}\rightarrow\left\\{\alpha^{-2}\hat{M},\alpha^{-3}\hat{E},\alpha^{-2}\hat{L}\right\\}.$
(51)

and the combination of the above transformations leads to

$\left\\{t,x,\Phi,\psi,\rho\right\\}\rightarrow\left\\{\Lambda^{-2}\alpha\hat{t},\Lambda^{-1}\hat{x},\Lambda^{2}\alpha^{-2}\hat{\Phi},\Lambda^{2}\alpha^{-1}\hat{\Psi},\Lambda^{4}\alpha^{-2}\hat{\rho}\right\\}$
(52)
$\left\\{M,E,L\right\\}\rightarrow\left\\{\Lambda\alpha^{-2}\hat{M},\Lambda^{3}\alpha^{-3}\hat{E},\Lambda\alpha^{-2}\hat{L}\right\\}.$
(53)

These scaling properties would allow one to use an individual simulation for
describing configurations with different physical parameters by appropriately
scaling the reference units.

##  Appendix B Numerical simulation details

###  B.1 Numerical methods