# Unified description of corpuscular and fuzzy bosonic dark matter II:
Dissipation and stochastic forces

Nick P. Proukakis1<EMAIL_ADDRESS>Gerasimos Rigopoulos1
<EMAIL_ADDRESS>Alex Soto1<EMAIL_ADDRESS>1
School of Mathematics, Statistics and Physics, Newcastle University, Newcastle
upon Tyne, NE1 7RU, UK

(January 2024)

###### Abstract

We extend our previous work (Proukakis et al., Phys. Rev. D 108, 083513 (2023)
[1]) on the dynamics of bosonic, non-relativistic and self-interacting dark
matter that simultaneously contains both a “fuzzy” low-momentum component and
one with higher momenta that may be well approximated as a collection of
distinct particles and described by a corresponding phase-space distribution.
Starting from the non-relativistic Schwinger-Keldysh action and working beyond
leading-order in the Keldysh basis fields, encoding stochastic fluctuations of
the slow modes and all fluctuations of the fast modes, we obtain stochastic
self-consistently coupled Gross-Pitaevskii, collisional Boltzmann kinetic and
Poisson equations. Our final set of equations, which feature various
collisional (dissipative and scattering) contributions and two corresponding
independent stochastic force terms, are consistent with generalized
fluctuation-dissipation type relations in the limit of thermal equilibrium
between the particles.

##  I Introduction

Over the last decade or so an alternative idea for the microphysics of dark
matter has been rising in popularity among a growing number of cosmologists -
see [2, 3, 4, 5, 6] for recent reviews and references to the literature.
Already introduced more than 2 decades ago [7] and often called fuzzy dark
matter (FDM) or $\psi$DM, the model assumes that dark matter is composed of
ultra-light bosons that can exhibit wave phenomena on galactic scales
(indicatively $\lesssim$ 1 kpc). Interestingly, in addition to reproducing the
successful predictions of the corpuscular model of Cold Dark Matter (CDM) at
early times and on larger scales, such a model – and its extension to include
self-interactions – does in fact present distinct phenomenology on short
scales [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
26]. Mathematically, this model is normally described by the Schrödinger-
Poisson or Gross-Pitaevskii-Poisson systems of coupled equations, essentially
a wave equation that includes the (Newtonian) self-Gravity of the classical
field representing the bosons - see [27, 28, 29, 30, 31, 32, 23, 25, 26, 24,
33] for increasingly sophisticated numerical simulations of self-gravitating
classical fields. However, matter made up of bosonic particles can fall
anywhere between the limiting descriptions of a completely continuous field
(wavelike) state and a particle-like one, depending on the underlying velocity
distribution. Even if a wavelike description is ultimately justified for all
modes, due to high phase space occupancy
$\frac{\rho}{m}\lambda_{\mathrm{dB}}^{3}=\frac{h^{3}}{m^{4}}\frac{\langle\rho\rangle}{\langle|{v}|\rangle^{3}}\gtrsim\mathcal{O}(1)$,
a situation expected in cosmologically relevant configurations if the bosons
are light enough, it may nonetheless be convenient or desirable to describe a
sub-population of the particles using an N-body phase space description. For
example this would be the case if there is a large range of velocities and
therefore the related de Broglie wavelengths span a correspondingly large
range without a single peak around a dominating characteristic value - see
[34] for an application of such a configuration in the case of cosmological
perturbation theory, and similarly [35, 36] for situations where a mixed
configuration has been studied.

In our earlier work [1], to which we also refer the reader for a more extended
discussion of the context of these computations and the relation to other
known formalisms [2, 27, 28, 29, 30, 31, 32, 8, 11, 12, 14, 27, 37, 35, 13,
38, 10, 23, 9, 25, 36, 33, 24, 26, 34, 39, 40, 41, 42, 43, 44], we derived a
set of coupled equations that describe the slow modes via a Gross-Pitaevskii
equation and the fast modes via a phase space kinetic equation, also including
their mutual gravitational attraction. In this work we extend our earlier
computations to include fluctuations of the fast modes and obtain equations
that go beyond the leading semi-classical dynamics, finding a novel set of two
‘dissipative/collisional’ kernels and two associated stochastic noise forces,
one additive and one multiplicative, which are related to the fluctuations of
the high momentum particles. The respective noises and dissipative/collisional
kernels come in pairs, each exhibiting the expected structure in terms of the
fast modes’ phase space distribution. We do not assume thermal equilibrium and
the distribution of the higher momentum particles is dynamically evolved via a
collisional Boltzmann equation. However, if a thermal Bose-Einstein
distribution is assumed for the fast modes, fluctuation-dissipation-like
relations between the stochastic noises and the dissipative/collisional
kernels are retrieved.

The next section presents the derived stochastic-dissipative dynamical
equations, while section III describes the steps necessary for their
derivation, with various intermediate technical computations relegated to
Appendices. Section IV briefly presents the emergence of fluctuation-
dissipation relations in the thermal limit, with some final thoughts and
pointers to possible future work presented in section V.

##  II Stochastic - Dissipative dynamical equations for condensate and
particles

Let us start by considering the relativistic covariant action for an
ultralight bosonic particle coupled minimally to gravity, given by

$S=\int d^{4}x\sqrt{-\bar{g}}\bigg{(}\frac{1}{16\pi
G}R-\frac{1}{2}\bar{g}^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-m^{2}\phi^{2}-\frac{\lambda}{3}\phi^{4}\bigg{)}\;,$
(1)

where $G$ is the gravitational constant, $R$ is the Ricci scalar containing
the information of the system’s geometry, and in consequence the gravity,
being a function of the metric $\bar{g}_{\mu\nu}$. The value $\bar{g}$ that
appears in the action is just the determinant of this metric. We have used the
shorthand $d^{4}x=dtd^{3}\mathbf{r}$ which we will also employ throughout the
paper and we work in natural units ($c=\hbar=1$). The scalar field $\phi$
represents a bosonic particle of mass $m$ and we consider a self-interaction
with coupling $\lambda$. If we assume a flat spacetime metric with a small
scalar perturbation describing weak gravity in the Newtonian gauge, we have

$ds^{2}=-(1+2V)dt^{2}+(1-2V)\delta_{ij}dx^{i}dx^{j}\;,$ (2)

where $V<<1$ will describe the typical Newtonian gravitational potential in
the non-relativistic limit. Considering the action up to second order in $V$,
redefining $\lambda=\frac{m^{2}}{3}g$, writing

$\phi=\frac{1}{\sqrt{2m}}\left(\psi e^{-imt}+\psi^{\ast}e^{imt}\right)$ (3)

and going to the non-relativistic limit, we obtain

$S=\int
dtd^{3}\mathbf{r}\left(i\psi^{\ast}\dot{\psi}+\frac{1}{2m}\psi^{\ast}\nabla^{2}\psi-\frac{g}{2}|\psi|^{4}-mV|\psi|^{2}+\frac{1}{8\pi
G}V\nabla^{2}V\right)\;.$ (4)

For studying cosmic large scale structures it would be desirable to also
consider the expansion of the Universe. In such case, we can start from a FLRW
metric with a scalar perturbation $V$ in the Newtonian gauge, defined by

$ds^{2}=-(1+2V)dt^{2}+a^{2}(1-2V)\delta_{ij}dx^{i}dx^{j}\;,$ (5)

where $a(t)$ is the scale factor of the universe. Following the same steps as
before we can arrive at the action

$S=\int dtd^{3}\mathbf{r}\
a^{3}\bigg{(}i\psi^{\ast}\dot{\psi}+\frac{3}{2}iH|\psi|^{2}+\frac{1}{2ma^{2}}\psi^{\ast}\nabla^{2}\psi-\frac{g}{2}|\psi|^{4}-mV|\psi|^{2}-\frac{1}{8\pi
Ga^{2}}(\nabla V)^{2}\bigg{)}\;.$ (6)

where now $\mathbf{r}$ refers to comoving coordinates. Here we note the
presence of an extra term, $\frac{3}{2}iH|\psi|^{2}$, where $H=\dot{a}/a$ is
the Hubble parameter. It is convenient to make a redefinition to the _comoving
field_ $\psi\to a^{-3/2}\psi$ and thus obtain

$S=\int
dtd^{3}\mathbf{r}\left(i\psi^{\ast}\dot{\psi}+\frac{1}{2ma^{2}}\psi^{\ast}\nabla^{2}\psi-\frac{g}{2a^{3}}|\psi|^{4}-mV|\psi|^{2}-\frac{a}{8\pi
G}(\nabla V)^{2}\right)\;,$ (7)

which is the same action as that used in our previous work [1]. It is
important to remark that we must require a regime where the scale factor $a$
varies sufficiently slowly in order to be able to apply the same formalism
within this paper.

From now on we will work with action (7) with an external potential $V_{ext}$
added

$S=\int
dtd^{3}\mathbf{r}\left(i\psi^{\ast}\dot{\psi}+\frac{1}{2ma^{2}}\psi^{\ast}\nabla^{2}\psi-\frac{g}{2a^{3}}(\psi^{\ast}\psi)^{2}-V_{ext}\psi^{\ast}\psi-
mV\psi^{\ast}\psi+\frac{a}{8\pi G}V\nabla^{2}V\right)\>.$ (8)

While such an external potential is usually zero in typical cosmological
contexts we nonetheless include it here, as it is possible to add external
point particle potentials such as those stemming from black holes residing in
the centres of halos [45]. Note that the parameter $m$ attached to the second-
to-last term now plays the role of the strength of the coupling of the boson
with the field $V$. We will consider this parameter small, in the same way as
we will consider the parameter $g$, since we will work with both in a
perturbative treatment.

Note that by using a Hubbard-Stratonovich transformation in the term involving
the gravitational interaction, the non-relativistic action (8) is equivalent
to

$S=\int dt\bigg{(}i\int d^{3}\mathbf{r}\ \psi^{\ast}\dot{\psi}-H\bigg{)}$ (9)

where the the Hamiltonian $H$ has the following form

$H=\int
d^{3}\mathbf{r}\left(-\frac{1}{2ma^{2}}\psi^{\ast}\nabla^{2}\psi+V_{ext}\psi^{\ast}\psi\right)+\frac{1}{2}\int
d^{3}\mathbf{r}\int
d^{3}\mathbf{r}^{\prime}\psi^{\ast}(t,\mathbf{r})\psi^{\ast}(t,\mathbf{r}^{\prime})U(\mathbf{r},\mathbf{r}^{\prime})\psi(t,\mathbf{r}^{\prime})\psi(t,\mathbf{r})\;.$
(10)

The term $U(\mathbf{r},\mathbf{r}^{\prime})$ describes the bosonic interaction
potential, which can be written in terms of two parts, in the form

$U(\mathbf{r},\mathbf{r}^{\prime})=g\frac{\delta(\mathbf{r}-\mathbf{r}^{\prime})}{a^{3}}-\frac{Gm^{2}}{a|\mathbf{r}-\mathbf{r}^{\prime}|}\;,$
(11)

where the second term denotes the usual gravitational potential. The first
term denotes a contact interaction, stemming from from the $\phi^{4}$ self-
interaction term and parameterized by the coupling constant $g$. The
incorporation of such an interaction term is very well motivated. For example,
attractive self-interactions appear in theoretical considerations for light
CP-odd axions [46, 47, 48, 49], while particle models with repulsive self-
interactions have been also constructed [50]. The interest in self-interacting
models has increased and some constraints have been placed on the coupling
[15, 16, 19, 22, 9].

For the convenience of the readers not interested in following all the
technical steps, we present below the equations that consitute the main
results of this paper. The steps required to obtain them are explained in the
next section with further technical details relegated to the appendices.
Working in the Schwinger-Keldysh formalism with action (8), we split the field
into slow ($\Phi_{0}$) and fast ($\varphi$) components

$\psi=\Phi_{0}+\varphi$ (12)

and derive a phase space density $f(\mathbf{r},\mathbf{p},t)$ from the
(appropriate Keldysh component of the) two point correlator
$\langle\varphi^{\ast}(x)\varphi(x^{\prime})\rangle$. We will also be
referring to $\Phi_{0}$ and $\varphi$ as the “coherent”, or “condensate”, and
“incoherent”, or “particle”, parts respectively. Their defining difference is
that $\Phi_{0}$ and $\varphi$ will be associated with relatively slow and fast
moving particles respectively.

Within a set of approximations and assumptions detailed in the subsequent
section and appendices, the slow field $\Phi_{0}$ and $f$ are found to satisfy
a set of self-consistent equations which generalise our study in [1]:

$\displaystyle\bullet\quad$ $\displaystyle
i\frac{\partial\Phi_{0}(x)}{\partial
t}=\left(-\frac{1}{2ma^{2}}\nabla^{2}+V_{ext}(x)+V_{c}(x)\right)\Phi_{0}(x)$
(13c) $\displaystyle\hskip 42.67912pt-iR(x)\,\,\Phi_{0}(x)+\xi_{1}(x)$
$\displaystyle\hskip 42.67912pt-2g\int
d^{4}x^{\prime}\Pi^{R}(x^{\prime},x)V_{nc}(x^{\prime})\,\,\Phi_{0}(x)+g\Phi_{0}(x)\,\,\xi_{2}(x)$
$\displaystyle\bullet\quad$ $\displaystyle\frac{\partial f}{\partial
t}+\frac{\mathbf{p}}{ma^{2}}\cdot\nabla
f-\nabla\bigg{(}V_{ext}(x)+V_{nc}(x)\bigg{)}\cdot\nabla_{\mathbf{p}}f=\frac{1}{2}(I_{a}+I_{b})$
(14) $\displaystyle\bullet\quad$ $\displaystyle\frac{1}{4\pi
G}\nabla^{2}V^{cl}(x)=\frac{m}{a}(n_{c}(x)+\tilde{n}(x))-ma^{2}\int
d^{4}x^{\prime}\Pi^{R}(x^{\prime},x)V_{nc}(x^{\prime})+\frac{1}{2}ma^{2}\,\xi_{2}(x)$
(15)

where the mean field potentials for the coherent and non-coherent parts are
respectively

$\displaystyle V_{c}(x)$ $\displaystyle=$ $\displaystyle
mV^{cl}(x)+\frac{g}{a^{3}}(n_{c}(x)+2\tilde{n}(x))$ (16) $\displaystyle
V_{nc}(x)$ $\displaystyle=$ $\displaystyle
mV^{cl}(x)+2\frac{g}{a^{3}}(n_{c}(x)+\tilde{n}(x))$ (17)

and the number densities for the coherent and non-coherent parts are

$n_{c}=|\Phi_{0}|^{2},\quad\tilde{n}=\int\frac{d^{3}p}{(2\pi)^{3}}f(x,\mathbf{p})\,.$
(18)

The terms on the right hand side of (14) correspond to collisional terms of
this Boltzmann equation and are given by

$\displaystyle I_{a}$ $\displaystyle=$ $\displaystyle
4\frac{g^{2}}{a^{6}}n_{c}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{2}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{p}_{2}-\mathbf{p}_{1}-\mathbf{q}+\mathbf{p}_{3})$
$\displaystyle\hskip
28.45274pt\times(\delta(\mathbf{p}_{1}-\mathbf{p})-\delta(\mathbf{p}_{2}-\mathbf{p})-\delta(\mathbf{p}_{3}-\mathbf{p}))((1+f_{1})f_{2}f_{3}-f_{1}(1+f_{2})(1+f_{3}))$
$\displaystyle I_{b}$ $\displaystyle=$ $\displaystyle
4\frac{g^{2}}{a^{6}}\int\frac{d^{3}p_{2}d^{3}p_{3}d^{3}p_{4}}{(2\pi)^{5}}\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})$
(19) $\displaystyle\hskip
28.45274pt\times[f_{3}f_{4}(f+1)(f_{2}+1)-ff_{2}(f_{3}+1)(f_{4}+1)]\,.$

The two functions $R$ and $\Pi^{R}$ appearing within (13c)-(13c) are
respectively

$\displaystyle R$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{a^{6}}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{5}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{q}+\mathbf{p}_{1}-\mathbf{p}_{3}-\mathbf{p}_{2})$
(20) $\displaystyle\hskip
28.45274pt\times\bigg{[}f_{1}(1+f_{2})(1+f_{3})-(1+f_{1})f_{2}f_{3}\bigg{]}$
$\displaystyle=$
$\displaystyle\frac{1}{4n_{c}}\int\frac{d^{3}p}{(2\pi)^{3}}I_{a}$

and

$\Pi^{R}(x,\mathbf{k})=\frac{1}{a^{3}}\int\frac{d^{3}p_{1}d^{3}p_{2}}{(2\pi)^{3}}\frac{1}{\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{1}}+i\sigma}\delta(\mathbf{k}+\mathbf{p}_{2}-\mathbf{p}_{1})\bigg{[}f_{1}(1+f_{2})-f_{2}(1+f_{1})\bigg{]}$
(21)

where the $i\sigma$ terms arise from integration in the complex plane to
obtain the correct causal structure for the retarded object. The quantities
$\xi_{1}$ and $\xi_{2}$ are Gaussian stochastic forces with correlation
functions

$\displaystyle\langle\xi_{1}^{\ast}(x)\xi_{1}(x^{\prime})\rangle$
$\displaystyle=$
$\displaystyle\frac{i}{2}\Sigma_{(c)}^{K}(x)\delta(x-x^{\prime})$ (22)
$\displaystyle\langle\xi_{2}(x)\xi_{2}(x^{\prime})\rangle$ $\displaystyle=$
$\displaystyle-2i\Pi^{K}(x,x^{\prime})$ (23)

where

$\displaystyle\Sigma_{(c)}^{K}(x)$ $\displaystyle=$
$\displaystyle-2i\frac{g^{2}}{a^{6}}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{5}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{q}+\mathbf{p}_{1}-\mathbf{p}_{2}-\mathbf{p}_{3})$
(24) $\displaystyle\hskip
28.45274pt\times\bigg{[}f_{1}(1+f_{2})(1+f_{3})+(1+f_{1})f_{2}f_{3}\bigg{]}$

and

$\displaystyle\Pi^{K}(x,\mathbf{k})$ $\displaystyle=$
$\displaystyle\frac{i}{a^{3}}\int\frac{d^{3}p_{1}d^{3}p_{2}}{(2\pi)^{2}}\delta(\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{1}})\delta(\mathbf{k}+\mathbf{p}_{2}-\mathbf{p}_{1})\bigg{[}f_{1}(1+f_{2})+f_{2}(1+f_{1})\bigg{]}\;.$
(25)

An important point to highlight here is the presence of two independent noise
terms, $\xi_{1}$ and $\xi_{2}$, in the dynamical equation for the coherent
part (13c)-(13c), which are directly related to the two $g^{2}$-order
dynamical terms $-iR$ and $-2g\int
d^{4}x^{\prime}\Pi^{R}(x^{\prime},x)V_{nc}(x^{\prime})$ respectively. Looking
at the correlations of these noise terms [Eqs. (22)-(23)], we highlight that
the first one ($\xi_{1}$) is a complex noise term whose strength is set by the
local $\Sigma_{(c)}^{K}(x)$, while the second term ($\xi_{2}$) is a real noise
term whose strength is set by the non-local $\Pi^{K}(x,x^{\prime})$. A quick
inspection of their defining equations allows us to observe that $\Sigma^{K}$
is “dual” to $R$ while $\Pi^{K}$ is “dual” to $\Pi^{R}$, as expounded upon
below. In both these cases, the evolution of the local value of $\Phi_{0}$
emerges through the dynamical balancing of the rates of two competing
collisional/scattering processes involving a $\Phi_{0}$ particle.

Let us first consider the $R$ term: Its appearance in the form
$i\partial\Phi_{0}/\partial t=-iR\Phi_{0}$ implies that such term can directly
change the norm of the coherent part, $\Phi_{0}$. In fact, looking at (20), we
see that, at any given time, the local value of $\Phi_{0}$ is set by the
cumulative difference between the local growth terms leading to the creation,
and the competing decay terms leading to the loss of coherent $\Phi_{0}$
particles through energy- and momentum-conserving binary collisions involving
incoherent particles, via a collisional process of the schematic form
$\varphi+\varphi\leftrightarrow\Phi_{0}+\varphi$. The corresponding noise term
amplitude, $\xi_{1}$, representing fluctuations in such a collisional
processes, is thus set by the sum of such incoming and outgoing contributions:
this correspondence can be easily seen by a comparison of the integrands
appearing within (20) and (24). We highlight the fact that the presence of the
fluctuations $\xi_{1}$ in (13c) implies that the coherent component can grow
dynamically from initial fluctuations, even for an initial $\Phi_{0}=0$. As a
result, such equations can provide useful insight into the question of the
formation of coherent galactic solitonic cores [51, 10, 38], to be reported
elsewhere. Since this collisional process leads to a change in the number of
$\Phi_{0}$ – and given that the total particle number (i.e. the sum of
coherent and incoherent particles) is conserved within our formalism – this
term thus also contributes to a directed collisional integral, $I_{a}$ in the
Boltzmann equation (14), through the relation $R=(1/4n_{c})\int
d^{3}p/(2\pi)^{3}\,I_{a}$.

Contrary to the above, inspection of (13c) and the form of (21) shows that the
process generating the noise term $\xi_{2}$ only changes the local value of
$\Phi_{0}$, but in a manner which independently preserves $\int
d^{3}x\,n_{c}(x)$ and $\int d^{3}x\,\tilde{n}(x)$. As such, this $O(g^{2})$
term in (13c) does not in fact lead to an explicit additional collisional
integral in the Boltzmann equation (14), although locally $f$ is modified
through the $V_{nc}(x)$ term appearing on its left-hand-side. Such a process
does however lead to a change in the energy and momentum of a $\varphi$
particle, associated with the related
$f_{1}(1+f_{2})-(1+f_{1})f_{2}=(f_{1}-f_{2})$ difference factor in (21). The
strength of the noise term associated with such a process is set, through
(25), by the sum of the two factors, i.e. $f_{1}(1+f_{2})+(1+f_{1})f_{2}$. The
process giving rise to such terms can in fact be traced to a ‘scattering’
process involving one coherent and one incoherent particle which exchange
energy, but independently preserve their total number, schematically depicted
as $\varphi+\Phi_{0}\leftrightarrow\varphi+\Phi_{0}$. The $\Pi^{R}$
contribution appears in the formalism accompanied by the non-condensate mean
field potential, $V_{nc}(x)$ [Eq. (17)], as evident from (13c), (15).

Moreover, the $\Pi^{R}$ and $\Pi^{K}$ terms also contribute directly to the
Poisson equation, which now also acquires an explicit stochastic noise term
associated with such ‘scattering’ processes between coherent and incoherent
subsets of the total particle populations comprising the total gravitating
mass of the system.

Interestingly – and in addition to all the non-stochastic models arising as
limiting cases of our theory as previously discussed in [1] – the above set of
stochastic equations can be viewed as a formal methodological extension of
related stochastic models in different physical systems. We briefly mention
two relevant cases here: Firstly, in the context of ultracold atomic gases
[52] (where there is no time dependent scale factor, i.e. $a\rightarrow 1$),
and based on distinct formalisms, subsets of (13c), (13c), (13c) and (14) have
been previously derived: Specifically, a coupled stochastic Gross-Pitaevskii-
Boltzmann system of equations (13c), (13c), (14) was formulated in [53, 54,
55], but without the inclusion of the ‘scattering’ terms of (13c) related to
$\Pi^{R}$ and $\xi_{2}$. Such scattering terms were in fact explicitly
included in [56, 57, 58, 59] (from where the ‘scattering’ terminology has been
borrowed), but such equation appeared under the assumption of a thermal
equilibrium (and the further approximation of the Bose-Einstein distribution
by a Rayleigh-Jeans distribution for highly-occupied modes), thus dropping the
Boltzmann equation (14), while simultaneously treating (13c) in an approximate
manner. Such a relation is discussed in more detail in our accompanying work
focussing directly on ultracold gases [60]. Moreover, equations (14), (15)
have been discussed in the context of collisionless plasma physics, but only
in the presence of incoherent particles, i.e. setting $\Phi_{0}=n_{c}=0$ and
$I_{a}=I_{b}=0$ in them [61].

##  III Step-by-step Recipe for obtaining the dynamical equations

In this section we describe in more detail how equations (13c), (13c), (13c),
(14), (15) can be derived, starting from the action (8). To aid readability,
we focus on mapping out the basic steps and relegate some technical details to
the appendices. For convenience, we will work in flat space to make the
construction simpler, and at the end we will make use of appropriate
transformations to express our results for an expanding universe.

###  III.1 Step I: Schwinger-Keldysh action

We will work in the Schwinger-Keldysh formalism [62, 63, 64, 53, 65, 61],
where the evolution of a non-equilibrium density matrix involves a doubling of
the Hilbert space, with one half of the degrees of freedom evolving forward in
time and the other half evolving backwards. In practice, we double the field
degrees of freedom in our starting action, denoting them by $\psi^{+}$ and
$\psi^{-}$, depending on whether they are on the forward or backward contour
respectively, and further satisfying $\psi^{+}(t_{f})=\psi^{-}(t_{f})$ at some
final time $t_{f}$ in the far future. The same is done for the field $V$, and
we write the corresponding Schwinger-Keldysh action as
$S=S[\psi^{+},V^{+}]-S[\psi^{-},V^{-}]$, where the minus sign comes from the
fact that we are going backwards in our temporal integration. With this, our
Schwinger-Keldysh action explicitly reads

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
dtd^{3}\mathbf{r}\bigg{(}i\psi^{\ast+}\dot{\psi}^{+}+\frac{1}{2m}\psi^{\ast+}\nabla^{2}\psi^{+}-\frac{g}{2}(\psi^{+\ast}\psi^{+})^{2}-V_{ext}\psi^{+\ast}\psi^{+}+\frac{1}{8\pi
G}V^{+}\nabla^{2}V^{+}-mV^{+}\psi^{+\ast}\psi^{+}$ (26)
$\displaystyle-i\psi^{\ast-}\dot{\psi}^{-}-\frac{1}{2m}\psi^{\ast-}\nabla^{2}\psi^{-}+\frac{g}{2}(\psi^{-\ast}\psi^{-})^{2}+V_{ext}\psi^{-\ast}\psi^{-}-\frac{1}{8\pi
G}V^{-}\nabla^{2}V^{-}+mV^{-}\psi^{-\ast}\psi^{-}\bigg{)}$

We then use a Keldysh rotation, defining

$\displaystyle\psi^{cl}=\frac{\psi^{+}+\psi^{-}}{\sqrt{2}}\,,\quad\psi^{q}=\frac{\psi^{+}-\psi^{-}}{\sqrt{2}}$
(27) $\displaystyle V^{cl}=\frac{V^{+}+V^{-}}{2}\,,\quad
V^{q}=\frac{V^{+}-V^{-}}{2}$ (28)

where the indices $cl$ and $q$ stands for the so called “classical” and
“quantum” components in this formalism, the latter signifying fluctuations
around the on-shell solution. With these transformations the action reads

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
dtd^{3}\mathbf{r}\bigg{(}\psi^{q\ast}\left(i\partial_{t}+\frac{1}{2m}\nabla^{2}-V_{ext}\right)\psi^{cl}+\psi^{cl\ast}\left(i\partial_{t}+\frac{1}{2m}\nabla^{2}-V_{ext}\right)\psi^{q}$
(29)
$\displaystyle-\frac{g}{2}\bigg{(}\psi^{q\ast}(|\psi^{cl}|^{2}+|\psi^{q}|^{2})\psi^{cl}+\psi^{cl\ast}(|\psi^{cl}|^{2}+|\psi^{q}|^{2})\psi^{q}\bigg{)}+\frac{1}{4\pi
G}\left(V^{q}\nabla^{2}V^{cl}+V^{cl}\nabla^{2}V^{q}\right)$ $\displaystyle-
mV^{cl}(\psi^{q\ast}\psi^{cl}+\psi^{cl\ast}\psi^{q})-mV^{q}(\psi^{cl\ast}\psi^{cl}+\psi^{q\ast}\psi^{q})\bigg{)}$

###  III.2 Step II: Split the bosonic field into a slow and fast components

Following the basic premise of our work, i.e. that the we can separate fast-
moving from slow-moving modes, we split the classical and quantum fields into
a slow part (that can be taken to represent modes belonging to the coherent
part) and a fast one (the non-coherent part or the “particles”) as111At this
point it is important to remark that here the splitting in the bosonic field
is arbitrary, and means that the field $\Phi$ will contain the condensate
phase, defined through the Penrose-Onsager mode [66], plus some other
components which are non-condensed but exhibit a high degree of coherence,
similar to the condensate. This higher degree of coherence differentiates them
from the modes in $\varphi$, which we will define as highly incoherent.

$\phi^{cl}=\Phi_{0}+\varphi,\qquad\phi^{q}=\Phi^{q}+\varphi^{q}\,.$ (30)

We will further work in an approximation where we consider the gravitational
potential $V$ and the external potential $V_{ext}$ as slow varying quantities
only. The most important condition that we will highlight for our slow fields
$\Phi$ and $V$ is that we will keep terms only up to second order in their
quantum components $\Phi^{q}$ and $V^{q}$: the first order terms will provide
a semiclassical description of the coherent part $\Phi_{0}$ and the potential
$V$, while the next order will describe stochastic fluctuations. In contrast,
we will keep all higher orders in $\varphi^{q}$, the quantum component of the
field’s fast part, tracking all its fluctuations. This will be necessary as we
will be integrating out the the fast fields to obtain effective equations for
the slow fields.

With such considerations, the action (29) now reads

$S=S_{0(fast)}+S_{0(slow)}+S_{I}$ (31)

where

$\displaystyle S_{0(fast)}$ $\displaystyle=$ $\displaystyle\int
d^{4}x\left(\begin{array}[]{cc}\varphi^{\ast}\>,&\varphi^{q\ast}\end{array}\right)\left(\begin{array}[]{cc}0&i\partial_{t}+\frac{1}{2m}\nabla^{2}-V_{ext}\\\
i\partial_{t}+\frac{1}{2m}\nabla^{2}-V_{ext}&0\end{array}\right)\left(\begin{array}[]{c}\varphi\\\
\varphi^{q}\end{array}\right)$ (37) $\displaystyle S_{0(slow)}$
$\displaystyle=$ $\displaystyle\int
d^{4}x\left(\begin{array}[]{cc}\Phi_{0}^{\ast}\>,&\Phi^{q\ast}\end{array}\right)\left(\begin{array}[]{cc}-mV^{q}&i\partial_{t}+\frac{1}{2m}\nabla^{2}-V_{ext}-V_{c}\\\
i\partial_{t}+\frac{1}{2m}\nabla^{2}-V_{ext}-V_{c}&0\end{array}\right)\left(\begin{array}[]{c}\Phi_{0}\\\
\Phi^{q}\end{array}\right)$ (50) $\displaystyle+\frac{1}{4\pi G}\int
d^{4}x\left(\begin{array}[]{cc}V^{cl}\>,&V^{q}\end{array}\right)\left(\begin{array}[]{cc}0&\nabla^{2}\\\
\nabla^{2}&0\end{array}\right)\left(\begin{array}[]{c}V^{cl}\\\
V^{q}\end{array}\right)$ $\displaystyle-m\int
d^{4}xV^{q}\langle\varphi^{\ast}\varphi\rangle$

where we have defined a mean field potential $V_{c}$ for the slow component
$\Phi$ as

$V_{c}(x)=mV^{cl}(x)+g\bigg{(}\frac{1}{2}|\Phi_{0}(x)|^{2}+\langle\varphi^{\ast}(x)\varphi(x)\rangle\bigg{)}$
(51)

Also, the interaction action $S_{I}$ is defined by two pieces:

$S_{I}=S_{g}+S_{m}$ (52)

where

$S_{g}=S_{g}^{(2)}+S_{g}^{(3)}+S_{g}^{(4)}$ (53)

with

$\displaystyle S_{g}^{(2)}$ $\displaystyle=$ $\displaystyle-\frac{g}{2}\int
d^{4}x\bigg{(}2\Phi_{0}\Phi_{0}^{\ast}\varphi\varphi^{q\ast}+2\Phi_{0}\Phi^{q\ast}\varphi\varphi^{\ast}+2\Phi^{q\ast}\Phi_{0}\varphi^{q\ast}\varphi^{q}+2\Phi^{q\ast}\Phi^{q}\varphi^{q\ast}\varphi+2\Phi_{0}^{\ast}\Phi^{q}\varphi^{\ast}\varphi+2\Phi_{0}^{\ast}\Phi_{0}\varphi^{\ast}\varphi^{q}$
(54)
$\displaystyle+2\Phi^{q}\Phi^{q\ast}\varphi^{q}\varphi^{\ast}+2\Phi^{q}\Phi_{0}^{\ast}\varphi^{q}\varphi^{q\ast}-2\Phi^{q\ast}\langle\varphi^{\ast}\varphi\rangle\Phi_{0}-2\Phi_{0}^{\ast}\langle\varphi^{\ast}\varphi\rangle\Phi^{q}\bigg{)}$
$\displaystyle S_{g}^{(3)}$ $\displaystyle=$ $\displaystyle-\frac{g}{2}\int
d^{4}x\bigg{(}2\Phi_{0}\varphi\varphi^{q\ast}\varphi^{\ast}+\Phi_{0}^{\ast}\varphi^{q\ast}\varphi\varphi+\Phi^{q\ast}\varphi^{\ast}\varphi\varphi+2\Phi^{q\ast}\varphi^{q\ast}\varphi^{q}\varphi+\Phi_{0}\varphi^{q\ast}\varphi^{q\ast}\varphi^{q}+\Phi^{q}\varphi^{q\ast}\varphi^{q\ast}\varphi$
(55)
$\displaystyle+2\Phi_{0}^{\ast}\varphi^{\ast}\varphi\varphi^{q}+\Phi^{q}\varphi^{\ast}\varphi^{\ast}\varphi+\Phi_{0}\varphi^{\ast}\varphi^{\ast}\varphi^{q}+2\Phi^{q}\varphi^{q}\varphi^{\ast}\varphi^{q\ast}+\Phi^{q\ast}\varphi^{\ast}\varphi^{q}\varphi^{q}+\Phi_{0}^{\ast}\varphi^{q\ast}\varphi^{q}\varphi^{q}\bigg{)}$
$\displaystyle S_{g}^{(4)}$ $\displaystyle=$ $\displaystyle-\frac{g}{2}\int
d^{4}x\bigg{(}\varphi^{q\ast}\varphi^{\ast}\varphi\varphi+\varphi^{q\ast}\varphi^{q\ast}\varphi^{q}\varphi+\varphi^{\ast}\varphi^{\ast}\varphi\varphi^{q}+\varphi^{\ast}\varphi^{q\ast}\varphi^{q}\varphi^{q}\bigg{)}$
(56)

and where

$\displaystyle S_{m}=-m\int d^{4}x\
V^{cl}(\varphi^{q\ast}\varphi+\varphi^{\ast}\varphi^{q})-m\int d^{4}x\
V^{q}(\varphi^{\ast}\varphi-\langle\varphi^{\ast}\varphi\rangle+\varphi^{q\ast}\varphi^{q})$
(57)

Important comments are in order. First, after the splitting in slow and fast
fields using (30), we have discarded all terms that violate energy-momentum
conservation. Thus, terms with a single $\varphi$ field do not contribute.
This class of objects are: terms with one slow field $\Phi$ and one fast
$\varphi$, terms with three slow $\Phi$ and one fast $\varphi$, and terms with
one slow $\Phi$, one slow $V$ and one fast $\varphi$. Second, we have placed
the slow terms $V^{cl}$ and $|\Phi_{0}|^{2}$ inside the coherent field
hamiltonian as parts of the mean field potential $V_{c}$ since they are
composed just by classical fields. This means that combinations involving both
slow and fast fields are kept in $S_{I}$. This consideration departs from the
procedures followed in [61] and in our previous work [1]. Considering only up
to order one in the quantum fields leads to no difference in doing this.
However, this step becomes important in order to get correct physical results
at higher order in the quantum fields. Also, for later convenience that will
become apparent in the next section, we have added the term
$g\langle\varphi^{\ast}\varphi\rangle$ in the mean field potential and
subtracted it in the interacting part $S_{g}^{(2)}$. In the same way, we have
placed a term $mV^{q}\langle\varphi^{\ast}\varphi\rangle$ in the slow action
and subtracted it in $S_{g}^{(4)}$. These introduced bracketed objects are two
point functions, corresponding to the dressed propagator for the fast
$\varphi$-particles. At this point, we observe that the hamiltonian for the
incoherent (fast) $\varphi$-particles is just the free one plus the external
potential $V_{ext}$. However, we will see later that it will be modified when
we compute the corresponding self-energies.

###  III.3 Step III: Effective action for the slow fields via perturbation
theory in the generating functional

Now, we turn our attention to the generating functional

$Z=\int\mathcal{D}[\Phi V\varphi]e^{iS_{0(slow)}+iS_{0(fast)}+iS_{I}}$ (58)

which we expand in powers of the interacting part as

$\displaystyle Z\approx\int\mathcal{D}[\Phi
V]e^{iS_{0(slow)}}\left(1+\frac{\int\mathcal{D}[\varphi]e^{iS_{0(fast)}}(iS_{g}+iS_{m})}{\int\mathcal{D}[\varphi]e^{iS_{0(fast)}}}+\frac{1}{2}\frac{\int\mathcal{D}[\varphi]e^{iS_{0(fast)}}(iS_{g}+iS_{m})^{2}}{\int\mathcal{D}[\varphi]e^{iS_{0(fast)}}}\right)\int\mathcal{D}[\varphi]e^{iS_{0(fast)}}$
(59)

In the Schwinger-Keldysh formalism we have that

$\int\mathcal{D}[\varphi]e^{iS_{0(fast)}}=1$ (60)

and therefore,

$Z\approx\int\mathcal{D}[\Phi V]e^{iS_{0(slow)}}\left(1+i\langle
S_{g}\rangle+i\langle S_{m}\rangle-\frac{1}{2}\langle
S_{g}^{2}\rangle-\frac{1}{2}\langle S_{m}^{2}\rangle-\langle
S_{g}S_{m}\rangle\right)$ (61)

Due to our addition and subtraction of terms in the mean field potential and
interaction action $S_{I}$, a straightforward computation shows that $\langle
S_{g}\rangle=\langle S_{m}\rangle=0$ and as a consequence we can rewrite the
expansion as an exponential again, such that we have

$Z\approx\int\mathcal{D}[\Phi V]e^{iS_{eff}}$ (62)

where we have defined an effective action $S_{eff}$ as

$S_{eff}=S_{0(slow)}+\frac{1}{2}i\langle S_{g}^{2}\rangle+\frac{1}{2}i\langle
S_{m}^{2}\rangle+i\langle S_{g}S_{m}\rangle$ (63)

When computing the bracketed quantities, the below Schwinger-Keldysh _bare
propagators_ for the fast field make a repeated appearance

$\displaystyle\langle\varphi^{\ast}(x)\varphi(x^{\prime})\rangle$
$\displaystyle=$ $\displaystyle iG^{K}_{0(nc)}(x,x^{\prime})$ (64)
$\displaystyle\langle\varphi^{q\ast}(x)\varphi(x^{\prime})\rangle$
$\displaystyle=$ $\displaystyle iG^{R}_{0(nc)}(x,x^{\prime})$ (65)
$\displaystyle\langle\varphi^{\ast}(x)\varphi^{q}(x^{\prime})\rangle$
$\displaystyle=$ $\displaystyle iG^{A}_{0(nc)}(x,x^{\prime})$ (66)

where $G_{0(nc)}^{R(A)}$ stands for the retarded (advanced) propagator and
$G_{0(nc)}^{K}$ corresponds to the Keldysh propagator of the fast fields. On
the other hand, we can relate these bare propagators to the _dressed
propagator_

$G_{(nc)}=\left(\begin{array}[]{cc}G_{(nc)}^{K}&G_{(nc)}^{R}\\\
G_{(nc)}^{A}&0\end{array}\right)$ (67)

through the Schwinger-Dyson equation for the two point correlation function

$G_{(nc)}=G_{0(nc)}+G_{0(nc)}\otimes\Sigma_{(nc)}\otimes G_{(nc)}$ (68)

where the product sign $\otimes$ implies appropriate convolutions in spacetime
and also multiplication in the Keldysh matrix space. The object
$\Sigma_{(nc)}(x,x^{\prime})$, also a matrix with the Keldysh structure (67),
corresponds to the self-energy of the particles and is of order $g$ and
higher. This means that for the objects in equation (63) we can approximate

$G_{0(nc)}\approx G_{(nc)}$ (69)

in our expressions. In consequence, any propagator appearing in the action
will be taken to correspond to the dressed one.

We will be using the Popov approximation [67, 68, 69] where

$\langle\varphi\varphi\rangle\simeq\langle\varphi^{\ast}\varphi^{\ast}\rangle\simeq\langle\varphi^{q}\varphi^{q}\rangle\simeq\langle\varphi^{q\ast}\varphi^{q\ast}\rangle\simeq
0$ (70)

and the following identities for the Schwinger-Keldysh propagators [61]

$\displaystyle G^{R}(x,x)+G^{A}(x,x)=0$ (71) $\displaystyle
G^{R(A)}(x,x^{\prime})G^{R(A)}(x^{\prime},x)=0$ (72) $\displaystyle
G^{R}(x,x^{\prime})G^{A}(x,x^{\prime})=0$ (73)

With the above and after a long computation the effective action reads

$\displaystyle S_{eff}$ $\displaystyle=$ $\displaystyle\int
d^{4}x\left(\begin{array}[]{cc}\Phi_{0}^{\ast}(x)\>,&\Phi^{q\ast}(x)\end{array}\right)\left(\begin{array}[]{cc}0&i\partial_{t}-H_{c}\\\
i\partial_{t}-H_{c}&0\end{array}\right)\left(\begin{array}[]{c}\Phi_{0}(x)\\\
\Phi^{q}(x)\end{array}\right)$ (86) $\displaystyle+\int
d^{4}x^{\prime}d^{4}x\left(\begin{array}[]{cc}\Phi_{0}^{\ast}(x^{\prime})\>,&\Phi^{q\ast}(x^{\prime})\end{array}\right)\left(\begin{array}[]{cc}0&-\Sigma_{(c)}^{A}(x^{\prime},x)\\\
-\Sigma_{(c)}^{R}(x,x^{\prime})&0\end{array}\right)\left(\begin{array}[]{c}\Phi_{0}(x)\\\
\Phi^{q}(x)\end{array}\right)$ $\displaystyle+\frac{1}{2\pi G}\int
d^{4}xV^{q}(x)\nabla^{2}V^{cl}(x)-m\int
d^{4}xV^{q}(x)\bigg{(}|\Phi_{0}(x)|^{2}+\langle\varphi^{\ast}(x)\varphi(x)\rangle\bigg{)}$
$\displaystyle+2m\int
d^{4}x^{\prime}d^{4}xV^{q}(x)\Pi^{R}(x^{\prime},x)V_{nc}(x^{\prime})$
$\displaystyle+\int
d^{4}x^{\prime}d^{4}x\bigg{(}g(\Phi^{q\ast}\Phi_{0}+\Phi_{0}^{\ast}\Phi^{q})(x^{\prime})+mV^{q}(x^{\prime})\bigg{)}\Pi^{K}(x^{\prime},x)\bigg{(}g(\Phi^{q\ast}\Phi_{0}+\Phi_{0}^{\ast}\Phi^{q})(x)+mV^{q}(x)\bigg{)}$
$\displaystyle-\int
d^{4}x^{\prime}d^{4}x\Phi^{q\ast}(x^{\prime})\Sigma_{(c)}^{K}(x^{\prime},x)\Phi^{q}(x)$

where we have denoted by

$H_{c}=-\frac{1}{2m}\nabla^{2}+V_{ext}+V_{c}(x)-2g\int
d^{4}x^{\prime}\Pi^{R}(x^{\prime},x)V_{nc}(x^{\prime})$ (87)

with $V_{c}$ defined as in (51) and

$V_{nc}(x)=mV^{cl}(x)+g\bigg{(}|\Phi_{0}(x)|^{2}+\langle\varphi^{\ast}(x)\varphi(x)\rangle\bigg{)}\,.$
(88)

The sub-index $nc$ indicates here, as we will see soon, that it corresponds to
the mean-field potential for the non-coherent, fast particles. In the above
expressions we have defined

$\displaystyle\Pi^{K}(x^{\prime},x)$ $\displaystyle=$
$\displaystyle-\frac{i}{2}\bigg{(}G^{K}(x^{\prime},x)G^{K}(x,x^{\prime})+G^{A}(x^{\prime},x)G^{R}(x,x^{\prime})+G^{R}(x^{\prime},x)G^{A}(x,x^{\prime})\bigg{)}$
(89) $\displaystyle\Pi^{R}(x^{\prime},x)$ $\displaystyle=$
$\displaystyle-\frac{i}{2}\bigg{(}G^{R}(x^{\prime},x)G^{K}(x,x^{\prime})+G^{K}(x^{\prime},x)G^{A}(x,x^{\prime})\bigg{)}$
(90) $\displaystyle\Sigma_{(c)}^{K}(x^{\prime},x)$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{2}\bigg{(}G^{K}(x^{\prime},x)G^{K}(x,x^{\prime})G^{K}(x,x^{\prime})+G^{K}(x^{\prime},x)G^{R}(x,x^{\prime})G^{R}(x,x^{\prime})+G^{K}(x^{\prime},x)G^{A}(x,x^{\prime})G^{A}(x,x^{\prime})$
(91)
$\displaystyle+2G^{K}(x,x^{\prime})G^{A}(x^{\prime},x)G^{R}(x,x^{\prime})+2G^{K}(x,x^{\prime})G^{R}(x^{\prime},x)G^{A}(x,x^{\prime})\bigg{)}$
$\displaystyle\Sigma_{(c)}^{R}(x^{\prime},x)$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{2}\bigg{(}G^{R}(x^{\prime},x)G^{R}(x^{\prime},x)G^{A}(x,x^{\prime})+G^{K}(x^{\prime},x)G^{K}\left(x^{\prime},x\right)G^{A}(x,x^{\prime})$
(92)
$\displaystyle+2G^{R}(x^{\prime},x)G^{K}(x^{\prime},x)G^{K}(x,x^{\prime})\bigg{)}$
$\displaystyle\Sigma_{(c)}^{A}(x^{\prime},x)$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{2}\bigg{(}G^{R}(x^{\prime},x)G^{K}(x,x^{\prime})G^{K}(x,x^{\prime})+G^{A}(x,x^{\prime})G^{A}\left(x,x^{\prime}\right)G^{R}(x^{\prime},x)$
(93)
$\displaystyle+2G^{K}(x^{\prime},x)G^{A}(x,x^{\prime})G^{K}(x,x^{\prime})\bigg{)}$

To simplify the notation a little bit, we will use the fact that for any
Keldysh propagator $iG^{K}(x,x)=2n_{B}+1$, where $n_{B}$ is the bosonic
occupation number. Note that $|\Phi_{0}|^{2}$ could be seen as the Keldysh
propagator for the coherent part. This means that

$\int d\mathbf{r}|\Phi_{0}|^{2}=2\int d\mathbf{r}n_{c}=2N_{c}$ (94)

where $n_{c}$ corresponds to the number density of the coherent (slow) part
which, when integrated in space, gives the number of condensed particles
$N_{c}$. Here we have neglected the factor 1 in the Keldysh expression since
this term gives a vacuum energy which can be subtracted by means of a
renormalisation of the field. We also observe that the normalization of the
$\Phi_{0}$ field corresponds in the Schwinger-Keldysh formalism to twice the
number of constituent particles, To normalise to a particle number $N_{c}$, as
is usual, we will redefine $\Phi_{0}\rightarrow\sqrt{2}\Phi_{0}$ at the end of
the computations. Meanwhile, we will replace in our expressions

$n_{c}=\frac{1}{2}|\Phi_{0}|^{2}\,.$ (95)

For the non-coherent (fast) part, we use the usual way of writing any Keldysh
function as [61, 70, 71]

$G^{K}=G^{R}\circ F-F\circ G^{A}$ (96)

where the symbol $\circ$ represents convolution and $F(x,x^{\prime})$ is a
Hermitian matrix in the spacetime domain. Its Wigner transform (see (145)) is
$\tilde{F}(x,\mathbf{p})$ and here any energy dependence in $\tilde{F}$ is
disregarded assuming that it is a slower function than $G^{R}-G^{A}$ [61].
This term is related to the Keldysh propagator as

$\tilde{F}(x,\mathbf{p})=iG^{K}(x,\mathbf{p})$ (97)

and with it we can set the following relation

$\langle\varphi^{\ast}(x)\varphi(x)\rangle=iG_{(nc)}^{K}(x,x)=\underset{\mathbf{p}}{\sum}\tilde{F}(x,\mathbf{p})=\underset{\mathbf{p}}{\sum}(2f(x,\mathbf{p})+1)$
(98)

where $f$ is the distribution function associated with the occupation number
of the mode $\mathbf{p}$. Hence, we can define the non-coherent number density

$\tilde{n}=\underset{\mathbf{p}}{\sum}f$ (99)

and neglecting the last sum over $1$ in (98), since as in the coherent case it
gives a vacuum energy which can be subtracted by means of a renormalisation,
we can write

$2\tilde{n}=\langle\varphi^{\ast}(x)\varphi(x)\rangle$ (100)

In addition to this, to get rid of the non-locality of the terms containing
the coherent self-energies in equation (86), we can take a Wigner transform;
for details on that computation and the terms arising, we refer to appendix A.

The interpretation of the effective action and the resulting equations of
motion is further facilitated if the action is expressed in terms of a single
quantum field $\Phi^{q}$ or $V^{q}$. This is achieved by a Hubbard-
Stratonovich transformation of the last two lines. With these considerations,
the effective action becomes

$\displaystyle S_{eff}$ $\displaystyle=$ $\displaystyle\int
d^{4}x\left(\begin{array}[]{cc}\Phi_{0}^{\ast}&\Phi^{q\ast}\end{array}\right)\left(\begin{array}[]{cc}0&i\partial_{t}-H_{c}-iR\\\
i\partial_{t}-H_{c}+iR&0\end{array}\right)\left(\begin{array}[]{c}\Phi_{0}\\\
\Phi^{q}\end{array}\right)$ (107) $\displaystyle+\frac{1}{2\pi G}\int
d^{4}x\,\,V^{q}(x)\nabla^{2}V^{cl}(x)-2m\int
d^{4}x\,\,V^{q}(x)(n_{c}(x)+\tilde{n}(x))$ $\displaystyle+2m\int
d^{4}x^{\prime}d^{4}x\,\,V^{q}(x)\Pi^{R}(x^{\prime},x)V_{nc}(x^{\prime})$
$\displaystyle+\int
d^{4}x^{\prime}d^{4}x\,\,\xi_{1}^{\ast}(x^{\prime})(\Sigma_{(c)}^{K}(x,x^{\prime}))^{-1}\xi_{1}(x)$
$\displaystyle-\int d^{4}x\,\,(\xi_{1}^{\ast}\Phi^{q}+\Phi^{q\ast}\xi_{1})$
$\displaystyle-\frac{1}{4}\int
d^{4}x^{\prime}d^{4}x\,\,\xi_{2}(x^{\prime})(\Pi^{K}(x^{\prime},x))^{-1}\xi_{2}(x)$
$\displaystyle-g\int
d^{4}x\,\,\xi_{2}(\Phi^{q\ast}\Phi_{0}+\Phi_{0}^{\ast}\Phi^{q})-m\int
d^{4}x\,\,\xi_{2}V^{q}$

where we recall that

$H_{c}=-\frac{1}{2m}\nabla^{2}+V_{ext}(x)+V_{c}(x)-2g\int
d^{4}x^{\prime}\,\,\Pi^{R}(x^{\prime},x)V_{nc}(x^{\prime})$ (108)

and

$\displaystyle R$ $\displaystyle=$ $\displaystyle
g^{2}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{5}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{q}-\mathbf{p}_{3}+\mathbf{p}_{1}-\mathbf{p}_{2})$
(109)
$\displaystyle\times\bigg{[}f_{1}(1+f_{2})(1+f_{3})-(1+f_{1})f_{2}f_{3}\bigg{]}$
$\displaystyle\Sigma_{(c)}^{K}(x)$ $\displaystyle=$
$\displaystyle-2ig^{2}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{5}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{q}+\mathbf{p}_{1}-\mathbf{p}_{2}-\mathbf{p}_{3})$
(110)
$\displaystyle\bigg{[}f_{1}(1+f_{2})(1+f_{3})+(1+f_{1})f_{2}f_{3}\bigg{]}$
$\displaystyle\Pi^{R}(x,\mathbf{k})$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{1}d^{3}p_{2}}{(2\pi)^{3}}\frac{1}{\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{1}}+i\sigma}\delta(\mathbf{k}+\mathbf{p}_{2}-\mathbf{p}_{1})\bigg{[}f_{1}(1+f_{2})-f_{2}(1+f_{1})\bigg{]}$
(111) $\displaystyle\Pi^{K}(x,\mathbf{k})$ $\displaystyle=$ $\displaystyle
i\int\frac{d^{3}p_{1}d^{3}p_{2}}{(2\pi)^{2}}\delta(\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{1}})\delta(\mathbf{k}+\mathbf{p}_{2}-\mathbf{p}_{1})\bigg{[}f_{1}(1+f_{2})+f_{2}(1+f_{1})\bigg{]}$
(112)

where we have used the shorthand $f_{i}=f(x,\mathbf{p}_{i})$ and the details
of the computation of the last two terms $\Pi^{K}(x,\mathbf{k})$ and
$\Pi^{R}(x,\mathbf{k})$ are in the appendix B. Action (107) is the main result
of this subsection and corresponds to _stochastic and dissipative/collisional_
dynamics for the slow, coherent fields, arising from fluctuations of the fast
fields that have been integrated out.

###  III.4 Step IV: Fast mode (non-coherent) dynamics from the propagator
relations

To get the effective action (107) we integrated out the non-coherent (fast)
part, meaning that its dynamics is not available from the action anymore. In
order to get the equations that describe this component we will now focus on
the Schwinger-Dyson equation (68) and we will use the corresponding Feynman
diagrams to construct the self-energies. Multiplying this equation by
$G_{0(nc)}^{-1}$ from the left we get

$(G_{0(nc)}^{-1}-\Sigma_{(nc)})\otimes G_{(nc)}=I$ (113)

where $I$ is the unit matrix. Explicitly, we have that

$\left(\begin{array}[]{cc}0&i\partial_{t}-H_{0}-\Sigma_{0(nc)}^{A}-\Sigma_{(nc)}^{A}\\\
i\partial_{t}-H_{0}-\Sigma_{0(nc)}^{R}-\Sigma_{(nc)}^{R}&-\Sigma_{(nc)}^{K}\end{array}\right)\left(\begin{array}[]{cc}G_{(nc)}^{K}&G_{(nc)}^{R}\\\
G_{(nc)}^{A}&0\end{array}\right)=\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)$ (114)

where we have split the advanced and retarded self-energies in two pieces. The
first one, $\Sigma_{0(nc)}^{R}$ (analogously for $\Sigma_{0(nc)}^{A}$), will
be the corrections of order one in the perturbative parameters, which we
display in figure 1.

Figure 1: Non-coherent self-energy $\Sigma_{0(nc)}^{R}$ at order one in the
perturbative parameters, describing the corrections to the propagator
$G_{(nc)}^{R}$, symbolized by the continuous line in the left half and dashed
in the right half (the opposite case is $G_{(nc)}^{A}$). The left diagram is
the usual loop from the vertex $\varphi^{q\ast}\varphi^{\ast}\varphi\varphi$.
The completely continuous line in the loop corresponds to the non-condensed
dressed propagator $G_{(nc)}^{K}$. The diagrams in the center and right are
not closed since they are composed of classical slow fields. They come from
the vertices $\Phi_{0}\Phi_{0}^{\ast}\varphi\varphi^{q\ast}$ and
$V^{cl}\varphi^{q\ast}\varphi$ respectively. The double line represents the
$\Phi_{0}$ fields and the wavy one the $V^{cl}$ field.

The first element (left diagram in figure 1) is a typical loop coming from the
vertex $\varphi^{q\ast}\varphi^{\ast}\varphi\varphi$ and the other two
diagrams are not closed in loops since they are formed with classical slow
fields. The rule to deal with them is the following: the classical fields
don’t enter in loops, but they enter as corrections just as they were a source
through the expression defined by the corresponding vertex. Thus,

$\Sigma_{0(nc)}^{A}(x,x^{\prime})=\Sigma_{0(nc)}^{R}(x,x^{\prime})=\left[mV^{cl}+2g(n_{c}+\tilde{n})\right]\delta(x-x^{\prime})=V_{nc}\,\delta(x-x^{\prime})$
(115)

The second piece of self-energies, without subindex $0$, will be composed by
the higher order corrections.

The matrix equation (114) gives rise to the following relations:

$\displaystyle\int
d^{4}x^{\prime}\left[(G_{0(nc)}^{A})^{-1}(x,x^{\prime})-V_{nc}-\Sigma_{(nc)}^{A}(x,x^{\prime})\right]G_{(nc)}^{A}(x^{\prime},y)=\delta(x-y)$
(116) $\displaystyle\int
d^{4}x^{\prime}\left[(G_{0(nc)}^{R})^{-1}(x,x^{\prime})-V_{nc}-\Sigma_{(nc)}^{R}(x,x^{\prime})\right]G_{(nc)}^{R}(x^{\prime},y)=\delta(x-y)$
(117) $\displaystyle\int
d^{4}x^{\prime}\left[(G_{0(nc)}^{R})^{-1}(x,x^{\prime})-V_{nc}-\Sigma_{(nc)}^{R}(x,x^{\prime})\right]G_{(nc)}^{K}(x^{\prime},y)=\int
d^{4}x^{\prime}\Sigma_{(nc)}^{K}(x,x^{\prime})G_{(nc)}^{A}(x^{\prime},y)$
(118)

We use the relation (96), which relates $G_{(nc)}^{K}$ with the Hermitian
function $F$. Furthermore, from equation (98) we observed that this function
$F$ is related with the occupation number of the particles $f$, as

$F=2f+1.$ (119)

and with it and multiplying (118) from the right by
$((G_{0(nc)}^{A})^{-1}-V_{nc}-\Sigma_{(nc)}^{A})$, using (116) and (117) and
after some algebra, we arrive at

$-\left[i\partial_{t}+\frac{1}{2m}\nabla^{2}-V_{ext}(y)-V_{nc}(y),F\right]=\Sigma_{(nc)}^{K}-[\Sigma_{(nc)}^{R},F]-F(\Sigma_{(nc)}^{R}-\Sigma_{(nc)}^{A})$
(120)

After a Wigner transform in (120), and using (119), we get

$\displaystyle\left(1-\frac{\partial}{\partial\varepsilon}\Re(\Sigma_{(nc)}^{R})\right)\frac{\partial}{\partial
t}f(x,\mathbf{p})+\left(\frac{\mathbf{p}}{m}+\nabla_{\mathbf{p}}\Re(\Sigma_{(nc)}^{R})\right)\cdot\nabla
f(x,\mathbf{p})$ (121)
$\displaystyle-\nabla\bigg{(}V_{ext}+V_{nc}+\Re(\Sigma_{(nc)}^{R})\bigg{)}\cdot\nabla_{\mathbf{p}}f(x,\mathbf{p})=\frac{1}{2}I^{coll}[f]$

where

$I^{coll}[f]=i\Sigma_{(nc)}^{K}+2(2f+1)\Im(\Sigma_{(nc)}^{R})$ (122)

and here $x$ corresponds to the central point coordinates
$({x_{1}+x_{2}})/{2}$, $\mathbf{p}$ is the momentum related with the Fourier
transform of the relative coordinate $x_{1}-x_{2}$, and $\Re$, $\Im$
respectively denote real and imaginary parts. We have left the explicit
computation of the self-energies $\Sigma_{(nc)}^{R}$ and $\Sigma_{(nc)}^{K}$
for appendix C, where we have observed that the real part of
$\Sigma_{(nc)}^{R}$ is zero and as a consequence the non-coherent self-
energies only contribute to the collisional terms.

###  III.5 Step V: The full equations of the system

Now we can gather all the elements to present the equations for the system and
summarise the important equations. In the effective action, as we mentioned
before, we will make the redefinition $\Phi_{0}\to\sqrt{2}\Phi_{0}$ to get as
normalization of the field the coherent number of elements $N_{c}$. Also, for
convenience, we make the redefinition $\xi_{1}\to\sqrt{2}\xi_{1}$. After this,
we vary with respect to $\Phi^{q\ast}$ and $V^{q}$ to get the equations of
motion of the coherent field and the field $V^{cl}$ respectively. In addition,
we use the equation (121), considering the results of the self-energies
computed in appendix C. Thus, our equations are

$\displaystyle i\frac{\partial\Phi_{0}(x)}{\partial
t}=\left(-\frac{1}{2m}\nabla^{2}+V_{ext}(x)+V_{c}(x)\right)\Phi_{0}(x)$ (124)
$\displaystyle-iR(x)\Phi_{0}(x)+\xi_{1}(x)-2g\int
d^{4}x^{\prime}\Pi^{R}(x^{\prime},x)V_{nc}(x^{\prime})\Phi_{0}(x)+g\xi_{2}(x)\Phi_{0}(x)$
$\displaystyle\frac{1}{4\pi
G}\nabla^{2}V^{cl}(x)=m(n_{c}(x)+\tilde{n}(x))-m\int
d^{4}x^{\prime}\Pi^{R}(x^{\prime},x)V_{nc}(x^{\prime})+\frac{1}{2}m\xi_{2}(x)$
$\displaystyle\frac{\partial f}{\partial t}+\frac{\mathbf{p}}{m}\cdot\nabla
f-\nabla\bigg{(}V_{ext}(x)+V_{nc}(x)\bigg{)}\cdot\nabla_{\mathbf{p}}f=\frac{1}{2}(I_{a}+I_{b})$
(125)

where we recall that the mean field potentials for the coherent and non-
coherent parts are respectively

$\displaystyle V_{c}(x)$ $\displaystyle=$ $\displaystyle
mV^{cl}(x)+g(n_{c}(x)+2\tilde{n}(x))$ (126) $\displaystyle V_{nc}(x)$
$\displaystyle=$ $\displaystyle mV^{cl}(x)+2g(n_{c}(x)+\tilde{n}(x))$ (127)

and the number densities for the coherent and non-coherent part are

$n_{c}=|\Phi_{0}|^{2},\quad\tilde{n}=\int\frac{d^{3}p}{(2\pi)^{3}}f(x,\mathbf{p})\,.$
(128)

The terms in the right hand side of (125) corresponding to collisional terms
of this Boltzmann equation are given by

$\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}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{p}_{2}-\mathbf{p}_{1}-\mathbf{q}+\mathbf{p}_{3})$
$\displaystyle\times(\delta(\mathbf{p}_{1}-\mathbf{p})-\delta(\mathbf{p}_{2}-\mathbf{p})-\delta(\mathbf{p}_{3}-\mathbf{p}))((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}}\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})$
(129)
$\displaystyle\times[f_{3}f_{4}(f+1)(f_{2}+1)-ff_{2}(f_{3}+1)(f_{4}+1)]\,.$

From the explicit expression for $R$ in (109) we have that

$\displaystyle R$ $\displaystyle=$ $\displaystyle
g^{2}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{5}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{q}+\mathbf{p}_{1}-\mathbf{p}_{3}-\mathbf{p}_{2})$
(130)
$\displaystyle\times\bigg{[}f_{1}(1+f_{2})(1+f_{3})-(1+f_{1})f_{2}f_{3}\bigg{]}$
$\displaystyle=$
$\displaystyle\frac{1}{4n_{c}}\int\frac{d^{3}p}{(2\pi)^{3}}I_{a}\,.$

Also,

$\Pi^{R}(x,\mathbf{k})=\int\frac{d^{3}p_{1}d^{3}p_{2}}{(2\pi)^{3}}\frac{1}{\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{1}}+i\sigma}\delta(\mathbf{k}+\mathbf{p}_{2}-\mathbf{p}_{1})\bigg{[}f_{1}(1+f_{2})-f_{2}(1+f_{1})\bigg{]}$
(131)

and the noise terms $\xi_{1}$ and $\xi_{2}$ are Gaussian, satisfying the
correlations

$\displaystyle\langle\xi_{1}^{\ast}(x)\xi_{1}(x^{\prime})\rangle$
$\displaystyle=$
$\displaystyle\frac{i}{2}\Sigma_{(c)}^{K}(x)\delta(x-x^{\prime})$ (132)
$\displaystyle\langle\xi_{2}(x)\xi_{2}(x^{\prime})\rangle$ $\displaystyle=$
$\displaystyle-2i\Pi^{K}(x,x^{\prime})$ (133)

where we remember that

$\displaystyle\Sigma_{(c)}^{K}(x)$ $\displaystyle=$
$\displaystyle-2ig^{2}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{5}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{q}+\mathbf{p}_{1}-\mathbf{p}_{2}-\mathbf{p}_{3})$
(134)
$\displaystyle\times\bigg{[}f_{1}(1+f_{2})(1+f_{3})+(1+f_{1})f_{2}f_{3}\bigg{]}$
$\displaystyle\Pi^{K}(x,\mathbf{k})$ $\displaystyle=$ $\displaystyle
i\int\frac{d^{3}p_{1}d^{3}p_{2}}{(2\pi)^{2}}\delta(\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{1}})\delta(\mathbf{k}+\mathbf{p}_{2}-\mathbf{p}_{1})\bigg{[}f_{1}(1+f_{2})+f_{2}(1+f_{1})\bigg{]}\,.$
(135)

In order to obtain a generalization in expanding universe, we make the
following transformations: we pass the physical coordinates to comoving ones
$\mathbf{r}\to a\mathbf{r}$, which means that we transform the derivatives
$\nabla\to{a}^{-1}\nabla$. Also, we perform the change $\mathbf{p}\to
a^{-1}\mathbf{p}$ and the densities are transformed to the comoving ones
$n_{c}\to a^{-3}n_{c}$ and $\tilde{n}\to a^{-3}\tilde{n}$. The fields
appearing in the equations are then to be understood as comoving, relating to
the physical ones as $\Phi_{0}=\Phi^{\rm phys}_{0}\,a^{-3/2}$ and $f=f^{\rm
phys}a^{-3}$, consistent with the definitions used in [1]. The resulting
equations are those presented in section II and constitute the main equations
for a bosonic system composed by slow (coherent) and fast (non-coherent)
particles, interacting with an additional potential that we have encoded in
$V^{cl}$ in an expanding universe.

The set of equations generalises the result obtained in [1], where it was
shown that in the limit $g\to 0$ and when all modes are coherent, we recover
the commonly used Fuzzy Dark Matter equations, and when all mades are
incoherent, we recover the typical Vlasov-Poisson equations which are used to
describe CDM. Also, if we had worked only up to order one in the quantum
components $V^{q}$ and $\Phi^{q}$ and turning off the gravitational
interaction we would have recovered the ZNG models [41, 42], used in Cold Atom
physics to describe bosonic atoms. The incorporation up to second order in
quantum components allowed us to get the noises $\xi_{1}$ and $\xi_{2}$, which
can relate with the Stochastic Projected Gross-Pitaevskii equations [56, 57,
58, 59] in the absence of gravity, when we additionally ignore the Boltzmann
equation for the incoherent particles and we consider them as described by an
equilibrium Rayleigh-Jeans distribution, as a semi-classical approximation to
the full Bose-Einstein distribution. Alternatively, ignoring terms $\Pi^{R}$,
associated noise $\xi_{2}$ and gravity, our equations coincide with the
alternative cold atom stochastic formulation of Stoof [53, 54] and its
resulting near-equilibrium stochastic equation [55], with the relation between
such models in the cold atom context discussed in our companion paper [60].

##  IV Thermal state and fluctuation-dissipation relation

In the equations we derived above, the fast (incoherent) modes are described
via an N-body phase space distribution function $f$ which is to be solved in
conjunction with the wave-like equation for the slow (coherent) component -
see III.5. The equations involve dissipative/collisional terms and stochastic
forces where a characteristic pattern with similar factors of $f$ appears
between the kernels and the variances of the corresponding stochastic forces.
As a consistency check, we now confirm that when thermal equilibrium is
assumed for the fast modes described by $f$, standard fluctuation-dissipation
relations are recovered.

Indeed, if we take the equilibrium Bose-Einstein distribution for the fast
modes

$f(\varepsilon_{i},\mu)=\frac{1}{e^{\beta(\varepsilon_{i}-\mu)}-1}$ (136)

where $\beta=1/k_{B}T$ and $\epsilon_{i}$ is the energy of the particle $i$,
and we apply this distribution in the equations (134) and (135), we can show
that

$\displaystyle iR(x,\mathbf{q})$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\Sigma_{(c)}^{K}(x,\mathbf{q})\frac{1}{1+2f(\varepsilon_{\mathbf{q}},\mu)}$
(137) $\displaystyle\Pi^{R}(x,\mathbf{k})-\Pi^{A}(x,\mathbf{k})$
$\displaystyle=$
$\displaystyle\Pi^{K}(x,\mathbf{k})\frac{1}{1+2f(\varepsilon_{\mathbf{k}},0)}$
(138)

The function $\Pi^{A}(x,\mathbf{k})$ is obtained from $\Pi^{R}(x,\mathbf{k})$
knowing that $\Pi^{A}(x,x^{\prime})=\Pi^{R}(x^{\prime},x)$ and proceeding as
in Appendix B. Then, we use that $G^{R}(p)-G^{A}(p)=-2\pi
i\delta(\epsilon-\epsilon_{\mathbf{p}})$ to arrive to (138). In the above
expressions $\varepsilon_{\mathbf{q}}$ corresponds to the energy of the slow
particle involved in the collisional process with a non-coherent, fast
particle and $\varepsilon_{\mathbf{k}}$ is an energy with values in the
coherent (slow) band; as we describe in appendix B it doesn’t correspond to a
single particle energy but a combination of slow particles. Notice that since
we have used the Bose-Einstein distribution, we have recovered in our equation
(137) the results obtained in [53, 54] in the context of Stochastic Gross-
Pitaevskii equations.

We can simplify our expressions by the use of the function $F$, defined as
$F=2f+1$ in (98), which for the Bose-Einstein distribution reads:

$F(\varepsilon,\mu)=\coth\left(\frac{\beta}{2}(\varepsilon-\mu)\right)$ (139)

Also, we use that $\Sigma^{R}_{(c)}(x,q)=-iR$,
$\Sigma_{(c)}^{A}=\Sigma_{(c)}^{R\ast}$ (see Appendix A). With that, equations
(137) and (138) become

$\displaystyle\Sigma_{(c)}^{K}(x,\mathbf{q})$ $\displaystyle=$ $\displaystyle
F(\varepsilon_{\mathbf{q}},\mu)(\Sigma^{R}_{(c)}(x,\mathbf{q})-\Sigma^{A}_{(c)}(x,\mathbf{q}))$
(140) $\displaystyle\Pi^{K}(x,\mathbf{k})$ $\displaystyle=$ $\displaystyle
F(\varepsilon_{\mathbf{k}},0)(\Pi^{R}(x,\mathbf{k})-\Pi^{A}(x,\mathbf{k}))\,,$
(141)

recovering the standard fluctuation-dissipation theorem for bosonic particles
[61].

##  V Final remarks

In this work we have extended our earlier first-principles derivation [1] of a
set of dynamical equations for a generic dark matter model involving a scalar
particle with mass $m$ and quartic self-interaction with self-coupling
strength $g$. Starting from the non-relativistic, non-equilibrium partition
function, we describe the system in terms of two components, one involving
fast modes described by a kinetic equation for their phase space distribution
$f$, also dubbed the non-coherent part here, and one involving slower modes
with longer de-Broglie wavelengths, which are described by a wave equation,
all coupled via gravity and their local self-interaction. The new feature here
is the derivation of a set of two dissipation/collisional terms and two
corresponding stochastic noise terms that affect the coherent part’s dynamics.
These terms reflect the influence of the fluctuations of the faster modes and
arise when these faster modes are integrated out to provide effective
equations for the slow modes. The kinetic equation for $f$ is in turn obtained
from the Schwinger-Dyson Equations of the fast modes.

The new dissipation/collisional kernels and the variances of the stochastic
noise forces derived here are expressed in terms of the fast modes’ phase
space distribution and come in two pairs, depending on the number of $f$
factors appearing in the relevant expressions. Their forms ensure that the
expected Fluctuation-Dissipation relations are satisfied if the fast particles
are taken to obey the thermal Bose-Einstein distribution, confirming the
consistency of our computation. Note that while the $\xi_{2}$ noise is
multiplicative, being proportional to the coherent component, the $\xi_{1}$
noise is additive. Hence, the fast fluctuations described by the latter can
source the growth of $\Phi_{0}$ ab initio. Interestingly, $\xi_{2}$ and the
corresponding dissipation kernel also affect the gravitational potential via
additions to the Poisson equation. Note that all these terms are proportional
to the self-coupling $g$ and hence the fluctuations we capture with our
present approximations do not involve effects mediated by the gravitational
sector but are restricted to the self-coupling. Gravitational fluctuation
effects would require us to go to higher order in $V^{q}$ and, presumably,
$\Phi^{q}$, a task we leave for future investigations. It would nevertheless
be interesting to see how our dynamical equations may be used to study the
soliton nucleation process of scalar dark matter [51, 38, 10]. It would also
be interesting to compare our dynamical equations with other kinetic
approahces [72, 73] or the non-equilibrium QFT formalism as in [74], which
however does not use self-interactions, or in [75] which ustilises the 2PI
effective action. Finally, the mechanism explored in [76] may be relevant to
the existence of fast and slow components that we describe here. We leave such
interesting questions for future work.

## Acknowledgements

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

##  Appendix A Wigner transforms of the coherent self-energies in the
effective action

We observe that the terms containing the coherent self-energies
$\Sigma_{(c)}^{R}$,$\Sigma_{(c)}^{A}$ and $\Sigma_{(c)}^{K}$ in (86) have a
double integration ($d^{4}xd^{4}x^{\prime}$). However, it is possible to
reduce the double integration to a single one using a Wigner transform. We
will analyse in detail the term containing the retarded self-energy for the
coherent (slow) part. The same steps can be used for the advanced and Keldysh
self-energies, so we will pay more attention to this first case, briefly
sketching the results for the other self-energies.

###  A.1 Coherent Retarded Self-Energy $\Sigma_{(c)}^{R}$ and $R$

. The term containing the retarded self-energy for the coherent part is

$\displaystyle\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x^{\prime})(-\Sigma_{(c)}^{R}(x,x^{\prime}))\Phi_{0}(x)$
$\displaystyle=$ $\displaystyle\frac{g^{2}}{2}\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x)\bigg{(}2G_{(nc)}^{K}(x^{\prime},x)G_{(nc)}^{R}(x,x^{\prime})G_{(nc)}^{K}(x,x^{\prime})$
$\displaystyle+G_{(nc)}^{A}\left(x^{\prime},x\right)(G_{(nc)}^{K}(x,x^{\prime}))^{2}+G_{(nc)}^{A}(x^{\prime},x)(G_{(nc)}^{R}(x,x^{\prime}))^{2}\bigg{)}\Phi_{0}(x^{\prime})$

For convenience, we replace the elements $G_{(nc)}^{R}(x,x^{\prime})$ in the
last term inside the parenthesis in (A.1) by
$(G_{(nc)}^{R}(x,x^{\prime})-G_{(nc)}^{A}(x,x^{\prime}))$, since
$G_{(nc)}^{A}(x^{\prime},x)G_{(nc)}^{A}(x,x^{\prime})=0$. With this, we have
that

$\displaystyle\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x^{\prime})(-\Sigma_{(c)}^{R}(x,x^{\prime}))\Phi_{0}(x)$
$\displaystyle=$ $\displaystyle g^{2}\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x^{\prime})G_{(nc)}^{K}(x^{\prime},x)G_{(nc)}^{R}(x,x^{\prime})G_{(nc)}^{K}(x,x^{\prime})\Phi_{0}(x)$
(143) $\displaystyle+\frac{g^{2}}{2}\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x^{\prime})G_{(nc)}^{A}(x^{\prime},x)G_{(nc)}^{K}(x,x^{\prime})G_{(nc)}^{K}(x,x^{\prime})\Phi_{0}(x)$
$\displaystyle+\frac{g^{2}}{2}\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x^{\prime})G_{(nc)}^{A}(x^{\prime},x)(G_{(nc)}^{R}(x,x^{\prime})-G_{(nc)}^{A}(x,x^{\prime}))$
$\displaystyle\
\times(G_{(nc)}^{R}(x,x^{\prime})-G_{(nc)}^{A}(x,x^{\prime}))\Phi_{0}(x)$

We will only look in detail at the procedure for the integral in the first
line in (143), which is

$J=\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x^{\prime})G_{(nc)}^{K}(x^{\prime},x)G_{(nc)}^{R}(x,x^{\prime})G_{(nc)}^{K}(x,x^{\prime})\Phi_{0}(x)$
(144)

and we will repeat the same procedure for the next lines, skipping for brevity
the details. In this term we use the Wigner transform

$G(x,x^{\prime})=\underset{p}{\sum}e^{ip(x-x^{\prime})}G\left(\frac{x+x^{\prime}}{2},p\right)$
(145)

where $G$ corresponds to any function of two points $x$ and $x^{\prime}$
(remembering that $x$ and $x^{\prime}$ are space-time coordinates) and the sum
$\underset{p}{\sum}=\underset{\mathbf{p}}{\sum}\int\frac{d\epsilon}{2\pi}$,
with $\epsilon$ the energy and $\mathbf{p}$ the corresponding momentum.
Therefore,

$J=\int d^{4}x\int
d^{4}x^{\prime}\underset{p_{1},p_{2},p_{3},q}{\sum}e^{i(p_{3}+p_{2}-p_{1}-q)(x-x^{\prime})}G_{(nc)}^{K}\left(\frac{x+x^{\prime}}{2},p_{1}\right)G_{(nc)}^{R}\left(\frac{x+x^{\prime}}{2},p_{3}\right)G_{(nc)}^{K}\left(\frac{x+x^{\prime}}{2},p_{2}\right)D\left(\frac{x+x^{\prime}}{2},q\right)$
(146)

where we have defined

$D(x^{\prime},x)=\Phi^{q\ast}(x^{\prime})\Phi_{0}(x)$ (147)

noting that $\Phi^{q\ast}(x^{\prime})\Phi_{0}(x)$ is a correlation function
for the coherent part, and we have taken the Wigner transform of this object.

Now, we make the change of variables

$x_{a}=\frac{x+x^{\prime}}{2},\qquad x_{b}=x-x^{\prime}$ (148)

where we observe that $x_{a}$ corresponds to the central point coordinate and
$x_{b}$ is the relative coordinate. With these new variables:

$J=\int d^{4}x_{a}\int
d^{4}x_{b}\underset{p_{1},p_{2},p_{3},q}{\sum}e^{i(p_{3}+p_{2}-p_{1}-q)x_{b}}G_{(nc)}^{K}(x_{a},p_{1})G_{(nc)}^{R}(x_{a},p_{3})G_{(nc)}^{K}(x_{a},p_{2})D(x_{a},q)$
(149)

Doing the integration in $x_{b}$ gives a delta function that we use to
eliminate the sum in $r$. Hence,

$J=\int
d^{4}x_{a}\underset{p_{1},p_{2},q}{\sum}G_{(nc)}^{K}(x_{a},p_{1})G_{(nc)}^{R}(x_{a},p_{1}+q-p_{2})G_{(nc)}^{K}(x_{a},p_{2})D(x_{a},q)$
(150)

Now, we observe that

$\underset{q}{\sum}D(x_{a},q)=\underset{q}{\sum}e^{iq(x_{a}-x_{a})}D\left(\frac{x_{a}+x_{a}}{2},q\right)=D(x_{a},x_{a})=\Phi_{0}(x_{a})\Phi^{q\ast}(x_{a})$
(151)

With that and renaming $x_{a}\rightarrow x$ we have that

$J=\int
d^{4}x\underset{p_{1},p_{2}}{\sum}\Phi^{q\ast}(x)G_{(nc)}^{K}(x,p_{1})G_{(nc)}^{R}(x,p_{1}+q-p_{2})G_{(nc)}^{K}(x,p_{2})\Phi_{0}(x)$
(152)

Repeating the same steps for the next two terms in (143) we end up with

$\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x)(-\Sigma_{(c)}^{R}(x,x^{\prime}))\Phi_{0}(x^{\prime})=-\int
d^{4}x\Phi^{q\ast}(x)\Sigma^{R}_{(c)}(x,q)\Phi_{0}(x)$ (153)

where

$\displaystyle\Sigma^{R}_{(c)}(x,q)$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{2}\underset{p_{1},p_{2}}{\sum}\bigg{[}2G_{(nc)}^{K}(x,p_{1})G_{(nc)}^{R}(x,p_{1}+q-p_{2})G_{(nc)}^{K}(x,p_{2})+G_{(nc)}^{A}(x,p_{1})G_{(nc)}^{K}(x,p_{2})G_{(nc)}^{K}\left(x,p_{1}+q-p_{2}\right)$
(154)
$\displaystyle+G_{(nc)}^{A}(x,p_{1})(G_{(nc)}^{R}(x,p_{2})-G_{(nc)}^{A}(x,p_{2}))(G_{(nc)}^{R}(x,p_{1}+q-p_{2})-G_{(nc)}^{A}(x,p_{1}+q-p_{2}))\bigg{]}$

Now, we use [61]

$(G_{(nc)}^{R}(p)-G_{(nc)}^{A}(p))=-2\pi
i\delta(\varepsilon-\varepsilon_{\mathbf{p}}),\qquad G_{(nc)}^{K}(p)=-2\pi
i\tilde{F}(p)\delta(\varepsilon-\varepsilon_{\mathbf{p}})$ (155)

and we symmetrize the first term between $p_{2}$ and $p_{1}+q-p_{2}$. Using
these relations, we can split $\Sigma_{(c)}^{R}$ in a real and imaginary part,
and we will observe that every term in the real part contains $(G^{R}+G^{A})$,
which will be zero after the energy integration. Therefore, we are left with
just the imaginary part. Then, we do the energy integration in the imaginary
part, and using that $F=2f+1$ we get that

$\displaystyle\Im(\Sigma_{(c)}^{R})$ $\displaystyle=$ $\displaystyle-2\pi
g^{2}\underset{\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3}}{\sum}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{1}+\mathbf{q}-\mathbf{p}_{2}})\delta(\mathbf{q}+\mathbf{p}_{1}-\mathbf{p}_{2}-\mathbf{p}_{3})$
(156)
$\displaystyle\times\bigg{[}f_{1}(1+f_{2})(1+f_{3})-(1+f_{1})f_{2}f_{3}\bigg{]}$

where we used the shorthand $f_{i}=f(x,\mathbf{p}_{i})$. Naming $R$ to minus
this imaginary part and considering a thermodynamical limit, we can pass the
sum on the momenta to an integral, we get that

$\displaystyle R$ $\displaystyle=$
$\displaystyle-\Im(\Sigma_{(c)}^{R})=g^{2}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{5}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{q}+\mathbf{p}_{1}-\mathbf{p}_{2}-\mathbf{p}_{3})$
(157)
$\displaystyle\times\bigg{[}f_{1}(1+f_{2})(1+f_{3})-(1+f_{1})f_{2}f_{3}\bigg{]}$

and with that

$\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x)(-\Sigma_{(c)}^{R}(x,x^{\prime}))\Phi_{0}(x^{\prime})=\int
d^{4}x\Phi^{q\ast}(x)iR\Phi_{0}(x)$ (158)

Analogously we can do the same for the term involving
$\Sigma_{(c)}^{A}(x,x^{\prime})$, and we will get

$\int d^{4}x\int
d^{4}x^{\prime}\Phi_{0}^{\ast}(x)(-\Sigma_{(c)}^{A}(x,x^{\prime}))\Phi^{q}(x^{\prime})=-\int
d^{4}x\Phi_{0}^{\ast}(x)iR\Phi^{q}(x)$ (159)

notice that this result preserves that
$\Sigma_{(c)}^{A}=\Sigma_{(c)}^{R\ast}$.

###  A.2 Coherent Keldysh Self-Energy $\Sigma^{K}_{(c)}$

The Keldysh self-energy is given by

$\displaystyle-\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x^{\prime})\Sigma_{(c)}^{K}(x,x^{\prime})\Phi^{q}(x)$
$\displaystyle=$ $\displaystyle\frac{g^{2}}{2}\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x^{\prime})\bigg{(}G_{(nc)}^{K}(x^{\prime},x)(G_{(nc)}^{K}(x,x^{\prime}))^{2}+G_{(nc)}^{K}(x^{\prime},x)(G_{(nc)}^{R}(x,x^{\prime}))^{2}$
(160)
$\displaystyle+G_{(nc)}^{K}(x^{\prime},x)\left(G_{(nc)}^{A}\left(x,x^{\prime}\right)\right)^{2}+2G_{(nc)}^{A}(x^{\prime},x)G_{(nc)}^{K}(x,x^{\prime})G_{(nc)}^{R}(x,x^{\prime})$
$\displaystyle+2G_{(nc)}^{R}(x^{\prime},x)G_{(nc)}^{A}(x,x^{\prime})G_{(nc)}^{K}(x,x^{\prime})\bigg{)}\Phi^{q}(x)$

After using the identities

$\displaystyle(G_{(nc)}^{R}(x,x^{\prime})-G_{(nc)}^{A}(x,x^{\prime}))(G_{(nc)}^{R}(x,x^{\prime})-G_{(nc)}^{A}(x,x^{\prime}))$
$\displaystyle=$
$\displaystyle(G_{(nc)}^{R}(x,x^{\prime}))^{2}+(G_{(nc)}^{A}(x,x^{\prime}))^{2}$
(161)
$\displaystyle(G_{(nc)}^{R}(x^{\prime},x)-G_{(nc)}^{A}(x^{\prime},x))(G_{(nc)}^{R}(x,x^{\prime})-G_{(nc)}^{A}(x,x^{\prime}))$
$\displaystyle=$
$\displaystyle-G^{A}(x^{\prime},x)G_{(nc)}^{R}(x,x^{\prime})-G_{(nc)}^{R}(x^{\prime},x)G_{(nc)}^{A}(x,x^{\prime})$

and performing the Wigner transforms as we did for $\Sigma_{({c})}^{R}$ and
change of variable we can arrive to

$\int d^{4}x\int
d^{4}x^{\prime}\Phi^{q\ast}(x)(-\Sigma_{(c)}^{K}(x,x^{\prime}))\Phi^{q}(x^{\prime})=-\int
d^{4}x\Phi^{q\ast}(x)\Sigma^{K}_{(c)}(x,q)\Phi^{q}(x)$ (163)

where

$\displaystyle\Sigma_{(c)}^{K}(x,q)$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{2}\underset{p_{1},p_{2}}{\sum}\bigg{[}G_{(nc)}^{K}(x,p_{1})G_{(nc)}^{K}(x,p_{2})G_{(nc)}^{K}(x,p_{1}+q-p_{2})$
(164)
$\displaystyle+G_{(nc)}^{K}(x,p_{1})(G_{(nc)}^{R}(x,p_{2})-G_{(nc)}^{A}(x,p_{2}))\left(G_{(nc)}^{R}\left(x,p_{1}+q-p_{2}\right)-G_{(nc)}^{A}(x,p_{1}+q-p_{2})\right)$
$\displaystyle-2G_{(nc)}^{K}(x,p_{1}+q-p_{2})(G_{(nc)}^{R}(x,p_{1})-G_{(nc)}^{A}(x,p_{1}))(G_{(nc)}^{R}(x,p_{2})-G_{(nc)}^{A}(x,p_{2}))\bigg{]}$

Symmetrizing between $p_{2}$ and $p_{1}+q-p_{2}$ in the last term, using the
relations in equation (155) and the relation $F=2f+1$, doing the energy
integration and passing the sum on the momenta to integrals as before we end
up with

$\displaystyle\Sigma_{(c)}^{K}(x)$ $\displaystyle=$
$\displaystyle-2ig^{2}\int\frac{d^{3}p_{1}d^{3}p_{2}d^{3}p_{3}}{(2\pi)^{5}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{q}+\mathbf{p}_{1}-\mathbf{p}_{2}-\mathbf{p}_{3})$
(165)
$\displaystyle\times\bigg{[}f_{1}(1+f_{2})(1+f_{3})+(1+f_{1})f_{2}f_{3}\bigg{]}$

##  Appendix B Self-energies $\Pi^{K}(x^{\prime},x)$ and
$\Pi^{R}(x^{\prime},x)$

###  B.1 Term $\Pi^{K}(x^{\prime},x)$

We have one term containing $\Pi^{K}(x^{\prime},x)$, made by the combination

$\int d^{4}x^{\prime}\int
d^{4}x(g(\Phi^{q\ast}\Phi_{0}+\Phi^{q\ast}\Phi_{0})(x^{\prime})+mV^{q}(x^{\prime}))\Pi^{K}(x^{\prime},x)(g(\Phi^{q\ast}\Phi_{0}+\Phi_{0}^{\ast}\Phi^{q})(x)+mV^{q}(x))$
(166)

We observe that the combination
$g(\Phi^{q\ast}\Phi_{0}+\Phi_{0}^{\ast}\Phi^{q})(x)+mV^{q}(x)$ is real, in
contrast with the terms involving the coherent self-energies in appendix A,
which are complex. As we did in the appendix A, using the identity

$-(G_{(nc)}^{R}(x,x^{\prime})-G_{(nc)}^{A}(x,x^{\prime}))(G_{(nc)}^{R}(x^{\prime},x)-G_{(nc)}^{A}(x^{\prime},x))=G_{(nc)}^{A}(x,x^{\prime})G_{(nc)}^{R}(x^{\prime},x)+G_{(nc)}^{R}(x,x^{\prime})G_{(nc)}^{A}(x^{\prime},x)$
(167)

and the explicit expression for $\Pi(x^{\prime},x)$ in equation (89) we obtain
that the Wigner transform of this object is

$\Pi^{K}(x,\mathbf{k})=-\frac{i}{2}\underset{p_{2}}{\sum}\bigg{(}G_{(nc)}^{K}(x,p_{2})G_{(nc)}^{K}(x,p_{2}+k)-(G_{(nc)}^{R}(x,p_{2})-G_{(nc)}^{A}(x,p_{2}))(G_{(nc)}^{R}(x,p_{2}+k)-G_{(nc)}^{A}(x,p_{2}+k))\bigg{)}$
(168)

where it is important to remark that the momentum $\mathbf{k}$ associated of
the real combination
$g(\Phi^{q\ast}\Phi_{0}+\Phi_{0}^{\ast}\Phi^{q})(x)+mV^{q}(x)$ doesn’t
correspond to that of a single particle but a combination of them in a slow
band. In consequence, it is not uniquely defined. In this way, the trick used
here just allows us to compute the Wigner transform of $\Pi(x,x^{\prime})$,
but it doesn’t make the noise correlators in equation (133) local. With the
use of (155) and the relation $F=2f+1$, doing the energy integration and
passing the sum on the momenta to integrals as before, we obtain after some
arrangements

$\Pi^{K}(x,\mathbf{k})=i\int\frac{d^{3}p_{1}d^{3}p_{2}}{(2\pi)^{2}}\delta(\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{1}})\delta(\mathbf{k}+\mathbf{p}_{2}-\mathbf{p}_{1})\bigg{[}f_{1}(1+f_{2})+(1+f_{1})f_{2}\bigg{]}$
(169)

###  B.2 Term $\Pi^{R}(x^{\prime},x)$

We focus on the term

$\int d^{4}x^{\prime}\int
d^{4}xV^{q}(x^{\prime})\Pi^{R}(x,x^{\prime})V^{cl}(x)$ (170)

Using the explicit form of $\Pi^{R}(x^{\prime},x)$, found in (90) and
following the same steps as before, it is possible to arrive at

$\Pi^{R}(x,\mathbf{k})=-\frac{i}{2}\underset{p_{2}}{\sum}\bigg{(}G_{(nc)}^{R}(x,p_{2}+k)G_{(nc)}^{K}(x,p_{2})+G_{(nc)}^{K}(x,p_{2}+k)G_{(nc)}^{A}(x,p_{2})\bigg{)}$
(171)

Then, we use the equation (155) for $G_{(nc)}^{K}$, the explicit relations for
the retarded and advanced propagators [61]

$G^{R}(p)=\frac{1}{\varepsilon-\varepsilon_{\mathbf{p}}+i\sigma},\quad
G^{A}(p)=\frac{1}{\varepsilon-\varepsilon_{\mathbf{p}}-i\sigma}$ (172)

and the relation $F=2f+1$, doing the energy integration and passing the sum on
the momenta to integrals we get finally that

$\Pi^{R}(x,\mathbf{k})=\int\frac{d^{3}p_{1}d^{3}p_{2}}{(2\pi)^{3}}\frac{1}{\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{1}}+i\sigma}\delta(\mathbf{k}+\mathbf{p}_{2}-\mathbf{p}_{1})\bigg{[}f_{1}(1+f_{2})-f_{2}(1+f_{1})\bigg{]}$
(173)

##  Appendix C Particle self-energies computation and collisional terms

The interaction action $S_{I}$ in equation (52), defines through the terms in
equations (54), (55), (56) and (57) the elements for constructing the Feynman
diagrams making up the self-energies for the non-coherent particles (see
figures 2 and 3). We will compute the retarded and Keldysh self-energies. The
advanced one is obtained from the retarded one via complex conjugation.

###  C.1 Retarded Particle Self-energy $\Sigma_{(nc)}^{R}$

Figure 2: Feynman diagrams for the particle retarded self-energy at order
$g^{2}$. The three upper diagrams correspond to $\Sigma^{R}_{(a)}$ (eq.
(175)), the remaining corresponds to $\Sigma^{R}_{(b)}$ (eq. (176)). The
double lines in the upper diagrams give origin to the $|\Phi_{0}|^{2}$ term in
eq. (175).

The contributions to the retarded self-energy $\Sigma_{(nc)}^{R}$, depicted in
figure 2 are written as

$-i\Sigma_{(nc)}^{R}=-i\Sigma_{(a)}^{R}-i\Sigma_{(b)}^{R}$ (174)

where we have split it in two different pieces

$\displaystyle-i\Sigma_{(a)}^{R}$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{4}|\Phi_{0}|^{2}\bigg{[}2G_{(nc)}^{A}(x^{\prime},x)G_{(nc)}^{K}(x,x^{\prime})+2G_{(nc)}^{R}(x,x^{\prime})G_{(nc)}^{K}(x^{\prime},x)$
(175)
$\displaystyle+2G_{(nc)}^{R}(x,x^{\prime})G_{(nc)}^{K}(x,x^{\prime})\bigg{]}$
$\displaystyle-i\Sigma_{(b)}^{R}$ $\displaystyle=$ $\displaystyle
i\frac{g^{2}}{8}\bigg{[}2G_{(nc)}^{R}(x^{\prime},x)G_{(nc)}^{A}(x,x^{\prime})G_{(nc)}^{R}(x,x^{\prime})+4G_{(nc)}^{K}(x^{\prime},x)G_{(nc)}^{K}(x,x^{\prime})G_{(nc)}^{R}(x,x^{\prime})$
(176)
$\displaystyle+2G_{(nc)}^{A}(x^{\prime},x)G_{(nc)}^{K}(x,x^{\prime})G_{(nc)}^{K}(x,x^{\prime})+2G_{(nc)}^{A}(x^{\prime},x)G_{(nc)}^{A}(x,x^{\prime})G_{(nc)}^{A}(x,x^{\prime})$
$\displaystyle+2G_{(nc)}^{A}(x^{\prime},x)G_{(nc)}^{R}(x,x^{\prime})G_{(nc)}^{R}(x,x^{\prime})\bigg{]}$

We use the identities

$\displaystyle G^{A}(x,x^{\prime})^{2}+G^{R}(x,x^{\prime})^{2}$
$\displaystyle=$ $\displaystyle(G^{R}(x,x^{\prime})-G^{A}(x,x^{\prime}))^{2}$
(177) $\displaystyle G^{R}(x^{\prime},x)G^{R}(x,x^{\prime})$ $\displaystyle=$
$\displaystyle 0$ (178)

and we symmetrize the last term of $\Sigma_{(a)}^{R}$ in (175) and the second
one of $\Sigma_{(b)}^{R}$ in (176), we take a Wigner transform, similarly to
the coherent part in appendix A and we use the equation (155). Doing this, as
we did in the coherent retarded self-energy, we can split in a real and
imaginary part. Every term in the real part contains $G^{A}+G^{R}$, which
becomes zero after the energy integration. Therefore, only the imaginary part
survives and we proceed with it doing the energy integrations, using $F=2f+1$
and the fact that in the Schwinger-Keldysh $|\Phi_{0}|^{2}=2n_{c}$. Thus, we
get that the imaginary part of the two pieces reads

$\displaystyle\Im(\Sigma_{(a)}^{R})$ $\displaystyle=$ $\displaystyle 2\pi
g^{2}n_{c}\underset{\mathbf{p}}{\sum}\bigg{[}2\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}}-\varepsilon_{\mathbf{k}}-\varepsilon_{\mathbf{p}+\mathbf{q}-\mathbf{k}})(f_{p}-f_{p+q-k})$
(179)
$\displaystyle-\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{k}}-\varepsilon_{\mathbf{p}}-\varepsilon_{\mathbf{k}+\mathbf{q}-\mathbf{p}})(f_{p}+f_{k+q-p}+1)\bigg{]}$
$\displaystyle\Im(\Sigma_{(b)}^{R})$ $\displaystyle=$ $\displaystyle 4\pi
g^{2}\underset{\mathbf{p},\mathbf{q}}{\sum}\delta(\varepsilon_{\mathbf{p}}+\varepsilon_{\mathbf{k}-\mathbf{q}}-\varepsilon_{\mathbf{k}}-\varepsilon_{\mathbf{p}-\mathbf{q}})[f_{k-q}f_{p}-f_{p-q}f_{k-q}-f_{p-q}f_{p}-f_{p-q}]$
(180)

We have considered only the imaginary parts, since they are part of the
collisional term in (122). We ignored the real parts since they only will
contribute to small shifts in the energies.

###  C.2 Particle Keldysh Self-energy $\Sigma_{(nc)}^{K}$

Figure 3: Feynman diagrams for the particle Keldysh self-energy at order
$g^{2}$. The six upper diagrams correspond to $\Sigma^{K}_{(a)}$ (eq. (182)),
the remaining corresponds to $\Sigma^{K}_{(b)}$ (eq. (183)).

The contributions to the Keldysh self-energy $\Sigma_{(nc)}^{K}$ are given by

$-i\Sigma_{(nc)}^{K}=-i\Sigma_{(a)}^{K}-i\Sigma_{(b)}^{K}$ (181)

where as before, we have split them in two parts, defined as

$\displaystyle-i\Sigma_{(a)}^{K}$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{4}|\Phi_{0}|^{2}\bigg{[}2G_{(nc)}^{K}(x,x^{\prime})^{2}+4G_{(nc)}^{A}(x^{\prime},x)G_{(nc)}^{R}(x,x^{\prime})+4G_{(nc)}^{R}(x^{\prime},x)G_{(nc)}^{A}(x,x^{\prime})$
(182)
$\displaystyle+4G_{(nc)}^{K}(x,x^{\prime})G_{(nc)}^{K}(x^{\prime},x)+2G_{(nc)}^{R}(x,x^{\prime})^{2}+2G_{(nc)}^{A}(x,x^{\prime})^{2}\bigg{]}$
$\displaystyle-i\Sigma_{(b)}^{K}$ $\displaystyle=$ $\displaystyle
i\frac{g^{2}}{4}\bigg{[}2G_{(nc)}^{K}(x^{\prime},x)G^{R}(x,x^{\prime})^{2}+2G_{(nc)}^{K}(x^{\prime},x)G_{(nc)}^{K}(x,x^{\prime})^{2}$
(183)
$\displaystyle+4G_{(nc)}^{A}(x^{\prime},x)G_{(nc)}^{K}(x,x^{\prime})G_{(nc)}^{R}(x,x^{\prime})+2G_{(nc)}^{K}(x^{\prime},x)G_{(nc)}^{A}(x,x^{\prime})^{2}$
$\displaystyle+4G_{(nc)}^{R}(x^{\prime},x)G_{(nc)}^{A}(x,x^{\prime})G_{(nc)}^{K}(x,x^{\prime})\bigg{]}$

We use the identities

$\displaystyle
G^{R}(x^{\prime},x)G^{A}(x,x^{\prime})+G^{A}(x^{\prime},x)G^{R}(x,x^{\prime})$
$\displaystyle=$
$\displaystyle-(G^{R}(x^{\prime},x)-G^{A}(x^{\prime},x))(G^{R}(x,x^{\prime})-G^{A}(x,x^{\prime}))$
(184) $\displaystyle G^{A}(x,x^{\prime})^{2}+G^{R}(x,x^{\prime})^{2}$
$\displaystyle=$ $\displaystyle(G^{R}(x,x^{\prime})-G^{A}(x,x^{\prime}))^{2}$
(185)

We symmetrize the last term in $\Sigma_{(b)}^{K}$ and following the same steps
as in the previous case, we get that

$\displaystyle\Sigma_{(a)}^{K}$ $\displaystyle=$ $\displaystyle-2i\pi
g^{2}n_{c}\underset{\mathbf{p}}{\sum}\bigg{[}4\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}}-\varepsilon_{\mathbf{k}}-\varepsilon_{\mathbf{p}+\mathbf{q}-\mathbf{k}})(2f_{p+q-k}f_{p}+f_{p}+f_{p+q-k})$
(186)
$\displaystyle+2\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{k}}-\varepsilon_{\mathbf{p}}-\varepsilon_{\mathbf{k}+\mathbf{q}-\mathbf{p}})(2f_{k+q-p}f_{p}+f_{p}+f_{k+q-p}+1)\bigg{]}$
$\displaystyle\Sigma_{(b)}^{K}$ $\displaystyle=$ $\displaystyle-4i\pi
g^{2}\underset{\mathbf{p},\mathbf{q}}{\sum}\delta(\varepsilon_{\mathbf{p}}+\varepsilon_{\mathbf{k}-\mathbf{q}}-\varepsilon_{\mathbf{k}}-\varepsilon_{\mathbf{p}-\mathbf{q}})[2f_{p-q}f_{k-q}f_{p}+f_{p-q}f_{p}+f_{p-q}f_{k-q}+f_{p-q}+f_{k-q}f_{p}]$
(187)

###  C.3 Collisional terms

Using the expression in equation (122) we can generate two kind of collisional
terms, constructed with the computed retarded and Keldysh self-energies:

$\displaystyle I_{a}[f(k)]$ $\displaystyle=$ $\displaystyle
i\Sigma_{(a)}^{K}+2(2f(k)+1)\Im(\Sigma_{(a)}^{R})$ (188) $\displaystyle
I_{b}[f(k)]$ $\displaystyle=$ $\displaystyle
i\Sigma_{(b)}^{K}+2(2f(k)+1)\Im(\Sigma_{(b)}^{R})$ (189)

Explicitly, after some algebra and passing the sums to integrals, we can get
that

$\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}}\delta(\varepsilon_{\mathbf{q}}+\varepsilon_{\mathbf{p}_{1}}-\varepsilon_{\mathbf{p}_{2}}-\varepsilon_{\mathbf{p}_{3}})\delta(\mathbf{p}_{2}-\mathbf{p}_{1}-\mathbf{q}+\mathbf{p}_{3})\bigg{(}\delta(\mathbf{p}_{1}-\mathbf{p})-\delta(\mathbf{p}_{2}-\mathbf{p})$
(190)
$\displaystyle-\delta(\mathbf{p}_{3}-\mathbf{p})\bigg{)}\bigg{[}(1+f_{1})f_{2}f_{3}-f_{1}(1+f_{2})(1+f_{3})\bigg{]}$
$\displaystyle I_{b}$ $\displaystyle=$ $\displaystyle
4g^{2}\int\frac{d^{3}p_{2}d^{3}p_{3}d^{3}p_{4}}{(2\pi)^{5}}\delta(\varepsilon_{\mathbf{p}_{3}}+\varepsilon_{\mathbf{p}_{4}}-\varepsilon_{\mathbf{p}}-\varepsilon_{\mathbf{p}_{2}})\delta(\mathbf{p}_{2}-\mathbf{p}_{3}-\mathbf{p}_{4}+\mathbf{p})$
(191)
$\displaystyle\times\bigg{[}f_{3}f_{4}(f+1)(f_{2}+1)-ff_{2}(f_{3}+1)(f_{4}+1)\bigg{]}$

## References

  * [1] Nick P. Proukakis, Gerasimos Rigopoulos, and Alex Soto.  Unified description of corpuscular and fuzzy bosonic dark matter.  Phys. Rev. D, 108:083513, Oct 2023.
  * [2] David J. E. Marsh.  Axion cosmology.  Physics Reports, 643:1–79, July 2016.
  * [3] Lam Hui.  Wave Dark Matter.  Annual Review of Astronomy and Astrophysics, 59:247–289, September 2021.
  * [4] Elisa G. M. Ferreira.  Ultra-light dark matter.  The Astronomy and Astrophysics Review, 29(1):7, December 2021.
  * [5] Ciaran A. J. O’Hare.  Cosmology of axion dark matter.  PoS, COSMICWISPers:040, 2024.
  * [6] Tonatiuh Matos, Luis A. Ureña López, and Jae-Weon Lee.  Short review of the main achievements of the scalar field, fuzzy, ultralight, wave, BEC dark matter model.  Front. Astron. Space Sci., 11:1347518, 2024.
  * [7] Wayne Hu, Rennan Barkana, and Andrei Gruzinov.  Cold and fuzzy dark matter.  Physical Review Letters, 85:1158–1161, 2000.
  * [8] C. G. Boehmer and T. Harko.  Can dark matter be a Bose-Einstein condensate?  JCAP, 06:025, 2007.
  * [9] Vicente Delgado and Antonio Muñoz Mateo.  Self-interacting superfluid dark matter droplets.  Mon. Not. Roy. Astron. Soc., 518(3):4064–4072, 2022.
  * [10] Kay Kirkpatrick, Anthony E. Mirasola, and Chanda Prescod-Weinstein.  Analysis of Bose-Einstein condensation times for self-interacting scalar dark matter.  Physical Review D, 106(4):043512, 2022.
  * [11] 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), August 2011.
  * [12] 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:043532, August 2011.
  * [13] Pierre-Henri Chavanis.  Jeans Instability of Dissipative Self-Gravitating Bose–Einstein Condensates with Repulsive or Attractive Self-Interaction: Application to Dark Matter.  Universe, 6(12):226, 2020.
  * [14] Tanja Rindler-Daller and Paul R. Shapiro.  Angular Momentum and Vortex Formation in Bose-Einstein-Condensed Cold Dark Matter Haloes.  Mon. Not. Roy. Astron. Soc., 422:135–161, 2012.
  * [15] Bohua Li, Tanja Rindler-Daller, and Paul R. Shapiro.  Cosmological constraints on bose-einstein-condensed scalar field dark matter.  Phys. Rev. D, 89:083536, Apr 2014.
  * [16] Bohua Li, Paul R. Shapiro, and Tanja Rindler-Daller.  Bose-einstein-condensed scalar field dark matter and the gravitational wave background from inflation: New cosmological constraints and its detectability by ligo.  Phys. Rev. D, 96:063505, Sep 2017.
  * [17] 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:023529, Jan 2018.
  * [18] Pierre-Henri Chavanis.  Jeans mass-radius relation of self-gravitating Bose-Einstein condensates and typical parameters of the dark matter particle.  Phys. Rev. D, 103(12):123551, 2021.
  * [19] Sayan Chakrabarti, Bihag Dave, Koushik Dutta, and Gaurav Goswami.  Constraints on the mass and self-coupling of ultra-light scalar field dark matter using observational limits on galactic central mass.  JCAP, 09:074, 2022.
  * [20] Bihag Dave and Gaurav Goswami.  Self-interactions of ULDM to the rescue?  JCAP, 07:015, 2023.
  * [21] Milos Indjin, I-Kang Liu, Nick P. Proukakis, and Gerasimos Rigopoulos.  Virialized profiles and oscillations of self-interacting fuzzy dark matter solitons.  Phys. Rev. D, 109(10):103518, 2024.
  * [22] S. T. H. Hartman, H. A. Winther, and D. F. Mota.  Constraints on self-interacting Bose-Einstein condensate dark matter using large-scale observables.  JCAP, 02(02):005, 2022.
  * [23] 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.
  * [24] Matthias Stallovits and Tanja Rindler-Daller.  Single and merger soliton dynamics in scalar field dark matter with and without self-interactions.  6 2024.  arXiv:2406.07419 [astro-ph.CO].
  * [25] 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.
  * [26] Connor A. Painter, Michael Boylan-Kolchin, Philip Mocz, and Mark Vogelsberger.  An attractive model: simulating fuzzy dark matter with attractive self-interactions.  2 2024.
  * [27] 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.
  * [28] Bodo Schwabe, Jens C. Niemeyer, and Jan F. Engels.  Simulations of solitonic core mergers in ultralight axion dark matter cosmologies.  Physical Review D, 94(4):043513, August 2016.
  * [29] 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.  Monthly Notices of the Royal Astronomical Society, 471(4):4559–4570, November 2017.
  * [30] Jan Veltmaat, Jens C. Niemeyer, and Bodo Schwabe.  Formation and structure of ultralight bosonic dark matter halos.  Physical Review D, 98(4):043509, 2018.
  * [31] Simon May and Volker Springel.  Structure formation in large-volume cosmological simulations of fuzzy dark matter: impact of the non-linear dynamics.  Monthly Notices of the Royal Astronomical Society, 506(2):2603–2618, June 2021.
  * [32] Matteo Nori, Andrea V. Macciò, and Marco Baldi.  Fuzzy Aquarius: evolution of a Milky-way like system in the Fuzzy Dark Matter scenario.  10 2022.  arXiv:2210.08022 [astro-ph.CO].
  * [33] Hsun-Yeong Yang, Barry T. Chiang, Guan-Ming Su, Hsi-Yu Schive, Tzihong Chiueh, and Jeremiah P. Ostriker.  Galactic disc heating by density granulation in fuzzy dark matter simulations.  Mon. Not. Roy. Astron. Soc., 530(1):129–148, 2024.
  * [34] Nick P. Proukakis, Gerasimos Rigopoulos, and Alex Soto.  Hybrid model of condensate and particle dark matter: Linear perturbations in the hydrodynamic limit.  Phys. Rev. D, 110(2):023504, 2024.
  * [35] 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.
  * [36] Alex Laguë, Bodo Schwabe, Renée Hložek, David J. E. Marsh, and Keir K. Rogers.  Cosmological simulations of mixed ultralight dark matter.  Phys. Rev. D, 109(4):043507, 2024.
  * [37] Alan H. Guth, Mark P. Hertzberg, and C. Prescod-Weinstein.  Do Dark Matter Axions Form a Condensate with Long-Range Correlation?  Physical Review D, 92(10):103513, 2015.
  * [38] Jiajun Chen, Xiaolong Du, Erik W. Lentz, and David J. E. Marsh.  Relaxation times for Bose-Einstein condensation by self-interaction and gravity.  Physical Review D, 106(2):023009, 2022.
  * [39] F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro.  Large-scale structure of the Universe and cosmological perturbation theory.  Physics Reports, 367(1-3):1–248, September 2002.
  * [40] Raul E. Angulo and Oliver Hahn.  Large-scale dark matter simulations.  Living Reviews in Computational Astrophysics, 8(1):1, December 2022\.
  * [41] Eugene Zaremba, Tetsuro Nikuni, and Allan Griffin.  Dynamics of trapped bose gases at finite temperatures.  Journal of Low Temperature Physics, 116:277–345, 1999.
  * [42] Allan Griffin, Tetsuro Nikuni, and Eugene Zaremba.  Bose-Condensed Gases at Finite Temperatures.  Cambridge University Press, 2009.
  * [43] Nick P Proukakis and Brian Jackson.  Finite-temperature models of bose–einstein condensation.  J. Phys. B: At. Mol. Opt. Phys, 41(20):203002, 2008.
  * [44] Kean Loon Lee and Nick P Proukakis.  Non-equilibrium atomic condensates and mixtures: collective modes, condensate growth and thermalisation.  Journal of Physics B: Atomic, Molecular and Optical Physics, 49(21):214003, October 2016.
  * [45] Elliot Yarnell Davies and Philip Mocz.  Fuzzy Dark Matter Soliton Cores around Supermassive Black Holes.  Mon. Not. Roy. Astron. Soc., 492(4):5721–5729, 2020.
  * [46] Steven Weinberg.  A new light boson?  Phys. Rev. Lett., 40:223–226, Jan 1978.
  * [47] F. Wilczek.  Problem of strong $p$ and $t$ invariance in the presence of instantons.  Phys. Rev. Lett., 40:279–282, Jan 1978.
  * [48] R. D. Peccei and Helen R. Quinn.  CP Conservation in the Presence of Instantons.  Phys. Rev. Lett., 38:1440–1443, 1977.
  * [49] R. D. Peccei and Helen R. Quinn.  Constraints imposed by $\mathrm{CP}$ conservation in the presence of pseudoparticles.  Phys. Rev. D, 16:1791–1797, Sep 1977.
  * [50] JiJi Fan.  Ultralight Repulsive Dark Matter and BEC.  Phys. Dark Univ., 14:84–94, 2016.
  * [51] D. G. Levkov, A. G. Panin, and I. I. Tkachev.  Gravitational Bose-Einstein condensation in the kinetic regime.  Phys. Rev. Lett., 121(15):151301, 2018.
  * [52] Nick Proukakis, Simon Gardiner, Matthew Davis, and Marzena Szymańska, editors.  Quantum Gases: Finite Temperature and Non-Equilibrium Dynamics: 1 (Cold Atoms).  ICP, 2013.
  * [53] Henk TC Stoof.  Coherent versus incoherent dynamics during bose-einstein condensation in atomic gases.  Journal of Low Temperature Physics, 114(1):11–108, 1999.
  * [54] R. A. Duine and H. T. C. Stoof.  Stochastic dynamics of a trapped bose-einstein condensate.  Physical Review A, 65(1), dec 2001.
  * [55] H. T. C. Stoof and M. J. Bijlsma.  Dynamics of fluctuating Bose–Einstein condensates.  J. Low Temp. Phys., 124:431–442, 2001.
  * [56] C W Gardiner and M J Davis.  The stochastic gross–pitaevskii equation: Ii.  Journal of Physics B: Atomic, Molecular and Optical Physics, 36(23):4731, nov 2003.
  * [57] A. S. Bradley and C. W. Gardiner.  The stochastic gross-pitaevskii equation: III.  2006\.  arXiv:cond-mat/0602162.
  * [58] A. S. Bradley, C. W. Gardiner, and M. J. Davis.  Bose-einstein condensation from a rotating thermal cloud: Vortex nucleation and lattice formation.  Phys. Rev. A, 77:033616, Mar 2008.
  * [59] P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, and C. W. Gardiner.  Dynamics and statistical mechanics of ultra-cold Bose gases using c-field techniques.  Advances in Physics, 57(5):363–455, September 2008.
  * [60] Nick P. Proukakis, Gerasimos Rigopoulos, and Alex Soto.  (In Preparation).
  * [61] Alex Kamenev.  Field Theory of Non-Equilibrium Systems.  Cambridge University Press, 2011.
  * [62] R. D. Jordan.  Effective field equations for expectation values.  Physical Review D, 33(2):444–454, January 1986.
  * [63] E. Calzetta and B. L. Hu.  Nonequilibrium quantum fields: Closed-time-path effective action, Wigner function, and Boltzmann equation.  Physical Review D, 37(10):2878–2900, May 1988.
  * [64] H. T. C. Stoof.  Initial stages of bose-einstein condensation.  Phys. Rev. Lett., 78:768–771, February 1997.
  * [65] Jørgen Rammer.  Quantum Field Theory of Non-equilibrium States.  Cambridge University Press, 2007.
  * [66] Oliver Penrose and Lars Onsager.  Bose-einstein condensation and liquid helium.  Phys. Rev., 104:576–584, Nov 1956.
  * [67] V.N. Popov.  Functional Integrals and Collective Excitations.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 1987.
  * [68] A. Griffin.  Conserving and gapless approximations for an inhomogeneous bose gas at finite temperatures.  Phys. Rev. B, 53:9341–9347, Apr 1996.
  * [69] D. A. W. Hutchinson, E. Zaremba, and A. Griffin.  Finite Temperature Excitations of a Trapped Bose Gas.  Phys. Rev. Lett., 78:1842–1845, 1997.
  * [70] Alexander Altland and Ben D Simons.  Condensed matter field theory.  Cambridge University Press, 2010.
  * [71] L. M. Sieberer, M. Buchhold, and S. Diehl.  Keldysh Field Theory for Driven Open Quantum Systems.  Rept. Prog. Phys., 79(9):096001, 2016.
  * [72] Ben Bar-Or, Jean-Baptiste Fouvry, and Scott Tremaine.  Relaxation in a Fuzzy Dark Matter Halo.  Astrophys. J., 871(1):28, 2019.
  * [73] Pierre-Henri Chavanis.  Landau equation for self-gravitating classical and quantum particles: application to dark matter.  Eur. Phys. J. Plus, 136(6):703, 2021.
  * [74] Pavel Friedrich and Tomislav Prokopec.  Field-theoretic approach to large-scale structure formation.  Physical Review D, 100(10):103527, 2019.
  * [75] Wen-Yuan Ai, Ankit Beniwal, Angelo Maggi, and David J. E. Marsh.  From QFT to Boltzmann: freeze-in in the presence of oscillating condensates.  JHEP, 02:122, 2024.
  * [76] Robert Brandenberger, Vahid Kamali, and Rudnei O. Ramos.  Decay of ALP condensates via gravitation-induced resonance.  JCAP, 11:009, 2023.