Beyond Schr\"{o}dinger-Poisson: Nonrelativistic Effective Field Theory for Scalar Dark Matter

Massive scalar fields provide excellent dark matter candidates, whose dynamics are often explored analytically and numerically using nonrelativistic Schr\"{o}dinger-Poisson (SP) equations in a cosmological context. In this paper, starting from the nonlinear and fully relativistic Klein-Gordon-Einstein (KGE) equations in an expanding universe, we provide a systematic framework for deriving the SP equations, as well as relativistic corrections to them, by integrating out `fast modes' and including nonlinear metric and matter contributions. We provide explicit equations for the leading-order relativistic corrections, which provide insight into deviations from the SP equations as the system approaches the relativistic regime. Upon including the leading-order corrections, our equations are applicable beyond the domain of validity of the SP system, and are simpler to use than the full KGE case in some contexts. As a concrete application, we calculate the mass-radius relationship of solitons in scalar dark matter and accurately capture the deviations of this relationship from the SP system towards the KGE one.


Introduction
Axions and axion-like particles are well motivated candidates for dark matter [1][2][3][4][5][6][7][8][9][10]. In cosmological and astrophysical contexts, the typical occupation number of the axion-like fields is large (for masses m a O[10]eV), allowing for a classical field description of the dynamics [11,12]. The classical field/wave dynamics are manifest in effects on the scale of the de Broglie wavelength of the particles. Since the mass of such axion-like fields is typically taken to be m a 10 −5 eV, the effects of wave dynamics can occur on macroscopic or even astrophysical scales, giving rise to the possibility of distinguishing such dark matter from other dark matter candidates. Novel dynamics such as the formation of solitons [13][14][15][16], interference patterns [14], transient vortices [17], and suppression of power below the de Broglie length-scale of the particles [8,18,19] is quite generic in axion-like dark matter. For recent reviews, see Refs. [20][21][22].
The axion field oscillates rapidly on the time-scale of order m −1 a , whereas its spatial variations are on length-scales L ∼ (vm a ) −1 , where v 1 is the typical velocity of the axion particles. Moreover m a /H 1 (where H is the Hubble parameter) within a few e-folds of expansion after the axion field starts oscillating. Together, these considerations indicate that a nonrelativistic description of the field, obtained by integrating out the rapid variations (in time), might be possible and fruitful for cosmological and astrophysical applications. Such an effective nonrelativistic theory would be

Klein-Gordon-Einstein
Nonrelativistic EFT for 'slow' modes = Schrödinger-Poisson + corrections integrate out 'fast' modes extremely useful (both analytically and computationally), since one would no longer need to resolve the rapid oscillations of the field.
In the present paper, we start from the relativistic Lagrangian of a classical, real-valued scalar field within general relativity. By systematically integrating out relativistic degrees of freedom we obtain an effective nonrelativistic description for the system. Our specific approach was first incorporated in Ref. [23] to obtain an effective field theory (EFT) in Minkowski spacetime for a self-interacting scalar field. It was then generalized for curved spacetimes in Ref. [24], and more specifically applied to the case of a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) universe, with the analysis restricted to linearized perturbations. *1 However, one important feature of dark matter is its ability to form dense, nonlinear structures due to gravitational instability in an expanding universe. The focus of this work is therefore to develop an EFT without any assumptions regarding the amplitude of the density perturbations of dark matter within an expanding universe. In this sense we obtain an EFT for axion dark matter in the nonlinear regime. Although metric perturbations are expected to remain small (at least in typical cosmological contexts [29,30]), we systematically go beyond linear order in the metric perturbations as well.
The leading-order result in our EFT is consistent with the Schrödinger-Poisson (SP) system in an expanding universe, which is widely used in the literature [8]. For example, the SP system has enabled long-time-scale simulations of nonlinear structure formation of axion-like fields [14,31,32]. It has also been used to understand the cosmological formation, gravitational clustering, and scattering of solitons with strong self-interactions in the early and contemporary universe [16]. Mirroring the late-universe simulations, purely gravitational growth of structure in the very early universe was pursued in Ref. [33] with the help of the SP system. The SP system was used for numerically exploring mergers and collisions of solitons with and without self-interactions in axion-like dark matter [34,35], along with their non-gravitational consequences [36,37]. The SP system was at the heart of exploring dynamical friction [38], relaxation [39], turbulence [40], halo substructure [41,42], kinetic nucleation of solitons [15,43], and the dynamics of transient vortices in fuzzy dark matter scenarios [17]. A number of existing numerical algorithms and codes are being used to explore the nonlinear dynamics of the SP system. (See, e.g., Refs. [40,44,45].) Given its importance and widespread use, it is critical to understand the domain of validity of the SP system as well as expected deviations from it. With our systematic expansion, which *1 Note that the terminology of effective field theory refers to two different approaches. One approach is bottom-up, in which all relevant operators that are consistent with the symmetries are included and then the coefficients are fixed by matching with experiments. This approach is incorporated for example in the EFT of inflation [25] and large-scale structure formation [26]. In contrast, our approach here is top-down, in which an EFT is obtained by taking the low-energy limit of a more complete theory. In this case, the coefficients appearing in the EFT are fixed by the parameters given in the more complete theory. This approach has been used for axion dark matter, for example in Refs. [23,27,28]. Useful comparisons of the different top-down results are also provided in the same papers. relies upon integrating out the dynamics on short time-scales, we go beyond the leading-order SP system of equations and capture quantitative deviations expected due to relativistic corrections. See Fig. 1. These deviations are expected to be small in most cosmological contexts in the late universe, when the fields are essentially nonrelativistic. Nevertheless, explicit expressions for the relativistic corrections to the Schrödinger-Poisson system can pinpoint which particular physical attribute of the system dominates the corrections. For example, one can hope to understand the relative importance of large gradients in the field, deviations from Newtonian gravity, and selfinteractions of the scalar field, and at what order in the relativistic corrections vector and tensor perturbations of the metric become relevant as one moves beyond the Schrödinger-Poisson system. This understanding, in turn, can clarify the domain of validity of the Schrödinger-Poisson system, and provide insights into the most efficient way of including relativistic corrections in different physical contexts (such as those discussed in the previous paragraph). The corrections can also point the way towards exploring deviations from general relativity or characterizing the type of field content making up the dark matter. Furthermore, they might point to symmetries in the problem that are lost or restored as we go beyond the SP system.
As an explicit application of our nonrelativistic EFT equations, we explore the mass-radius relation for dense solitons in the axion field. We demonstrate that our EFT better approximates the fully relativistic solution within the mildly relativistic regime compared to the SP equations alone. Although our EFT with leading-order relativistic corrections is more complicated than the Schrödinger-Poisson system, it is still easier to use numerically and analytically than the fully relativistic Klein-Gordon-Einstein (KGE) equations.
The rest of the paper is organized as following. In Section 2, we provide the fundamental equations and model for a real-valued scalar field within general relativity. We define the small parameters and state the relevant approximations which allow us to simplify the general system. In Section 3 we define and apply our procedure to systematically remove the relativistic modes from our system, and arrive at the nonrelativistic EFT. The main results of our paper are included in this section. In Section 4, we use our EFT equations to improve upon the mass-radius relationship of solitons obtained from the SP system. We end with discussion and conclusions in Section 5. A number of appendices provide technical details of calculations to which we allude but do not explicitly provide within the main text.

Setting the stage
In this paper we consider the nonlinear and inhomogeneous dynamics of a scalar field. Having in mind cosmological applications, we consider an expanding universe which contains a perfect fluid in addition to the scalar field; the additional fluid contributes to the homogeneous and isotropic background. However, we neglect perturbations of the perfect fluid and assume that gravitational collapse is only sourced by the scalar field. Our theory therefore takes the following form: GeV is the reduced Planck mass. The scalar field and the background fluid are described by a Lagrangian density (L ϕ ) and an energy-momentum tensor (T f ), respectively, where m is the mass of the scalar field, λ and κ are dimensionless coupling constants, and Λ is a large mass scale (compared to m). Note that we do not necessarily assume that the scalar field is the axion. However, for the specific case of an axion field, nonperturbative effects generate a periodic potential which is usually approximated by where f a is the axion decay constant and the last approximate equality holds for small field values (compared to f a ), which is a good approximation in the nonrelativistic limit. As a result, for the axion field, the parameters defined in Eq. (2.2) are not independent; namely, we have λ = −κ = −m 2 /f 2 a and Λ = f a . In a more general scalar effective field theory, the parameters λ and κ are independent, but the κ term is expected to be suppressed compared to the λ term for a sufficiently large cutoff Λ. In what follows, we will assume such a hierarchy but our approach for obtaining the nonrelativistic EFT can easily be extended to a situation with no hierarchy.
As discussed in the introduction, in the nonrelativistic regime the dynamics of the scalar field are dominated by oscillations with frequency almost equal to its mass m. Thus it is reasonable to rewrite the scalar field in terms of a new, complex variable ψ (the so-called "nonrelativistic field") by The remaining time or space dependence, encoded in ψ, is expected to be dominated by low-energy physics (i.e., lower than the mass scale), so that ψ is a slowly varying function of time and space (compared to the dominant frequency of the system given by m). However, due to the nonlinearities involved in the system, high-frequency oscillations appear in ψ with small amplitudes. The task of the following section is to integrate out such high-frequency modes and obtain an effective theory for the slowly varying mode. It should be noted that the field redefinition of Eq. (2.5) is not a one-to-one correspondence between the real field ϕ and the complex nonrelativistic field ψ. This fact, which is usually overlooked in the literature, was first addressed in Ref. [23] for Minkowski spacetime. In that work, the authors assumed a relation similar to Eq. (2.5) as a transformation in phase space with an accompanying redefinition for the conjugate momentum, which together make the transformation canonical and invertible. A nonlocal operator was also introduced in Ref. [23], which simplifies the derivation of the EFT in Minkowski spacetime. However, as discussed in Ref. [24], this strategy is not very helpful for the case of curved spacetimes. An alternative approach is to remove the redundancy in Eq. (2.5) by adding a constraint on the nonrelativistic field ψ(t, x). One convenient choice for the constraint is [24] e −imtψ + e imtψ * = 0 , (2.6) where the overdot denotes a time derivative. This constraint implies an equation of motion for ψ that is first order in time derivatives. By using Eq. (2.6) one can show thaṫ It is possible to apply the above field redefinitions at the level of the action and write the Lagrangian of Eq. (2.2) in terms of ψ and ψ * (see Ref. [24]). However, we are most interested in the equations of motion in terms of the nonrelativistic field. Applying Eqs. (2.5) and (2.7) to the Klein-Gordon equation yields where D is a differential operator defined by and V int is given in Eq. (2.2), which here is written in terms of ψ and ψ * using Eq. (2.5). Eq. (2.8) is a generalized Schrödinger equation in an arbitrary spacetime. Notice that Eq. (2.8) is exact and there exists a one-to-one map from the complex field ψ to the real field ϕ and its conjugate momentum. Also note the appearance of rapidly oscillating factors. A common approximation is to neglect such terms, under the assumption that they average to zero. Whereas this is true at leading order (which leads to the SP equations for the scalar dark matter, as we discuss below), these terms will play a crucial role in the derivation of our EFT, as shall be discussed in Sec. 3. Note that the oscillatory terms are present even in a free theory with V int = 0, which then leads to a tower of higher spatial-derivative terms in the free EFT; these terms can be thought of as the expansion of the relativistic energy in the large-mass limit [23].
To fully describe the system, the Schrödinger equation of Eq. (2.8) must be accompanied by the Einstein field equations as well as the energy-momentum conservation for the fluid, Let us emphasize again that in what follows, for simplicity, we will ignore perturbations of the background fluid.

FLRW metric with perturbations
Since we have in mind the application of our EFT to cosmology, we consider a perturbed expanding universe as a specific case of the general discussion in the previous section. As noted above, we intend to study the inhomogeneities in the scalar field ψ(t, x) nonlinearly. This implies that, contrary to the case for linear perturbation theory, the vector and tensor modes of the spacetime metric may play a nontrivial role in the dynamics of the scalar degrees of freedom. As a result, here we start from a general metric, including all forms of the metric perturbations, and then estimate the contribution of each type of mode to the dynamics of the scalar field. It is convenient to work with the ADM metric decomposition, which is given by 11) where N , N i and γ ij are the lapse function, shift vector, and the first fundamental form, respectively. We remove the gauge redundancy by the following choice of the metric components: where we have and we lower and raise the Latin indices with δ ij and δ ij . The background geometry is FLRW spacetime and a(t) is its corresponding scale factor. We have fixed the gauge by requiring the scalar mode of g 0i and the vector and some of the scalar modes of g ij to vanish. This choice of gauge can be retained to all orders of perturbations and is a natural generalization of the Newtonian gauge, which is particularly convenient for the system in the nonrelativistic regime. Note that we think of the above metric as perturbative in Φ, Ψ, σ i and h ij (while we treat the scalar field ψ nonperturbatively), which we will justify shortly. From Eqs. (2.11)-(2.13), one has which may be used in Eq. (2.8) as well as for the Einstein field equations.

Power counting
As we take the nonrelativistic limit of a relativistic theory, several small parameters/operators appear, which allows us to organize different terms that arise in the EFT. In this section we shall identify these small parameters. Furthermore, in the nonrelativistic limit and as a result of the source of gravitational perturbations being a scalar field, there exists a hierarchy among the amplitudes of the scalar, vector, and tensor modes of the perturbed spacetime metric. As we will see below, the scalar modes dominate and the amplitude of the vector mode is larger than that of the tensor modes. Let us start by looking at the Klein-Gordon equation (ignoring metric perturbations): (2.15) As stated above, in the nonrelativistic regime the mass term is the dominant contribution to the time evolution of the scalar field, and all other terms are suppressed. We have written the equation of motion in Eq. (2.15) in two different lines to emphasize this hierarchy. Demanding that the terms on the second line are smaller than those on the first line, we identify the following small parameters in the nonrelativistic limit: The first parameter quantifies the smallness of the expansion rate compared to the mass scale. In the opposite regime, when the Hubble scale is larger than m, the scalar field does not oscillate and cannot mimic dark matter behavior. The second parameter quantifies the smallness of the typical momentum of the dark matter "particles" compared to m, while the last two parameters specify the smallness of self-interactions. Note that if we assume that λ and κ are of the same order and Λ is a very large mass scale, then one can see that κ is much smaller than λ . In fact for the special case of the axion-like field we have κ = 2 λ . Although we do not restrict ourselves to the axion, we do assume a similar hierarchy between these two parameters, with κ ∼ O( 2 λ ). Next we study the hierarchy among the dynamical variables, which follows as a consequence of the system being nonrelativistic. For these estimates, one can use the Einstein field equations. However, most of the approximate relations can also be estimated by considering symmetries and other simple relationships. First we note that in order for the system to remain nonrelativistic even amid the gravitational dynamics, the gravitational potentials must remain small, The fact that Φ and Ψ are expected to be of the same order in g is discussed below. From the 00 component of the Einstein field equations one can obtain the Poisson-like equation for Ψ, leading to Further, the Poisson equation implies that there is another small parameter related to the amplitude of the scalar field, from which we find 2 ϕ ∼ x g . In addition, by using the trace-free part of the ij-component of the Einstein field equations we can see that that is, the difference between the two gravitational potentials is typically one order smaller than the gravitational potentials themselves. Further, from the 0i-component of the Einstein field equations we find 1 a where we have usedφ ∼ mϕ. Note that the relation between the vector mode and the scalar field (which acts as its source) could also be identified from symmetries and dimensional analysis. Finally, by using the ij-component of the Einstein field equations (or, again, by symmetries), we find In the following, instead of keep tracking of these small parameters individually, we collectively denote all of them by = { H , x , λ , g , ϕ } and work up to appropriate order in . This effectively means that we assume all small parameters are of the same order (except for κ , which is one order smaller). Depending on the application, it is expected that a hierarchy among the small parameters exists, in which case our EFT would be simplified. By using our approach, any higher-order term in the EFT can be derived systematically. However, in the main text of this paper we focus on the leading relativistic corrections to the equations of motion. In other words, we are primarily interested in terms which are one order smaller in than the SP system of equations, which constitute the leading relativistic terms. In Appendix C we proceed one order beyond the leading-order corrections for the specific case of spherically symmetric solitonic solutions.

Scalar and vector equations
By using our power-counting arguments, one can obtain a set of equations for the gravitating scalar dark matter in an expanding background at the requisite order in . Because our primary interest is the evolution of the scalar modes, we will be dealing with scalar equations of motion; hence the tensor modes h ij cannot appear by themselves, but will always enter the equations of motion accompanied by at least two spatial derivatives. This implies that the tensor modes would only appear at O 3 , according to Eq. (2.22). Similarly, the vector mode σ i appears with at least one spatial derivative which, upon using Eq. (2.21), implies that it appears at O 2 in the scalar equations of motion. Based on the above considerations, the generalized Schrödinger equation of Eq. (2.8), the Einstein field equations that reduce to the Poisson equation in the nonrelativistic limit, and the equation for the vector mode take the following approximate form: where we have defined Within the expressions in Eqs. (2.27)-(2.29), the fields ϕ andφ need to be replaced by ψ and ψ * using Eqs. (2.5) and (2.7) (see Appendix A). Note that we have replaced Ψ with Φ in several terms, because the difference between the two gravitational potentials is one order smaller than Φ and Ψ themselves. (See Appendix A for equations that are nonperturbative in Φ and Ψ, though still linear in the vector and tensor modes.) To avoid clutter we did not expand the exponential factors, but one needs to keep in mind that they are only relevant to appropriate order in their Taylor expansion. Notice that since -at this stage -different variables may contain highly oscillating contributions, it is not necessarily the case that the time derivative operator is small. Moreover, as a result of the assumed hierarchy between the self-interaction terms (i.e., κ ∼ 2 λ ), the κ term only appears in the Schrödinger equation at the current working order.
In general, the order of terms that are neglected must be compared with the leading-order terms. For example, in the Schrödinger equation of Eq. (2.23), the leading-order terms (such as ∇ 2 ψ/m 2 a 2 ) are already of O 2 , according to our power counting. This implies that we are neglecting some terms that are at least two orders higher in . It is thus evident that we have only kept the leading-order nontrivial corrections, which is indeed the case for all other equations. To go beyond that, one needs a more accurate set of equations. Toward that end, in Appendix A we consider the same set of equations but treat the scalar gravitational potentials nonperturbatively. Then we may apply our EFT beyond leading-order corrections for specific scenarios, such as spherically symmetric solitonic solutions (see Appendix C). In such cases, in which one may appeal to additional symmetries, deriving the next-higher-order corrections can be simplified. Throughout Sec. 3, however, we develop an EFT applicable to more general situations.
We shall see in Sec. 3 that one can obtain the SP system of equations from Eqs. (2.23) and (2.24) to leading order in the EFT, while corrections appear at the next order. As a final remark, note that at the background level the above set of equations reduce to where an overbar indicates that the quantity is evaluated at the spatially homogeneous background level (see Ref. [24]). In Eq. (2.30) we have also included the continuity equation for the additional perfect fluid, which is assumed to be spatially homogeneous. So far, by removing unnecessary terms, we have already taken the first step toward identifying the leading-order corrections in the nonrelativistic EFT. In principle, one can solve Eqs. (2.23)-(2.26) numerically to obtain the dynamics of the scalar field. However, such numerical computation is, in general, a difficult task due to the rapidly oscillating factors appearing in those equations. In Sec. 3 we shall remove such factors in a systematic way (instead of naively neglecting them) and obtain their corresponding corrections to the slowly varying variables. We will see that this procedure leads to nontrivial corrections at a given working order in , and hence cannot be neglected.

Effective field theory in the nonrelativistic limit
As stated in the previous section, we are interested in the slowly varying modes of dynamical variables. However, due to the appearance of oscillatory factors in the equations of motion, the dynamics of slowly varying quantities will be affected by rapidly oscillating terms. The situation is illustrated in Fig. 2, which shows the typical behavior of the variables in the frequency domain (i.e., the Fourier transform of time-dependent variables). As shown in Fig. 2, the zero mode (which translates to the slowly varying mode in the time domain) dominates, although modes with nonzero frequencies -close to integer multiples of the mass m -also exist in the spectrum, albeit with subdominant amplitudes.
To obtain a theory entirely in terms of slowly varying quantities, we must integrate out modes  associated with these rapid oscillations. This is a nontrivial task because the modes associated with rapid oscillations are sourced by the slowly varying mode, and they in turn backreact on the evolution of the slowly varying mode, affecting its dynamics.
Working in the time domain, we may apply a smearing operator to each variable in order to extract the slowly varying part of the variables. That is, we may take a time average of each variable with a suitable choice of window function [24]: in which W (t) is the window function and X s is the "slow mode" of the variable X. In Ref. [24] it has been shown that the top-hat window function in the frequency domain, which becomes sinc(t) in the time domain, is a suitable choice. Besides the slow mode X s , each variable also contains a tower of modes associated with rapid oscillations. Quite generally, one has where the coefficients X ν depend on both time and space. We define the "slow mode" as X s = X ν=0 , and refer to the modes associated with rapid oscillations as "nonzero modes," that is, modes X ν with ν = 0. According to the definition in Eq. (3.2), we have (X * ) ν = X * −ν , and therefore if X is real-valued then the modes X ν obey the constraint X * ν = X −ν . Note first that the expansion in Eq. (3.2) is exact, as a result of the appropriate choice of the window function. (See Ref. [24] for an outline of the proof.) Second, we emphasize that the coefficients X ν (even with ν = 0) are themselves slowly varying functions of time (compared to the frequency of the oscillations) since, as noted in the second expression in Eq. (3.2), the X ν may be represented by the smearing operator acting on X (weighted by an appropriate phase). In fact, as shall be made explicit below, the ν = 0 modes can be written in terms of the slow mode X s (since, to reiterate, they are sourced by the slow mode). We may therefore identify a new small operator, namely, the time-derivative operator acting on the slow modes, where X can be any of our variables after the field redefinition, including ψ, Φ, Ψ, σ i , h ij , and a; the subscript indicates that the slow mode (with ν = 0) is considered. *2 Thus, we can include the time-derivative operator within the set of small parameters/operators identified below Eq. (2.22), and our EFT will be an expansion in powers of = { t , H , x , λ , g , ϕ }. Note that these small parameters are not all independent; one may derive relations among them by using the equations of motion. Let us emphasize that t is defined as the operator that acts on modes X ν , rather than on the full fields. In the latter case -according to the definition in Eq. (3.2) -the time derivative would act on the oscillatory factors, which are not slowly varying, which is why, for example,Ψ in Eq. (2.24) should not be considered as contributing to O 2 in the power counting. We will be interested in the effective equations for the slow modes X s = X ν=0 ; therefore we will systematically remove nonzero modes X ν with ν = 0 from the theory. For a nonrelativistic system, all ν = 0 modes are suppressed by factors of the small parameters/operators denoted collectively by . Using power counting to estimate the size of each term that appears in the equations of motion for the nonzero modes, we may solve for them perturbatively in . To achieve this, we expand the nonzero modes as a power series in : The superscript (n) denotes the order of magnitude relative to the slow mode: This expansion allows us to solve for the ν = 0 modes perturbatively and substitute the solutions back into the equations for the slow mode, resulting in an EFT for the slow modes. This procedure has been done explicitly for an interacting theory in Minkowski spacetime in Ref. [23] and extended to the case of a linearly perturbed FLRW universe in Ref. [24]. Here we outline essential steps in the derivation, focusing mainly on the Schrödinger equation, and discuss additional details in Appendix B.
Applying the mode expansion of Eq. (3.2) to Eq. (2.23) yields the following equation for each mode: in which repeated indices are summed over. *4 This equation makes it evident that modes of different ν couple to each other; in particular, nonzero modes (ν = 0) affect the dynamics of the slow mode (ν = 0) and the nonzero modes are sourced by the slow mode. Similar results follow from Eqs. (2.24) for Φ and (2.25) for Ψ (see Appendix B). We can then solve for the nonzero modes perturbatively. At leading order we find

7)
*2 Because the functions Xν with ν = 0 are also slowly varying and can be expressed in terms of Xs, we also have t ∼ |Ẋν /mXν | -at least for ν not too large. *3 Note that the power series we use here is slightly different from the one introduced in Ref. [24]. X where ψ (1) ν is derived from Eq. (3.5) while Ψ (1) ν and H (1) ν can be obtained by solving equations after the mode decomposition of Eqs. (2.24) and (2.30), respectively. Note that δ ν,i is the Kronecker delta function, and the superscript (1) denotes that each term on the right hand side is suppressed by O( ) compared to ψ s , Ψ s or H s . An expression can also be derived for the ν = 0 modes of Φ, with the additional complication that the solution would be nonlocal. Fortunately, among the leading-order corrections (on which we focus in this section), nonzero modes of Φ do not contribute. Note also that the leading-order nonzero modes of the scale factor vanish.
We can now use these solutions to replace nonzero modes that appear in the equations for the slow modes. Furthermore, based on the power counting discussed in Sec. 2.2, we can neglect terms that are at higher order compared to the leading-order corrections. After significant algebraic simplification, for the Schrödinger equation we find Interestingly, notice that the other gravitational potential, Ψ s , decouples from Φ s and ψ s to this order. However, to close the system of equations, we must add one for the vector modes, which at this order is simply given by Eq. (2.26) with all variables replaced by their corresponding slow modes: Equations (3.8)-(3.10) are sufficient for obtaining the leading-order corrections to the Schrödinger-Poisson system of equations. However, the gravitational potential Ψ s might also be of interest for some purposes, such as lensing effects of compact objects. We therefore present its corresponding effective equation, which can be obtained by a similar procedure, starting from Eq. (2.25): (3.12) We do not need Eqs. (3.10) and (3.11) for the evolution of ψ s at this order. At the FLRW background level, it is easy to check that at leading order there are no corrections to the Friedmann equation or the continuity equation. That is, we have while the background Schrödinger equation receives corrections; these can be obtained by the replacement ψ s →ψ s in Eq. (3.8) and setting metric perturbations to zero. See Ref. [24], in which higher-order corrections are considered, leading to an interesting effective-fluid description for the scalar field with nontrivial pressure and viscosity. The effect of the new terms beyond the SP equations can be explored by numerical simulations. To test the EFT in a simple (but still nontrivial) example, in the next section we study solitonic solutions to the above equations under the assumption of spherical symmetry, neglecting the expansion of universe and the additional fluid component. Before we do that, we end this section with a few remarks regarding our results.
Note that from the bottom-up EFT point of view one can expect the appearance of all correction terms in Eq. (3.8) but with unknown coefficients. However, it is not the case that all terms consistent with the symmetries appear: for example, terms like Φ s λ|ψ s | 2 ψ s or H s ∇ 2 ψ s did not appear. This can be thought as the consequence of the original theory with which we started, namely general relativity with a scalar field minimally coupled to gravity. One may expect new terms to appear if one considers a modified theory of gravity, which may also change the coefficients of terms already identified in Eqs. (3.8)- (3.11). It would be interesting to explore such a possibility, which is beyond the scope of this paper. Furthermore, note that a term proportional to λ 2 has appeared in Eq. (3.8) with the same structure as the κ term. Therefore, we see that a single term in the original theory gives rise a tower of terms in the low-energy EFT as a result of integrating out high-energy modes. See Refs. [23,28] for more details.
Finally, we note that the SP equations (3.12) and (3.13) in an expanding universe possess a scaling symmetry [46] in the following sense. If ψ s (t, x) and Φ s (t, x) are a set of solutions, then ηψ s (ηt, √ ηx) and ηΦ s (ηt, √ ηx) are also solutions if we make the replacement λ → λ/η, ρ s → ρ s η 2 and H s → ηH s for any constant η. And the small parameters ( x , λ , g , ϕ ) get multiplied by η.
In general, this particular scaling symmetry of the solutions is lost as we include corrections to the SP system.

Approximate solitonic solutions
Equations (3.8)-(3.11) of our EFT, which include relativistic corrections to the SP system, can be incorporated in many different contexts and the solutions will take different forms. In this section we study one of the simplest solutions: spherically symmetric, stationary solutions of the form with Φ s = Φ s (r) in Eqs. (3.8) and (3.9) and Ψ s = Ψ s (r) in Eq. (3.11). As we will see below, this form corresponds to approximate solitonic solutions. Under spherical symmetry, the vector and tensor modes vanish identically. Moreover, by definition ψ s is a slowly varying function of time, so we must have µ/m ∼ t 1. In this section we ignore FLRW expansion, setting a(t) = 1, and also ignore contributions from the background fluid.
The specific form of the field in Eq. (4.1) resembles the wavefunction of stationary states in quantum mechanics. Although strictly speaking the field ψ s is not a wavefunction, its time evolution is governed by Eq. (3.8) which, at leading order, resembles the conventional Schrödinger equation. In the quantum-mechanical context, stationary solutions correspond to states with vanishing probability current and hence no actual time evolution. Similarly, the ansatz of Eq. (4.1) corresponds to a time-independent energy density. *5 By using the stationary, spherically symmetric form of ψ s in Eq. (4.1) and the time independence of Φ s in Eqs. (3.8) and (3.9), we find *5 One crucial difference between our system and conventional quantum mechanics is that our system is nonlinear (even without relativistic corrections), so that a superposition of solutions fails to be a solution. As a result, unlike in quantum mechanics, one cannot express the general time evolution in terms of a superposition of various stationary states. and where for simplicity we have set κ = 0. The above equations are the time-independent nonrelativistic EFT system of equations; in both Eq. (4.2) and (4.3), terms on the second line are smaller than those on the first by O( ) (while the first lines are already O 2 according to our power counting). Note the explicit appearance of the parameter µ in these equations. Since the stationary ansatz of Eq. (4.1) removes all time derivatives from Eqs. (3.8) and (3.9), we have used the leading-order equations to remove all spatial derivatives in subleading terms, which yields multiple terms proportional to µ in the final result. *6 Moreover, the ∇ 4 term in Eq. (3.8), which appears due to integrating out nonzero modes X ν with ν = 0, can also be removed by a similar manipulation. This removal of higher spatial derivatives makes the system more suitable for numerical calculations. We look for spatially localized, nodeless and regular solutions. That is, we demand that f (r) and Φ s (r) vanish fast enough at infinity; that f (0) = 0; and that the solutions are monotonic. Such solutions are expected to describe long-lived solitonic solutions.
We wish to compare solutions of our EFT, Eqs. To answer the first question we try to reconstruct the scalar field ϕ, or equivalently ψ, from the knowledge of ψ s and nonzero modes ψ ν . From Eq. (3.6) for the stationary solution we have where we have used the leading-order Schrödinger equation to simplify terms. Using ψ = ψ s + ν =0 ψ (1) ν e iνmt + O 3 and Eq. (2.5), we see that the relativistic scalar field will take the form where the higher-order terms also include higher multiples of the frequency ω ≡ m − µ. We have defined time-independent coefficients ϕ 1 = (1 + µ/2m)f and ϕ 3 = λf 3 /96m 3 at this working order. As a result, the specific form of Eq. (4.1) implies a periodic solution for the scalar field with a period 2π/ω, and we must look for this type of solution in the relativistic theory; keeping in mind that the true relativistic solutions also include radiating modes (leading to deviations from periodicity) [47][48][49], which are not captured here. Such field configurations correspond to *6 Some appropriate field redefinitions can remove µ completely from the equations, but change the asymptotic behavior of Φ to a nonzero constant. In this case Ψ can no longer be replaced by Φ to the working order, so that both gravitational potentials must be solved simultaneously. Such a system is easier to solve numerically, and the plots in this section take advantage of this procedure. We discuss this further in Appendix B and show the resulting simplified set of equations in Appendix C, in which we also go one order higher in the EFT expansion.
As for the second question, one important and reliable observable we have for solitonic solutions is their mass. We use the ADM definition of mass [76] (see also Ref. [77]), which is the Schwarzschild mass for an observer at infinity, where G N is Newton's gravitational constant, related to the reduced Planck mass by G N = 1/(8πM 2 Pl ). In the second equality we have used the Einstein field equations and the expression for the G 0 0 component of the Einstein tensor. Since the ADM mass is time independent, in the language of the mode expansion of Eq.
to working order. This means that, although the KGE and SP (plus suitable corrections) systems are two different theories, we can compare their solutions by demanding that they both yield the same solitonic mass. For a given mass, the radius is a good measure with which to compare SP (plus corrections) and KGE solutions, since it involves comparing the density of the two solutions and roughly shows how fast the fields decay with distance from the origin. The radius of the soliton can be defined as the distance enclosing most of the mass (the integrand of Eq. (4.6)). However, since the integrand of Eq. (4.6) depends on time, this definition results in an oscillating radius for the soliton. Whereas in most applications such oscillations may not be an issue, here we must treat the oscillations carefully. To avoid any confusion, we define the radius R 95 as the distance enclosing 95 percent of the mass after taking the time average of the mass density: *7 where the term in square brackets is given by the corresponding term in the integrand of Eq. (4.7). In practice, it is easier to fix the value of the field at the origin, rather than the mass, and find the corresponding solution in both theories. The mass and radius of the solution can then be obtained from Eq. (4.7). This can be done for our EFT system at lowest order corresponding to the usual SP equations with no corrections, as well as incorporating corrections at both O( ) and O( 2 ) beyond SP. That is, we consider SP, SP×[1 + O( )] and SP×[1 + O( ) + O( 2 ) ] equations respectively. The latter theory has been made explicit in Appendix C. *7 If we first find the 95 percent radius and then take its time average we get a result which is the same as above at leading order but has higher order corrections as well. Similarly, results for the KGE equations are obtained by following Ref. [74] and expanding the fields and gravitational potentials in terms of Fourier cosine series up to the frequency 4ω. We note that it is much simpler to do the shooting for solutions in the KGE system in spherical coordinates, whereas results in the effective theory are given in isotropic coordinates. The procedure for relating the two coordinate systems is discussed in Appendix D. We emphasize that, instead of comparing theories at some fixed central field amplitudes, we compare them using the mass-radius curve.
In Fig. 3 we compare the mass-radius relationship for solitons in the free theory (λ = 0) obtained with the KGE and with SP equations (with and without corrections). In the figure, "SP" refers to results from the lowest-order SP equations, neglecting all corrections, as given by the first lines of Eqs. (4.2) and (4.3). The EFT including O( ) corrections to the SP equations is given by Eqs. (4.2) and (4.3) (first and second lines), and the O 2 corrections are given in Appendix C. Compared with the SP equations, our effective equations with just the O( ) corrections improve the massradius relation significantly in the mildly relativistic regime, for R 95 10 m −1 . Note that at this radius, the mass calculated using the SP equations differs from that obtained from the relativistic KGE calculation by > 50%, whereas including O( ) corrections in our EFT leads to a discrepancy of < 10% from the relativistic KGE solutions. We can improve the results further by including O 2 corrections, which match the relativistic KGE calculation to within ∼ 1% around R 95 10 m −1 .
We have also confirmed the improvement in the mass-radius relationship obtained from our equations compared to the SP system in theories with repulsive (λ > 0) and attractive (λ < 0) selfinteractions in Fig. 4. To check the validity of our effective theory, we plot various small parameters for the theory with repulsive self-interactions in Fig. 5. The left panel shows the profile of in The results for much larger |λM 2 Pl /m 2 | (as would be the case for QCD and ultra-light axions) can also be obtained within our EFT. In this case the deviations from the SP system with attractive self-interactions still appear at mR 95 50. For the repulsive case, the deviation from the SP system becomes significant at larger and larger mR 95 as |λM 2 Pl /m 2 | increases.
terms of the radius, while the right panel shows the maximum value of in terms of the 95% radius. Our perturbative scheme fails when ∼ 1. Also note that the value of x is a measure of particle momentum, and we see that particles are indeed mildly relativistic when R 95 10 m −1 . We note that our reasons for choosing λ = ±12m 2 /M 2 Pl in Fig. 4 are that (i) by this choice all small parameters become the same order of magnitude in the mildly relativistic regime and (ii) it can make the comparison with the λ = 0 easier. With more canonical parameter choices, the value of |λM 2 Pl /m 2 | can be very large, for example, for QCD axions λM 2 where f a ∼ 10 11 GeV. A natural question is at what mR 95 does the mass-radius relation of SP equations start to deviate significantly from that of KGE equations when |λM 2 Pl /m 2 | 10? We have confirmed that few percent level differences from the SP system always start appearing at mR 95 50 as long as λM 2 Pl /m 2 −10. *8 For the repulsive case, and with λM 2 Pl /m 2 10, there is a minimum radius for the soliton (in the SP system). Around the minimum, the soliton is formed by a balance between gravity and self-interactions which implies that mR min ∝ λM 2 Pl /m 2 . The deviation between results from SP and KGE equations is large at mR min . As a result, relativistic corrections become important at larger mR 95 as we increase the value of λM 2 Pl /m 2 -which shifts the minimum (and therefore the whole curve) to the right in the mass-radius plot. The difference between attractive and repulsive cases is because the former has a balance between the gradient term and attractive self-interactions. For the repulsive case, also see Ref. [78].

Summary and conclusion
In this article we have provided a nonrelativistic effective field theory (EFT) for scalar dark matter in a cosmological context. We obtained the EFT by systematically integrating out rapidly oscil- lating modes from our system. We have treated matter nonperturbatively while keeping metric perturbations to the desired order in the EFT (including nonlinear metric perturbations). Our approach recovers the familiar Schrödinger-Poisson system at the lowest order in the EFT and provides explicit expressions for the corrections due to scalar field self-interactions and general-relativistic corrections, as identified in Eqs. (3.8) -(3.11).
Our framework allows for a systematic assessment of the relative importance of the various corrections to the Schrödinger-Poisson (SP) system in the mildly relativistic regime. For example, we showed that for leading-order corrections to the SP system, it is sufficient to include scalar modes at sub-leading order and vector modes at leading order, while neglecting tensor modes. As another example, if self-interactions are sufficiently strong it might be possible to ignore gravitational corrections [16]. Furthermore, a straightforward consequence of our results is that we see the loss of a particular scaling symmetry present in the SP system when we include leading-order corrections to the SP system (even with the quartic self-coupling set to zero). Our expressions also reveal that not all possible correction terms that might be expected from a bottom-up effective field theory consideration are present in general-relativistic systems. This result might be relevant for constraining departures from general relativity.
On the one hand our work provides a robust justification for using the Schrödinger-Poisson equations in many contexts relevant for structure formation. On the other hand, it provides guidance about how to improve upon the Schrödinger-Poisson system when relativistic corrections begin to become important without resorting to a fully relativistic analysis. Even when incorporating relativistic corrections, our EFT approach does not require tracking (or resolving) short time-scale dynamics of the system. The equations for the slowly varying, nonrelativistic modes acquire modifications upon integrating out the rapidly oscillating modes.
Implications of the sub-leading terms in the EFT can become important in the mildly relativistic regime. This could arise because one or all of the small parameters (defined in Section 2.2) are not sufficiently small. For example, we used our EFT equations to derive the mass-radius relationship of spherically symmetric dense solitons. We were able to obtain a significant improvement in the accuracy of this curve by using our EFT to go beyond the Schrödinger-Poisson system. Our EFT with corrections is simpler to use than the fully relativistic treatment in the mildly relativistic regime. However, as a limitation, our EFT cannot capture the decay of solitons by number-changing processes [27,49,79,80].
Whereas we only studied the stationary soliton in this paper, the time-dependent problem of mildly relativistic collapse is the next logical step to be studied with our EFT (see Ref. [81] for non-relativistic collapse). This study is possible within our EFT primarily because there is a hierarchy between the scale associated with nonrelativistic gravitational collapse, given by the Jeans scale (λ J ), and the Compton length-scale (λ C ) at which relativistic corrections become important: λ J λ C . Equivalently, one can see that many scalar-field oscillations with period similar to the Compton time (τ C ) fit within the Jeans time (τ J ) associated with matter collapse. *9 In certain regimes, the collapse can be protected from becoming too relativistic by gradients and anisotropies in the field.
Beyond spherical collapse, we can also employ our EFT to study transient vortices in fuzzy dark matter [17]. Our results might be useful for studying linear and nonlinear structure formation in axion-like fields in the early universe, and for exploring the gravitational and self-interaction driven clustering dynamics of a weakly coupled inflaton field after the end of inflation [16,33] (or, more generally, moduli fields [82]). Whereas SP systems are prominently used in late-universe applications, most studies of the early universe use a nonlinear Klein-Gordon equation in an expanding universe, with occasional inclusion of gravitational perturbations [83,84]. Our work can also be seen as a link connecting the nonrelativistic SP system with the partially relativistic system used in the early-universe context. Our EFT incorporates nonlinear aspects of gravity, and can thus provide a stepping stone towards exploring systems (such as scalar field collapse to black holes [85][86][87]), for which it may be necessary to deploy the full power of general relativity.
On the more formal side, it would be interesting to consider an effective fluid description of our corrections to the SP system using the Madelung-like transforms [24,[88][89][90]. Finally, repeating our work for more general matter content (higher-spin fields, such as vector dark matter [91,92]), or different relativistic theories of gravity [93], could change the nature of the corrections to the SP system (or even the SP system itself). This could further allow for constraining the nature of the underlying fields that constitute dark matter.

A Equations for a perturbed FLRW universe
In this Appendix we list the set of equations for a scalar field and an additional homogeneous and isotropic perfect fluid in a perturbed FLRW universe. Compared to what appeared in Sec. 2.3, we keep all scalar modes nonperturbatively while we only consider linear corrections to the equations due to the vector mode. We derive these equations with two limits in mind, for which such a framework is self-consistent: In Sec. 3 we proceed one order beyond the SP system, for which our equations below reduce to Eqs. (2.23)-(2.26), and in Sec. 4 and Appendix C we work under the assumption of spherical symmetry (in which case the vector mode would not be excited). For the same reason, we also neglect tensor modes, as a result of estimates discussed in Sec. 2.2. With these considerations in mind, from the Einstein field equations we obtain where we have defined with For the homogeneous and isotropic background we have Based on the power counting discussed in Sec. 2.2, one can check that the above relations reduce to Eqs. (2.23)-(2.30) as far as the leading-order corrections to the SP system are concerned.

B Details of the EFT derivation
In this Appendix we provide some details for obtaining the effective equations in the nonrelativistic limit. We derive the leading-order corrections here while in Appendix C we obtain next-to-leadingorder corrections for the spherically symmetric solitonic solutions.

B.1 EFT for the background
As a warm up, let us start with the equations for the spatially homogeneous background quantities which, as far as the leading-order corrections are concerned, are given by Eq. (2.30). Note that Eq. (2.30) is slightly simplified compared to Eq. (A.10) in that the self-interaction κ term is neglected except in the Schrödinger equation.
To obtain the effective equations for the slowly varying part of the variables, we start by applying the mode expansion described in Sec. 3 to the background equations and obtain equations for each mode ν as follows: where summation over repeated indices is understood, and we have (ψ n+1 ) ν = α 1 ,α 2 ,...,αnψ and we recall that overbars denote spatially homogeneous quantities,X →X(t). Similarly, for the other two equations we have We are interested in the system of equations for the slow modes (ν = 0). However, the slow modes are coupled to infinitely many nonzero modes (i.e., modes with ν = 0). Fortunately, due to the nonrelativistic nature of the system there is a hierarchy among different quantities, which enables us to solve for the nonzero modes perturbatively. The formal solution for the nonzero modes (ν = 0) is and Note that in the formal solutions above we have used the fact that all modes,ψ ν and ρ ν , are slowly varying and their time derivatives are suppressed by at least one factor of . In fact, all terms on the right-hand side are formally suppressed by at least one factor of compared to the slow mode of the variable appearing on the left-hand side. By using the power counting presented in Sec. 2.2 one can find solutions for the nonzero modes at leading order as follows: Note that at leading order, there are only a finite number of nonzero modes that remain nonvanishing. Next we use the above results to remove the nonzero modes that appear in the equations for the slow modes,ψ s = ψ , H s = H , and ρ s = ρ , which results in where dots stand for higher-order terms. It is instructive to have a comparison with the naive equations in which one neglects all the oscillating terms in Eq. (A.10) and obtains In comparison, there are several terms that appear in the systematic EFT expansion of the Schrödinger equation. The effect of these terms becomes more and more important as the system becomes more and more relativistic, that is, when small parameters, denoted by , grow. The two other equations do not get any explicit corrections at this order. This does not mean that the behavior of H s and ρ s are not affected at this order. Rather, due to the coupling toψ s they deviate from the naive solution. Furthermore, by considering next-to-leading-order corrections, one does find new terms appearing in the equations for the Hubble parameter and for the fluid (see Ref. [24]).
One can obtain an effective equation for the slow mode of the scale factor by a similar procedure. The easiest way to do this is to start from the definition of the Hubble parameter and employ the mode decomposition, which results iṅ It is then easy to see that a

B.2 Inhomogeneities
We now turn to the EFT for inhomogeneities. The procedure is similar to what we have done for the background theory. In contrast to what is done in Sec. 3, here we keep Φ and Ψ distinguished for future convenience of field redefinitions. This allows us to remove one parameter from the theory governing spherically symmetric solitons. This can be done by a redefinition of variables, as shall be discussed shortly.
We start by a formal mode decomposition of the equations for inhomogeneities, Eqs. (A.1)-(A.4), which results in Note that the mode decomposition of exponential functions must be understood as the mode decomposition of its Taylor expansion. For example, we have e 2Φ ν = δ ν,0 + 2Φ ν + 2 Φ 2 ν + . . . . The other equation we need for a closed system is the one for the vector mode. However, at the working order, vector modes are treated only at leading order and no correction is needed for them. As a result, the effective equation for the vector slow mode is simply Eq. (A.4) with all variables replaced by their slow modes.
As the next step, we need to solve for the nonzero modes perturbatively. For inhomogeneous variables, however, the relevant equations are not necessarily algebraic and can be differential equations including spatial derivatives. Whereas one can write a formal solution for a differential equation, for example by using inverse Laplacian operators, the solution might appear to be nonlocal. As was the case for the background equations, the nonzero modes appear as the corrections to the SP system. As a result, "integrating out" nonzero modes might in general result in a set of equations which are nonlocal in space. The appearance of nonlocal terms is a general feature of gauge theories (including general relativity), and is not necessarily a sign of the breakdown of any fundamental physical principle. Therefore, we expect such nonlocal contributions to show up at some level in a perturbative expansion. Fortunately, up to leading-order corrections to the SP system, these nonlocal solutions do not appear in our EFT, though they do at higher orders; see Appendix C.
To leading order, the resulting solutions for nonzero modes that will be needed are Next we must substitute these expressions into to the equations for the slow modes. From the Schrödinger equation we obtain where for the operators we havẽ

(B.24)
Note that in the expression forD s , terms that appear on the second line are one order smaller than those on the first line. Substituting the solutions for the nonzero modes to the above expression s ∇ 2 ψ * s + 2|ψ s | 2 ∇ 2 ψ s + ψ * s (∇ψ s ) 2 + · · · = 0 . (B.25) Likewise, the effective equations for the gravitational potential take the following form: All the equations obtained here will be reduced to our main results in Sec. 3 if one uses |Φ−Ψ| ∼ O 2 to replace Ψ with Φ in several places, as long as the difference does not contribute to the EFT at the current working order. However, keeping the two gravitational potentials distinct helps us to remove one parameter from the theory in the case of spherically symmetric solitons. In this case, as discussed in Sec. 4, we use the ansatz ψ = f (r)e iµt . After neglecting the expansion of the universe, we take the following steps to remove µ from all equations: We first remove higher-order spatial derivatives that appear in the EFT equations by using lower-order equations. The result will contain higher-order time derivatives instead. However, dealing with time derivatives is not a problem as a result of the simple ansatz we consider. Then, we use the following set of redefinitions to completely remove µ from the effective equations: After this redefinition we haveΦ − Ψ ∼ µ/m ∼ O( ), which justifies our attempt to keep the two gravitational potentials distinct even at this order. In Appendix C we present the results after this redefinition, where we also go one order higher in our perturbative expansion.

C Effective equations with O( 2 ) corrections
In this Appendix we provide the time-independent effective equations including O( 2 ) corrections for a spherically symmetric spacetime, which can be used to further justify the effective field theory through the mass-radius relation. As stated in Sec. 4, here we ignore the vector and tensor components of metric perturbations as well as the expansion of the universe. Also for simplicity, we set κ = 0. To achieve these results, we again remove µ from all equations using the procedure outlined at the end of Appendix B, including the redefinitions of Eq. (B.28). The effective equations that are shown here are the ones that we use for our numerical solutions. Since we have already discussed the steps required to obtain the EFT, we only present the final results here. The effective equation for the profile of the fieldf (r) is s λ 2f 5 768m 5 + λ 2f 3 (∇f ) 2 1024m 7 + λ 3f 7 73728m 8 + · · · = 0 , (C.1) and the equations for the gravitational potentials are given by
Here we provide the leading-order coordinate transformations, Φ 0 (r) = α 0 (r) + O( 2 ) , (D.8) where the first two equations are given by the coefficients of the cos(0ωt) and cos(2ωt) terms in Eq. (D.3) and the third equation is given by the coefficient of the cos(0ωt) term in Eq. (D.5). The derivation of higher-order transformations is straightforward.