\section{Introduction}\label{sec:intro}
Ultra-light particles have been considered a promising candidate for dark matter \cite{Hu:2000ke,Hui:2016ltb,Ferreira:2020fam}. Axions, axion-like particles, and fuzzy dark matter are examples of this class of dark matter models. They arise naturally in UV-complete theories like string theory \cite{Arvanitaki:2009fg}. Presumably any such dark-matter candidates would arise in a fully relativistic setting, though their astrophysical implications would become most relevant deep in the nonrelativistic regime. Hence it is imperative to develop a self-consistent, effective description of the low-energy dynamics of such particles in an expanding spacetime, which incorporates corrections from subleading contributions. Building on Refs.~\cite{Guth:2014hsa,Namjoo:2017nia}, we develop such a formalism here. Our analysis applies to any massive scalar field oscillating in an expanding universe, for which interactions with relativistic fields can be neglected, including cases of physical interest such as ultra-light dark matter.

In this paper we study the dynamics of a massive scalar field $\phi$ minimally coupled to gravity. Our main interest concerns the case in which such a field plays the role of dark matter, although our analysis applies to other situations as well. The particles being ultra-light implies that, to account for the observed density of dark matter, the occupation numbers of $\phi$ particles must be very large. This, in turn, suggests that the system should be describable by a classical scalar field to a good approximation \cite{Guth:2014hsa,Hertzberg:2016tal}. On the other hand, the dark matter is supposed to be cold, and therefore nonrelativistic. Putting these two statements together, we deduce that a {\it classical field theory in the nonrelativistic limit} should provide an appropriate description of such a system.

The dynamics of a classical scalar field of mass $m$ in the nonrelativistic limit should be dominated by oscillations with nearly constant frequency $m$. That is, the energy of $\phi$ should be dominated by its rest mass. Although the scalar field oscillates very fast (compared to the Hubble scale, for example), we develop a formalism with which to study the dynamics of the system over time-scales much longer than $1/m$. If we introduce the ansatz \cite{Guth:2014hsa}
\eq{
	\label{eq:anphi}
	\phi (t, {\bf x} )= \dfrac{1}{\sqrt{2m}}\left( e^{-imt}\psi (t, {\bf x} ) +e^{imt}\psi^* (t, {\bf x} ) \right)\,,	
}   
where $\psi$ is a slowly varying function of time, then (as we will demonstrate below) the energy density of $\phi$ in a Friedmann-Lem\^{a}itre-Robertson-Walker (FLRW) universe is given (to leading order) by
\eq{
	\label{eq:endens}
	\rho= m|\psi|^2\propto a^{-3}\,,
}
which is slowly varying compared to the frequency of oscillations $m$. On the other hand, to leading order the pressure only contains rapidly oscillating terms which average out to zero on time scales larger than $1/m$. 

This motivates an approach that focuses on the slowly varying variable $\psi$.  We therefore take the relation \eqref{eq:anphi} as a field redefinition, with which to define the {\it nonrelativistic field} $\psi$. We will consider this field redefinition in more detail in the next section. Note that the $\psi$ field will not remain entirely as a slowly varying function of time; subdominant but highly oscillating contributions will develop as a result of nonlinear dynamics of the system. Because we are interested in the slowly varying part of the field, we develop techniques with which to remove the rapidly oscillating terms systematically. The result is an effective field theory (EFT) for the slowly varying part. The main aim of this paper is to derive such an EFT. 

Our nonrelativistic EFT will be an expansion in powers of several small parameters. To identify the relevant small parameters, 
 consider a massive, possibly self-interacting, real-valued scalar field $\phi$ coupled minimally to gravity, with action given by
\eq{
	\label{eq:action}
	S=\int\dd[4]{x}\left[\frac{1}{2}\Mpl^2\sqrt{-g}R+\mathcal{L}_\phi\right]\,,
}
where the matter Lagrangian is
\eq{
	\label{eq:lagphi}
	\mathcal{L}_\phi=-\frac{\sqrt{-g}}{2}\left[g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi+m^2 \phi^2+\dfrac{\lambda}{4!}\phi^4\right]\,.
}
The resulting Klein-Gordon equation in an expanding universe reads
\ba 
\ddot \phi +3H\dot \phi -\dfrac{\nabla^2}{a^2}\phi+m^2 \phi +\dfrac{\lambda}{6}\phi^3=0 .
\label{eq:KleinGordon}
\ea 
 In the nonrelativistic limit, any dimensionful parameter with dimension mass should be smaller than the mass of the scalar field. Given the form of Eq.~(\ref{eq:KleinGordon}), this suggests three relevant, small quantities: 
\ba 
\label{eq:eps}
\epsilon_{\scriptscriptstyle H} \sim \dfrac{H}{m}\, , \qquad \epsilon_k \sim \dfrac{k_p^2}{m^2} \, , \qquad \epsilon_\lambda \sim \lambda \dfrac{\phi^2}{m^2},
\ea 
where $k_{p}$ is a typical physical wavenumber. Furthermore, after removing $\phi$ in favor of the new field $\psi$ as in Eq.~\eqref{eq:anphi}, all remaining variables are expected to be slowly varying compared to the time-scale $1/m$. This suggests a fourth small parameter, 
\ba 
\epsilon_t \sim \dfrac{1}{mX}\dfrac{\partial X}{\partial t},
\label{eq:ept}
\ea 
where $X$ is a variable in the theory. We will make our expansion in terms of such small parameters more explicit in the following sections. We will see that this perturbative expansion also allows us to construct the full solution, including oscillations, by solving well-behaved non-oscillatory equations. This can be considered a reliable method for finding solutions of an otherwise generically stiff system of equations with rapidly oscillating terms, with much broader application than what shall be discussed here.

In this paper we mainly focus on the case of a free, massive scalar field, setting $\lambda = 0$. As we will see, even in the absence of tree-level self-interactions for $\phi$, the backreaction from rapidly oscillating modes on the evolution of the nonrelativistic system in an expanding FLRW spacetime induces subdominant self-interaction terms in the evolution of $\psi$. The induced self-coupling could be of interest in the light of astrophysical observations \cite{Bullock:2017xww} which suggest that cold dark matter particles might depart from a purely collisionless regime \cite{Fan:2016rda,Tulin:2017ara,Ferreira:2020fam}, though a direct comparison of effects arising from the induced self-interaction with recent astrophysical observations remains beyond the scope of the present paper.
In Appendix~\ref{app:interaction} we consider the case of tree-level self-interactions as well.

Besides deriving an EFT for the nonrelativistic field, we also develop an equivalent description of the system in terms of an imperfect fluid. We show that in the latter description, one needs to identify a nonzero pressure, a nontrivial sound speed, and an effective bulk viscosity to fully express the EFT, up to the working order, in the fluid language. 

While the effect of oscillations of ultra-light particles as dark matter are neglected in a majority of related analyses, they are studied in a rather restricted literature.  The effect of oscillations on the dynamics of binaries is studied in Ref.~\cite{Blas:2019hxz}, while Ref.~\cite{Boskovic:2018rub} explores their impact on the orbital motion of stars in galactic halos. The resulting parametric resonance are also studied e.g. in Ref.~\cite{Tkachev:1986tr} and recently in Ref.~\cite{Hertzberg:2020dbk}. Ref.~\cite{Cookmeyer:2019rna} also compares the predictions for effects on the cosmic microwave background (CMB) from the naive theory, in which oscillations are neglected, with predictions from the exact theory.


 
Let us make some remarks on terminology used in this paper. We work up to linear order in spatially varying fluctuations around a homogeneous and isotropic background. We are interested in nonlinearities in the time evolution of various quantities which occur even at this order; these nonlinearities arise because of the coupling between different time-dependent variables, which yield a {\it mode coupling} once we decompose the relevant quantities into series of modes with distinct characteristic time dependencies. To avoid confusion and simplify the terminology, we refer to quantities that vary in space (around a homogeneous, isotropic background) as ``fluctuations,'' and restrict ``(non)linearities'' to features of the time evolution of various quantities. Furthermore, as noted above, we are interested in an EFT for the slowly varying parts of the fields, which we denote the ``slow modes,'' while we systematically remove the rapidly oscillating contributions, which correspond to nonzero modes in the mode decomposition (in time) of various quantities. Finally, we adopt a perturbative scheme, accounting for various contributions to the effective equations of motion in terms of different powers of small parameters. We therefore distinguish between ``perturbations'' (as powers of these small expansion parameters) and ``fluctuations'' (as spatially varying quantities around the FLRW background).

The rest of the paper is organized as follows. In Sec.~\ref{sec:redef} we introduce a suitable field redefinition for the scalar field in a general spacetime geometry as well as in a perturbed FLRW universe. In Sec.~\ref{sec:EFT}, we derive our EFT for the background evolution as well as fluctuations. In Sec.~\ref{sec:fluid} we re-express our EFT results in terms of fluid dynamics. In Sec.~\ref{sec:num} we confirm the validity of our results by numerical analysis, and in Sec.~\ref{sec:con} we summarize our main results and discuss possible future directions. Furthermore, in a number of appendices, we explore different corners of the subject and extend the analyses presented in the main body of the paper. In Appendix~\ref{app:Hamiltonian} we derive the Hamiltonian of the system after the field redefinition discussed in Sec.~\ref{sec:redef}. In Appendix~\ref{app:linhigher} we present the equations for higher-order nonzero modes, ignored in the results of Sec.~\ref{sec:EFT}. In Appendix~\ref{app:vis} we discuss how an imperfect fluid can be described in the general relativistic framework, the results of which are used for the fluid interpretation of our EFT in Sec.~\ref{sec:fluid}. In Appendices~\ref{app:interaction} and \ref{app:multi-component}, we extend our results to the self-interacting case and to a multicomponent universe, respectively. In Appendix~\ref{app:nonlocal} we employ a nonlocal field redefinition in Minkowski spacetime to confirm the higher-order, momentum-dependent terms we obtained in the effective sound speed (which has a different numerical factor compared to previous analyses). The resulting sound speed is valid for arbitrary momentum and reduces to the well-known results in the limits of small and large momenta. Finally, in Appendix~\ref{app:gauge} we derive the gauge transformation that can be used to express our EFT (written in Newtonian gauge) in other gauges.

\section{A suitable field redefinition}\label{sec:redef}
Taking the nonrelativistic limit of a scalar field theory is usually done by starting with an appropriate field redefinition, which is the subject of this section. Similar to the approach in Ref.~\cite{Namjoo:2017nia}, we aim to take the nonrelativistic limit and obtain higher-order corrections to a low-energy EFT {\it after} performing the field redefinition. 

We are interested in studying the system at low energies, for which the  mass term is the dominant contribution in the system's evolution. That is, we expect the main time dependence of the system to be rapid oscillations with frequency approximately equal to the mass $m$. This motivates a natural choice for a complex-valued nonrelativistic field, $\psi$, related to $\phi$ as in Eq.~(\ref{eq:anphi}). The prefactor $1/\sqrt{2m}$ in Eq.~(\ref{eq:anphi}) has been chosen so that $|\psi|^2$ corresponds to the number density of particles, as will become evident below. We will also see that this field transformation, to leading order, yields a Schr\"odinger-like equation for $\psi$ in Minkowski spacetime. 


As discussed above, $\psi$ is expected to be a slowly varying function of time in the nonrelativistic limit. Note that since, by Einstein's equations, the metric components are sourced by the energy-momentum tensor, which is quadratic in the original field $\phi$ --- and hence in the new field $\psi$ --- the metric components are dominated by slowly varying functions of time, although they also contain subdominant oscillatory contributions. For this reason, we do not make a similar field redefinition for the spacetime metric. 

As mentioned in Ref.~\cite{Namjoo:2017nia}, the field redefinition of Eq.~\eqref{eq:anphi} may not be considered complete if the real and imaginary parts of $\psi$ are treated as two independent variables, since that would consist of replacing one degree of freedom with two degrees of freedom. In Ref.~\cite{Namjoo:2017nia}, this is addressed by adding to Eq.~\eqref{eq:anphi} a comparable field redefinition for $\pi$, the conjugate momentum of $\phi$, such that the transformation from $(\phi,\pi)$ to $(\psi,\psi^*)$ becomes a canonical transformation. In addition, in Ref.~\cite{Namjoo:2017nia} a nontrivial nonlocal operator is used within the transformation, which simplifies the equation of motion for $\psi$ in Minkowski spacetime. In principle, a similar procedure can be employed when studying the behavior of $\psi$ in a curved spacetime. However, as we discuss in Appendix~\ref{app:nonlocal}, a nonlocal field redefinition does not in general result in any simplification in curved spacetime, even for the simple case of an exact FLRW universe. Notice that the nonlocal field redefinition is not essential for obtaining the nonrelativistic limit and we will not exploit it in a curved geometry.\footnote{See, however, Ref.~\cite{Friedrich:2019zic} where a generalization of the nonlocal operator introduced in Ref.~\cite{Namjoo:2017nia} to the case of expanding background has been employed.} Furthermore, making a transformation for $\pi$ similar to what was used in Ref.~\cite{Namjoo:2017nia} results in equations of motion different from the Schr\"odinger equation written in an FLRW background \cite{Guth:2014hsa}. Finally, we note that in terms of the original field $\phi$, the conjugate momentum differs from $\dot \phi$ in a curved spacetime, leaving more freedom for an appropriate field redefinition. 

Given these considerations, in this paper we follow a different approach compared to Ref.~\cite{Namjoo:2017nia}. From Eq.~\eqref{eq:anphi}, it is evident that there is a redundancy in the definition $\psi$ of the form
\eq{
	\label{eq:gaugetr}
	\psi (t, {\bf x}) \to\psi (t, {\bf x}) +ie^{imt}\eta (t, {\bf x}) \,,
}   
with $\eta$ any real function of space and time. It is easy to show that this transformation leaves $\phi$ --- and, as a result, the Lagrangian --- invariant. To fix the gauge and remove the redundancy, the following transformation from an arbitrary gauge seem to be appropriate:
\eq{
	\label{eq:gaugefix}
	\eta=\frac{1}{2m}\left(e^{-imt}\dot{\psi}+e^{imt}\dot{\psi}^*\right)\,.
} 
Note that this choice of the gauge guarantees that the following constraint is satisfied in the new gauge:
\eq{
\label{eq:gaugecon}
e^{-imt}\dot{\psi}+e^{imt}\dot{\psi}^*=0\,,
}
which allows us to have a first order, Schr\"odinger-like equation for $\psi$. Therefore, we may obtain the equations of motion for the new fields after gauge fixing by varying the following Lagrangian with respect to $\psi$ and $\psi^*$:
\eq{
\label{eq:lagmultiplier}
\mathcal{L}_\phi(\psi,\psi^*)+\xi (e^{-imt}\dot{\psi}+e^{imt}\dot{\psi}^*)\,,
} 
where $\mathcal{L}_\phi(\psi,\psi^*)$ is the Lagrangian given in Eq.~\eqref{eq:lagphi} with $\phi$ replaced by the field redefinition of Eq.~\eqref{eq:anphi}, and $\xi$ is a Lagrange multiplier, which ensures that the constraint of Eq.~(\ref{eq:gaugecon}) is satisfied as a result of gauge fixing. Note that in this gauge, we have
\eq{
	\label{eq:phidot}
	\dot{\phi}=-i\sqrt{\frac{m}{2}}\left(e^{-imt}\psi-e^{imt}\psi^*\right)\,,
}
or, by using Eqs.~\eqref{eq:anphi} and \eqref{eq:phidot}, 
\eq{
\label{eq:psiphi}
	\psi=\sqrt{\frac{m}{2}}e^{imt}\left(\phi+i\frac{\dot{\phi}}{m}\right)\,,
}
which shows that our field redefinition is invertible. Notice that the equation of motion for $\xi$ is a differential equation, rather than an algebraic constraint, and it is not possible to substitute the solution of $\xi$ back into the Lagrangian. In other words, the gauge condition cannot be imposed at the level of the Lagrangian. This makes the Lagrangian rather complex.
It is, however, possible to derive a different, and simpler, Lagrangian which yields the same equations of motion via the following procedure. First, in the original theory in terms of $\phi$, we introduce a new field by the relation $\chi = \dot \phi$ (at the level of equations of motion); next we replace $\dot \phi$ with $\chi$ in the Lagrangian and add a Lagrange multiplier to guarantee the equivalence of the two theories: 
\eq{
\label{eq:barL}
\bar{\mathcal{L}}_\phi(\phi,\chi)=\mathcal{L}_\phi(\phi,\chi,\partial_i\phi)+\pi(\dot{\phi}-\chi)\,.
}
Note that at this stage, the field $\pi$ is only a Lagrange multiplier and \emph{not} the momentum conjugate of $\phi$. The equations of motion for this Lagrangian are as follows:
\eq{
	\label{eq:chiphi}
	\dot{\phi}-\chi=0\,,\qquad\pdv{\mathcal{L}_\phi}{\chi}-\pi=0\,,\qquad\pdv{\mathcal{L}_\phi}{\phi}-\partial_i\pdv{\mathcal{L}_\phi}{\partial_i\phi}-\dot{\pi}=0\,.
}
The first equation is the constraint. The second equation, together with $\chi=\dot{\phi}$, shows that the Lagrange multiplier is nothing but the momentum conjugate of the scalar field in the original Lagrangian (hence our motivation for the notation, $\pi$). Finally, using the first two equations, the last one gives the correct equation of motion for $\phi$. Therefore, the two theories $\mathcal{L}_\phi$ and $\bar{\mathcal{L}}_\phi$ are equivalent. 

In the theory governed by $\bar{\mathcal{L}}_\phi$, the parameter $\pi$ is nondynamical, so one can replace it by its solution to simplify the Lagrangian. From Eq.~\eqref{eq:chiphi} we have:
\eq{
\pi=	-\sqrt{-g}(g^{00}\chi+g^{0i}\partial_i\phi)\,.
}
Plugging this solution back into the Lagrangian results in
\eq{
	\label{eq:phichilag}
\mathcal{L}(\phi,\chi)=-\frac{1}{2}\sqrt{-g}\left[-g^{00}\chi^2+2(g^{00}\chi+g^{0i}\partial_i\phi)\dot{\phi}+g^{ij}\partial_i\phi\partial_j\phi+m^2\phi^2\right]\, ,
} 
where, from now on, we consider a non-self-interacting theory, i.e. we set $\lambda=0$ in Eq.~\eqref{eq:lagphi}. (We will consider tree-level self-interactions in Appendix~\ref{app:interaction}.)
Again, one can show that the Lagrangian in Eq.~(\ref{eq:phichilag}) yields consistent equations of motion. The next step is to introduce the nonrelativistic fields $\psi$ and $\psi^*$, which we define by the following expressions:
\eqa{
\label{eq:cantr1}
\phi=\frac{1}{\sqrt{2m}}\left(e^{-imt}\psi+e^{imt}\psi^*\right) \, , \qquad 
\chi=-i\sqrt{\frac{m}{2}}\left(e^{-imt}\psi-e^{imt}\psi^*\right)\,,
}  
which are indeed suggested by Eqs.~\eqref{eq:anphi} and \eqref{eq:phidot}. With these definitions, the Lagrangian in Eq.~\eqref{eq:phichilag} may be written
\eq{
\spl{
\mathcal{L}(\psi,\psi^*)=-\frac{\sqrt{-g}}{2}\Bigg[ &ig^{0\mu}(\psi^*\partial_\mu\psi-\psi\partial_\mu\psi^*)+\frac{1}{m}g^{\mu'\nu'}\partial_{\mu'}\psi\partial_{\nu'}\psi^*+m(g^{00}+1)\psi^*\psi\\
&+\left\{\frac{e^{-2imt}}{2m}\bigg(g^{\mu'\nu'}\partial_{\mu'}\psi\partial_{\nu'}\psi+\left(m^2(g^{00}+1)+\frac{im}{\sqrt{-g}}\partial_\mu(\sqrt{-g}g^{0\mu})\right)\psi^2\bigg)+\text{c.c.}\right\}\Bigg]\, .
}
\label{eq:Lpsipsistar}
}
Here ``$\text{c.c.}$'' denotes the complex conjugate of the immediately preceding expression, and the primes on indices in $g^{\mu'\nu'}$ indicate that $\mu$ and $\nu$ cannot both be time components, i.e., $\mu\times \nu\neq0$. (In Appendix~\ref{app:Hamiltonian}, we consider the Hamiltonian of this theory in terms of $\psi$ and $\psi^*$.) From Eq.~(\ref{eq:Lpsipsistar}), the equation of motion for $\psi$ can be written as
\eq{
ig^{00}\dot{\psi}+\mathcal{D}\psi+e^{2imt}\mathcal{D}^*\psi^*=0\,,
\label{eq:psidotD}
}
where we have defined the operator
\eq{
\label{eq:D}
\mathcal{D}=\frac{m}{2}\left(g^{00}+1\right)+ig^{0i}\partial_i+\frac{i}{2\sqrt{-g}}\partial_\mu(\sqrt{-g}g^{0\mu})-\frac{1}{2m\sqrt{-g}}\partial_\mu(\sqrt{-g}g^{\mu i})\partial_i-\frac{1}{2m}g^{ij}\partial_i\partial_j\,.
}
As promised, the equation of motion for $\psi$ is first order in time derivatives and is identical to the equation of motion from the constrained Lagrangian of  Eq.~\eqref{eq:lagmultiplier}. In Minkowski spacetime the Lagrangian reduces to 
\ba 
\mathcal{L}(\psi,\psi^*)= \dfrac{i}{2}(\psi\dot \psi^*-\psi^*\dot \psi)-\frac{1}{2m}\partial_{i}\psi\partial_{j}\psi^*
-\left\{\frac{e^{-2imt}}{4m}\partial_{i}\psi\partial_{j}\psi+\text{c.c.}\right\}
 \, \, \text{for Minkowski}.\quad 
\ea 
Neglecting the last term, which is rapidly oscillatory, the resulting theory has a global $U(1)$ symmetry associated to which there is a conserved charge, $\int d^3 x|\psi|^2$, which is the number of particles, a quantity that is conserved in a nonrelativistic theory. The equation of motion reads
\eq{
i\dot{\psi}+\dfrac{\nabla^2}{2m}(\psi+e^{2imt}\psi^*)=0  \qquad \text{for Minkowski},
}
which, again after neglecting the oscillating term, is just the Schr\"odinger equation with vanishing potential. Note, however, that so far our equations are exact and only a field redefinition has been performed. We will see how one can remove the oscillatory terms in a systematic way, rather than just neglecting them, in Sec.~\ref{sec:EFT}.

In a curved background, the complete set of equations also includes the Einstein equations, $\Mpl^2G_{\mu\nu}=T_{\mu\nu}$, where the energy-momentum tensor for the scalar field is
\eq{
	T_{\mu\nu}=\partial_\mu\phi\partial_\nu\phi-\frac{1}{2}g_{\mu\nu}\left(g^{\alpha\beta}\partial_\alpha\phi\partial_\beta\phi+m^2\phi^2\right)\,,
	\label{eq:Tmn}
}
in which we must insert definitions of Eqs.~\eqref{eq:anphi} and \eqref{eq:phidot}. To avoid clutter we do not write out $T_{\mu\nu} (\psi, \psi^*)$ explicitly here.

\subsection{FLRW background}
What has been derived so far is valid for a general curved geometry. For the remainder of this paper, we consider the case of a spatially flat FLRW spacetime with small fluctuations around it. Neglecting the tensor and vector degrees of freedom, in Newtonian gauge, the line element takes the form
\eq{
	\label{eq:ds2}
	\dd{s}^2=-(1+2\Phi)\dd{t}^2+a^2(t)(1-2\Phi)\delta_{ij}\dd{x}^i\dd{x}^j\,,
}
where $\Phi (t, {\bf x})$ is the Newtonian potential. Since the anisotropic stress vanishes (to leading order in fluctuations) for an FLRW spacetime filled with a scalar field, the fluctuation
in the $00$-component of the metric is identical to the
fluctuation in the spatial components. We split the field into a background component and small fluctuations, $\psi (t, {\bf x}) =\bar{\psi} (t) +\delta\psi (t, {\bf x} )$, where the splitting becomes unambiguous by requiring that the spatial average of $\delta\psi (t, {\bf x})$ vanishes. The operator $\mathcal{D}$ defined in Eq.~\eqref{eq:D} up to linear order in $\Phi$ is
\eq{
\label{eq:DFRW}
\mathcal{D}=-\frac{3iH}{2}-\frac{\nabla^2}{2ma^2}+m\Phi+3iH\Phi+2i\dot{\Phi}\, .
} 
Hence, from Eqs.~(\ref{eq:psidotD}) and (\ref{eq:Tmn}), the equations governing the background evolution are
\eqa{
\label{eq:bckg1}
i\dot{\bar{\psi}}+i\frac{3}{2}H\left(\bar{\psi}-e^{2imt}\bar{\psi}^*\right)&=0\\
\label{eq:bckg2}
3\Mpl^2H^2-m\bar{\psi}^*\bar{\psi}&=0\, .
}  
To linear order in spatially varying quantities, we find from the Klein-Gordon equation and the $0i$-component of Einstein's equations (after some algebra) the coupled equations
\eqa{
\label{eq:linpert1}
i\dot{\delta\psi}+[\mathfrak{D}\delta\psi-(m-2iH)\bar{\psi}\Phi]
	+e^{2imt}[\mathfrak{D}^*\delta\psi^*-(m+2iH)\bar{\psi}^*\Phi]&=0 ,\\
\label{eq:linpert2}
\dot{\Phi}+H\Phi+\frac{i}{4\Mpl^2}\left(\bar{\psi}\delta\psi^*-\bar{\psi}^*\delta\psi+e^{-2imt}\bar{\psi}\delta\psi-e^{2imt}\bar{\psi}^*\delta\psi^*\right)&=0\, ,
}
in which we have defined a new operator
\eq{
\label{eq:Df}
\mathfrak{D}=\frac{3iH}{2}+\frac{\nabla^2}{2ma^2}-\frac{1}{2\Mpl^2}\left(e^{-imt}\bar{\psi}-e^{imt}\bar{\psi}^*\right)^2\,.
}
Note the difference between the new operator $\mathfrak{D}$ and the operator $\mathcal{D}$ in Eq.~\eqref{eq:D} for a general metric and its form for FLRW universe in Eq.~\eqref{eq:DFRW}. In deriving Eq.~\eqref{eq:linpert1} for $\delta\psi$ we have used Eq.~\eqref{eq:linpert2} to eliminate the $\dot{\Phi}$ term appearing in Eq.~\eqref{eq:DFRW}. Further, note that $\Phi$ is a non-dynamical degree of freedom and can be eliminated by solving the $00$-component of Einstein's field equations,
\eq{
\label{eq:pois}
\frac{\nabla^2\Phi}{a^2}=\frac{m}{2\Mpl^2}(\bar{\psi}\delta\psi^*+\bar{\psi}^*\delta\psi)+3H\dot{\Phi}+
\frac{m}{4\Mpl^2}(e^{-imt}\bar{\psi}+e^{imt}\bar{\psi}^*)^2\Phi\, ,
}
which is the relativistic form of the Poisson equation. Here we have not explicitly eliminated the $\dot{\Phi}$ term from Eq.~\eqref{eq:linpert2} to avoid unnecessary complications, but we will do so in the calculations of the next section.  
Although it seems redundant to consider Eq.~\eqref{eq:linpert2} for $\Phi$, it would be much easier to obtain an effective description using Eq.~\eqref{eq:linpert2} rather than \eqref{eq:pois}. But we will keep track of Eq.~\eqref{eq:pois} and its effective form in the following sections since it helps us to interpret the results in the form of an effective imperfect fluid.   

Note again that the equations are so far exact, up to ${\cal O} (\Phi,\delta\psi)$. However, in the nonrelativistic  (i.e. large mass) limit, we may neglect terms that are rapidly oscillating, since they average out to zero, as well as terms that are suppressed by factors of $(H/m)$. The above system of equations then reduces to the well-known Schr\"odinger-Friedmann and Schr\"odinger-Poisson equations for background and fluctuations, respectively
 \cite{Guth:2014hsa}:
\eqa{
\label{eq:leadbck}
	i\dot{\bar{\psi}}+i\frac{3}{2}H\ \bar{\psi} \simeq 0 \, , \qquad 
	3\Mpl^2H^2 =m\bar{\psi}^*\bar{\psi} \, ,
	\\
	i\dot{\delta\psi}+\dfrac32 iH\delta \psi+\dfrac{\nabla^2}{2ma^2}\delta \psi -m\bar{\psi}\Phi \simeq 0\, ,
\\
		\left(\frac{\nabla^2}{a^2} +\dfrac{3}{2}H^2\right)\Phi \simeq \frac{m}{2\Mpl^2}(\bar{\psi}\delta\psi^*+\bar{\psi}^*\delta\psi)\, ,
	}
where the second term on the left-hand side of the last equation is usually neglected, as these equations are usually considered for short (sub-horizon) scales. Note that the above equations are equivalent to a system of equations governing the evolution of a matter-dominated universe, showing that a nonrelativistic scalar field can play the role of cold dark matter (to leading order; see Sec.~\ref{sec:EFT} and Sec.~\ref{sec:fluid} for subleading deviations from this statement). 

We expect the above simplified equations to be valid to a good approximation. However we are interested in how subleading contributions change the dynamics of the approximated system. In particular, we want to take into account the backreaction of rapidly oscillating modes on the evolution of the slowly varying quantities. We therefore keep the subleading terms and later on will remove the fast oscillating modes in a systematic way to obtain an effective theory for slow modes. Note that, as discussed in the introduction, we are considering a system in which terms that vary on different characteristic time-scales couple to each other, which gives rise to the backreaction.
 
\section{The effective field theory in the nonrelativistic limit}
\label{sec:EFT}
In this section we outline the derivation and present the results of our effective field theory. We first introduce an appropriate mode decomposition (in time) to disentangle the part of the fields in which we are most interested (which varies slowly in time) from the parts that need to be systematically removed (which oscillate rapidly). We then discuss a procedure for removing the latter part, resulting in an effective theory for the former. We mainly follow the procedure outlined in Ref.~\cite{Namjoo:2017nia} (but see also Refs.~\cite{Mukaida:2016hwd,Braaten:2018lmj} for different, but equivalent, approaches). Here we consider the case in which the single field $\phi$ dominates the FLRW universe, and neglect any tree-level self-interactions. The effect of self-interaction shall be detailed in Appendix \ref{app:interaction} while the EFT for the case of a multicomponent universe is studied in Appendix \ref{app:multi-component}.  

\subsection{Smearing and mode expansion}\label{sec:mexp}

By applying the field redefinition introduced in the last section, in which we isolated the main time dependence of the field as an oscillatory factor (to leading order), we expect that apart from explicit factors like $\exp(2imt)$ appearing in the equations of motion, all remaining functions to be slowly varying in time. Let us denote the variables appearing in the equations of the last section collectively by $X(t)$, where we suppress any possible spatial dependence. By the above reasoning we expect $\dot{X}(t)\ll m X(t)$ with $m$ the mass of the scalar field. In other words, the spectrum of $X(t)$ in the frequency domain, $\hat X(\omega)$, would be localized around $\omega=0$. However, this cannot be the whole story because the explicit rapidly oscillating factors in the equation of motion, proportional to $\exp(2imt)$, will induce high-frequency components in the spectrum of $X(t)$. This happens through nonlinear multiplicative terms in Eqs.~\eqref{eq:bckg1}-\eqref{eq:pois}, which cause different frequencies to couple. Since the spectrum of the oscillatory factors like $e^{i \nu m t}$ (with $\nu$ an integer) would be a delta function at frequencies $\omega=\nu m$, we expect that at higher orders a localization happens around the frequencies $\omega=\nu m$. These subleading contributions would then backreact, again due to nonlinearities in the equations, on the slowly varying modes, whose spectra are localized around $\omega=0$. Our aim in this section is to study this backreaction to see how it affects the dynamics of the slow modes of different variables. 


Since we are interested in the slow mode of functions of time (i.e. the portion of functions that are slow compared to oscillations with frequency $m$ or higher), as discussed in Sec.~\ref{sec:intro}, we can define a {\it smearing} operator acting on each function, which can be thought of as a time-average of the function, subject to a window function $W(t)$:
\eq{
\label{eq:avg}
	\ev{X(t)}=\int_{-\infty}^\infty \dd{t'}X(t')W(t-t')\, .
}  
This is the temporal counterpart of the spatial average considered for example in Ref.~\cite{Baumann:2010tm}. To adopt a suitable window function, note that in Fourier space, the right-hand side of Eq.~\eqref{eq:avg} takes the form ${\hat X}(\omega)\hat{W}(\omega)$.  As mentioned above, since the slow mode is localized around $\omega=0$ in the frequency domain, a natural choice for the window function would be a top-hat in Fourier space:
\eq{
{\hat W}(\omega)=
\begin{cases}
1\quad|\omega|<m/2\\
0\quad\text{otherwise}\,.
\end{cases}	
}   
Note that we have chosen the half-width of the square function to be $m/2$, which means that the window function removes the portion of the field that is spread away more than $m/2$ from each side of $\omega=0$. Choosing such a width is, first, consistent with requiring that the resulting smeared function has to be slowly varying compared to frequency $m$ and, second, leads to a convenient (and exact) mode expansion of $X(t)$, as we will explain below.
The window function in the time domain is then
\eq{
\label{eq:w}
W(t)=\int_{-\infty}^\infty \frac{\dd{\omega}}{2\pi}{\hat W}(\omega)e^{i\omega t}=\frac{\sin(mt/2)}{\pi t}\,.
}
Figure~\ref{fig:window} visualizes the window function in time and frequency domains. Note that the peculiar form of the window function in the time domain, with many zeros, is essential for capturing the low-frequency part (i.e. the slow mode) of the function. 
\begin{figure}
\begin{tikzpicture}[yscale=1.5]
\pgfmathdeclarefunction{sinc}{1}{%
        \pgfmathparse{abs(#1)<0.01 ? int(1) : int(0)}%
        \ifnum\pgfmathresult>0 \pgfmathparse{1}\else\pgfmathparse{sin(#1 r)/#1}\fi%
        }
\draw [thick , ->] (-4,0) -- (4,0);
\draw [thick , ->] (0,-0.5) -- (0,1.8);
\draw[domain=-3.8:3.8, samples=250,thick,blue] plot (\x,{sinc(3.5*\x)});

\node [above] at (0,1.8) {$W(t)$};
\node [right] at (4,0) {$t$};
\node at (0.3,1.15) {\small $\frac{m}{2\pi}$};
\node at (0.8,-0.18) {\small $\frac{2\pi}{m}$};

\end{tikzpicture}
\hspace{0.5 cm}
\begin{tikzpicture}[yscale=1.5]
\pgfmathdeclarefunction{rect}{1}{%
        \pgfmathparse{abs(#1)<1.5 ? 0.8 : int(0)}%
        }
\draw [thick ,->] (-3,0) -- (3,0);
\draw [thick ,->] (0,-0.5) -- (0,1.8);
\draw[domain=-3:3, samples=300,thick, blue] plot (\x,{rect(\x)});

\node [above] at (0,1.8) {$\hat{W}(\omega)$};
\node [right] at (3,0) {$\omega$};
\node at (1.5,-0.18) {\small $\frac{m}{2}$};
\node at (0.2,0.95) {\small $1$};
\end{tikzpicture}
\caption{The window function $W(t)$ (left) and its Fourier transform (right).}
\label{fig:window}
\end{figure}

As mentioned earlier, in our nonlinear system, even if we start with purely slow functions of time, subdominant nonzero modes (which mainly vary with frequencies $\nu m$ with nonzero integer $\nu$) would develop. To capture these nonzero modes we can apply the same smearing operator and define the mode $\nu$ of the function by 
\eq{
	\label{eq:fnu_def}
	X_\nu(t)=\ev{X(t)e^{-i\nu mt}}=\int_{-m/2}^{m/2}\frac{\dd{\omega}}{2\pi}\hat{X}(\omega+\nu m)e^{i\omega t}\, ,
}
where in the last equality we explicitly used the adopted window function in the frequency domain. 
This result allows us to write any function of time by a mode expansion as follows: 
\eq{
	\label{eq:modexp}
	X(t)=\sum_{\nu=-\infty}^{\infty}X_\nu(t)e^{i\nu mt}\,.
}  
This is the key relation that will be used in the derivation of our EFT and from now on, our discussion will not explicitly rely on the smearing operator. Note that the specific choice of the window function makes the smearing operator a projection operator, i.e. $\ev{\ev{X(t)}}=\ev{X(t)}$. Note also that Eq.~\eqref{eq:modexp} is indeed exact, as a result of the appropriate choice of window function and its width; as can be checked by replacing the definition of $X_\nu$ from \eqref{eq:fnu_def} into \eqref{eq:modexp}. For later usage, note that the following relations hold for arbitrary functions of time $X(t)$ and $Y(t)$:
\eq{
	(XY)_\nu=\sum_{\alpha}X_\alpha Y_{\nu-\alpha}\,, \qquad 
	(X e^{i\mu m t})_\nu=X_{\nu-\mu}\, , \qquad (X^*)_\nu = X_{-\nu}^{*},
}
where the last equality should be understood as a relation between the modes of the complex conjugate of the field and the complex conjugate of the modes of the field. 
It is also easy to show that for a real function $R(t)$ we have $R_{-\nu}=R^*_\nu$. 

In the nonrelativistic limit, we will be interested in functions that are mainly localized around integer multiples of the characteristic frequency of the system $m$. In this regime, all coefficients $X_\nu(t)$ are slowly varying in time. Furthermore, in the same limit, the fields are mainly concentrated around $\omega=0$, i.e. the slow mode dominates over all other modes. This implies that the functions of time that we will be dealing with are slowly varying, so that a dimensionless parameter like  $\epsilon_t \sim |\dot{X}/mX|$ would always be small. Our EFT will then be an expansion in powers of such a  parameter, although other small parameters will come in as well, as discussed in the introduction and will be made explicit in the next section.

Now that we have a procedure to expand any function of time in terms of its different modes, we can derive differential equations for each mode. We are specially interested in the equation for the slow mode which, for later convenience, we denote by $X_s(t)\equiv X_{\nu=0}(t)$. The
evolution of each mode involves all other modes due to nonlinearities in the equations. In the nonrelativistic limit, in which our time-dependent functions should be dominated by the slow mode, we require $|X_s| \gg |X_\nu|$ for $\nu \neq0$. In this regime, we can solve equations for nonzero modes perturbatively in the small quantities discussed in the Introduction. To achieve this, we further decompose each nonzero mode by a perturbative expansion
\eq{
\label{eq:pertexp}
X_\nu=\sum_{n=1}^{\infty}X_\nu^{(n)}= X_\nu^{(1)}+X_\nu^{(2)}+\dots\,, \qquad {\text{for}} \, \nu \neq 0,
} 
where for each $n$ we have $X_\nu^{(n)}=\order{\epsilon^n}$ and $\epsilon$ is a small parameter, to be determined by the equations governing the dynamics of the system, as in Eqs.~\eqref{eq:eps} and \eqref{eq:ept}. We always work up to the same order in all small parameters so that we can schematically show the orders in a perturbative expansion by just one parameter $\epsilon$ as above. After obtaining the solutions of the nonzero modes (in terms of the slow modes) order by order in perturbation theory, we plug the solutions back into the equations for the slow mode. This procedure leads to an effective equation for the slow mode. The tools that we have developed here (mostly implicitly) will allow us to obtain EFTs for the background evolution as well as for the small fluctuations, as we shall show explicitly below.    

\subsection{Effective field theory for the background}\label{sec:eftbkg}
In this section we apply the procedure we outlined in the last section to the equations of motion for the spatially homogeneous scalar field in an FLRW background, Eqs.~\eqref{eq:bckg1} and \eqref{eq:bckg2}. For convenience, we introduce the rescaled, dimensionless quantities 
\eq{
\label{eq:rescale}
\tilde{\psi} (t) \equiv\frac{\bar{\psi} (t) }{\sqrt{m}\Mpl}\,,\qquad\tilde{H} (t) \equiv\frac{H (t)}{m}\,,
}
which we use in the course of calculations. We also rescale the time $\tilde{t}=mt$. Application of the mode expansion introduced in the previous section to Eqs.~\eqref{eq:bckg1} and \eqref{eq:bckg2}, with rescaled variables, results in the following equations for each mode $\nu$:
\eqa{
	\label{eq:psinu}
	\tilde{\psi}'_\nu+
	i\nu \tilde{\psi}_\nu+\frac{3}{2}\tilde{H}_\alpha(\tilde{\psi}_{\nu-\alpha}-\tilde{\psi}^*_{2+\alpha-\nu})&=0\\
	\label{eq:Hnu}
	3\tilde{H}_\alpha \tilde{H}_{\nu-\alpha}-\tilde{\psi}_\alpha \tilde{\psi}^*_{\alpha-\nu}&=0\,,
}
where the summation over repeated indices is understood and  the prime denotes the derivative with respect to $\tilde{t}$. We are, in particular, interested in the equations for the slow mode quantities $\bar{\psi}_s=\ev{\bar{\psi}}$ and $H_s=\ev{H}$, with the corresponding equations after rescaling the variables:
\eq{
\label{eq:psis0}
\tilde{\psi}'_s+\frac{3}{2}\bigg(\tilde{H}_s\tilde{\psi}_s+\sum_{\alpha\neq0}\tilde{H}_\alpha\tilde{\psi}_\alpha\bigg)-\frac{3}{2}\bigg(\tilde{H}_s\tilde{\psi}_2^*+\tilde{H}_{-2}\tilde{\psi}_s^*+\sum_{\alpha\neq0,-2}\tilde{H}_\alpha\tilde{\psi}^*_{\alpha+2}\bigg)=0\,,
}
and
\eq{
\label{eq:Hs0}
3\bigg(\tilde{H}_s^2+\sum_{\alpha\neq0}|\tilde{H}_\alpha|^2\bigg)-\bigg(|\tilde{\psi}_s|^2+\sum_{\alpha\neq0}|\tilde{\psi}_\alpha|^2\bigg)=0\,,
}
which are obtained simply by setting $\nu=0$ in Eqs.~\eqref{eq:psinu} and \eqref{eq:Hnu}. We see clearly from the last set of equations that the slow  mode of each function is coupled to all nonzero modes due to nonlinearities of the system. We can solve the equations for the nonzero modes perturbatively in small parameters. As discussed in the introduction, we identify a new small parameter, in addition to $\epsilon_t$, from the Hubble scale
\eq{
\epsilon_{\scriptscriptstyle H}\sim\frac{H}{m} =\tilde H\,.
}
It is evident that such a quantity has to be small for a nonrelativistic field rapidly oscillating around the minimum of its potential. 
Thus, at the background level, we have two small parameters in the problem, which we collectively denote by $\epsilon=\{\epsilon_t,\epsilon_{\scriptscriptstyle H}\}$. In this section we will keep terms up to third order in $\epsilon$ but our formalism is general and can be continued to higher orders. Note that the rescaled variables introduced in Eq.~\eqref{eq:rescale} are first order in $\epsilon$, i.e. $\tilde{\psi}_s=\order{\epsilon}$ and $\tilde{H}_s=\order{\epsilon}$, which makes the power counting in the derivation of EFT fairly easy. We also have $\tilde{\psi}_s'\sim\epsilon\tilde{\psi}_s=\order{\epsilon^2}$. It is important to notice that working up to order $n$ in solving for nonzero modes would result in an EFT for the slow modes up to $\order{\epsilon^{n-1}}$ corrections once the corrections are compared to leading order terms. This is because even the zeroth-order equations for the slow modes are suppressed by a factor of $\order{\epsilon}$ (due to, for example, the time derivative operator or the Hubble parameter). The same statement will be valid for fluctuations as well.

To proceed systematically, we apply the perturbative expansion of \eqref{eq:pertexp} to the nonzero modes. Note that we do not apply the perturbative expansion to the slow modes since we are interested in the effective equations for those quantities, rather than their corresponding solutions. From \eqref{eq:psinu} and for $\nu\neq0$ we have 
\eq{
\label{eq:psinun}
\spl{
\tilde{\psi}_\nu^{(n)}=&-\frac{\tilde{\psi}_\nu^{(n-1)'}}{i\nu}-\frac{3}{2i\nu }\bigg[\tilde{H}_s\tilde{\psi}_\nu^{(n-1)}+\tilde{H}_\nu^{(n-1)}\tilde{\psi}_s+\sum_{\ell=1}^{n-2}\sum_{\alpha\neq0,\nu}\tilde{H}_{\alpha}^{(\ell)}\tilde{\psi}_{\nu-\alpha}^{(n-\ell-1)}\bigg]\\
&+\frac{3}{2i\nu}\bigg[\tilde{H}_s\tilde{\psi}_s^*\delta_{\nu,2}\delta_{n,2}+(1-\delta_{\nu,2})\big(\tilde{H}_s\tilde{\psi}_{2-\nu}^{(n-1)}{}^*+\tilde{H}_{\nu-2}^{(n-1)}\tilde{\psi}_s^*\big) +\sum_{\ell=1}^{n-2}\sum_{\alpha\neq0,\nu-2}\tilde{H}_{\alpha}^{(\ell)}\tilde{\psi}_{2+\alpha-\nu}^{(n-\ell-1){}^*}\bigg]\,,
}
}
where $n\geq1$. Similarly From Eq.~\eqref{eq:Hnu} for $\nu\neq0$ we get
\eq{
\label{eq:Hnun}
\tilde{H}_\nu^{(n)}=\frac{1}{6\tilde{H}_s}\bigg[\tilde{\psi}_s\tilde{\psi}^{(n)}_{-\nu}{}^*+\tilde{\psi}_s^*\tilde{\psi}^{(n)}_{\nu}+\sum_{\ell=1}^{n}\sum_{\alpha\neq0,\nu}\tilde{\psi}_\alpha^{(\ell)}\tilde{\psi}_{\alpha-\nu}^{(n+1-\ell)}{}^*\bigg]-\frac{1}{2\tilde{H}_s}\sum_{\ell=1}^{n}\sum_{\alpha\neq0,\nu}\tilde{H}_\alpha^{(\ell)}\tilde{H}_{\nu-\alpha}^{(n+1-\ell)}\,.
} 
Note that since the rescaled variables $\tilde{\psi}$ and $\tilde{H}$ are first order in $\epsilon$ we expect the nonzero modes $\tilde{\psi}_\nu$ and $\tilde{H}_\nu$ to start from the second order. We will see below that it is in fact the case. 
By using Eqs.~\eqref{eq:Hnun} and \eqref{eq:psinun} we see that 
\eq{
\tilde{\psi}_\nu^{(1)}=\tilde{H}_\nu^{(1)}=0\,,
}
and for the leading nonzero modes we find
\begin{eBox}
\eq{
\label{eq:psinu1}
\tilde{\psi}_\nu^{(2)}=
\frac{3}{4i}\tilde{H}_s\tilde{\psi}_s^* \, \delta_{\nu,2}\, , \qquad 
\tilde{H}_\nu^{(2)}=\frac{1}{8i}\tilde{\psi}_s^*{}^2 \, \delta_{\nu,2}-\frac{1}{8i}\tilde{\psi}_s^2 \, \delta_{\nu,-2}\,.
}
\end{eBox}
In obtaining this result there was no need to solve a differential equation, since the time derivative term is  subleading, contributing to higher orders. The relation between nonzero modes and the zero modes is one of the main results of our paper.
We can easily proceed and compute higher-order terms in $\epsilon$. After some algebra we find that for $n=3$, 
\eq{
\spl{
\tilde{\psi}_\nu^{(3)}=&\left(\frac{3}{8}(\tilde{H}_s\tilde{\psi}_s^*)'+\frac{9}{16}\tilde{H}_s^2\tilde{\psi}_s^*+\frac{3}{32}|\tilde{\psi}_s|^2\tilde{\psi}_s^*\right) \, \delta_{\nu,2}+\frac{3}{32}\tilde{\psi}_s^3 \, \delta_{\nu,-2}-\frac{3}{64}\tilde{\psi}_s^*{}^3 \, \delta_{\nu,4}\,,
}
}
and similarly for the Hubble parameter,
\eq{
\spl{
\tilde{H}_\nu^{(3)}=&\frac{1}{6\tilde{H}_s}\left(\frac{3}{8}(\tilde{H}_s\tilde{\psi}_s^*)'\tilde{\psi}_s^*+\frac{9}{16}\tilde{H}_s^2\tilde{\psi}_s^*{}^2+\frac{3}{32}|\tilde{\psi}_s|^4+\frac{3}{32}\tilde{\psi}_s^4\right) \, \delta_{\nu,2}\\
&+\frac{1}{6\tilde{H}_s}\left(\frac{3}{8}(\tilde{H}_s\tilde{\psi}_s)'\tilde{\psi}_s+\frac{9}{16}\tilde{H}_s^2\tilde{\psi}_s{}^2+\frac{3}{32}|\tilde{\psi}_s|^4+\frac{3}{32}\tilde{\psi}_s^*{}^4\right) \, \delta_{-\nu,2}\,.
}
}
 Since in this section we will only consider corrections up to order $n=3$ for the background equations, we can use the leading-order equations \eqref{eq:leadbck} for the slow modes  to simplify the above relations to the same order, which yields
\begin{eBox}
\eq{
\label{eq:psinu2}
\tilde{\psi}_\nu^{(3)}=-\frac{3}{32}|\tilde{\psi}_s|^2\tilde{\psi}_s^* \, \delta_{\nu,2}+\frac{3}{32}\tilde{\psi}_s^3 \, \delta_{\nu,-2}-\frac{3}{64}\tilde{\psi}_s^*{}^3 \, \delta_{\nu,4}\,,\qquad\tilde{H}_\nu^{(3)}=0\,.
}
\end{eBox}
Now that we have obtained the solutions for nonzero modes, we may substitute into the equations of motion for the slow modes, Eqs.~\eqref{eq:psis0} and \eqref{eq:Hs0}. Considering terms up to order $n=3$ we have
\begin{equation}
\label{eq:psisorders}
\begin{split}
&\tilde{\psi}_s'
+\frac{3}{2}\tilde{H}_s\tilde{\psi}_s
\\
&-\frac{3}{2}\left(\tilde{H}_s\tilde{\psi}_2^{(2)}{}^*+\tilde{H}_{-2}^{(2)}\tilde{\psi}_s^*\right)\\
&-\frac{3}{2}\bigg(\tilde{H}_s\tilde{\psi}_2^{(3)}{}^*+\tilde{H}_{-2}^{(3)}\tilde{\psi}_s^*-\tilde{H}_2^{(2)}\tilde{\psi}_2^{(2)}\bigg)+\order{\epsilon^5}=0\,,
\end{split}
\end{equation}
and
\begin{equation}
\begin{split}
&3\tilde{H}_s^2-|\tilde{\psi}_s|^2\\
&+3\left(\tilde{H}_{2}^{(2)}\right)^2+3\left(\tilde{H}_{-2}^{(2)}\right)^2-\vert \tilde{\psi}_2^{(2)}\vert^2+\order{\epsilon^5}=0\,,
\end{split}
\end{equation}
where we have written the terms that are at the same order in the same line and we have excluded the terms that contain nonzero modes at order $n=1$ as they vanish identically.  After substitution of nonzero modes we finally obtain the effective equations of motion for the slow modes. For convenience we write the equations in this last step using the original, dimensionful variables:
\begin{eBox}
\eqa{
\label{eq:EFTbck1}
&i\dot{\bar{\psi}}_s+i\frac{3}{2}H_s\bar{\psi}_s+\frac{9}{16}\frac{|\bar{\psi}_s|^2}{\Mpl^2} \bar{\psi}_s+i\frac{9}{32}\frac{H_s}{m}\frac{|\bar{ \psi}_s|^2}{\Mpl^2}\bar{\psi}_s+\order{\epsilon^3}H_s\bar{\psi}_s=0 \, ,\\
\label{eq:EFTbck2}
&3\Mpl^2H_s^2=m|\bar{\psi}_s|^2+\frac{3}{32\Mpl^2}|\bar{\psi}_s|^4+\order{\epsilon^4}m|\bar{\psi}_s|^2\,.
} 
\end{eBox}
This is one of our main results in this paper. Note that, interestingly, the backreaction of the nonzero modes on the slow mode effectively induces a kind of self-interaction for the slow mode and a new contribution to the energy density, which sources the Hubble parameter. (We will consider characteristics of the induced self-interaction below, after deriving our EFT for the fluctuations.) We expect these nontrivial corrections to the Schr\"odinger and Friedmann equations to help improve the validity of the simplified nonrelativistic equations when compared to the exact equations. In Sec.~\ref{sec:num}, where we present the numerical solutions, we will confirm  that this is indeed the case, and illustrate that the solutions of our EFT follow the exact solutions closely. We can also provide further insight into the above results by describing the system as a fluid in an expanding background, which will be done in Sec.~\ref{sec:fluid}. 

\subsubsection*{The scale factor}
Since in a flat FLRW universe, the background equations do not depend on the scale factor explicitly, we did not need to solve for it. However, since the scale factor would show up in the equations governing small fluctuations, we would need to have a perturbative expansion of it as well. A simple way to achieve this is to use its differential equation $\dot{a}=aH$, which, after the mode expansion and using the rescaled variables, results in   
\begin{equation}
a'_\nu+i\nu a_\nu=a_\alpha \tilde{H}_{\nu-\alpha}\,.
\end{equation}
As before, we perform a perurbative expansion for nonzero modes, and it is easy to see that $a_\nu^{(1)}=0$ for all $\nu$. For $n>1$ we have
\begin{equation}
a_\nu^{(n)}=-\frac{a_\nu^{(n-1)}{'}}{i\nu}+\frac{a_s\tilde{H}_\nu^{(n)}}{i\nu}+\frac{a_\nu^{(n-1)}\tilde{H}_s}{i\nu}+\sum_{\ell=1}^{n-1}\sum_{\alpha\neq\{0,\nu\}}\frac{a_\alpha^{(\ell)}\tilde{H}_{\nu-\alpha}^{(n-\ell)}}{i\nu}\,.
\end{equation}
Specifically, for $n=2$ the we have
\begin{eBox}
\begin{equation}
a_\nu^{(2)}=\frac{1}{2i}\tilde{H}_\nu^{(2)}a_s=\left(-\frac{1}{16}\tilde{\psi}_s^*{}^2 \, \delta_{\nu,2}-\frac{1}{16}\tilde{\psi}_s^2 \, \delta_{\nu,-2}\right)a_s\,.
\end{equation}
\end{eBox}
Since, for the fluctuations, we work up to $n=2$ in perturbative expansion, we only need the scale factor up to this order. Therefore, we neglect higher-order corrections to the scale factor.  
Finally, the effective equation for the slow-mode scale factor $a_s=\ev{a}$ to this order becomes
\begin{equation}
\dot{a}_s=a_s H_s+\order{\epsilon^2}a_sH_s\, ,
\end{equation}
that is, to this order, the slow mode of the scale factor does not receive any correction from nonzero modes. Before concluding this section, we define variables for the inverse of the scale factor and its square for later use:  
\eq{
\label{eq:r}
q(t)\equiv\frac{1}{a(t)}\,,\qquad r(t)\equiv\frac{1}{a(t)^2}\,.
}
Then it is straightforward to show that $q_\nu^{(1)}=r_\nu^{(1)}=0$ and for $n=2$ we have $q_\nu^{(2)}=-a_\nu^{(2)}/a_s^2$ and $r_\nu^{(2)}=-2a_\nu^{(2)}/a_s^3$.

\subsection{Effective field theory for fluctuations}\label{sec:eftpert}
We now proceed to find the effective equations of motion for the slow modes of the linear fluctuations around FLRW background. The procedure is similar to the last section. We start by applying the mode expansion to Eqs.~\eqref{eq:linpert1} and \eqref{eq:linpert2}. Note that we need to use background quantities and their mode expansion as well. We again use the rescaled variables defined in Eq.~\eqref{eq:rescale} as well as
\eq{
\tilde{\delta\psi} (t, {\bf x} )\equiv\frac{\delta\psi (t, {\bf x}) }{\sqrt{m}\Mpl}\,.
}
Furthermore, to make all the expressions dimensionless (which makes the power counting more straightforward) we rescale the spatial coordinates by $\tilde{x}=mx$ and use the notation $\tilde{\nabla}^2$ for the corresponding comoving spatial Laplacian operator. Our aim is to compute the effective equations of motion for $\delta\psi_s=\ev{\delta\psi}$ and $\Phi_s=\ev{\Phi}$. The relevant equations for general mode $\nu$ from Eqs.~\eqref{eq:linpert1} and \eqref{eq:linpert2} in rescaled variables read
\eqa{
\label{eq:dpsinu}
&i\tilde{\delta\psi}'_\nu-\nu\tilde{\delta\psi}_\nu+(\tilde{\mathfrak{D}}_\alpha\tilde{\delta\psi}_{\nu-\alpha}+\tilde{\mathfrak{D}}^*_{-\alpha}\tilde{\delta\psi}^*_{2+\alpha-\nu}) \nonumber \\ 
&\quad\quad\quad  -\Phi_{\alpha}(\tilde{\psi}_{\nu-\alpha}+\tilde{\psi}^*_{2+\alpha-\nu}-2i\tilde{H}_\beta(\tilde{\psi}_{\nu-\alpha-\beta}-\tilde{\psi}^*_{2+\alpha+\beta-\nu}))=0\,,\\
\label{eq:Phinu}
&\Phi'_\nu+i\nu \Phi_\nu+\tilde{H}_\alpha\Phi_{\nu-\alpha}-\frac{1}{4i}(\tilde{\psi}_\alpha\tilde{\delta\psi}_{\alpha-\nu}^*-\tilde{\psi}_{-\alpha}^*\tilde{\delta\psi}_{\nu-\alpha}+\tilde{\psi}_\alpha\tilde{\delta\psi}_{\nu+2-\alpha}-\tilde{\psi}^*_{-\alpha}\tilde{\delta\psi}^*_{2+\alpha-\nu})=0\,,
}
where the rescaled operator is
\eq{
\tilde{\mathfrak{D}}_\nu=\frac{1}{2}\Big(3i\tilde{H}_\nu-r_\nu\tilde{\nabla}^2+2\tilde{\psi}^*_{-\beta}\tilde{\psi}_{\nu-\beta}-\tilde{\psi}_\beta\tilde{\psi}_{\nu+2-\beta}-\tilde{\psi}^*_{-\beta}\tilde{\psi}_{2+\beta-\nu}\Big)\,,
}
and we have used the notation defined in Eq.~\eqref{eq:r} for the mode expansion of $1/a^2$.
The Poisson equation, Eq.~\eqref{eq:pois}, leads to
\begin{equation}
\spl{
r_\alpha\tilde{\nabla}^2\Phi_{\nu-\alpha}=&\frac{1}{2}(\tilde{\psi}^*_{-\alpha}\tilde{\delta\psi}_{\nu-\alpha}+\tilde{\psi}_{\alpha}\tilde{\delta\psi}^*_{\alpha-\nu})+3\tilde{H}_{\nu-\alpha}(\Phi_{\alpha}'+i\alpha\Phi_{\alpha})\\
&+\frac{1}{4}(\tilde{\psi}_\alpha\tilde{\psi}_{\nu+2-\alpha-\beta}+\tilde{\psi}^*_{-\alpha}\tilde{\psi}^*_{2+\alpha+\beta-\nu}+2\tilde{\psi}^*_{-\alpha}\tilde{\psi}_{\nu-\alpha-\beta})\Phi_{\beta}\,.
}
\end{equation}
Note that we can eliminate the $\Phi'_\alpha$ term by using Eq.~\eqref{eq:Phinu}.

Once again we have an infinite set of coupled equations and to proceed we have to identify small parameters in the problem. Besides the previously introduced parameters $\epsilon_{\scriptscriptstyle H}$ and $\epsilon_t$, there is another one due to the fact that we are now dealing with field fluctuations with nonzero momenta. In the nonrelativistic limit we expect $\epsilon_k \sim (k^2/m^2a^2)$ to be small, where $k$ is the comoving momentum of the nonrelativistic field. Furthermore, in order to make the power counting in the perturbative expansion tractable, we take into account the fact that the rescaled fluctuations are small  (compared to the background quantities), which we quantify by yet another small parameter $\epsilon_g$, so that we have $\tilde{\delta\psi}=\order{\epsilon_g}$. It is important to note that since we are working only to linear order in spatially varying quantities, the new parameter $\epsilon_g$ becomes irrelevant. Nevertheless, introducing this parameter is helpful for the perturbative expansion of nonzero modes, when we must evaluate the size of each term. Therefore, our EFT for the field fluctuations would be an expansion in four small parameters, $\epsilon=\{\epsilon_t,\epsilon_{\scriptscriptstyle H},\epsilon_k,\epsilon_g\}$. Note that we have not made any assumption about the possible hierarchy between these parameters. In particular, fluctuations can be both at sub-horizon and super-horizon scales without affecting the following analysis, as long as the Hubble parameter and the physical wavenumber are small compared to the mass of the field. Note that in terms of rescaled variables, we have a useful order of magnitude relation $\tilde{\delta\psi}_s\sim\Phi_s=\order{\epsilon}$.   

Similar to what has been done in the previous section, by using the perturbative expansion in Eq.~\eqref{eq:pertexp}, we can solve for the relevant nonzero modes and find the effective equations of motion for the slow modes. Again, similar to the background case, we can deduce that at order $n=1$ we have
\eq{
\tilde{\delta\psi}_\nu^{(1)}=\Phi_\nu^{(1)}=0\,,
}
and the leading nonzero modes corresponding to $n=2$ can be obtained by using Eqs.~\eqref{eq:dpsinu} and \eqref{eq:Phinu}
\begin{eBox}
\eq{
\label{eq:Phinu2}
\tilde{\delta\psi}_{\nu}^{(2)}=\left(\frac{3}{4i}\tilde{H}_s\tilde{\delta\psi}_s^*+\frac{\tilde{\nabla}^2\tilde{\delta\psi}_s^*}{4a_s^2}-\frac{1}{2}\tilde{\psi}_s^*\Phi_s\right)\delta_{\nu,2}\,,\quad
\Phi_\nu^{(2)}=\frac{1}{8}\bigg(\tilde{\psi}_s^*\tilde{\delta\psi}_s^*\delta_{\nu,2}+\tilde{\psi}_s\tilde{\delta\psi}_s\delta_{\nu,-2}\bigg)\,.
}
\end{eBox}
We will present the equations for higher orders ($n>2$) in Appendix~\ref{app:linhigher} and limit ourselves to order $n=2$ (and hence up to $\order{\epsilon}$ corrections to the equations for the slow-modes) in this section. As we will see, even at this order we find nontrivial terms in the effective equations of motion.
Having found the solutions for the leading-order nonzero modes, the effective equations governing the slow mode fluctuations $\delta\psi_s$ and $\Phi_s$ can be obtained by setting $\nu=0$ in Eqs.~\eqref{eq:dpsinu} and \eqref{eq:Phinu} and replacing the nonzero modes of the background and fluctuation variables by their corresponding solutions. Returning to the original variables at the last step, this procedure results in
\begin{eBox}
\eq{
\spl{
	\label{eq:eftsch}
	&i\dot{\delta\psi_s}+i\frac{3}{2}H_s\delta\psi_s+\frac{\nabla^2\delta\psi_s}{2ma_s^2}-m\bar{\psi}_s\Phi_s\\
	&+\frac{9}{8}\frac{|\bar{ \psi}_s|^2}{\Mpl^2} \delta\psi_s-\frac{7}{16}\frac{\bar{\psi}_s^2}{\Mpl^2}\delta\psi_s^*+\frac{\nabla^4\delta\psi_s}{8m^3a_s^4}+2iH_s \bar{\psi}_s \Phi_s+\order{\epsilon^2}H_s\delta\psi_s=0\,,
}
}
\end{eBox}
as the effective Schr\"odinger equation, and 
\begin{eBox}
\eq{
\spl{
	\label{eq:eftphidot}
	&\dot{\Phi}_s+H_s\Phi_s+\frac{i}{4\Mpl^2}(\bar{\psi}_s \delta\psi_s^*-\bar{\psi}_s^*\delta\psi_s)\\
	&+\frac{3H_s}{8m\Mpl^2}(\bar{\psi}_s\delta\psi_s^* + \bar{\psi}_s^*\delta\psi_s)+\frac{i}{8m\Mpl^2}\frac{\nabla^2}{2ma_s^2}(\bar{\psi}_s \delta\psi_s^*- \bar{\psi}_s^*\delta\psi_s)+\order{\epsilon^2}H_s\Phi_s=0\,,
}
}
\end{eBox}
as the equation for the slow-mode Newtonian potential $\Phi_s$. We also have the effective Poisson equation,
\begin{eBox}
\eq{
\spl{
\label{eq:eftpoisson}
\frac{\nabla^2\Phi_s}{a_s^2}=&\frac{m}{2\Mpl^2}(\bar{\psi}^*_s \delta\psi_s+ \bar{\psi}_s\delta\psi^*_s)\\
&+\frac{3i}{4\Mpl^2}H_s(\bar{\psi}^*_s \delta\psi_s- \bar{\psi}_s\delta\psi^*_s)-\frac{3}{2}H_s^2\Phi_s+\order{\epsilon^2}H_s^2\frac{\delta\psi_s}{\bar{\psi}_s}\,.
}
}
\end{eBox}
Once again we see that the nonzero modes induce nontrivial corrections to all equations for the slow modes. We stress that all these corrections have been mostly neglected in the literature to date. In Sec.~\ref{sec:fluid}, we provide an alternative description of the system in terms of an imperfect fluid, to extract useful effective quantities for the system, such as the sound speed and viscosity. Also, in Sec.~\ref{sec:num}, we confirm the validity of our EFT by numerical analysis. In Appendix~\ref{app:multi-component} we study the case of a multicomponent universe and in Appendix~\ref{app:interaction} we take into account the self-interaction of the scalar field, which has been neglected so far. 

In both Eq.~\eqref{eq:EFTbck1} for $\bar{\psi}_s$ and Eq.~\eqref{eq:eftsch} for $\delta \psi_s$, we find nonlinear terms in the effective equations of motion that take the form of self-interaction terms, even though we began with a free scalar field. To further characterize the induced self-interaction, we may consider Eqs.~\eqref{eq:eftbckslf} and \eqref{eq:fluctslambda} for the evolution of the slow modes $\bar{\psi}_s$ and $\delta \psi_s$ in the presence of a tree-level self-interaction of the form $V (\phi) = \lambda \phi^4 / (4!)$ for the original (relativistic) scalar field $\phi$. If we consider $\psi_s (t, {\bf x}) = \bar{\psi}_s (t) + \delta \psi_s (t, {\bf x})$ and work to first order in fluctuations, we may combine Eqs.~\eqref{eq:eftbckslf} and \eqref{eq:fluctslambda} to write
\begin{equation}
\begin{split}
i \dot{\psi}_s &+ \frac{3i}{2} H_s \psi_s + \frac{\nabla^2 \psi_s}{2m a_s^2}  - \frac{\lambda}{8 m^2} \vert \psi_s \vert^2 \psi_s + \frac{ 9}{16 \Mpl^2} \vert \psi_s \vert^2 \psi_s   \\
&\quad + \frac{\nabla^4 \psi_s}{8 m^3 a_s^4}  - m \bar{\psi}_s \left[ 1 - 2i \left( \frac{ H_s}{m} \right) \right] \Phi_s + {\cal O} (\epsilon^3 ) m \psi_s = 0.
\end{split}
\label{eq:eompsicombined}
\end{equation}
(In the limit of a rigid spacetime, with $a_s (t) \rightarrow 1$, $(m / \Mpl) \rightarrow 0$, and $\vert \Psi_s (t, {\bf x}) \vert \rightarrow 0$, Eq.~\eqref{eq:eompsicombined} matches the equation of motion found in Ref.~\cite{Namjoo:2017nia} for the slow mode $\psi_s (t, {\bf x})$ in Minkowski spacetime, to the appropriate order.) The form of Eq.~\eqref{eq:eompsicombined} suggests that the {\it induced self-interaction strength}, even in the absence of a tree-level self-coupling, takes the form
\begin{equation}
    \lambda_{\rm eff} = - \frac{9}{2} \left( \frac{m}{\Mpl} \right)^2 + {\cal O} (\epsilon^2) .
    \label{eq:lambdaeff}
\end{equation}
We see immediately that the induced self-interaction is attractive $(\lambda_{\rm eff} < 0$) and that it arises gravitationally (being proportional to $G \sim 1/\Mpl^2$).\footnote{The right-hand side of the effective Friedmann equation in Eq.~\eqref{eq:EFTbck2} includes a contribution to $\rho_{\rm eff}$ proportional to $\lambda_{\rm eff} \vert \bar{\psi}_s \vert^4$, scaling with $\lambda_{\rm eff}$ and $\bar{\psi}_s$ in the way one would naively expect from the form of the induced self-coupling, though with a different sign and overall coefficient. The term proportional to $\vert \bar{\psi}_s \vert^4$ in Eq.~\eqref{eq:EFTbck2} arises from backreaction of rapid oscillations on $H_2^{(2)}$ and $H_{-2}^{(2)}$. How best to incorporate global gravitational effects arising from the gravitationally induced self-interactions remains the subject of further study.} 

The emergence of a gravitationally induced self-coupling could be of interest in the context of self-interacting dark matter \cite{Bullock:2017xww,Fan:2016rda,Tulin:2017ara}, though estimating the magnitude of the effects remains subtle. To address various astrophysical observations, models of self-interacting dark matter typically introduce interactions that yield $\sigma / m \lesssim 1 \, {\rm cm}^2 / g = 4.6 \times 10^3 \, {\rm GeV}^{-3}$, where $\sigma$ is the total scattering cross section. Given the form of $\lambda_{\rm eff}$ in Eq. (3.41), the induced self-interactions we have identified appear, at least naively, to be much too weak to account for such interactions. However, a direct comparison between the cross section deduced from our EFT and the cross section required in self-interacting dark matter models is complicated for at least two reasons. First, our EFT calculation is only valid up to linear order in spatially varying fluctuations, whereas the self-interactions within dark matter halos require a fully nonlinear analysis. Second, the very high occupation numbers of particles in ultra-light dark matter models implies that classical field equations yield a more reliable description of the system \cite{Deng:2018jjz} than naive scattering amplitudes calculated within quantum field theory. We leave a more systematic study of this interesting topic for future work. As we will see in the next section, however, despite the small magnitude of the gravitationally induced self-coupling, the backreaction effects captured in our effective description lead to interesting -- and in principle observable -- phenomenological features, including nontrivial pressure, sound speed, and bulk viscosity.

As a final --- but important --- remark, notice that a by-product of our EFT is the ability to construct the full solution: Once the slow modes are found by (perhaps numerically) solving the effective equations, we are able to find the full solution, including its oscillatory behavior. For example, the scalar field $\psi$ can be constructed via $\psi = \sum_\nu \psi_\nu\, e^{i\nu mt}$, where $\psi_\nu$ for $\nu\neq 0$ are obtained order by order by our perturbative prescription and are given once the solution for $\psi_s$ is obtained from the EFT.  In Sec.~\ref{sec:num} we show that this procedure leads to results that match to the exact solution with very good accuracy (and one can obtain even more accurate results by going to higher orders in the perturbation theory). In situations where the oscillatory behavior is of interest, this procedure is expected to be much more efficient compared to solving the exact equations. This is because the exact equations of motion, involving rapidly oscillating factors, are generically expected to suffer from stiffness and instabilities. On the other hand, the EFT equations are well behaved and are expected to be solved easily by standard numerical algorithms. Such a theoretical framework would be appropriate for a variety of situations, relevant to cosmology and astrophysics. Examples of such situations are the change in the orbits of planets and stars in the dark matter halo \cite{Boskovic:2018rub} or the resonances in binary pulsars \cite{Blas:2019hxz} as results of the  oscillations of dark matter.

\section{Effective fluid description} 
\label{sec:fluid}
Having derived an EFT for the system, it is instructive to find a more intuitive way to interpret the results. In this section, we obtain an equivalent description of the system in terms of an imperfect fluid and argue that the higher-order terms in our EFT can be understood as new contributions to the variables describing the fluid. 

To start, we may interpret the equations \eqref{eq:EFTbck1} and \eqref{eq:EFTbck2} as a universe (effectively) expanding with rate $H_s$ and filled with some effective fluid. As a result, we may consider the right-hand side of the Friedmann equation, Eq.~\eqref{eq:EFTbck2}, as an effective energy density of the fluid: 
\begin{eBox}
\eq{
	\rho_{\text{eff}} = m \vert \bar{\psi}_s \vert^2 +\dfrac{3}{32\Mpl^2}\vert \bar{\psi}_s \vert^4\,.
	\label{eq:rho_eff}
}
\end{eBox}
Using the Schr\"odinger equation, Eq.~\eqref{eq:EFTbck1}, we can then derive a continuity equation for the fluid in an FLRW background,
\eq{
	\dot \rho_{\text{eff}} +3H_s(\rho_{\text{eff}}+p_{\text{eff}})=0\,,
}
where the effective pressure turns out to be
\begin{eBox}
\begin{equation}
p_{\text{eff}}=\frac{9}{32 \Mpl^2}|\bar{\psi}_s|^4\,.
\label{peff}
\end{equation}
\end{eBox}
Notice that the effective pressure is of order $\sim\order{\epsilon^2}\rho_{\text{eff}}$, so to $\order{\epsilon}\rho_{\text{eff}}$ the effective fluid behaves like pressureless matter, consistent with cold dark matter.\footnote{ The effective pressure in Eq.~(\ref{peff}) should not be confused with the smeared pressure of the scalar field appearing in the energy-momentum tensor at the level of the background, $\ev{p_\phi}= \langle \frac{1}{2} \dot \phi^2 -\frac{1}{2}m^2 \phi^2 \rangle$, as is done for example in Ref.~\cite{Blas:2019hxz}. The reason is that in the former, we absorbed another contribution from the left-hand side of the Friedmann equation, which appears as a result of integrating out the nonzero modes of the Hubble parameter, such as $H_2$ and $H_{-2}$. } Nonetheless, the rapidly oscillating modes induce a small effective pressure, which yields an effective equation of state $w_{\rm eff} = p_{\rm eff} / \rho_{\rm eff}$ of the form
\begin{equation}
    w_{\rm eff} = \frac{ 27}{32} \left( \frac{ H_s }{m} \right)^2 + {\cal O} (\epsilon^3) ,
\end{equation}
upon using Eq.~(\ref{eq:EFTbck2}) for $H_s$. It is evident that the effective pressure is a purely gravitational effect which induces a sort of interaction in the fluid, even though the original theory involves only a free scalar field.  


We aim to make a similar analogy between field fluctuations in our EFT and fluctuations of the fluid. Once we include spatially varying field fluctuations, our corresponding fluid description will feature an imperfect fluid. For such a description, we must incorporate the bulk viscosity, parameterized by the coefficient $\zeta$.\footnote{In principle, we expect other variables in the imperfect fluid, such as the shear viscosity and the effect of heat transfer, to appear as well. This is because such contributions are consistent with --- and hence allowed by --- the symmetries of the problem. However, as we will see, these additional variables are not required to fully describe the low-energy system under study to working order, although they may show up at higher orders, neglected in this paper.}
In Appendix~\ref{app:vis} we present some equations governing an imperfect fluid with bulk viscosity. As can be seen from the results of Appendix~\ref{app:vis}, at the background level, the pressure and bulk viscosity are degenerate and always appear in the form $p-3H\zeta$. It is this combination that we denoted by $p_{\text{eff}}$. However, when we incorporate fluctuations, the degeneracy will be broken. Note that the effective equation of motion for fluctuations, Eq.~(\ref{eq:eftsch}), has been obtained to $\order{\epsilon}$. At this order, the background fluid is effectively pressureless. As a result, in the equations for fluctuations of a viscous fluid we set $p_{\text{eff}}=0$ or $p=3H\zeta$ wherever these background quantities appear. 

We may define the {\it comoving overdensity} of the fluid, $\delta_{\rm eff}$, in terms of the right-hand side of the effective Poisson equation, Eq.~\eqref{eq:eftpoisson}:
\eq{
	\frac{\nabla^2\Phi_s}{a_s^2}=\frac{1}{2\Mpl^2}\rho_{\text{eff}} \, \delta_{\text{eff}}\,,	
}
where $\rho_{\text{eff}}$ is the effective energy density of the background in Eq.~\eqref{eq:rho_eff}, which, to ${\cal O}(\epsilon )$, is $\rho_{\text{eff}}=m|\bar{\psi}_s|^2$. This yields
\eq{
\delta_{\text{eff}}=\left(\frac{\delta\psi_s}{\bar{\psi}_s} + \frac{\delta\psi_s^*}{\bar{\psi}_s^*}\right) + i\frac{3H_s}{2m}\left(\frac{\delta\psi_s}{\bar{\psi}_s} - \frac{\delta\psi_s^*}{\bar{\psi}_s^*}\right)-\Phi_s\,.
\label{eq:deltaeffdef}
}
Note that although $\delta_{\rm eff}$ is constructed from complex quantities, the combination in Eq.~(\ref{eq:deltaeffdef}) remains real. Other fluid fluctuations, such as the effective fluctuations in density, velocity, and pressure, can also be derived by comparing the effective equations in Sec.~\ref{sec:EFT} with the fluid equations outlined in Appendix~\ref{app:vis}. This results in 
\eqa{
	\label{eq:drhoeff}
&\delta\rho_{\text{eff}}=m|\bar{\psi}_s|^2 \left( \frac{\delta\psi_s}{\bar{\psi}_s} + \frac{\delta\psi_s^*}{\bar{\psi} _s^*}-\Phi_s\right)\,,\\
\label{eq:dueff}
&\delta u_{\text{eff}}=\frac{1}{2mi}\left(\frac{\delta\psi_s}{\bar{\psi}_s} - \frac{\delta\psi_s^*}{\bar{\psi}_s^*}\right) + \frac{3H_s}{4m^2}\left(\frac{\delta\psi_s}{\bar{\psi}_s} + \frac{\delta\psi_s^*}{\bar{\psi}_s^*}\right) + \frac{1}{2mi}\frac{\nabla^2}{4m^2a_s^2}\left(\frac{\delta\psi_s}{\bar{\psi}_s} - \frac{\delta\psi_s^*}{\bar{\psi}_s^*}\right)\,,\\
&\delta
\label{eq:dpeff} p_{\text{eff}}=m|\bar{\psi} _s|^2 \left[ \left( -\frac{\nabla^2}{4m^2a_s^2} - \frac{\nabla^4}{8m^4a_s^4} + \frac{21H_s^2}{16m^2}\right)\left(\frac{\delta\psi_s}{\bar{\psi}_s}+\frac{\delta\psi_s^*}{\bar{\psi}_s^*} \right) - \frac{3iH_s\nabla^2}{8m^3a_s^2}\left(\frac{\delta\psi_s}{\bar{\psi}_s} - \frac{\delta\psi_s^*}{\bar{\psi}_s^*}\right)\right]
.}
Note that the above gauge-dependent variables are written in the Newtonian gauge. (See Appendix \ref{app:gauge} for relevant gauge transformations and corresponding expressions in the {\it time-averaged comoving} gauge.)

Using Eq.~\eqref{eq:deltaeffdef} and Eqs.~\eqref{eq:eftsch} and \eqref{eq:eftphidot} we can derive a second-order differential equation for $\delta_{\text{eff}}$, from which we can read new variables for the fluid, by making the analogy between the resulting equation and the standard second-order equation for an imperfect fluid outlined in  Appendix~\ref{app:vis} (see Eq.~\eqref{eq:ddot_fluidp0}). In this way, we obtain
\eq{
	\label{eq:eftdelta}
	\ddot{\delta}_{\text{eff}}+2H_s\dot{\delta}_{\text{eff}}-c^2_{\text{eff}}\frac{\nabla^2\delta_{\text{eff}}}{a^2}-\frac{\zeta_{\text{eff}}}{\rho_{\text{eff}}}\frac{\nabla^2\dot{\delta}_{\text{eff}}}{a_s^2}=\frac{\rho_{\text{eff}}}{2\Mpl^2}\delta_{\text{eff}}\,,
}
where we have defined the effective speed of sound by
\begin{eBox}
\eq{
	\label{eq:eftcs}
	c^2_{\text{eff}}=\frac{k^2}{4m^2a_s^2}-\frac{k^4}{8m^4a_s^4}+\frac{15}{16}\frac{H_s^2}{m^2}\,,
}
\end{eBox}
and the effective coefficient of bulk viscosity by
\begin{eBox}
\eq{
	\label{eq:eftzeta}
	\zeta_{\text{eff}}=-\frac{H_s}{2m^2}\rho_{\text{eff}}\,.
}
\end{eBox}
The bulk viscosity coefficient may be re-expressed in terms of another quantity with dimension of velocity \cite{Carrasco:2012cv}:
\ba 
c_{bv}^2 \equiv  \dfrac{H_s}{\rho_{\text{eff}}} \zeta_{\text{eff}}=-\dfrac{H_s^2}{2m^2},
\ea 
which may be compared with $c^2_{\text{eff}}$ in Eq.~\eqref{eq:eftcs}.
	
Some remarks are in order regarding these results. The leading term in the effective speed of sound in Eq.~\eqref{eq:eftcs} is well known \cite{Hwang:2009js,Poulin:2018dzj}. However, we find two additional contributions. The second term is a higher-order, momentum-dependent contribution. A similar term also appears in the analyses of Refs.~\cite{Hwang:2009js,Poulin:2018dzj} but with a different coefficient. (Our term is larger than the result of Refs.~\cite{Hwang:2009js,Poulin:2018dzj} by a factor of 2.) We trace this discrepancy to the fact that in Refs.~\cite{Hwang:2009js,Poulin:2018dzj}, the backreaction of the nonzero modes has been neglected, though they contribute at the same order. Note also that, despite the claims of Ref.~\cite{Hwang:2009js} (and also recently in Ref.~\cite{Lague:2020htq}), we do not expect these results to hold for arbitrary momentum, since the whole formalism breaks down as soon as the momentum of the field becomes comparable to its mass. However, in Minkowski spacetime, one can make a nonlocal field redefinition which yields an EFT that is nonperturbative in $\epsilon_k$ and holds for arbitrary momentum. Upon doing so, we confirm our coefficient, which is different from that obtained in Refs.~\cite{Hwang:2009js,Poulin:2018dzj}. (See Appendix~\ref{app:nonlocal} for details, where we also show that, thanks to the nonlocal field redefinition, the resulting sound speed approaches unity in the large-momentum limit, consistent with the expectation for a scalar field with canonical kinetic energy.) 

Finally, the last term in Eq.~\eqref{eq:eftcs}, which is a nontrivial, momentum-independent contribution from the background evolution, is new and has not been identified in previous analyses. This term shows that all field fluctuations experience an effective sound speed, even modes that are well beyond the horizon. Note that since the effective sound speed is always positive, it suppresses structure formation. Explicit investigation of this effect would be interesting but is beyond the scope of this paper.

Besides the sound speed, we have also derived the coefficient of the bulk viscosity in Eq.~\eqref{eq:eftzeta}, which has also been omitted in previous analyses. In principle one could have anticipated that such a term would appear from an EFT perspective, since it is consistent with the symmetries of the system. Note that this coefficient is \emph{negative}. From the second law of thermodynamics it can be shown that the coefficient of bulk viscosity of an isolated, imperfect fluid in thermal equilibrium must be positive \cite{Landau:1959}. However, our nonrelativistic system is {\it not} an isolated system: it exchanges energy with the relativistic sector of the field/fluid, and our effective description remains ignorant about the latter. Therefore, the sign in Eq.~(\ref{eq:eftzeta}) is not inconsistent with the second law of thermodynamics.\footnote{Taking both sectors into account simultaneously must result in a non-negative viscosity coefficient, and, indeed, a canonically normalized (relativistic) scalar field, re-expressed in terms of fluid dynamics, shows no viscocity. In a similar way, the usual energy conditions, such as the strong, weak, and dominant energy conditions, need not hold in a low-energy EFT, even if the underlying (relativistic) theory obeys them \cite{Creminelli:2006xe,Azhar:2018nol}.} A negative bulk viscosity would lead to an enhancement in the development of structures, competing with the positive sound speed. It would be interesting to investigate which contribution dominates, and whether the balance depends on time or length-scale. 

\section{Comparison with numerical solutions}\label{sec:num}
In this section we show that the effective description for the slow modes, as well as the nontrivial expressions for the nonzero modes obtained in Sec.~\ref{sec:EFT}, are consistent with the smearing of the exact solutions by the method discussed in Sec.~\ref{sec:mexp}. We also show that the EFT is able to construct the full oscillating solution with good accuracy, while our effective equations are well-behaved with no oscillatory terms. 


In Sec.~\ref{sec:mexp} we adopted the window function $W(t)$ of Eq.~\eqref{eq:w} for obtaining different modes of each variable. It was an appropriate choice for a mathematically rigorous formulation of the smearing procedure and for the exact mode decomposition. However, it is impractical for numerical purposes. The amplitude of the sinc function does not fall off sufficiently rapidly and, further, the naive truncation of the function $W (t)$ beyond the range of the simulation time results in the Gibbs phenomenon at discontinuities after a (numerical) Fourier transformation \cite{Proakis}. To efficiently circumvent these issues in our numerical smearing (as is well-known in the field of signal processing) we replace $W(t)$ with $W(t) {\cal B}(t/T)$, where ${\cal B}(t)$ is the {\it Blackman} window function \cite{Proakis} and $T$ is the total time of the simulation.
Fig.~\ref{fig:modes} depicts the results of this smearing procedure for the slow mode as well as the leading nonzero modes for the exact numeric solutions of $\psi$ and $\delta \psi$ (i.e. the solutions of Eqs.~\eqref{eq:bckg1}, \eqref{eq:bckg2}, \eqref{eq:linpert1} and \eqref{eq:linpert2}). In the same figure, we compare these results with the predictions of our EFT, the solutions of Eqs.~\eqref{eq:EFTbck1}, \eqref{eq:EFTbck2}, \eqref{eq:eftsch}, \eqref{eq:eftphidot}, and \eqref{eq:eftpoisson}, along with the nonzero modes obtained in Eqs.~\eqref{eq:psinu1}, \eqref{eq:psinu2}, and \eqref{eq:Phinu2} in terms of the slow modes. This confirms that our EFT predictions match the numerical smearing with good accuracy. Fig.~\ref{fig:modes} also shows the hierarchy between nonzero modes among themselves and with the slow mode as predicted by the EFT. 


As we saw in Sec.~\ref{sec:EFT}, all the variables can be represented by a mode expansion. For instance for the background variable $\bar{\psi} (t)$, to working order, we can write
\eq{
\label{eq:psiexp}
\bar{\psi} (t)= \bar{\psi} _s + \bar{\psi}_2 e^{2imt} + \bar{\psi}_4 e^{4imt}+ \bar{\psi}_{-2} e^{-2imt}+\dots\,,
}
while similar expressions hold for the other variables. 
\fg{
\centering
\includegraphics[width=0.45\textwidth]{nonzero_mode}
\hskip 1cm
\includegraphics[width=0.45\textwidth]{nonzero_mode_pert}
\caption{Absolute values of the slow mode as well as the leading nonzero modes for $\tilde{\psi}(t)$ (left) and $\tilde{\delta\psi}(t)$ (right) obtained by smearing the exact numerical solutions (solid, blue curve) compared with the predictions of our EFT (dashed, orange curve).  We have set $H(0)/m=0.06$, $k/m=0.02$ and  $a_s(0)=1$ which, from the EFT, implies $a(0)=0.998$. The initial conditions for the slow mode are $\tilde{\psi}_s(0)=0.1+0.025i$ and $\tilde{\delta\psi}_s(0)=(1-1.75i)\times10^{-6}$  which, from the EFT, imply $\tilde{\psi}(0)=0.098+0.020i$ and  $\tilde{\delta\psi}(0)=(0.326-1.61i)\times10^{-6}$ for the exact equations. The horizon crossing for the (randomly chosen) mode occurs around $\tilde{t}=25$.}
\label{fig:modes}
}
As mentioned earlier, since the nonzero modes are also expressed in terms of the slow mode, solving the obtained effective equations for the slow mode allows us to construct the full solution by relations like Eq.~\eqref{eq:psiexp} order by order in our perturbation theory. In Figs.~\ref{fig:ReImpsi} and \ref{fig:ReImdpsi} we compare the exact solutions of $\bar{\psi}$ and $\delta \psi$ with the solution for the slow mode in our EFT (labeled ``EFT, slow mode''), the full solution constructed out of our EFT (labeled``EFT, constructed''), and the ``naive'' theory (where all oscillatory terms are simply neglected in the equations of motion, as is typically done in analyses of axion-like dark matter models). It is evident from these plots that the naive theory can cause notable error, while the EFT solution follows the exact solution quite reliably. We, however, warn that the amount of deviation is sensitive to the small parameters based on which our EFT is constructed. (Here, for example, we take $H/m=0.06$ and $k/am=0.02$ at the initial time of the simulation, a somewhat random choice.) An interesting investigation, which we will leave for future work, would be to study the error that one  encounters by using the naive theory in a realistic situation where, e.g., the other components of matter in the universe are taken into account and the parameter space is chosen based on observational constraints. 


\fg{
\includegraphics[width=0.47\textwidth]{Re_psi}
\hspace{5mm}
\includegraphics[width=0.48\textwidth]{Im_psi}
\caption{Comparison of the real (left) and imaginary part (right) of $\tilde{\psi}$ between the exact numerical results, the full solution constructed out of EFT, the slow mode, and the naive solutions. The choice of parameters and initial conditions are the same as Fig.~\ref{fig:modes}.}
\label{fig:ReImpsi}
}
\fg{
\includegraphics[width=0.48\textwidth]{Re_dpsi}
\hspace{5mm}
\includegraphics[width=0.47\textwidth]{Im_dpsi}
\caption{Comparison of the real (left) and imaginary part (right) of $\tilde{\delta\psi}$ between the exact numerical results, the full solution constructed out of EFT, the slow mode, and the naive solutions. The choice of parameters and initial conditions are the same as Fig.~\ref{fig:modes}.}
\label{fig:ReImdpsi}
}

\fg{
\includegraphics[width=0.48\textwidth]{errbckg}
\hspace{5mm}
\includegraphics[width=0.47\textwidth]{errpert}
\caption{Relative error in variables $H$, $a$, $\delta_c$, and $\delta u$ for the naive theory compared with the EFT. The choice of parameters and initial conditions are the same as Fig.~\ref{fig:modes}.}
\label{fig:H}
}

A remark regarding the choice of initial conditions for comparison is in order. For the naive theory (which simply neglects any rapidly oscillating contribution) we have no option other than choosing the initial conditions to be the same as the exact theory. In the EFT, on the other hand, we have perturbative access to the full solution, as just described above, so that we can match the initial conditions of the full solution constructed out of EFT to the ones chosen for the exact solution.{\footnote{In fact, we found it more convenient to fix the initial conditions for the EFT first and then obtain the corresponding initial conditions for the exact solution. Note, however, that the inverse procedure (which is more realistic) is also possible.}} That is why in Figs.~\ref{fig:ReImpsi} and \ref{fig:ReImdpsi} the initial condition for the slow mode is different from the exact solution (since nonzero modes contribute to the initial conditions as well). Note that this is the natural choice of the initial conditions for our EFT, suggested by the EFT itself.  

Finally, having confirmed the validity of our EFT, we can consider the results of the EFT to be a sufficiently accurate description of the dynamics of the system, and then compare it with the naive theory in variables that are of observational interest. In Fig.~\ref{fig:H} we make this comparison for the Hubble parameter, the scale factor, the density contrast, and the velocity potential. Depending on the parameters and initial conditions, the error seems to be in fact observable.  The quantification of that statement in realistic situations, however, is beyond the scope of this paper and will be studied elsewhere. 


 
\section{Summary and Outlook}\label{sec:con}
In this paper we have obtained an EFT for a massive, nonrelativistic scalar field in an expanding background by systematically integrating out rapidly oscillating modes. We applied our formalism to spatially homogeneous quantities as well as to spatially varying fluctuations (working to linear order in the fluctuations), but the same methods can be employed for studying any system in which there are fast oscillations while the physically interesting variables are slowly varying in time. For the sake of convenient access, we summarize our main results here. The effective equations governing the dynamics of the system (in an expanding massive-field-dominated universe) are 
\eqa{
i\dot{\bar{\psi}}_s +\frac{3i}{2}H_s \bar{\psi}_s + \frac{9}{16}\left[1+\frac{iH_s}{2m}\right] \frac{|\bar{\psi}_s|^2}{\Mpl^2} \bar{\psi}_s \simeq 0\, , \qquad 
3\Mpl^2H_s^2\simeq m|\bar{\psi}_s|^2 + \frac{3}{32\Mpl^2}|\bar{\psi}_s|^4\, ,
} 
for the background variables (up to ${\cal O}(\epsilon^2)$), and

\eqa{
&i\dot{\delta\psi_s}+\frac{3i}{2}H_s\delta\psi_s + \left( 1+\frac{\nabla^2}{4m^2a_s^2}\right) \frac{\nabla^2\delta\psi_s}{2ma_s^2} - \left(1-\frac{2iH_s}{m}\right) m\bar{\psi}_s \Phi_s \nonumber \\
&\quad\quad\quad\quad + \frac{9|\bar{\psi}_s|^2}{8\Mpl^2}\delta\psi_s - \frac{7\bar{\psi}_s^2}{16\Mpl^2}\delta\psi_s^* \simeq 0\, ,\\
&\dot{\Phi}_s + H_s\Phi_s + \frac{1}{4\Mpl^2} \left[i\left( 1+\frac{\nabla^2}{4m^2a_s^2}\right) (\bar{\psi}_s \delta\psi_s^* - \bar{\psi}_s^* \delta\psi_s) +\frac{3H_s}{2m}(\bar{\psi}_s \delta\psi_s^* + \bar{\psi}_s^* \delta\psi_s) \right]\simeq 0\,,\\ \nonumber \\
&\left[\frac{\nabla^2}{a_s^2}+\frac{3}{2}H_s^2\right]\Phi_s\simeq  \frac{m}{2\Mpl^2} (\bar{\psi}^*_s\delta\psi_s + \bar{\psi}_s \delta\psi^*_s) + \frac{3iH_s}{4\Mpl^2} (\bar{\psi}^*_s \delta\psi_s - \bar{\psi}_s\delta\psi^*_s)\,,
} 
\\
for fluctuations (up to ${\cal O}(\epsilon)$), where, as a reminder, we parameterize the field as $\psi (t, {\bf x}) = \bar{\psi} (t) + \delta \psi (t, {\bf x})$. To leading order, these equations correspond to the already well-known Schr\"odinger and Schr\"odinger-Poisson equations for the background evolution and for the fluctuations, respectively. However, in our EFT, we have obtained nontrivial corrections as a result of integrating out the nonzero, rapidly oscillating modes (rather than neglecting them).  Furthermore, we have also interpreted the results more intuitively by describing the system as an effective (imperfect) fluid. To fully describe the system, to working order, we identify a nonzero effective pressure as well as  an effective sound speed and a bulk viscosity as follows:
\eq{
p_{\text{eff}}=\frac{9}{32 \Mpl^2}|\bar{\psi}_s|^4\,,\qquad c^2_{\text{eff}}=\frac{k^2}{4m^2a_s^2}-\frac{k^4}{8m^4a_s^4}+\frac{15}{16}\frac{H_s^2}{m^2}\,,\qquad \zeta_{\text{eff}}=-\frac{H_s}{2m^2}\rho_{\text{eff}}\,.
}
Note that the pressure and the bulk viscosity were missing in all previous analyses. Furthermore, the second term in the sound speed has a different numerical prefactor, compared to other results in the literature (see e.g., Refs.~\cite{Hwang:2009js,Poulin:2018dzj}). The discrepancy seems to be due to an extra contribution as a result of nontrivial effects of oscillatory (nonzero) modes, neglected in other studies. The last term in the sound speed was also missing in previous analyses. The size of the error arising by neglecting these terms requires further investigation in realistic situations, which is beyond the scope of this paper.

The derived EFT is interesting from both a theoretical and a practical point of view. On the theoretical side we can see that gravity induces an effective, attractive self-interaction in a free scalar field theory, which manifests as nontrivial pressure, sound speed, and viscosity. These effects can be important for the background evolution as well as for the growth of overdensities. Note that the effective sound speed is positive while the viscosity is negative, so that they act in opposite directions: the former suppresses the growth of overdensities while the latter tends to enhance it. It would be interesting to investigate in which situations the various variables win, and how their incorporation changes the results compared to the naive theory. We leave such a study for future work. 

From a practical point of view, as we saw in Sec.~\ref{sec:num}, whereas the naive treatment can deviate substantially from the exact results and hence may cause error in interpreting observations, our EFT remains quantitatively reliable. In addition, our method paves the way for efficiently obtaining accurate solutions (including oscillatory behavior) without having to solve the exact equations (which are expected to be stiff due to rapid oscillations). Note that simulating the exact theory numerically requires time increments $\Delta t<1/m$ in order to accurately capture effects of the oscillating terms, whereas the corresponding equations within our EFT remain well-behaved, with no rapidly oscillating terms, so that it is sufficient to use $\Delta t<{\text{min}}(1/H,a/k)$; roughly speaking, this yields an $\order{1/\epsilon}$ increase in efficiency. The full solution can then be constructed order by order using the mode decomposition outlined in Sec.~\ref{sec:EFT} with no difficulty. Such a method for solving differential equations containing rapidly oscillating terms can have interesting applications in much broader situations of scientific interest. 

There are a number of different directions --- besides the ones already mentioned --- that we would like to explore in future work. Such studies include the application of our EFT to predictions of possible impacts of ultra-light dark matter models on the CMB and other observations; the EFT in the nonlinear regime and the corresponding corrections to structure formation (cf. Ref.~\cite{Musoke:2019ima}) as well as the dynamics of celestial objects as they move through dark matter halos; and the extension of our EFT to the situation in which the nonrelativistic field is coupled to another (perhaps relativistic) dynamical field. 


\acknowledgments

We are grateful to Xingang Chen, Mohammad Ali Gorji, Mahdiyar Noorbala, Katelin Schutz, and Tracy Slatyer for helpful discussions. B.S. thanks YITP at Kyoto University for hospitality during the time this work was in progress. Portions of this work were conducted in MIT's Center for Theoretical Physics and supported in part by the U.S. Department of Energy under Contract No.~DE-SC0012567.


