Effective theories for a nonrelativistic field in an expanding universe: Induced self-interaction, pressure, sound speed, and viscosity

A massive, nonrelativistic scalar field in an expanding spacetime is usually approximated by a pressureless perfect fluid, which leads to the standard conclusion that such a field can play the role of cold dark matter. In this paper, we systematically study these approximations, incorporating subleading corrections. We provide two equivalent effective descriptions of the system, each of which offers its own advantages and insights: (i) A nonrelativistic effective field theory (EFT) with which we show that the relativistic corrections induce an effective self-interaction for the nonrelativistic field. As a byproduct, our EFT also allows one to construct the exact solution, including oscillatory behavior, which is often difficult to achieve from the exact equations. (ii) An effective (imperfect) fluid description, with which we demonstrate that, for a perturbed Friedmann-Lema\^{i}tre-Robertson-Walker (FLRW) universe: (a) The pressure is small but nonzero (and positive), even for a free theory with no tree-level self-interactions. (b) The sound speed of small fluctuations is also nonzero (and positive), reproducing already known leading-order results, correcting a subdominant term, and identifying a new contribution that had been omitted in previous analyses. (c) The fluctuations experience a negative effective bulk viscosity. The positive sound speed and the negative bulk viscosity act in favor of and against the growth of overdensities, respectively. The net effect may be considered a smoking gun for ultra-light dark matter.

1 Introduction Ultra-light particles have been considered a promising candidate for dark matter [1][2][3]. 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 [4]. 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. [5,6], 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 φ 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 φ particles must be very large. This, in turn, suggests that the system should be describable by a classical scalar field to a good approximation [5,7]. On the other hand, the dark matter is supposed to be cold, and therefore nonrelativistic. Putting these two statements together, we deduce that a 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 φ 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 [5] φ(t, x) = 1 √ 2m e −imt ψ(t, x) + e imt ψ * (t, x) , (1.1) where ψ is a slowly varying function of time, then (as we will demonstrate below) the energy density of φ in a Friedmann-Lemâitre-Robertson-Walker (FLRW) universe is given (to leading order) by 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 ψ. We therefore take the relation (1.1) as a field redefinition, with which to define the nonrelativistic field ψ. We will consider this field redefinition in more detail in the next section. Note that the ψ 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 φ coupled minimally to gravity, with action given by where the matter Lagrangian is The resulting Klein-Gordon equation in an expanding universe reads (1.5) 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. (1.5), this suggests three relevant, small quantities: where k p is a typical physical wavenumber. Furthermore, after removing φ in favor of the new field ψ as in Eq. (1.1), all remaining variables are expected to be slowly varying compared to the time-scale 1/m. This suggests a fourth small parameter, 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 λ = 0. As we will see, even in the absence of tree-level self-interactions for φ, 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 ψ. The induced self-coupling could be of interest in the light of astrophysical observations [8] which suggest that cold dark matter particles might depart from a purely collisionless regime [3,9,10], 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 D 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. [11], while Ref. [12] explores their impact on the orbital motion of stars in galactic halos. The resulting parametric resonance are also studied e.g. in Ref. [13] and recently in Ref. [14]. Ref. [15] 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 timedependent variables, which yield a 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. 2 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. 3, we derive our EFT for the background evolution as well as fluctuations. In Sec. 4 we re-express our EFT results in terms of fluid dynamics. In Sec. 5 we confirm the validity of our results by numerical analysis, and in Sec. 6 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 A we derive the Hamiltonian of the system after the field redefinition discussed in Sec. 2. In Appendix B we present the equations for higher-order nonzero modes, ignored in the results of Sec. 3. In Appendix C 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. 4. In Appendices D and E, we extend our results to the self-interacting case and to a multicomponent universe, respectively. In Appendix F 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 G we derive the gauge transformation that can be used to express our EFT (written in Newtonian gauge) in other gauges.

A suitable field redefinition
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. [6], we aim to take the nonrelativistic limit and obtain higher-order corrections to a low-energy EFT 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, ψ, related to φ as in Eq. (1.1). The prefactor 1/ √ 2m in Eq. (1.1) has been chosen so that |ψ| 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ödinger-like equation for ψ in Minkowski spacetime.
As discussed above, ψ 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 φ -and hence in the new field ψ -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. [6], the field redefinition of Eq. (1.1) may not be considered complete if the real and imaginary parts of ψ are treated as two independent variables, since that would consist of replacing one degree of freedom with two degrees of freedom. In Ref. [6], this is addressed by adding to Eq. (1.1) a comparable field redefinition for π, the conjugate momentum of φ, such that the transformation from (φ, π) to (ψ, ψ * ) becomes a canonical transformation. In addition, in Ref. [6] a nontrivial nonlocal operator is used within the transformation, which simplifies the equation of motion for ψ in Minkowski spacetime. In principle, a similar procedure can be employed when studying the behavior of ψ in a curved spacetime. However, as we discuss in Appendix F, 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. 1 Furthermore, making a transformation for π similar to what was used in Ref. [6] results in equations of motion different from the Schrödinger equation written in an FLRW background [5]. Finally, we note that in terms of the original field φ, the conjugate momentum differs fromφ 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. [6]. From Eq. (1.1), it is evident that there is a redundancy in the definition ψ of the form with η any real function of space and time. It is easy to show that this transformation leaves φ -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: Note that this choice of the gauge guarantees that the following constraint is satisfied in the new gauge: which allows us to have a first order, Schrödinger-like equation for ψ. Therefore, we may obtain the equations of motion for the new fields after gauge fixing by varying the following Lagrangian with respect to ψ and ψ * :  which shows that our field redefinition is invertible. Notice that the equation of motion for ξ is a differential equation, rather than an algebraic constraint, and it is not possible to substitute the solution of ξ 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 φ, we introduce a new field by the relation χ =φ (at the level of equations of motion); next we replaceφ with χ in the Lagrangian and add a Lagrange multiplier to guarantee the equivalence of the two theories: Note that at this stage, the field π is only a Lagrange multiplier and not the momentum conjugate of φ. The equations of motion for this Lagrangian are as follows: The first equation is the constraint. The second equation, together with χ =φ, 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, π). Finally, using the first two equations, the last one gives the correct equation of motion for φ. Therefore, the two theories L φ andL φ are equivalent.
In the theory governed byL φ , the parameter π is nondynamical, so one can replace it by its solution to simplify the Lagrangian. From Eq. (2.8) we have: Plugging this solution back into the Lagrangian results in where, from now on, we consider a non-self-interacting theory, i.e. we set λ = 0 in Eq. (1.4).
(We will consider tree-level self-interactions in Appendix D.) Again, one can show that the Lagrangian in Eq. (2.10) yields consistent equations of motion. The next step is to introduce the nonrelativistic fields ψ and ψ * , which we define by the following expressions: Here "c.c." denotes the complex conjugate of the immediately preceding expression, and the primes on indices in g µ ′ ν ′ indicate that µ and ν cannot both be time components, i.e., µ × ν = 0. (In Appendix A, we consider the Hamiltonian of this theory in terms of ψ and ψ * .) From Eq. (2.12), the equation of motion for ψ can be written as where we have defined the operator (2.14) As promised, the equation of motion for ψ is first order in time derivatives and is identical to the equation of motion from the constrained Lagrangian of Eq. (2.4). In Minkowski spacetime the Lagrangian reduces to 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, d 3 x|ψ| 2 , which is the number of particles, a quantity that is conserved in a nonrelativistic theory. The equation of motion reads iψ + ∇ 2 2m (ψ + e 2imt ψ * ) = 0 for Minkowski, (2.16) which, again after neglecting the oscillating term, is just the Schrödinger 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. 3. In a curved background, the complete set of equations also includes the Einstein equations, M Pl 2 G µν = T µν , where the energy-momentum tensor for the scalar field is in which we must insert definitions of Eqs. (1.1) and (2.5). To avoid clutter we do not write out T µν (ψ, ψ * ) explicitly here.

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 where Φ(t, 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, ψ(t, x) =ψ(t)+δψ(t, x), where the splitting becomes unambiguous by requiring that the spatial average of δψ(t, x) vanishes. The operator D defined in Eq. (2.14) up to linear order in Φ is  ) 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 O(Φ, δψ). 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ödinger-Friedmann and Schrödinger-Poisson equations for background and fluctuations, respectively [5]: 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 matterdominated universe, showing that a nonrelativistic scalar field can play the role of cold dark matter (to leading order; see Sec. 3 and Sec. 4 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. 3 The effective field theory in the nonrelativistic limit 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. [6] (but see also Refs. [17,18] for different, but equivalent, approaches). Here we consider the case in which the single field φ dominates the FLRW universe, and neglect any tree-level self-interactions. The effect of self-interaction shall be detailed in Appendix D while the EFT for the case of a multicomponent universe is studied in Appendix E.

Smearing and mode expansion
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Ẋ(t) ≪ mX(t) with m the mass of the scalar field. In other words, the spectrum of X(t) in the frequency domain,X(ω), would be localized around ω = 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. (2.20)-(2.25), which cause different frequencies to couple. Since the spectrum of the oscillatory factors like e iνmt (with ν an integer) would be a delta function at frequencies ω = νm, we expect that at higher orders a localization happens around the frequencies ω = νm. These subleading contributions would then backreact, again due to nonlinearities in the equations, on the slowly varying modes, whose spectra are localized around ω = 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. 1, we can define a 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): This is the temporal counterpart of the spatial average considered for example in Ref. [19].
To adopt a suitable window function, note that in Fourier space, the right-hand side of Figure 1. The window function W (t) (left) and its Fourier transform (right).
Eq. (3.1) takes the formX(ω)Ŵ (ω). As mentioned above, since the slow mode is localized around ω = 0 in the frequency domain, a natural choice for the window function would be a top-hat in Fourier space:Ŵ 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 ω = 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 (ω)e iωt = sin(mt/2) πt .
(3.3) Figure 1 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. 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 νm with nonzero integer ν) would develop. To capture these nonzero modes we can apply the same smearing operator and define the mode ν of the function by 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: 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. X(t) = X(t) . Note also that Eq. (3.5) 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 ν from (3.4) into (3.5). For later usage, note that the following relations hold for arbitrary functions of time X(t) and Y (t): 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 −ν = R * ν . 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 ν (t) are slowly varying in time. Furthermore, in the same limit, the fields are mainly concentrated around ω = 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 ǫ t ∼ |Ẋ/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) ≡ X ν=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 | ≫ |X ν | for ν = 0. 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 where for each n we have X (n) ν = O(ǫ n ) and ǫ is a small parameter, to be determined by the equations governing the dynamics of the system, as in Eqs. (1.6) and (1.7). 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 ǫ 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.

Effective field theory for the background
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. (2.20) and (2.21). For convenience, we introduce the rescaled, dimensionless quantities which we use in the course of calculations. We also rescale the timet = mt. Application of the mode expansion introduced in the previous section to Eqs. (2.20) and (2.21), with rescaled variables, results in the following equations for each mode ν: where the summation over repeated indices is understood and the prime denotes the derivative with respect tot. We are, in particular, interested in the equations for the slow mode quantitiesψ s = ψ and H s = H , with the corresponding equations after rescaling the variables: and which are obtained simply by setting ν = 0 in Eqs. (3.9) and (3.10). 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 ǫ t , from the Hubble scale 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 ǫ = {ǫ t , ǫ H }. In this section we will keep terms up to third order in ǫ but our formalism is general and can be continued to higher orders. Note that the rescaled variables introduced in Eq. (3.8) are first order in ǫ, i.e.ψ s = O(ǫ) andH s = O(ǫ), which makes the power counting in the derivation of EFT fairly easy. We also haveψ ′ s ∼ ǫψ s = O ǫ 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 O ǫ 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 O(ǫ) (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 (3.7) 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 (3.9) and for ν = 0 we havẽ where n ≥ 1. Similarly From Eq. (3.10) for ν = 0 we get (3.15) Note that since the rescaled variablesψ andH are first order in ǫ we expect the nonzero modesψ ν andH ν to start from the second order. We will see below that it is in fact the case. By using Eqs. (3.15) and (3.14) we see that and for the leading nonzero modes we find 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 ǫ. After some algebra we find that for n = 3, and similarly for the Hubble parameter, 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 (2.26) for the slow modes to simplify the above relations to the same order, which yields Now that we have obtained the solutions for nonzero modes, we may substitute into the equations of motion for the slow modes, Eqs. (3.11) and (3.12). Considering terms up to order n = 3 we haveψ and 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: 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ödinger and Friedmann equations to help improve the validity of the simplified nonrelativistic equations when compared to the exact equations. In Sec. 5, 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. 4.

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ȧ = aH, which, after the mode expansion and using the rescaled variables, results in As before, we perform a perurbative expansion for nonzero modes, and it is easy to see that a (1) Specifically, for n = 2 the we have 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 = a to this order becomesȧ 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: Then it is straightforward to show that q

Effective field theory for fluctuations
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. (2.22) and (2.23). Note that we need to use background quantities and their mode expansion as well. We again use the rescaled variables defined in Eq. (3.8) as well as Furthermore, to make all the expressions dimensionless (which makes the power counting more straightforward) we rescale the spatial coordinates byx = mx and use the notatioñ ∇ 2 for the corresponding comoving spatial Laplacian operator. Our aim is to compute the effective equations of motion for δψ s = δψ and Φ s = Φ . The relevant equations for general mode ν from Eqs. (2.22) and (2.23) in rescaled variables read where the rescaled operator is Note that we can eliminate the Φ ′ α term by using Eq. (3.32). 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 ǫ H and ǫ 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 ǫ k ∼ (k 2 /m 2 a 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 ǫ g , so that we haveδψ = O(ǫ g ). It is important to note that since we are working only to linear order in spatially varying quantities, the new parameter ǫ 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, ǫ = {ǫ t , ǫ H , ǫ k , ǫ 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 relatioñ Similar to what has been done in the previous section, by using the perturbative expansion in Eq. (3.7), 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δ and the leading nonzero modes corresponding to n = 2 can be obtained by using Eqs. (3.31) and (3.32) (3.36) We will present the equations for higher orders (n > 2) in Appendix B and limit ourselves to order n = 2 (and hence up to O(ǫ) 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 δψ s and Φ s can be obtained by setting ν = 0 in Eqs. (3.31) and (3.32) 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 as the effective Schrödinger equation, anḋ as the equation for the slow-mode Newtonian potential Φ s . We also have the effective Poisson equation, 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. 4, 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. 5, we confirm the validity of our EFT by numerical analysis. In Appendix E we study the case of a multicomponent universe and in Appendix D we take into account the self-interaction of the scalar field, which has been neglected so far.
In both Eq. (3.23) forψ s and Eq. (3.37) for δψ 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. (D.3) and (D.5) for the evolution of the slow modesψ s and δψ s in the presence of a tree-level self-interaction of the form V (φ) = λφ 4 /(4!) for the original (relativistic) scalar field φ. If we consider ψ s (t, x) =ψ s (t) + δψ s (t, x) and work to first order in fluctuations, we may combine Eqs. (D.3) and (D.5) to write We see immediately that the induced self-interaction is attractive (λ eff < 0) and that it arises gravitationally (being proportional to G ∼ 1/M Pl 2 ). 2 The emergence of a gravitationally induced self-coupling could be of interest in the context of self-interacting dark matter [8][9][10], 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 σ/m 1 cm 2 /g = 4.6 × 10 3 GeV −3 , where σ is the total scattering cross section. Given the form of λ 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 [20] 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) 2 The right-hand side of the effective Friedmann equation in Eq. (3.24) includes a contribution to ρ eff proportional to λ eff |ψs| 4 , scaling with λ eff andψ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 |ψs| 4 in Eq. (3.24) arises from backreaction of rapid oscillations on H solving the effective equations, we are able to find the full solution, including its oscillatory behavior. For example, the scalar field ψ can be constructed via ψ = ν ψ ν e iνmt , where ψ ν for ν = 0 are obtained order by order by our perturbative prescription and are given once the solution for ψ s is obtained from the EFT. In Sec. 5 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 [12] or the resonances in binary pulsars [11] as results of the oscillations of dark matter.

Effective fluid description
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 (3.23) and (3.24) 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. (3.24), as an effective energy density of the fluid: where the effective pressure turns out to be Notice that the effective pressure is of order ∼ O ǫ 2 ρ eff , so to O(ǫ)ρ eff the effective fluid behaves like pressureless matter, consistent with cold dark matter. 3 Nonetheless, the rapidly 3 The effective pressure in Eq. (4.3) should not be confused with the smeared pressure of the scalar field appearing in the energy-momentum tensor at the level of the background, p φ = 1 2φ 2 − 1 2 m 2 φ 2 , as is done for example in Ref. [11]. 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 H2 and H−2.
oscillating modes induce a small effective pressure, which yields an effective equation of state w eff = p eff /ρ eff of the form upon using Eq. (3.24) 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 ζ. 4 In Appendix C we present some equations governing an imperfect fluid with bulk viscosity. As can be seen from the results of Appendix C, at the background level, the pressure and bulk viscosity are degenerate and always appear in the form p−3Hζ. It is this combination that we denoted by p eff . However, when we incorporate fluctuations, the degeneracy will be broken. Note that the effective equation of motion for fluctuations, Eq. (3.37), has been obtained to O(ǫ). At this order, the background fluid is effectively pressureless. As a result, in the equations for fluctuations of a viscous fluid we set p eff = 0 or p = 3Hζ wherever these background quantities appear.
We may define the comoving overdensity of the fluid, δ eff , in terms of the right-hand side of the effective Poisson equation, Eq. (3.39): where ρ eff is the effective energy density of the background in Eq. (4.1), which, to O(ǫ), is ρ eff = m|ψ s | 2 . This yields Note that although δ eff is constructed from complex quantities, the combination in Eq. (4.6) 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. 3 with the fluid equations outlined in Appendix C. This results in 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.
Note that the above gauge-dependent variables are written in the Newtonian gauge. (See Appendix G for relevant gauge transformations and corresponding expressions in the timeaveraged comoving gauge.) Using Eq. (4.6) and Eqs. (3.37) and (3.38) we can derive a second-order differential equation for δ 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 C (see Eq. (C.14)). In this way, we obtain where we have defined the effective speed of sound by 11) and the effective coefficient of bulk viscosity by The bulk viscosity coefficient may be re-expressed in terms of another quantity with dimension of velocity [21]: which may be compared with c 2 eff in Eq. (4.11). Some remarks are in order regarding these results. The leading term in the effective speed of sound in Eq. (4.11) is well known [22,23]. 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. [22,23] but with a different coefficient. (Our term is larger than the result of Refs. [22,23] by a factor of 2.) We trace this discrepancy to the fact that in Refs. [22,23], the backreaction of the nonzero modes has been neglected, though they contribute at the same order. Note also that, despite the claims of Ref. [22] (and also recently in Ref. [24]), 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 ǫ k and holds for arbitrary momentum. Upon doing so, we confirm our coefficient, which is different from that obtained in Refs. [22,23]. (See Appendix F 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. (4.11), 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. (4.12), 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 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 [25]. However, our nonrelativistic system is 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. (4.12) is not inconsistent with the second law of thermodynamics. 5 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.

Comparison with numerical solutions
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. 3, are consistent with the smearing of the exact solutions by the method discussed in Sec. 3.1. 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. 3.1 we adopted the window function W (t) of Eq. (3.3) 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 [28]. 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)B(t/T ), where B(t) is the Blackman window function [28] and T is the total time of the simulation.  5 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 [26,27]. slow modes. This confirms that our EFT predictions match the numerical smearing with good accuracy. Fig. 2 also shows the hierarchy between nonzero modes among themselves and with the slow mode as predicted by the EFT. As we saw in Sec. 3, all the variables can be represented by a mode expansion. For instance for the background variableψ(t), to working order, we can writē ψ(t) =ψ s +ψ 2 e 2imt +ψ 4 e 4imt +ψ −2 e −2imt + . . . , while similar expressions hold for the other variables. 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. (5.1) order by order in our perturbation theory. In Figs. 3 and 4 we compare the exact solutions ofψ and δψ 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. 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. 6 That is why in Figs. 3 and 4 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. 5 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.

Summary and Outlook
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 for the background variables (up to O(ǫ 2 )), and for fluctuations (up to O(ǫ)), where, as a reminder, we parameterize the field as ψ(t, x) = ψ(t) + δψ(t, x). To leading order, these equations correspond to the already well-known corresponding initial conditions for the exact solution. Note, however, that the inverse procedure (which is more realistic) is also possible.
Schrödinger and Schrödinger-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: 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. [22,23]). 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. 5, 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 ∆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 ∆t < min(1/H, a/k); roughly speaking, this yields an O(1/ǫ) increase in efficiency. The full solution can then be constructed order by order using the mode decomposition outlined in Sec. 3 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 mentionedthat 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. [29]) 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.

A The Hamiltonian
In this appendix we will find the Hamiltonian for the evolution of ψ in a general curved background. Since the gravity sector is standard and we have not performed any redefinition for the metric we will assume the gravity sector as a fixed background with given time evolution. Standard treatment of the Hamiltonian for the metric can be found for example in Ref. [30]. Here, instead of deriving the Hamiltonian from the rather involved Lagrangian for the ψ field in Eq. (2.12), we first work out the Hamiltonian of the Lagrangian in Eq. (2.10) in terms of φ and the auxiliary field χ and then perform the field redefinition, via a canonical transformation, to derive a Hamiltonian for ψ and ψ * . Recall that after removingφ in favor of a new variable χ and integrating out a nondynamical field, we obtained the Lagrangian in Eq. (2.10), which we repeat here: Here, we have again written the potential term in its general form V (φ) for the scalar field which can contain both the mass term and the interaction. Before computing the Hamiltonian it is instructive to count the number of degrees of freedom in our theory. Since the metric is assumed to be fixed, in the original scalar field theory we had only two degrees of freedom, (φ, π φ ). In Eq. (A.1), we introduced a different Lagrangian which is completely equivalent, in the sense that it leads to the same equations of motion. Naively, the Lagrangian in Eq. (A.1) appears to involve two fields, φ and χ, so in principle it may have four degrees of freedom, (φ, χ, π φ , π χ ). However, the new theory is a constrained system [31,32]. To see this we compute the conjugate momenta We therefore see that we cannot solve forφ andχ in terms of phase space variables, that is, the Lagrangian is degenerate. The dynamics in phase space are constrained to a part of the phase space specified by primary constraints C i = 0 for i ∈ {1, 2}, where These two constraints reduce the number of degrees of freedom to 4 − 2 = 2, which is consistent with the original theory. The Hamiltonian density can then be computed from the standard procedure for a constrained system: and u 1 and u 2 are unspecified functions. Note that the dynamics is controlled by the total Hamiltonian H T , and we are not allowed to impose constraints on the total Hamiltonian before computating the relevant Poisson brackets. The constraints must be preserved in time; we therefore have the following set of equations describing the system: where Q 1 = φ, Q 2 = χ, P i are their corresponding conjugate momenta, and H T = d 3 x H T is the total Hamiltonian. The first term on the right-hand side of the equation forĊ i is due to the fact that the metric components (and hence constraints) can explicitly depend on time. Note that the Poisson brackets must be understood as operators in terms of functional derivatives, i.e.
From the set of equations in Eq. (A.8), we find u 1 =χ and u 2 =φ as well as the Klein-Gordon equation for the original field φ. An alternative approach which can lead to a simplified Hamiltonian formulation in a constrained system is to use the Dirac formalism [31]. To employ this method, we first note that the Poisson bracket of constraints is nonzero. Let us then define It can be shown that where we have defined N ≡ √ −gg 00 and N i ≡ √ −gg 0i . A prime over N i′ denotes that the argument is evaluated at x ′ , while ∂ ′ i denotes differentiation with respect to x ′ i . Note that the bracket of C 2 with itself is nonzero due to the presence of the spatial derivative terms.
With the help of the inverse matrix ∆ ij , we can construct Dirac brackets defined by (A.12) The inverse can be computed to be where by inverse we mean (A.14) Using the Dirac brackets allows us to use the simplified Hamiltonian, in which the constraint terms are removed. In other words, we simply use the Hamiltonian H in Eq. (A.7), with the price paid that the Poisson brackets are replaced with Dirac ones. Two relevant brackets are and One can then see thatQ i = {Q i , H} D with H = Hd 3 x gives the correct equations of motion. The next step is to express everything in terms of ψ and ψ * with the help of the field redefinition of Eq. (2.11). We must add to Eq. (2.11) a suitable transformation for their conjugate momenta so the transformation is guaranteed to be a canonical transformation. We use the standard procedure of constructing a generating function [33]. Let us consider a generating function of old momenta π φ and π χ and new variables ψ and ψ * of the form F 3 (π φ , π χ , ψ, ψ * ). Then we must have which, by using Eq. (2.11), we can find F 3 to be The Hamiltonian for ψ and ψ * is then where the first term is the Hamiltonian of Eq. (A.7) expressed in terms of ψ and ψ * , which takes form g 00 m 2 ψ 2 + g ij ∂ i ψ∂ j ψ + c.c. , .20) and the second term is a partial time derivative of the density F 3 , defined by F 3 = d 3 x F 3 , which takes the form (A.21) where we have imposed the constraints of Eqs. (A.4) and (A.5) to eliminate old momenta at the level of the Hamiltonian. Therefore, we have The final step is to find consistent Dirac brackets for ψ and ψ * . From the inverse transformation we have It is then simple to check that the relationψ = {ψ, H ψ } D (and its complex conjugate) leads to the correct dynamics for the system, where H ψ = d 3 x H ψ is the Hamiltonian with the Hamiltonian density given in Eq. (A.22).

B Higher-order terms for fluctuations
In Sec. 3.3, we have only shown nonzero modes up to order n = 2 for the fluctuations δψ(t, x) and Φ(t, x). For completeness, we present the general form to arbitrary order and also derive equations up to order n = 3. From Eq. (3.32) we can deduce for nonzero modes and similarly from Eq. (3.31) we get H sδ ψ s δ ν,2 δ n,2 + (1 − δ ν,2 )(H sδ ψ and In the above expressions, a prime over the summation means ℓ +  must be less than the upper limit. Most of the above terms are zero at leading order, but at higher orders more and more terms contribute. For instance it can be shown that for n = 3 we havẽ and (B.8) This procedure can be continued straightforwardly to obtain higher orders. These terms can then be used to get an effective equation for the slow modes of different variables as done in Sec. 3. After lengthy algebra we obtain the effective equation for δψ s , up to this order, as follows Likewise, for Φ s we havė The effective Poisson equation takes form (B.11)

C Viscous fluid
In this appendix we review the equations for the evolution for an imperfect fluid, which are discussed in various references; see, e.g., Refs. [25,[34][35][36]. When describing an imperfect fluid, one typically characterizes the deviation from a perfect fluid by ∆T µν : where p, ρ and u µ are the pressure, density, and 4-velocity of the fluid. Here we adopt the convention of Refs. [25,35] that u i is the velocity of energy transport. In general, for firstorder hydrodynamics, one may consider effects of shear and bulk viscosity as well as heat transfer. However, for the purpose of our effective fluid description, up to the working order of Sec. 4 we only need to consider the bulk viscosity to fully describe the system. At higher order, on the hand, we expect the effective shear viscosity and the effective heat transfer to appear (since these contributions are consistent with the symmetries of the low-energy effective theory). In this simplified case ∆T µν takes the form where ζ is the coefficient of bulk viscosity (which can be both time and position dependent the latter of which will be ignored here); and a semicolon denotes covariant differentiation.
Note that the form of ∆T µν suggests modification of the pressure p with bulk viscosity pressure Π bv of the form: p → p + Π bv where Π bv = −ζu κ ;κ is proportional to the divergence of the velocity of the fluid. For the background equations this causes a trivial modification in the continuity equation of the forṁ where we have defined p eff = p − 3Hζ. As a result, the effect of bulk viscosity at the background level is degenerate with the effect of pressure. In the language of our effective fluid description for the background evolution, the effective pressure in Eq. (4.3) must include the bulk viscosity term, i.e. it is the expression for p eff . On the other hand, for fluctuations of the fluid, the bulk viscosity has nontrivial effects in the equations of motion, and the aforementioned degeneracy beteween pressure and the bulk viscosity will be broken. The energy conservation for fluctuations results iṅ where δρ and δp are fluctuations of density and pressure, respectively, and δu is the fluctuation in the velocity potential. Note that for fluctuations of the bulk viscosity pressure we set δΠ bv = −ζδ(u κ ;κ ), i.e., we will not consider fluctuations of ζ. We will see that this is sufficient to fully describe the system to the working order. Other equations for the fluctuations read We must add to these equations a relation between pressure and the density fluctuations. This is usually done in the comoving gauge, defined by the condition δu c = 0, where we define δp c = c 2 s δρ c . (C.8) Here the subscript 'c' denotes that Eq. (C.8) holds in the comoving gauge, and c s is the sound speed [37]. By a gauge transformation, we can find a similar relation in a general gauge as follows: where we have used Eq. (C.3) for the background energy conservation. Defining the comoving overdensity by and using other equations (and their time derivatives) to remove all variables other than δ, we obtain the following second-order differential equation for δ: in which we have defined For the special case of p eff = 0 we have a simpler form This is the equation that is of our interest in Sec. 4 since there we are working to O(ǫ) for the field fluctuations, and to that order, p eff vanishes; its nonzero value only appears at O(ǫ 2 ) and higher, which will be ignored. However, when we consider self-interaction in Appendix D the effective pressure is nonzero at leading order and we must use Eq. (C.11).
Other fluid fluctuations can be written in terms of δ as Before concluding this section we report similar equations in comoving gauge, i.e. δu c = 0 (which will be used in Appendix G). We write the line element in this gauge by The equations are then found to take the following form: and, for the continuity and Euler equations we havė

D EFT for the self-interacting field
In this appendix we consider a self-interaction term in the potential, Such self-interactions, even with a very weak coupling λ, might play an important role in the context of cold dark matter scenarios [5,6,9,10]. We follow the same procedure as in Sec. 3 to construct a low-energy EFT. We perform the field redefinition introduced in Sec. 2, expand the new field ψ as an infinite series of different modes, and find the equation of motion for each mode. Then we identify appropriate small parameters, perform a perturbative expansion of the nonzero modes, and incorporate the backreaction from (rapidly oscillating) nonzero modes into the effective equations of motion for the slow modes.
To identify appropriate small parameters, we consider the equation of motion in terms of φ as discussed in the introduction. Since we are interested in the regime of oscillating solutions with dominant frequency m, we must have where we have expressed what we had in Eq. (1.6) in terms of ψ. This new parameter, together with ǫ t , ǫ H , and ǫ k , forms the set of small parameters for our problem. As discussed in Sec. 3, for fluctuations, we also consider the small parameter ǫ g , to indicate that we are working only to linear order in spatially varying fluctuations and make the power counting more straightforward. Here we present the results of our analysis up to O ǫ 2 for background and O(ǫ) for fluctuations. For background variables we have as the effective Schrödinger and Friedmann equations, respectively. Note that the term proportional to λ 2 on the first line in Eq. (D.3) has also been obtained in Ref. [6] We can find the effective fluid density and pressure by the similar procedure outlined in Sec. 4, which results in The effective comoving overdensity can be read from Eq. (D.6): We then obtain a second-order differential equation for δ eff from which one can read the effective sound speed and the effective viscosity. Comparing with Eq. (C.11) we obtain where the leading term proportional to λ is consistent with the results of Refs. [9,38]. And, finally, the bulk viscosity coefficient reads We are not making any assumption about the relative fraction of the energy densities. However, we expect that by decreasing the share of the scalar field, the backreaction effects of rapid oscillations will be suppressed. We, of course, still need to assume that H ≪ m, such that the scalar field is oscillating, so our formalism is applicable. We present the results considering terms up to order n = 3 (inclusive) and hence up to O ǫ 2 in the effective equations for the slow modes: where ρ s and p s are the slow-mode parts of the density and pressure of the additional fluid. The effective Friedmann equation reads The effective equation of energy conservation for ρ s , at this order, is the same as Eq. (E.1) and does not get corrections, i.e.ρ s + 3H s (ρ s + p s ) = 0.
We now proceed to consider fluctuations. The fluctuations for a perfect fluid are characterized by δρ, δp and δu corresponding to fluctuations in energy density, pressure, and the velocity potential. The relevant, exact equations are where we have defined and D is defined in Eq. (2.24). We also need the relation δp c =c 2 s δρ c , in comoving gauge δu c = 0, wherec s is the speed of sound for the new component (not to be confused with the sound speed deduced from the EFT for the scalar field), which can be written in Newtonian gauge via a gauge transformation. Again following the same set of steps as in Sec. 3, we obtain the effective equations for slow-mode variables δψ s , Φ s , as well as δρ s , δp s , and δu s . Here we only work up to order n = 2 (which implies that in the effective equations for the slow modes, we neglect terms at the order of ǫ 2 and higher). At this order, equations of energy and momentum conservation for the fluid remain unchanged, i.e. Eqs. (E.6) and (E.7) hold with all variables replaced by their slow-mode counterparts. The equation for δψ s at this order becomes and for Φ s with only a few changeṡ and As for the effective fluid description, it is rather tricky to derive the corresponding sound speed and viscosity in the case of a multicomponent universe, as the definition of the density contrast for the individual components is ambiguous (partly because we are dealing with two systems that are not well separated and indeed coupled through gravity). We will leave this analyses for a future work and only present here the effective energy density and the effective pressure for the background evolution: F Nonlocal field redefinition and the sound speed at arbitrary scale In Ref. [6] a nonlocal field redefinition is introduced, which yielded a dramatic simplification when deriving an EFT in Minkowski spacetime. In this appendix, we use the same field redefinition -still in Minkowski spacetime -to derive a sound speed for the fluctuations, which we expect to hold for arbitrary momentum. We show that this leads to a sound speed, consistent with Eq. (4.11) in low-momentum limit, thereby confirming our results, including the coefficient of the subdominant term. Further, this sound speed converges to c s ≃ 1 in the relativistic regime, as expected. We also add to the analysis of Ref. [6] by deriving the nonlocal operator (rather than simply postulating it), by requiring the resulting theory to have certain properties. We will see that these requirements do not fix the field redefinition uniquely, but they do suggest that the redefinition proposed in Ref. [6] is the simplest possibility. Taking some steps more generally to also include the case of an unperturbed FLRW background, we will see that a similar set of requirements fail in an FLRW universe to give sufficiently simple relations, which justifies our approach in this paper of returning to a local field redefinition.
We start from the Lagrangian in terms of φ and χ given in Eq. (2.10) and consider a more general field redefinition to introduce the ψ field. Note that φ and χ are both real fields, therefore a general field redefinition may take the following form: This ensures that we obtain the Schrödinger equation in the same limit. Finally, we want the ψ field to be a 3-scalar, and to transform under spatial rotations in the standard way, which restricts the spatial dependence of O To further restrict the forms of O 1 and O 2 , we consider the Lagrangian in terms of ψ and ψ * . Plugging the definitions of Eq. (F.1) into the Lagrangian of Eq. (2.10), assuming an unperturbed FLRW spacetime -that is, neglecting backreaction of the field on the spacetime geometry -and performing some integrations by parts we obtain where we have assumed that O 1 and O 2 can be expressed as infinite series in powers of ∇ 2 , so that spatial integrations by parts do not introduce minus signs. We have also defined Before considering an FLRW background, let us first consider the Minkowski case with a(t) = 1 and H(t) = 0. In this case we can find a solution for Γ 2 = 0 if we assume that Q is real (or pure imaginary), which results in two equations, This reduces the general field redefinition to the nonlocal field redefinition introduced in Ref. [6]. The Lagrangian then takes the form L = i 2 ψ ψ * −ψ * ψ − mψ * (P − 1)ψ for Minkowski , (F.14) which has an explicit U (1) symmetry. Note that the U (1) symmetry will be broken in the presence of a self-interaction, indicating the violation of particle-number conservation. Note also that there is no oscillatory factor in this free theory, which is a major simplification, since mode coupling does not occur. This simplification, however, does not extend to the case of a self-interacting theory.
Returning to an FLRW background, we need to solve a differential equation for Q by requiring the terms in the brackets of Eq. (F.9) to vanish, if we insist that Γ 2 = 0. However, we were unable to find a simple solution for that equation, and even requiring Γ 2 to be rather simple does not seem to simplify the derivation of an EFT. We therefore reverted to the local field redefinition, as introduced in Sec. 2 (see, however, Ref. [16] for a generalization of the nonlocal operator to the case of curved geometry).
The nonlocal field redefinition of Eqs. (F.1) and (F.13) in a Minkowski background enables us to derive the sound speed for density fluctuations, applicable for a wide range of momenta. To see this, first note that the equation of motion in this case is given by iψ = m(P − 1)ψ.
(F. 15) In the low-momentum limit, we may expand the nonlocal operator P to obtain the Schrödinger equation. However, Eq. (F.15) is exact (in Minkowski spacetime) and holds for arbitrary momentum. The Hamiltonian density in terms of the original field φ is given by which, in terms of ψ, yields the Hamiltonian (representing the total energy of the system): where in the second equality we have performed an integration by parts. This implies that the energy density of the system can be written as ρ = m|P 1/2 ψ| 2 . Inverting this relation suggests the relation ψ = P −1/2 ρ/m e iθ , where θ is an arbitrary, but real, function of time and space. We assume that our Minkowski spacetime is filled with a homogeneous and isotropic condensate of particles, and then study small fluctuations around it (note also that we are ignoring the gravitational effects). Using the equation of motion for ψ, Eq. (F.15), at the background level we havė ρ = 0,θ = 0. and going to Fourier space, we conclude that the sound speed must take the form As emphasized above, we expect this relation to hold for arbitrary momentum (in Minkowski spacetime). In the small-momentum limit we obtain

G Gauge transformation
In this appendix we investigate gauge transformations for fluctuations in our EFT. We mostly follow the conventions of Ref. [35]. The most general perturbed FLRW metric, considering only scalar perturbations, may be written For a general coordinate transformation of the form x µ → x µ + δx µ , the metric fluctuations transform as where we have defined δx 0 = δt and δx i = ∂ i δx. The scalar field φ also transforms as ∆δφ =φδt. However, we need a gauge transformation for δψ. From Eq. (2.6) it is evident that δψ is not a scalar under general coordinate transformations. We have where we have defined x ′µ = x µ + δx µ . As a precaution, note that in this appendix we use the prime to denote the variables in the new gauge, not to be confused with the rescaled time derivative in Sec. 3. In the new coordinates we can writė As a result, from Eq. (2.6) which defines ψ we can show that Note that since, at this stage, δt and δx are arbitrary functions it is not necessarily the case that the slow modes (with ν = 0) will dominate over the others. However, we will restrict ourselves to transformations for which such a hierarchy exists, as otherwise the gauge transformation may spoil the EFT construction for fluctuations. This is as a result of the fact that the field redefinition (2.11) breaks the general covariance and ψ is not a scalar under general coordinate transformation. Substituting into the gauge transformation equations of Eq. (G.2), we find ∆E ν = 2 δ t ν + iνmδt ν , ∆A ν = 2H α δt ν−α , ∆B ν = −2r α δx ν−α , (G As expected, this coincides with the definition δρ ′ eff = ρ eff δ eff when we re-express the gauge invariant comoving overdensity in Eq. (4.6) in terms of variables in the new gauge. By further investigation and using the Euler equation in comoving gauge given in Eq. (C.19), we can obtain the effective sound speed and coefficient of bulk viscosity which, not surprisingly, coincide with what we have already obtained in Newtonian gauge. This confirms the consistency of our results and, as a working example, shows the way that one can obtain the EFT and the fluid description in other gauges, making use of our results in the Newtonian gauge.