# Hybrid model of condensate and particle Dark Matter: linear perturbations in
the hydrodynamic limit

Nick P. Proukakis<EMAIL_ADDRESS>Gerasimos Rigopoulos
<EMAIL_ADDRESS>Alex Soto<EMAIL_ADDRESS>School of Mathematics, Statistics and Physics,
Newcastle University, Newcastle upon Tyne, NE1 7RU, UK

(November 2023)

###### Abstract

We analyse perturbations of self-interacting, scalar field dark matter that
contains modes both in a coherent condensate state and an incoherent particle-
like state. Starting from the coupled equations for the condensate, the
particles and their mutual gravitational potential, first derived from first
principles in earlier work by the authors, we derive a hydrodynamic limit of
two coupled fluids and study their linearized density perturbations in an
expanding universe. We find that away from the condensate-only or particle-
only limits, and for certain ranges of the parameters, such self-interacting
mixtures can significantly enhance the density power spectrum above the
standard linear $\Lambda$CDM value at localised wavenumbers, even pushing
structure formation into the non-linear regime earlier than expected in
$\Lambda$CDM for these scales. We also note that such mixtures can lead to
degeneracies between models with different boson masses and self-coupling
strengths, in particular between self-coupled models and non-coupled Fuzzy
Dark Matter made up of heavier bosons. These findings open up the possibility
of a richer phenomenology in scalar field dark matter models and could further
inform efforts to place observational limits on their parameters.

##  I Introduction

For many decades now, Cold Dark Matter (CDM) has been the prevailing paradigm
for explaining structure formation in the universe. Indeed, the assumption
that the matter budget of the universe is dominated by a non-relativistic
particle with very little velocity dispersion and negligible coupling to
itself and the particles of the standard model has been able to explain a vast
swathe of observations, spanning the whole history of the observable universe
[1]. Despite these successes however, more detailed investigations of the
small scale structures predicted by numerical simulations of CDM have led to a
surge in the investigations of an alternative dark matter model called Fuzzy
Dark Matter (FDM) [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] \- see [15, 16,
17] for detailed reviews and comprehensive bibliography - which promises to
solve CDM’s small scale problems such as the Cusp-Core problem [18], the too-
big-to-fail problem [19, 20] or the missing satellites problem [21]. Although
the extent to which these problems pose fundamental challenges to CDM is still
unclear, their investigation, coupled with a lack of any compelling particle
candidate for dark matter, has made FDM and its generalizations an attractive
subject of study.

Unlike the corpuscular nature of CDM, with N-body simulations probing its
phase space distribution, FDM is modelled as a continuous field exhibiting
wave interference effects. Although the resulting large scale structure is
very similar to CDM, there are differences on smaller scales: halos now
exhibit a “granular” appearance and FDM predicts the formation of compact
objects from the balance of gravity and “quantum pressure”. Such objects are
the ground states of FDM and in a cosmological context reside in the centres
of more extended halos. These features present a natural cutoff in the
overdensity power spectrum of FDM, precisely because the quantum pressure acts
to prevent short scale overdensities from forming.

Although sometimes associated to, or motivated by, speculative axion-like
particles, FDM is in its basis a dark matter model composed of bosons, in
contrast to CDM which is presumably composed of some unknown type of fermion.
Bosons can exist in both condensed and non-condensed states, and FDM is
usually associated with the former state - see however [22] for a quantitative
assessment of the level of coherence, a related concept and a qualifier of
condensation, inside an FDM halo and [23] on how coherence of light
condensates formed in the early universe is prone to decay. Depending on how
the dark matter bosonic field came into existence in the early universe, it
may be that a particle description is also appropriate or at least offers a
sufficiently accurate dynamical picture for some of the field’s modes. In this
work we utilize a formalism developed by us in [24] to describe hybrid states,
where the bosonic field can have both condensed and non-condensed components,
and study linearized perturbations of such a model.111A similar dual particle-
condensate description for homogeneous states in the early universe was
recently put forward in [25] using the 2PI formalism. It would be interesting
to see how that approach relates to the one adopted in [24]. We work in the
non-relativistic limit using a hydrodynamical truncation of the Boltzmann
equation describing the non-condensed “particles” to leading order in the
self-coupling, ignoring collisional terms.

The outline of the paper is as follows: In the next section we describe our
basic equations and detail our approximations, leading to a hydrodynamical,
two-fluid description. Linearized density perturbation equations for these two
fluids are then given in section III and their solutions for the linear
density power spectra are shown in section IV. We observe that coupled
mixtures of condensed and non-condensed fluids, components that are
essentially differentiated by the existence or not of a quantum pressure term
$\propto k^{4}$, exhibit oscillatory features beyond the scale where FDM
models present a quantum-pressure-induced cutoff and can exhibit an increase
of power at certain localized wavenumbers, possibly pushing such scales into
the non-linear regime earlier at higher redshift than CDM. Of course, this
invalidates perturbation theory and calls for fully non-linear simulations
before any safe conclusions can be drawn on the effects of such models on
structure formation. We also suggest that such mixtures may show degeneracies
between models with different boson masses and self-coupling strengths, in
particular between self-coupled models and non-self-interacting FDM ($g=0$)
made up of _heavier_ bosons. Such degeneracies should be kept in mind when
placing observational limits on models whose only parameter is the boson mass.

##  II Hybrid Model and Hydrodynamic equations

First we discuss our hybrid condensate-particle model and derive the
corresponding hydrodynamic equations.

###  II.1 Hybrid Condensate-Particle Model

A generalised model encompassing both CDM and FDM as limiting cases was
previously derived by us in [24], based on the coexistence of two phases for a
single species of bosonic particle of mass $m$ under the effect of gravity. In
flat spacetime the equations that describe this system are

$\displaystyle i\hbar\frac{\partial\Phi_{0}}{\partial
t}=\left[-\frac{\hbar^{2}}{2m}\nabla^{2}+\left(mV+g(n_{c}+2\tilde{n})\right)-i\hbar
R\right]\Phi_{0}$ (1) $\displaystyle\frac{\partial f}{\partial
t}+\frac{\mathbf{p}}{m}\cdot\nabla
f-\nabla\left(mV+2g(n_{c}+\tilde{n})\right)\nabla_{\mathbf{p}}f=\frac{1}{2}(I_{a}+I_{b})$
(2) $\displaystyle\nabla^{2}V=4\pi Gm(n_{c}+\tilde{n})$ (3)

where equation (1) describes the part of the bosons that are in a condensate
state, (2) is a collisional Boltzmann equation describing the incoherent
particles, i.e. the non-condensed part, and (3) determines the gravitational
potential generated by both of these components, whose number densities are
defined as $n_{c}=|\Phi_{0}|^{2}$ for the condensate and $\tilde{n}$ for the
incoherent particles; the latter’s particle density is related with the
distribution function $f$ which appears in (2) via

$\tilde{n}=\int\frac{d^{3}p}{(2\pi\hbar)^{3}}f\;.$ (4)

The bosons are further allowed to self interact via a self-interaction
coupling $g$. The right-hand-side (rhs) of the Boltzmann equation involves the
collision integrals

$\displaystyle I_{a}$ $\displaystyle=$ $\displaystyle
4g^{2}n_{c}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{2}\hbar^{4}}\delta(\varepsilon_{c}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{p}_{2}-\mathbf{p}_{1}-\mathbf{q}+\mathbf{p}_{3})$
(5)
$\displaystyle\times(\delta(\mathbf{p}_{1}-\mathbf{p})-\delta(\mathbf{p}_{2}-\mathbf{p})-\delta(\mathbf{p}_{3}-\mathbf{p}))$
$\displaystyle\times((1+f_{1})f_{2}f_{3}-f_{1}(1+f_{2})(1+f_{3}))$
$\displaystyle I_{b}$ $\displaystyle=$ $\displaystyle
4g^{2}\int\frac{d^{3}p_{2}d^{3}p_{3}d^{3}p_{4}}{(2\pi)^{5}\hbar^{7}}\delta(\varepsilon_{\mathbf{p}_{3}}+\varepsilon_{\mathbf{p}_{4}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}})\delta(\mathbf{p}+\mathbf{p}_{2}-\mathbf{p}_{3}-\mathbf{p}_{4})$
(6) $\displaystyle\times[f_{3}f_{4}(f+1)(f_{2}+1)-ff_{2}(f_{3}+1)(f_{4}+1)]$

where the subscripts in the distribution function are a shorthand
$f_{i}=f(\mathbf{r},\mathbf{p}_{i},t)$; here $\varepsilon_{c}$ corresponds to
the energy of a condensate element with momentum $\mathbf{q}$, and the
energies of the non-condensed particles, $\varepsilon_{\mathbf{p}}$, are
respectively defined by

$\displaystyle\varepsilon_{c}=\frac{q^{2}}{2m}+mV+g(n_{c}+2\tilde{n})\;,$ (7)
$\displaystyle\varepsilon_{\mathbf{p}}=\frac{p^{2}}{2m}+mV+2g(n_{c}+\tilde{n})\;.$
(8)

Meanwhile, the condensate equation involves a damping term expressed in terms
of

$R=\frac{1}{4n_{c}}\int\frac{d^{3}p}{(2\pi\hbar)^{3}}I_{a}\;.$ (9)

Notice that in the limit $g\to 0$ and considering only the condensed component
with $f=0$, we recover the limit of Fuzzy Dark Matter (FDM), while if
$\Phi_{0}=0$ and we only have non-condensed particles, we recover the limit of
the equations that are typically used for Cold Dark Matter (CDM). In this
hybrid model, we can see that at order $g^{2}$ it is possible to have
interchange of particles and condensate reflected in the existence of $-i\hbar
R$ in (1) and $I_{a}$ in (2). In this way, the self-interaction is a
significant part of the phenomenology of the model. Even if we ignore these
particle interchanges and only consider phenomena up to order $g$, the
coupling between condensed and non-condensed components is still non-trivial.

We next map our equations (1)-(3) into a system of hydrodynamic equations.

###  II.2 Hydrodynamic Equations

The hydrodynamic equations for the condensate are routinely obtained via the
Madelung transformation $\Phi_{0}=\sqrt{n_{c}}e^{i\theta}$ [26, 27]: Splitting
the resulting equation in real and imaginary parts and defining the condensate
velocity $\mathbf{v}_{c}=\frac{\hbar\nabla\theta}{m}$, and with
$\rho_{c}=mn_{c}$ and $\tilde{\rho}=m\tilde{n}$ being the mass densities for
the condensate and particles respectively, we get for the condensate

$\displaystyle\frac{\partial\rho_{c}}{\partial
t}+\nabla\cdot(\rho_{c}\mathbf{v}_{c})$ $\displaystyle=$
$\displaystyle-2\rho_{c}R\;,$ (10)
$\displaystyle\left(\frac{\partial}{\partial
t}+\mathbf{v}_{c}\cdot\nabla\right)\mathbf{v}_{c}$ $\displaystyle=$
$\displaystyle-\nabla\left(-\frac{\hbar^{2}}{2m^{2}}\frac{\nabla^{2}\sqrt{\rho_{c}}}{\sqrt{\rho_{c}}}+V+\frac{g}{m^{2}}(\rho_{c}+2\tilde{\rho})\right)\;.$
(11)

For the particles on the other hand, we can compute from equation (2) the
first three moments, multiplying it by $1$, $(\mathbf{p}-m\mathbf{v}_{c})$ and
$(\mathbf{p}-m\mathbf{v}_{c})^{2}$ respectively, and integrating with respect
to $\mathbf{p}$. It can be verified that

$\displaystyle\int\frac{d^{3}p}{(2\pi\hbar)^{3}}I_{b}=0,\quad\int\frac{d^{3}p}{(2\pi\hbar)^{3}}(\mathbf{p}-m\mathbf{v}_{c})I_{a}=0,\quad\int\frac{d^{3}p}{(2\pi\hbar)^{3}}(\mathbf{p}-m\mathbf{v}_{c})I_{b}=0$
(12)
$\displaystyle\int\frac{d^{3}p}{(2\pi\hbar)^{3}}(\varepsilon_{\mathbf{p}}-\varepsilon_{c})I_{a}=0,\qquad\int\frac{d^{3}p}{(2\pi\hbar)^{3}}\varepsilon_{\mathbf{p}}I_{b}=0$
(13)

and with these identities and Eq. (9), we obtain

$\displaystyle\frac{\partial\tilde{\rho}}{\partial
t}+\nabla\cdot(\tilde{\rho}\widetilde{\mathbf{v}})$ $\displaystyle=$
$\displaystyle 2\rho_{c}R$ (14) $\displaystyle\left(\frac{\partial}{\partial
t}+\widetilde{\mathbf{v}}\cdot\nabla\right)\widetilde{\mathbf{v}}$
$\displaystyle=$
$\displaystyle-\nabla\left(V+\frac{2g}{m^{2}}(\rho_{c}+\tilde{\rho})\right)$
(15)
$\displaystyle-\frac{1}{\tilde{\rho}}\nabla_{i}P_{ij}-2\frac{\rho_{c}}{\tilde{\rho}}R(\widetilde{\mathbf{v}}-\mathbf{v}_{c})$
$\displaystyle\frac{\partial E}{\partial
t}+\nabla\cdot(E\widetilde{\mathbf{v}})$ $\displaystyle=$
$\displaystyle-\nabla\cdot\mathbf{Q}-D_{ij}P_{ij}$ (16)
$\displaystyle+2\rho_{c}R\left(\frac{1}{2}(\widetilde{\mathbf{v}}-\mathbf{v}_{c})^{2}-\frac{\hbar^{2}\nabla^{2}\sqrt{\rho_{c}}}{2m^{2}\sqrt{\rho_{c}}}-\frac{g}{m^{2}}\rho_{c}\right)$

with $D_{ij}=\frac{1}{2}(\nabla_{i}\tilde{v}_{j}+\nabla_{j}\tilde{v}_{i})$ and
where we have defined as $\widetilde{\mathbf{v}}$ the bulk velocity of the
particles as

$\widetilde{\mathbf{v}}=\int\frac{d^{3}p}{(2\pi\hbar)^{3}}\frac{\mathbf{p}}{m}\frac{f}{\tilde{n}}$
(17)

and

$\displaystyle P_{ij}$ $\displaystyle=$ $\displaystyle
m\int\frac{d^{3}p}{(2\pi\hbar)^{3}}\left(\frac{p_{i}}{m}-\tilde{v}_{i}\right)\left(\frac{p_{j}}{m}-\tilde{v}_{j}\right)f$
(18) $\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{1}{2m}\int\frac{d^{3}p}{(2\pi\hbar)^{3}}(\mathbf{p}-m\widetilde{\mathbf{v}})^{2}f$
(19) $\displaystyle\mathbf{Q}$ $\displaystyle=$
$\displaystyle\frac{1}{2m}\int\frac{d^{3}p}{(2\pi\hbar)^{3}}(\mathbf{p}-m\widetilde{\mathbf{v}})^{2}\left(\frac{\mathbf{p}}{m}-\widetilde{\mathbf{v}}\right)f$
(20)

It is important to remark that, as is usually the case, the Boltzmann equation
leads to an infinite hierarchy of moment equations and here we have only
displayed equations for the first two moments. Accordingly, the existence of
terms such as $\mathbf{Q}$ and $E$ shows that the above equations are not
closed, requiring the computation of the next, cubic moments. Those in turn,
would require quartic moments, etc, in a never-ending hierarchy of equations.
To truncate it, further assumptions are required.

In order to close the equations we will assume that $f$ is an even function of
$(\mathbf{p}-m\widetilde{\mathbf{v}})$, meaning that we are considering the
distribution of velocities to be symmetric with respect to the average bulk
velocity. With this symmetry assumption, $\mathbf{Q}=0$ and, for $P_{ij}$, the
only non-zero entries are the diagonal ones:

$P_{ii}=m\int\frac{d^{3}p}{(2\pi\hbar)^{3}}\left(\frac{p_{i}}{m}-\tilde{v}_{i}\right)^{2}f\;.$
(21)

Furthermore, we will assume that the three diagonal entries are equal:
$P_{11}=P_{22}=P_{33}\equiv P$, i.e. isotropy in the $P_{ij}$ term. With this,

$E=\frac{3}{2}P\;.$ (22)

Hence, equation (16) becomes

$\frac{\partial P}{\partial
t}+\nabla\cdot(P\widetilde{\mathbf{v}})=-\frac{2}{3}P\nabla\cdot\widetilde{\mathbf{v}}+\frac{4}{3}\rho_{c}R\left(\frac{1}{2}(\widetilde{\mathbf{v}}-\mathbf{v}_{c})^{2}-\frac{\hbar^{2}\nabla^{2}\sqrt{\rho_{c}}}{2m^{2}\sqrt{\rho_{c}}}-\frac{g}{m^{2}}\rho_{c}\right)\;.$
(23)

We observe that now the system of equations are closed.

We highlight here that in order to construct the quantum pressure term
$-\frac{\hbar^{2}}{2m^{2}}\frac{\nabla^{2}\sqrt{\rho_{c}}}{\sqrt{\rho_{c}}}$
in (11) we needed to divide by $\sqrt{\rho_{c}}$ in some point, meaning that
we can never take $\rho_{c}$ as zero. Also, to arrange equation (15) in a
similar form to (11), we divided by $\tilde{\rho}$. In this way, we can never
take $\tilde{\rho}$ as zero either in the above coupled system of equations
for the two components. If we want to consider the hydrodynamical equations
for only condensate or only particles we should start from Eqs. (1) - (3),
setting to zero the absent component there and then repeat the steps for the
remaining component.

Next, we formulate the above equations in the setting of an expanding
universe.

###  II.3 Formulation in an Expanding Universe

At this point, in the interest of simplicity, we will keep the next
computations at order $g$, meaning that henceforth we set $R=0$ and,
correspondingly, $I_{a}=I_{b}=0$. Including cosmic expansion we have that
$\mathbf{r}=a\mathbf{x}$, where $\mathbf{r}$ is the physical coordinate,
$\mathbf{x}$ the comoving coordinate and $a(t)$ the scale factor. Also, we
have that the velocity is described by $\mathbf{u}=H\mathbf{r}+\mathbf{v}$,
where $\mathbf{u}$ is the physical velocity, $\mathbf{v}$ the peculiar
velocity and $H=\dot{a}/a$. Considering the expansion, the time derivatives
for a function can be written as $\left.\frac{\partial}{\partial
t}\right|_{\mathbf{r}}=\left.\frac{\partial}{\partial
t}\right|_{\mathbf{x}}-H\mathbf{x}\cdot\nabla$, where the nabla symbol stands
for derivatives on the comoving coordinate $\mathbf{x}$. With this, the
hydrodynamical equations in an expanding universe up to order $g$ are

$\displaystyle\frac{\partial\rho_{c}}{\partial
t}+3H\rho_{c}+\frac{1}{a}\nabla\cdot(\rho_{c}\mathbf{v})=0$ (24)
$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\frac{1}{a}\mathbf{v}\cdot\nabla\mathbf{u}=-\nabla\left(-\frac{\hbar^{2}}{2m^{2}a^{3}}\frac{\nabla^{2}\sqrt{\rho_{c}}}{\sqrt{\rho_{c}}}+\frac{1}{a}V+\frac{g}{m^{2}a}(\rho_{c}+2\tilde{\rho})\right)$
(25) $\displaystyle\frac{\partial\tilde{\rho}}{\partial
t}+3H\tilde{\rho}+\frac{1}{a}\nabla\cdot(\tilde{\rho}\widetilde{\mathbf{v}})=0$
(26) $\displaystyle\frac{\partial\widetilde{\mathbf{u}}}{\partial
t}+\frac{1}{a}\widetilde{\mathbf{v}}\cdot\nabla\widetilde{\mathbf{u}}=-\nabla\left(\frac{1}{a}V+\frac{2g}{m^{2}a}(\rho_{c}+\tilde{\rho})\right)-\frac{1}{a\tilde{\rho}}\nabla
P$ (27) $\displaystyle\frac{\partial P}{\partial
t}+5HP+\frac{1}{a}\nabla\cdot(\widetilde{\mathbf{v}}P)=-\frac{2}{3a}P\nabla\cdot\widetilde{\mathbf{v}}$
(28) $\displaystyle\nabla^{2}V=4\pi Ga^{2}(\rho_{c}+\tilde{\rho})$ (29)

The above set of equations is closed and corresponds to a possible starting
point for e.g. a perturbative analysis. We notice however that if we choose a
particular equation of state

$P=\kappa\tilde{\rho}^{5/3}$ (30)

then equation (28) becomes equivalent to (26) and therefore redundant. We will
therefore use this equation of state for the non-condensed particles with
$\kappa$ a free parameter in the ensuing analysis. Note that the equation of
state (30) has been considered in [28], also assuming the same symmetries in
$f$ and $P_{ij}$.

We observe that according to our equations, the main distinction between the
condensate and the particles lies in that the condensate possesses a quantum
pressure term (first term on the rhs of (25)) while the particles exhibit a
different pressure given by (30), which will appear on the rhs of (27). Also,
the effect of self-interaction is slightly different in both cases: notice the
difference between the factors of $2$ in those terms on the rhs of (25) and
(27) which are $\propto g$. This originates in the different mean-field
interaction potentials exhibited by condensed and non-condensed particles, and
is well-known in condensed matter systems [29, 30].

##  III Linear Regime equations

We will now study the perturbations around background values for the equations
(24)-(27) and (29) along with (30). Thus, considering the perturbations at
first order

$\rho_{c}=\bar{\rho}_{c}(1+\delta_{c}),\hskip
28.45274pt\tilde{\rho}=\bar{\rho}_{nc}(1+\delta_{nc}),\hskip
28.45274ptV=V_{0}+\delta V,\hskip 28.45274ptP=P_{0}+\delta P$ (31)

and the peculiar velocities around zero, assuming the background follows the
standard cosmological expansion of CDM, using the order zero equations in the
first order ones and rearranging them, we have that the equations for the
condensate and particles density contrasts are

$\displaystyle\ddot{\delta}_{c}+2H\dot{\delta}_{c}+\left(\frac{\hbar^{2}k^{4}}{4m^{2}a^{4}}-4\pi
G\bar{\rho}f+\frac{g\bar{\rho}fk^{2}}{m^{2}a^{2}}\right)\delta_{c}$
$\displaystyle\hskip 28.45274pt-(1-f)\left(4\pi
G\bar{\rho}-\frac{2g\bar{\rho}k^{2}}{m^{2}a^{2}}\right)\delta_{nc}=0$ (32)
$\displaystyle\ddot{\delta}_{nc}+2H\dot{\delta}_{nc}-\left(4\pi
G\bar{\rho}(1-f)-\left(\frac{2g\bar{\rho}(1-f)}{m^{2}}+\frac{5\kappa\bar{\rho}^{2/3}(1-f)^{2/3}}{3}\right)\frac{k^{2}}{a^{2}}\right)\delta_{nc}$
$\displaystyle\hskip 28.45274pt-f\left(4\pi
G\bar{\rho}-\frac{2g\bar{\rho}k^{2}}{m^{2}a^{2}}\right)\delta_{c}=0$ (33)

where we defined the condensate fraction $f$ as222Not to be confused with the
phase space density of the non-condensed component used in section II.

$f=\frac{\bar{\rho}_{c}}{\bar{\rho}}\,.$ (34)

and

$\bar{\rho}=\bar{\rho}_{c}+\bar{\rho}_{nc}\,.$ (35)

Taking the limit $g\rightarrow 0$ and $\kappa\to 0$ we obtain the following
reduced equations

$\displaystyle\ddot{\delta}_{c}+2H\dot{\delta}_{c}+\left(\frac{\hbar^{2}k^{4}}{4m^{2}a^{4}}-4\pi
G\bar{\rho}f\right)\delta_{c}-4\pi G\bar{\rho}(1-f)\delta_{nc}=0\;,$ (36)
$\displaystyle\ddot{\delta}_{nc}+2H\dot{\delta}_{nc}-\bar{\rho}(1-f)4\pi
G\delta_{nc}-4\pi G\bar{\rho}f\delta_{c}=0\;.$ (37)

We note that these equations are formally identical to the equations used in
[31], which however considered a different physical scenario of two different
particle species making up the dark matter, one of which is light enough to
exhibit wave-like properties on the relevant scales and the other being
sufficiently massive such that it can be described by N-body simulations. This
is physically distinct to our setting where both coherent and incoherent
particles stem from an internal partition of the same field on the basis of
phase-space density, as previously discussed in [24, 22]. In our general
equations, setting $\kappa\to 0$ means that we are considering the non-
condensate overdensity as a CDM-like pressureless fluid. Thus, a hybrid model
with $\kappa=0$ describes the coexistence of a condensate with a bath of CDM-
like particles. A non-zero $\kappa$ value generalises this situation,
modelling the pressure arising from the velocity dispersion of the self-
interacting particles, and allowing new possible effects in the non-condensate
sector.

The equations can be manipulated a little bit further by making use of the
Friedmann equation

$H^{2}=\frac{8\pi
G}{3}\left(\bar{\rho}+\rho_{\gamma}\right)+\frac{\Lambda}{3}$ (38)

in the form

$1=\bar{\Omega}+\Omega_{\gamma}+\Omega_{\Lambda}$ (39)

where $\bar{\Omega}=\frac{8\pi G(\bar{\rho}_{c}+\bar{\rho}_{nc})}{3H^{2}}$ is
the total matter density parameter, composed by the mixture of the condensed
and the non-condensed particles defined by our equations,
$\Omega_{\gamma}=\frac{8\pi G\bar{\rho}_{\gamma}}{3H^{2}}$ is the radiation
density parameter, and $\Omega_{\Lambda}=\frac{\Lambda}{3H^{2}}$ the
cosmological constant density parameter. We have assumed the baryonic matter
to be negligible compared with dark matter. With these definitions we replace
the density $\bar{\rho}=\frac{3H^{2}\bar{\Omega}}{8\pi G}$ in (32) and (33)
and change the time variable to be the scale factor, obtaining

$\displaystyle\frac{d^{2}\delta_{c}}{da^{2}}+\left(\frac{3}{a}+\frac{d\ln
H}{da}\right)\frac{d\delta_{c}}{da}-\frac{3\bar{\Omega}}{2a^{2}}\left(f-\frac{\hbar^{2}k^{4}}{6H^{2}\bar{\Omega}m^{2}a^{4}}-\frac{gfk^{2}}{4\pi
Gm^{2}a^{2}}\right)\delta_{c}$ $\displaystyle\hskip
28.45274pt-\frac{3\bar{\Omega}}{2a^{2}}(1-f)\left(1-\frac{gk^{2}}{2\pi
Gm^{2}a^{2}}\right)\delta_{nc}=0\;,$ (40)
$\displaystyle\frac{d^{2}\delta_{nc}}{da^{2}}+\left(\frac{3}{a}+\frac{d\ln
H}{da}\right)\frac{d\delta_{nc}}{da}$ $\displaystyle\hskip
28.45274pt-\frac{1}{a^{2}}\left(\frac{3\bar{\Omega}}{2}(1-f)-\frac{1}{8\pi
Ga^{2}}\left(\frac{6\bar{\Omega}g(1-f)}{m^{2}}+\frac{5\kappa(1-f)^{2/3}\bar{\Omega}^{2/3}(8\pi
G)^{1/3}}{(3H^{2})^{1/3}}\right)k^{2}\right)\delta_{nc}$ $\displaystyle\hskip
28.45274pt-\frac{3\bar{\Omega}}{2a^{2}}f\left(1-\frac{gk^{2}}{2\pi
Gm^{2}a^{2}}\right)\delta_{c}=0\;.$ (41)

Note that the wavenumber enters (41) as $k^{2}$ only, i.e. like a typical
pressure term, while (40) contains the $k^{4}$ term characteristic of the
quantum pressure felt by the condensate but not by the incoherent particles.
The two density contrasts $\delta_{c}$ and $\delta_{nc}$ are coupled via terms
proportional to the condensate fraction $f$ and $(1-f)$. Equations (40) and
(41) will form the basis of the analysis that follows.

##  IV Solution of the density contrasts in a matter dominated universe

Consider first the case of a radiation dominated universe with a cosmological
constant, where $\bar{\Omega}\approx 0$. Our equations (40) and (41) read

$\displaystyle\frac{d^{2}\delta_{c}}{da^{2}}+\left(\frac{3}{a}+\frac{d\ln
H}{da}\right)\frac{d\delta_{c}}{da}+\frac{\hbar^{2}k^{4}}{4H^{2}m^{2}a^{6}}\delta_{c}$
$\displaystyle=$ $\displaystyle 0\;,$ (42)
$\displaystyle\frac{d^{2}\delta_{nc}}{da^{2}}+\left(\frac{3}{a}+\frac{d\ln
H}{da}\right)\frac{d\delta_{nc}}{da}$ $\displaystyle=$ $\displaystyle 0\;.$
(43)

We observe that $\delta_{c}$ and $\delta_{nc}$ are decoupled and the self-
interaction and particle pressure don’t play any role. The fact that both
components are decoupled and both solutions in this stage are known will be
important below when we consider the transfer functions used in the literature
for both components.

Now, we move to the matter dominated case, where both components are coupled
and the behaviour is more involved. Here, $\bar{\Omega}\approx 1$ and our
equations (40) and (41) read

$\displaystyle\frac{d^{2}\delta_{c}}{da^{2}}+\frac{3}{2a}\frac{d\delta_{c}}{da}-\frac{3}{2a^{2}}\left(f-\frac{\hbar^{2}k^{4}}{6H_{0}^{2}\bar{\Omega}_{0}m^{2}a}-\frac{g}{m^{2}}\frac{fk^{2}}{4\pi
Ga^{2}}\right)\delta_{c}$ $\displaystyle\hskip
28.45274pt-\frac{3}{2a^{2}}(1-f)\left(1-\frac{g}{m^{2}}\frac{k^{2}}{2\pi
Ga^{2}}\right)\delta_{nc}=0\;,$ (44)
$\displaystyle\frac{d^{2}\delta_{nc}}{da^{2}}+\frac{3}{2a}\frac{d\delta_{nc}}{da}-\frac{3}{2a^{2}}\left((1-f)-\left(\frac{g}{m^{2}}\frac{(1-f)}{2\pi
Ga^{2}}+\frac{5\kappa}{6a}\left(\frac{(1-f)^{2}}{3\pi^{2}G^{2}H_{0}^{2}\bar{\Omega}_{0}}\right)^{\frac{1}{3}}\right)k^{2}\right)\delta_{nc}$
$\displaystyle\hskip
28.45274pt-\frac{3}{2a^{2}}f\left(1-\frac{g}{m^{2}}\frac{k^{2}}{2\pi
Ga^{2}}\right)\delta_{c}=0\;,$ (45)

where we have used that in this case

$H^{2}=\frac{H_{0}^{2}\bar{\Omega}_{0}}{a^{3}}$ (46)

and the values $H_{0}$ and $\bar{\Omega}_{0}$ are the values measured today.

We will solve these equations numerically. To set the initial conditions we
can use that since condensate and particles are decoupled before the time of
radiation-matter equality, and they each obey known equations individually, we
can use the same transfer functions for $\delta_{c}$ and $\delta_{nc}$ as have
been used to study them separately in the literature. Thus, for the non-
condensed part $\delta_{nc}$ we use the typical power spectrum for
$\Lambda$CDM

$P(k)=(T(k))^{2}A_{s}k^{n_{s}}$ (47)

where $A_{s}k^{n_{s}}$ is the primordial power spectrum with the values for
$A_{s}$ and $n_{s}$ obtained from the PLANCK data [32], and we will use the
following fitted transfer function $T(k)$ [33] defined by

$T(k)=\frac{\ln\left(1+2.34\frac{k}{k_{eq}}\right)}{2.34\frac{k}{k_{eq}}}\left(1+3.89\frac{k}{k_{eq}}+\left(16.1\frac{k}{k_{eq}}\right)^{2}+\left(5.46\frac{k}{k_{eq}}\right)^{3}+\left(6.71\frac{k}{k_{eq}}\right)^{4}\right)^{-\frac{1}{4}}$
(48)

where $k_{eq}$ is the wavenumber associated to the mode crossing the horizon
in the radiation-matter equality. The values for the Hubble parameter and
matter density parameter will be those obtained by PLANCK [32].

For the condensate $\delta_{c}$ we will use

$P_{FDM}(k)=(T_{FDM})^{2}P(k)=T_{FDM}^{2}(T(k))^{2}A_{s}k^{n_{s}}$ (49)

where the transfer function $T_{FDM}$ is defined as

$T_{FDM}=\frac{\cos(x^{3})}{1+x^{8}}$ (50)

with $x=1.61\left(\frac{m}{10^{-22}eV}\right)^{1/18}\frac{k}{k_{J\ast}}$, and
$k_{J\ast}=9\sqrt{\frac{m}{10^{-22}eV}}$ 1/Mpc, according to [2], and used in
[34] and [35].

To extract the behaviour of our hybrid model we numerically solve equations
(44) and (45) and construct the total overdensity

$\delta=\frac{\bar{\rho}_{c}\delta_{c}+\bar{\rho}_{nc}\delta_{nc}}{\bar{\rho}_{c}+\bar{\rho}_{nc}}=f\delta_{c}+(1-f)\delta_{nc}$
(51)

for which we compute the power spectrum $P(k)$ and dimensionless power
spectrum $\Delta^{2}(k)$, defined as

$P(k)=|\delta(k)|^{2},\quad\Delta^{2}(k)=\frac{k^{3}}{2\pi^{2}}P(k)\;.$ (52)

Based on these, we examine the effect of a variable condensate-fraction $f$,
the role of interactions through appropriate choices of $g$ values and how
such findings change when we additionally consider the pressure term,
parametrised by $\kappa$ of Eq. (30).

Figure 1:  Power spectrum in our hybrid condensate-particle dark matter model
for different condensate fractions, plotted as a function of $k$ in
logarithmic scale for a redshift $z=3$ and a reference boson mass $m=2\cdot
10^{-22}\ eV/c^{2}$ in the absence of self-interactions ($g=0$) and pressure
($\kappa=0$). The established limiting cases of CDM (solid black line, top)
and FDM (dot-dashed black line, bottom) are shown for comparison. It is clear
that increasing the condensate fraction gradually shifts the spectrum
downwards, as explicitly shown for the hybrid cases $f=0.2$ (blue), $f=0.5$
(orange) and $f=0.9$ (green).

###  IV.1 Role of Condensate Fraction in Simplest Non-Interacting ($g=0$),
Pressure-Less Particle ($\kappa=0$) limit.

First, we analyse the case of $g=0$ and $\kappa=0$. Linearised density
perturbations in this context are well known for the cases of CDM, for which
the condensate fraction $f=0$ by construction, and have also been explored in
FDM (e.g. [2, 3, 7] which refer to the $g=0$ case), for which $f=1$: these
arise naturally as limiting cases of our generalised hybrid equations and are
included here as a means of tabulating our model, and for comparison. Indeed,
as shown in Fig. 1, the hybrid linearised model including both such types of
particles through a value of $f\neq\left\\{0,1\right\\}$, does indeed smoothly
interpolate – as expected – between CDM (solid black line) and FDM (dot-dashed
black line) values with increasing condensate fraction. This demonstrates the
power of our hybrid approach in probing the parameter space between these two
established models in the literature. In this figure, we have used a specific
reference value of $m=2\cdot 10^{-22}\ eV/c^{2}$ for the boson mass, but the
conclusions are valid for any boson mass with higher (lower) masses moving the
FDM side of the envelope to the right (left).

###  IV.2 Combined Effect of Repulsive Interactions ($g>0$) and Variable
Condensate Fraction $f$ for Pressure-Less Particle Fluid ($\kappa=0$).

As a natural next step, we consider the effect of self-interactions in our
system, an established generalization of pure FDM [27, 36, 37, 38, 39, 40]. We
focus here on the case of repulsive self-interactions ($g>0$) and, to better
understand these, we limit the discussion in this section to $\kappa=0$.

####  IV.2.1 Study of Limiting Cases: $f\to 0$ (all particles) and $f\to 1$
(all condensate)

We first consider the limit $f\to 0$, which is the situation where all the
system is non-condensed and, for reference, we use the same mass $m=2\cdot
10^{-22}\ eV/c^{2}$. As we increase the interaction strength $g$ to positive
values, we find that the power spectrum starts deviating, and dropping below
(i.e. bounded upwards by) the CDM one (solid black line) at large $k$ values,
while also exhibiting oscillations: such behaviour is shown for various
interaction strengths $g$ in the left plot of Fig. 2. The deviation point from
the CDM curve moves to lower $k$ and the spectrum exhibits more frequent
oscillations with increasing $g$, as expected from a perusal of Eq. (45) in
this limit. An oscillatory decaying behaviour of the power spectrum, but with
a much sharper cutoff due to the quantum pressure contribution, has been
previously noted in the case of (non-interacting) FDM [7, 15] (which would
amount here to the $f\to 1$ limit): this is shown here for comparison by the
dot-dashed line. Interestingly, although for small $g>0$ values the power
spectrum lies above that of the FDM (for which $g=0$), further increasing the
repulsive interaction strength can lead to the deviation from CDM occurring at
wavenumbers even smaller than those set by the FDM’s ($g=0$) quantum pressure
– an example of this for a value of $g=3\cdot 10^{-87}Jm^{3}$ is shown by the
orange curve in the left plot in figure 2.

Figure 2:  Role of interaction strength $g$ on the power spectrum plotted as a
function of $k$ in logarithmic scale for a redshift $z=3$ and a reference
boson mass $m=2\cdot 10^{-22}\ eV/c^{2}$ in the absence of pressure
($\kappa=0$) for the cases of purely particle ($f\to 0$, left plot) and purely
condensate ($f\to 1$, right plot) . The established limiting cases of CDM
(solid black line, top) – relevant in the $f\to 0$ case – and FDM (dot-dashed
black line, bottom) – correspondingly relevant in the $f\to 1$ limit – are
shown for comparison in both plots. In each case we see a shift of the power
spectrum to lower values of $k$ without any identified lower bounds, and the
emergence of an oscillatory behaviour, with the spectrum bounded from above by
the CDM ($f\to 0$ case) and FDM ($f\to 1$ case) respectively. Left: For $f\to
0$ the power spectrum of the interacting case is oscillatory at large $k$ with
a relatively slow decrease in amplitude. Right: for $f\to 1$, the limiting
behaviour at the interacting model at high wavenumbers is that of FDM with a
relatively sharp cutoff and no short scale oscillations.  Figure 3:  Left:
Typical power spectrum example, similar to Fig. 2, but now in the simultaneous
presence of both condensed and non-condensed components, shown here for the
specific case of $f=0.6$ (and, still, $\kappa=0$), $m=2\cdot 10^{-22}eV/c^{2}$
and $z=3$. While the generic features discussed previously at smaller $k$
still persist, we see here that for sufficiently high $g$ (orange curve), the
hybrid power spectrum can even exceed, at high $k$, the CDM spectrum. Right:
Corresponding dimensionless power spectrum as function of $k$ in logarithmic
scale. The green line defines the limit between the linear and non-linear
regimes. Interestingly for sufficiently high $g$ (orange curve), the spectrum
even enters into the non-linear regime where the presented perturbation theory
breaks down, indicating that the co-existence of condensed and non-condensed
particles could potentially enhance smaller-scale structure (although our
present analysis cannot describe such regime which would require fully non-
linear numerical approaches).

Moving now to the opposite, $f\to 1$ limit, which corresponds to having all
bosons in the condensate, we see that it is now the FDM case that provides an
upper bound for the power spectrum of the interacting system. In other words,
in both limiting cases considered in the appropriate limits ($f\to 0$, left
image vs. $f\to 1$, right image), it is actually the relevant model (CDM, FDM)
which sets an upper bound to the power spectrum, with no lower bound provided.
For example, also in the $f\to 1$ case we see that for values of $g>3\cdot
10^{-89}Jm^{3}$, the power spectrum again deviates from the CDM line at a
wavenumber to the left of the FDM scale (orange curve in the right plot in
figure 2).

To summarize, these findings imply that (at least in the $\kappa=0$ limit) the
power spectrum predictions in either CDM or FDM provide an _upper bound_ to
the $g\neq 0$ power spectrum, with increasing self interaction driving the
power below the CDM/FDM power, adding further oscillations with increasing $k$
and moving the point of deviation from the CDM/FDM spectra to lower values of
$k$.

####  IV.2.2 Case $0<f<1$

Figure 4:  Left: Characterisation of the parameter range in $g-m$ space for
which the power spectra of hybrid models lie, for a given condensate fraction
$f$, consistently below the CDM curve for all (probed) $k$ values: the blue
shaded region identifies such regime of CDM-bounded spectra for the $f=0.6$
case of Fig. 3, with the solid blue line showcasing the locus of all those
$(m,\,g)$ values for which the power spectrum touches the CDM curve for at
least one value of $k$ (at $z=3$). The $(m,\,g)$ values corresponding to the
red and orange curves shown in Fig. 3 are indicated by the like-coloured stars
in this plot, which are respectively located below and above the blue line.
The value of $f=0.6$ has been chosen as a conservative estimate, with higher
and lower condensate fractions leading to even broader parameter range,
identified as the regions below the corresponding lines indicated by $f=0.8$
(grey dashed) and $f=0.3$ (purple dot-dashed). Right: Corresponding $g-f$
parameter region for which the hybrid model power spectrum always lies below
the CDM one, shown here by the cyan-shaded region for a fixed mass $m=2\cdot
10^{-22}eV/c^{2}$; the existence of a local minimum around a condensate
fraction of $f\simeq 0.6$, explains such choice of $f$ on the left (and other
related) plots, as for such values the enhancement in power occurs for the
smallest values of $g$. In particular, we note that the curve increases
monotonically in the limits $f\to 0$ and $f\to 1$, indicating that all power
spectra lie below the CDM ones in such limits, as expected. The blue dots
indicate the values of $f$ that we have examined and the blue line separating
the two regions of the parameter space is constructed via interpolation.

We now turn our attention to condensate fractions between the two extreme
cases. Firstly, we note that the shift of the deviation of the power spectrum
from the CDM one to smaller $k$ values with increasing $g$, and its
oscillatory features, continue in this regime, again being unbounded from
below even in the presence of coupling between the two components (condensate,
particles). Our findings indicate that, for a given $g$ (and $m$), the value
of $k$ at which the spectrum starts bending away from the corresponding CDM
one, is not significantly modified by a variation in the condensate fraction
$f$. Instead, the main effect of $f$ is to modify the initial inclination of
the curve after it deviates from the CDM one, pushing it towards the CDM case
for $f\to 0$ or towards the FDM curve for $f\to 1$. More interesting, however,
is the observation that, for a fixed mass, and for big enough values of $g$, a
condensate fraction between $0$ and $1$ leads to an enhancement in the power
spectrum for a range of small scales such that it grows beyond the CDM case,
entering quickly into the non-linear regime. To see this effect we refer the
reader to the orange curve in figure 3, where the curves correspond to a
redshift $z=3$. For illustration we demonstrate the case $f=0.6$ which we also
found to be the approximate value which makes the effect maximal (see
subsequent figure 4). Within this linearized model, there is an enhancement of
power which exceeds the CDM result in a band of wavenumbers within the
oscillatory regime. To highlight this more clearly, in the right plot of
figure 3 we show the dimensionless power spectrum $\Delta^{2}(k)$ where the
green line defines the limit of the linear regime: for linear perturbation
theory to be trustworthy the power spectrum must lie below the green line for
which $\Delta^{2}=1$. Clearly, perturbation theory cannot be valid here and a
fully non-linear treatment is required. However, the fact that the power
spectrum is enhanced into the non-linear regime for situations where standard
CDM is not would imply that non-linear structures could be formed early in the
universe at such scales for such hybrid models.

This enhancement is negligible if the condensate fractions are closer to
almost all-condensate or all-particles, showing that it is generated by the
dynamical interplay of the two coupled components. The shaded area in the left
plot of figure 4 designates the region of values for $g$ where the hybrid
curve is always bounded by the CDM curve for a range of boson masses between
$10^{-22}\ eV/c^{2}$ to $10^{-21}\ eV/c^{2}$, considering $f=0.6$, close to
which we have found the enhancement to be strongest. For comparison, the
purple and gray dashed lines indicate where the limit of such a region would
be for $f=0.3$ and $f=0.8$ respectively, showing that any possible bounds on
$g$ are always above the ones found for $f=0.6$, as we point out in the right
plot of figure 4.

###  IV.3 Role of Particle-Component Pressure ($\kappa\neq 0$) with/without
Self-Interactions

The particle component may have a non-zero pressure, stemming from velocity
dispersion, that we have modelled here with the equation of state (30). We now
turn to the effects of such a non-zero pressure, discussing it in the non-
interacting limit and for both repulsive and attractive self-interactions.

####  IV.3.1 Non-Interacting Limit ($g=0$) with $\kappa\neq 0$

To isolate the effect of the particle pressure, we first discuss the power
spectra arising in the non-interacting limit ($g=0$). We first highlight the
limiting cases $f\to 0$ and $f\to 1$, which are shown in Fig. 5 for a mass of
$m=2\cdot 10^{-22}\ eV/c^{2}$. The left plot shows that increasing particle
pressure in the $f\to 0$ limit has a qualitatively similar effect as was
previously found (Fig. 2) for increasing interaction, i.e. a decrease of the
power spectrum compared to the CDM limit (solid black line) at high $k$, and
an oscillatory behaviour, with higher $\kappa$ amplifying the effect and
enhancing the presence of oscillations. The similarity between the effects of
varying $g$ and varying $\kappa$ in the $f\to 0$ limit is expected from a
perusal of Eq. (45).333The terms proportional to these two parameters are not
identical however: since they carry different factors of $a$, the $g$ term is
more dominant than the $\kappa$ term at earlier times. We leave a more
detailed study of the redshift dependence of perturbations in our model for
future work. In particular, for larger values of $\kappa$, the spectrum
deviates from the corresponding CDM one at a lower $k$ value than the
corresponding FDM ($g=0$, $\kappa=0$, dot-dashed line) prediction – this can
be seen by the green curve in Fig. 5 for a value of $\kappa=3\cdot 10^{21}\ \
{m^{4}}{s^{-2}kg^{-2/3}}$. Since the $\kappa$ terms become more dominant at
later times, they will eventually lead to a steeper decay of the oscillations
at a high $k$, something we do observe for our probed values.

Figure 5:  Role of particle pressure $\kappa$ on the power spectrum plotted
for three different values of $\kappa$ (see figure key) as a function of $k$
in logarithmic scale for a redshift $z=3$ and a reference boson mass $m=2\cdot
10^{-22}\ eV/c^{2}$ for a non-interacting system ($g=0$) for the cases of
purely particle ($f\to 0$, left plot) and purely condensate ($f\to 1$, right
plot). Left: In the $f\to 0$ limit, there is a shift of the power spectrum to
lower values of $k$ without any identified lower bounds, and the emergence of
an oscillatory behaviour, with the spectrum bounded from above by the CDM
case. Right: In the $f\to 1$ limit all $\kappa\neq 0$ spectra overlap exactly
with the FDM prediction.

We stress that the inclusion of $\kappa$ has been on purely phenomenological
grounds444Such a term also simplifies the system of hydrodynamical equations
by reducing the number of independent equations that have to be solved, see
section II.3. and we do not have an a-priori value for this parameter. In
order to estimate the order of magnitude of $\kappa$ for which the relevant
terms can have an effect at the redshifts we consider here, we can directly
compare with the values of $g$ used through the text from equation (27).
Imposing that the particle-component pressure is comparable with that coming
from the self-interaction around $z=0$, we must have that

$\frac{g}{m^{2}}\tilde{\rho}\sim\kappa\tilde{\rho}^{2/3}$ (53)

Considering our reference values of $m=2\cdot 10^{-22}eV/c^{2}$, $g$ of the
order of $10^{-87}Jm^{3}$ and as reference density the local dark matter
density $0.0122\,{M_{\odot}}/{pc^{3}}$ [41], the order of magnitude of
$\kappa$ is around $10^{21}\ {m^{4}}{s^{-2}kg^{-2/3}}$, which explains the
apparently high numerical values of $\kappa$ used in our analysis. A first
principles calculation would be needed to assess whether such a value is
reasonable and, consequently, whether such a term is relevant at all.

For the case $f\to 1$ (all condensate and absence of particles) a non-zero
$\kappa$ plays no role since the particle pressure is not present in the
condensate equation. This is shown in the right plot of figure 5, giving the
same FDM curve for any choice of $\kappa$, as expected.

####  IV.3.2 Repulsive Interactions ($g>0$) with $\kappa\neq 0$

Given the qualitative similarity of an increasing $g>0$ (with $\kappa=0$)
[Fig. 3] and increasing $\kappa>0$ (with $g=0$) [Fig. 5] in the $f\to 0$
limit, it is natural to see how the two compare and – more importantly – to
examine what would be the combined effect of $g>0$ and $\kappa>0$ for a
general intermediate $0<f<1$ value. We discuss this in the context of our
previously-considered $f=0.6$ representative case.

Interestingly, the previously-discussed enhancement of the power spectrum at
large $k$ for large $g>0$ and $f=0.6$ [Fig. 3, orange curve] seems to be
largely alleviated by the addition of the particle pressure ($\kappa\neq 0$):
this is evident in Fig. 6, comparing the hybrid power spectra for $f=0.6$ with
$g=0$, $\kappa\neq 0$ (green), $g\neq 0$, $\kappa=0$ (orange) and $g\neq 0$,
$\kappa\neq 0$ (blue). We observe that the existence of pressure in the
particle component (at least modelled as in Eq. (30)) dampens the previous
excess and the spectrum is now again bounded by the CDM curve. Comparing the
blue curve to the green curve, representing a case with pressure but no self-
coupling for the particles, we see that the introduction of a non-zero $g>0$
is to enhance the perturbations on small scales, somewhat countering the
damping effects of the particle pressure.

Figure 6:  Interplay of repulsive interactions ($g>0$) and particle pressure
$\kappa>0$ on the power spectrum for $f=0.6$ for a redshift $z=3$ and a
reference boson mass $m=2\cdot 10^{-22}\ eV/c^{2}$. The addition of a non-zero
$\kappa$ term can alleviate the large $k$ interaction-induced increase of the
hybrid power spectrum beyond the corresponding CDM one (orange plot), as
evident from the blue curve.

####  IV.3.3 Attractive Interactions ($g<0$) with $\kappa\neq 0$

So far we have considered a positive value of $g$. However, an attractive
self-interaction ($g<0$) can appear in axion contexts and already studied for
example in [27, 36, 42, 43, 44]. In particular, it was shown in [43] that even
a tiny attractive self-interaction could be important in the stability of
structures, and numerical studies have been carried out recently [44].

Figure 7:  Interplay of particle mass, attractive interactions ($g<0$) and
particle pressure $\kappa>0$ on the power spectrum for a small condensate
fraction of $f=0.2$ at a redshift $z=3$. The presence of attractive
interactions without particle pressure (grey line) sends the system into the
nonlinear regime for sufficiently large $|g|$, significantly exceeding the CDM
line (black). The addition of increasing values of particle pressure $\kappa$
counterbalances this effect (blue and orange lines with increasing $|g|$ and
$\kappa\neq 0$). Remarkably, the combined effect of modest attractive
interaction ($g=-1.57\cdot 10^{-87}Jm^{3}$) and particle pressure
($\kappa=8.6\cdot 10^{21}\ {m^{4}}{s^{-2}kg^{-2/3}}$) for a particle of mass
$m=5\cdot 10^{-22}\ eV/c^{2}$ (orange curve) leads to effectively the same
power spectrum as that of FDM (for which $g=\kappa=0$) but for a particle mass
one order of magnitude higher ($m=5\cdot 10^{-21}\ eV/c^{2}$, purple dashed
line).

As an example, we consider a value of $g=-1.57\cdot 10^{-87}Jm^{3}$,
corresponding to one of the values for the scattering length used in [44] and
a mass of $m=5\cdot 10^{-22}\ eV/c^{2}$, along with a condensate fraction of
$f=0.2$. As pointed out in [43], negative values of $g$ for a pure condensate
case ($f=1$) can lead to instabilities and we observe this in our analysis
even in the hybrid case: see the rapid growth of power at small scales in fig.
7 (grey line) which implies a quick transition to the non-linear regime at
high $k$. This would of course call for a full non-linear treatment before any
valid conclusions could be reached as the presented linearized analysis is
obviously no longer valid. With the addition of the particle pressure however,
values of $\kappa$ can be chosen so that the destabilising effect of a
negative $g$ can be counteracted, keeping the spectrum in the linear regime,
as indicated by the orange curve in fig 7.

Interestingly, adjusting the particle pressure can result in a curve that is
very close to a typical FDM curve but for a _heavier_ mass of $m=5\cdot
10^{-21}\ eV/c^{2}$. With this in mind, the existence of these two competing
effects can lead to a _parameter degeneracy_ that could bypass mass
constraints as those presented in [45], where the boson mass is estimated as
$m>10^{-21}\ eV/c^{2}$ with Ly-$\alpha$ forest data, but without the
consideration of self-interactions and the existence of a particle component,
or the higher bounds using Ly-$\alpha$ obtained in [46], where $m>2\cdot
10^{-20}\ eV/c^{2}$. In any case, our analysis here shows that the particle
pressure would allow to have bosonic masses below such bounds but with similar
behaviour to FDM models of higher masses. Finally, it is interesting to note
that the balance between a negative self-coupling and the pressure in the
particle sector can also lead to localised increases of power that could push
some scales into the non-linear regime earlier - see the blue line in fig. 7 .

##  V Conclusions

In this work we have examined the behaviour of linear density perturbations in
a model of dark matter that consists of both a condensed, coherent part,
governed by the Gross-Pitaevskii equation, as well as a non-condensed,
incoherent part described by a Boltzmann equation. The system is self-
gravitating in an expanding universe with Hubble and density parameters set at
the Planck satellite values, and is also self-interacting with a quartic
interaction with a self-coupling constant $g$. The full equations for such a
hybrid model were derived from first principles from the non-relativistic,
bosonic action in [24], and our present work analysed the hydrodynamical
limit, including mean-field effects in a self-consistent manner to order $g$,
while ignoring collision terms of order $g^{2}$. Moreover, in order to close
the Boltzmann moment equations, we made some simplifying assumptions. We
further assumed that the pressure due to velocity dispersion in the particle
component can be modeled via (30) which allowed us to further drop one of the
resulting hydrodynamical equations.

This hybrid model reproduces CDM and fuzzy dark matter (FDM, with $g=0$) in
the respective limits $f\to 0$ and $f\to 1$ of the condensate fraction $f$.
Inclusion of a repulsive self-interaction ($g>0$) introduces damped
oscillations on small scales in the particle-only limit. Similar oscillations
exist for $f\to 1$ for sufficiently large $g$ but the behaviour tends to that
of FDM, with a sharp cutoff for small $g$. An interesting behaviour arises for
mixtures where both components are present and comparable: a repulsive self-
coupling can couple the perturbations in the two components leading to
localized enhancement in the power spectrum that can even push perturbations
into the non-linear regime at times earlier than those expected from $\Lambda
CDM$. However, such enhancements can be moderated by the inclusion of pressure
in the non-condensed, particle component, at least if this pressure is modeled
by (30). The strength of such effects naturally depends on the choice of model
parameters and a better understanding of any such closure conditions for the
Boltzmann hierarchy is required before any definite conclusion can be drawn on
the existence, importance and impact of $\kappa\neq 0$ terms.

We also considered an attractive self-coupling ($g<0$) and observed that it
induces instability and perturbation growth at small scales. Such
instabilities for pure condensate states were already pointed out in [43] and
we observe similar behaviour even at a relatively low condensate fraction of
$f=0.2$. The addition of pressure on the incoherent particles, modelled here
by $\kappa\neq 0$, can counteract such instability. An interesting consequence
is that the combined effect of the self interaction and the pressure can
result in a spectrum that resembles that of pure FDM ($f=1$, $g=0$,
$\kappa=0$) but of a _higher boson mass_ which, for the values of interaction
and pressure presumed here, could amount to an increase of approximately one
order of magnitude. Hence, mixtures of condensed and particle components, like
the ones described here, introduce degeneracies that should be kept in mind
when placing observational limits on the boson mass. Furthermore, particle
pressure can result in localized features of enhanced power at specific
wavenumbers, depending on the choice of parameters. The result of such
enhancements can only be assessed via fully non-linear hybrid simulations
based on our full model of Eqs. (1)-(3), which is left to future work.

Bosons can behave both as incoherent collections of particles and also as
coherent waves. This work examined hybrid states where both behaviours can
coexist and showed that bosonic dark matter could have a rather rich
phenomenology in terms of the matter power spectrum, not necessarily limited
to the wavy effects of what is usually referred to as fuzzy dark natter. This
phenomenology depends in turn on whether such states can be realized in the
universe and, ultimately, the dark matter production mechanism for which we
can only be agnostic at the moment. It is worth noting here that our model
bears some similarities to models of mixed dark matter that have started being
investigated in the literature [31, 47] but with the important difference that
these works consider different particle species, with correspondingly
different masses. In our case, self-interactions play an important role and
they are naturally inherited by both components which are now coupled
directly, beyond the mere effect of gravity. It is such coupling that gives
rise to the phenomenology we discuss here. Formally however, the equations are
very similar and one might further envisage couplings even between different
species. Such models would have comparable phenomenology, depending on the
form and values of these couplings. Finally, when the phenomenology we
examined here is taken into account, it would affect any parameter estimation
from observations, not least due to the existence of degeneracies as we have
alluded to above. Further investigation is therefore required, especially in
understanding the corpuscular, collisional component, before observations can
safely place limits on generic models of bosonic dark matter.

## Acknowledgements

This work is supported by the Leverhulme Trust, Grant no. RPG-2021-010.

## References

  * [1] C. S. Frenk and S. D. M. White.  Dark matter and cosmic structure.  Annalen der Physik, 524(9-10):507–534, October 2012.
  * [2] Wayne Hu, Rennan Barkana, and Andrei Gruzinov.  Cold and fuzzy dark matter.  Phys. Rev. Lett., 85:1158–1161, 2000.
  * [3] David J. E. Marsh and Pedro G. Ferreira.  Ultra-Light Scalar Fields and the Growth of Structure in the Universe.  Phys. Rev. D, 82:103528, 2010.
  * [4] David J. E. Marsh and Joe Silk.  A Model For Halo Formation With Axion Mixed Dark Matter.  Mon. Not. Roy. Astron. Soc., 437(3):2652–2663, 2014.
  * [5] Hsi-Yu Schive, Tzihong Chiueh, and Tom Broadhurst.  Cosmic structure as the quantum interference of a coherent dark wave.  Nature Physics, 10(7):496–499, July 2014.
  * [6] Hsi-Yu Schive, Ming-Hsuan Liao, Tak-Pong Woo, Shing-Kwong Wong, Tzihong Chiueh, Tom Broadhurst, and W. Y. Pauchy Hwang.  Understanding the Core-Halo Relation of Quantum Wave Dark Matter from 3D Simulations.  Phys. Rev. Lett., 113(26):261302, 2014.
  * [7] Renée Hlozek, Daniel Grin, David J. E. Marsh, and Pedro G. Ferreira.  A search for ultralight axions using precision cosmological data.  Phys. Rev. D, 91(10):103512, 2015.
  * [8] Hsi-Yu Schive, Tzihong Chiueh, Tom Broadhurst, and Kuan-Wei Huang.  Contrasting Galaxy Formation from Quantum Wave Dark Matter, $\psi$DM, with $\Lambda$CDM, using Planck and Hubble Data.  Astrophys. J., 818(1):89, 2016.
  * [9] Philip Mocz, Mark Vogelsberger, Victor H. Robles, Jesús Zavala, Michael Boylan-Kolchin, Anastasia Fialkov, and Lars Hernquist.  Galaxy formation with BECDM – I. Turbulence and relaxation of idealized haloes.  Mon. Not. Roy. Astron. Soc., 471(4):4559–4570, 2017.
  * [10] Jan Veltmaat, Bodo Schwabe, and Jens C. Niemeyer.  Baryon-driven growth of solitonic cores in fuzzy dark matter halos.  Phys. Rev. D, 101(8):083518, 2020.
  * [11] Xiaolong Du, Bodo Schwabe, Jens C. Niemeyer, and David Bürger.  Tidal disruption of fuzzy dark matter subhalo cores.  Phys. Rev. D, 97(6):063507, 2018.
  * [12] Mattia Mina, David F. Mota, and Hans A. Winther.  Solitons in the dark: First approach to non-linear structure formation with fuzzy dark matter.  Astron. Astrophys., 662:A29, 2022.
  * [13] Alex Laguë, J. Richard Bond, Renée Hložek, David J. E. Marsh, and Laurin Söding.  Evolving ultralight scalars into non-linearity with Lagrangian perturbation theory.  Mon. Not. Roy. Astron. Soc., 504(2):2391–2404, 2021.
  * [14] S. T. H. Hartman, H. A. Winther, and D. F. Mota.  Cosmological simulations of self-interacting Bose-Einstein condensate dark matter.  Astron. Astrophys., 666:A95, 2022.
  * [15] David J. E. Marsh.  Axion Cosmology.  Phys. Rept., 643:1–79, 2016.
  * [16] Lam Hui.  Wave Dark Matter.  Ann. Rev. Astron. Astrophys., 59:247–289, 2021.
  * [17] Elisa G. M. Ferreira.  Ultra-light dark matter.  Astron. Astrophys. Rev., 29(1):7, 2021.
  * [18] B. Moore.  Evidence against dissipationless dark matter from observations of galaxy haloes.  Nature, 370:629, 1994.
  * [19] Michael Boylan-Kolchin, James S. Bullock, and Manoj Kaplinghat.  Too big to fail? The puzzling darkness of massive Milky Way subhaloes.  Mon. Not. Roy. Astron. Soc., 415:L40, 2011.
  * [20] Michael Boylan-Kolchin, James S. Bullock, and Manoj Kaplinghat.  The Milky Way’s bright satellites as an apparent failure of LCDM.  Mon. Not. Roy. Astron. Soc., 422:1203–1218, 2012.
  * [21] G Kauffmann, Simon D. M. White, and B. Guiderdoni.  The Formation and Evolution of Galaxies Within Merging Dark Matter Haloes.  Mon. Not. Roy. Astron. Soc., 264:201, 1993.
  * [22] I-Kang Liu, Nick P. Proukakis, and Gerasimos Rigopoulos.  Coherent and incoherent structures in fuzzy dark matter haloes.  Mon. Not. Roy. Astron. Soc., 521(3):3625–3647, 2023.
  * [23] Robert Brandenberger, Vahid Kamali, and Rudnei O. Ramos.  Decay of ALP condensates via gravitation-induced resonance.  JCAP, 11:009, 2023.
  * [24] Nick P. Proukakis, Gerasimos Rigopoulos, and Alex Soto.  Unified description of corpuscular and fuzzy bosonic dark matter.  Phys. Rev. D, 108(8):083513, 2023.
  * [25] Wen-Yuan Ai, Ankit Beniwal, Angelo Maggi, and David J. E. Marsh.  From QFT to Boltzmann: Freeze-in in the presence of oscillating condensates.  10 2023.
  * [26] Abril Suarez and Tonatiuh Matos.  Structure Formation with Scalar Field Dark Matter: The Fluid Approach.  Mon. Not. Roy. Astron. Soc., 416:87, 2011.
  * [27] Pierre-Henri Chavanis.  Mass-radius relation of newtonian self-gravitating bose-einstein condensates with short-range interactions. i. analytical results.  Physical Review D, 84(4), aug 2011.
  * [28] Kyung-Jin Ahn and Paul R. Shapiro.  Formation and evolution of the self-interacting dark matter halos.  Mon. Not. Roy. Astron. Soc., 363:1092–1124, 2005.
  * [29] Nick P Proukakis and Brian Jackson.  Finite-temperature models of bose–einstein condensation.  J. Phys. B: At. Mol. Opt. Phys, 41(20):203002, 2008.
  * [30] Allan Griffin, Tetsuro Nikuni, and Eugene Zaremba.  Bose-Condensed Gases at Finite Temperatures.  Cambridge University Press, 2009.
  * [31] Bodo Schwabe, Mateja Gosenca, Christoph Behrens, Jens C. Niemeyer, and Richard Easther.  Simulating mixed fuzzy and cold dark matter.  Phys. Rev. D, 102(8):083518, 2020.
  * [32] N. Aghanim et al.  Planck 2018 results. VI. Cosmological parameters.  Astron. Astrophys., 641:A6, 2020.  [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [33] James M. Bardeen, J. R. Bond, Nick Kaiser, and A. S. Szalay.  The Statistics of Peaks of Gaussian Random Fields.  Astrophys. J., 304:15–61, 1986.
  * [34] Jiajun Zhang, Hantao Liu, and Ming-Chung Chu.  Cosmological Simulation for Fuzzy Dark Matter Model.  Front. Astron. Space Sci., 5:48, 2019.
  * [35] Mihir Kulkarni and Jeremiah P. Ostriker.  What is the halo mass function in a fuzzy dark matter cosmology?  Mon. Not. Roy. Astron. Soc., 510(1):1425–1430, 2021.
  * [36] Pierre-Henri Chavanis and Luca Delfini.  Mass-radius relation of newtonian self-gravitating bose-einstein condensates with short-range interactions. II. numerical results.  Physical Review D, 84(4), aug 2011.
  * [37] Pierre-Henri Chavanis.  Growth of perturbations in an expanding universe with Bose-Einstein condensate dark matter.  Astron. Astrophys., 537:A127, 2012.
  * [38] JiJi Fan.  Ultralight Repulsive Dark Matter and BEC.  Phys. Dark Univ., 14:84–94, 2016.
  * [39] Abril Suárez and Pierre-Henri Chavanis.  Cosmological evolution of a complex scalar field with repulsive or attractive self-interaction.  Phys. Rev. D, 95(6):063515, 2017.
  * [40] Vicente Delgado and Antonio Muñoz Mateo.  Self-interacting superfluid dark matter droplets.  Mon. Not. Roy. Astron. Soc., 518(3):4064–4072, 2022.
  * [41] S. Sivertsson, H. Silverwood, J. I. Read, G. Bertone, and P. Steger.  The localdark matter density from SDSS-SEGUE G-dwarfs.  Mon. Not. Roy. Astron. Soc., 478(2):1677–1693, 2018.
  * [42] Pierre-Henri Chavanis.  Collapse of a self-gravitating Bose-Einstein condensate with attractive self-interaction.  Phys. Rev. D, 94(8):083007, 2016.
  * [43] Vincent Desjacques, Alex Kehagias, and Antonio Riotto.  Impact of ultralight axion self-interactions on the large scale structure of the Universe.  Phys. Rev. D, 97(2):023529, 2018.
  * [44] Philip Mocz et al.  Cosmological structure formation and soliton phase transition in fuzzy dark matter with axion self-interactions.  Mon. Not. Roy. Astron. Soc., 521(2):2608–2615, 2023.
  * [45] Vid Iršič, Matteo Viel, Martin G. Haehnelt, James S. Bolton, and George D. Becker.  First constraints on fuzzy dark matter from Lyman-$\alpha$ forest data and hydrodynamical simulations.  Phys. Rev. Lett., 119(3):031302, 2017.
  * [46] Keir K. Rogers and Hiranya V. Peiris.  Strong Bound on Canonical Ultralight Axion Dark Matter from the Lyman-Alpha Forest.  Phys. Rev. Lett., 126(7):071302, 2021.
  * [47] Sophie M. L. Vogt, David J. E. Marsh, and Alex Laguë.  Improved mixed dark matter halo model for ultralight axions.  Phys. Rev. D, 107(6):063526, 2023.