Classical Nonrelativistic Effective Field Theory and the Role of Gravitational Interactions

Coherent oscillation of axions or axion-like particles may give rise to long-lived clumps, called axion stars, because of the attractive gravitational force or its self-interaction. Such a kind of configuration has been extensively studied in the context of oscillons without the effect of gravity, and its stability can be understood by an approximate conservation of particle number in a non-relativistic effective field theory (EFT). We extend this analysis to the case with gravity to discuss the longevity of axion stars and clarify the EFT expansion scheme in terms of gradient energy and Newton's constant. Our EFT is useful to calculate the axion star configuration and its classical lifetime without any ad hoc assumption. In addition, we derive a simple stability condition against small perturbations. Finally, we discuss the consistency of other non-relativistic effective field theories proposed in the literature.


Introduction
Light scalar particles arise in numerous theories of physics beyond the Standard Model. A prime example is the axion, typically realized as a pseudo Nambu-Goldstone boson (pNGB) of a broken U (1) symmetry at some high cosmological scale [1][2][3]. Such particles are also natural candidates for the particle nature of dark matter [4,5].
Axions are produced early in the universe, and naturally occupy states with very high occupation numbers, that coherently oscillate with dominant frequency ω ∼ m φ , where m φ is the axion mass. Such states are well-described by a classical field. Owing to the smallness of m φ , higher-order contributions to the oscillation frequency ω n ∼ n m φ are extremely rapid and are typically neglected at leading order. There now exists a standard procedure for deriving the low energy limit of the Klein-Gordon equation for a real scalar field: one expands the relativistic field φ in a nonrelativistic wavefunction ψ using the relation (1.1) and drops any term beyond leading order in rapidly oscillating factors e ±i m φ t . The resulting equation of motion, generically a nonlinear Schrödinger equation, is both classical and non-relativistic.
Localized quasi-stable solutions to the classical, non-relativistic equations of motion are known as oscillons [6][7][8] or boson stars. 1 Such objects can be supported by a balance of attractive and repulsive forces, sometimes including gravity. The original solutions for gravitationally bound (but otherwise non-interacting) boson stars were found by [11,12], though self-interactions have been included in recent years either in generic φ 4 scalar theory [13][14][15], and also in the specific case of the axion potential [16,17]. In the latter case then configurations are referred to as axion condensates, or more often, axion stars.
There are certain applications in which the non-relativistic limit may be insufficient. Over the last few years, a number of procedures have appeared for organizing corrections to the non-relativistic limit, which can generically be referred to as non-relativistic effective field theories (NREFTs). One such method was presented by some of the present authors (hereafter MTY) [8], in which the scalar field was decomposed into non-relativistic and rapidly oscillating parts, φ = φ NR +δφ. MTY presented a scheme for integrating out δφ perturbatively in the self-interaction coupling. This gave rise to corrections to the self-interaction couplings, as well as new terms in the equations of motion which were higher-order in spatial gradients as well as time derivatives. The method of MTY thereby accounted for all relevant relativistic corrections. Some relativistic corrections give imaginary terms in the EFT, which leads to the decay of oscillon states [7,[18][19][20]. The lifetime of an oscillon can be estimated from the EFT and the result is consistent with the classical numerical simulation.
A different method for constructing the NREFT of a real scalar field was performed by Braaten, Mohapatra, and Zhang (hereafter BMZ) [21], who matched relativistic scattering amplitudes at high energies to effective operators in the low energy theory, giving rise to modified self-interaction couplings. The original work of BMZ did not take into account higher-order gradients and time derivatives, which give important contributions at the same order in the EFT expansion as the self-interaction corrections they calculated. In a recent update [22], they have presented a more complete NREFT which includes these corrections as well.
Very recently, yet another method was presented by Namjoo, Guth, and Kaiser (hereafter NGK) [23]. In this work, NGK defined the relativistic field using a nonlocal operator which generalized Eq. (1.1), then expanded the wavefunction ψ and its conjugate ψ * as a tower of oscillating modes with frequencies ±n m φ . The equation of motion for the lowest order mode, taken to have frequency +m φ , was determined by integrating out the other modes. 2 It has been shown that, under appropriate field redefinitions, the methods of MTY, BMZ, and NGK give rise to matching elements in the low-energy S-matrix at ψ 6 order [23] and in the T-matrix at ψ 8 order [22].
A typical boson star is dilute and weakly bound, and as a result, relativistic corrections are extremely small. A theoretically ideal application in which relativistic corrections become important is in a dense axion star; in this configuration, the binding energy of the axions is large and the radius of the star is extremely small (close to the Compton wavelength of individual axions) [24]. This is a clear scenario in which the approximation of neglecting rapidly oscillating terms breaks down badly, and one must instead integrate out these modes. In [25], the leading nontrivial corrections to the wavefunctions of dense axion stars were calculated by generalizing the field operator for the scalar to include these fast-oscillating modes. The resulting self-interaction potential also agreed well with the effective interactions calculated in [8] at ψ 6 order. Importantly, in the dense axion star, gravity does not play an important role, as the self-interactions terms dominate [20,26], and so the application of a non-gravitational NREFT was appropriate. This is not the case in a dilute axion star, where gravity contributes at leading order. To analyze such states in an NREFT, gravity must be taken into account.
Thus, today there exist a plethora of methods of calculation in the NREFT of real scalar fields. It is relevant to understand the consistency of the results between methods, and also to point out any relevant advantages of one over another. In this note, we will clarify the general MTY procedure for calculating corrections in NREFT [8]. We will argue for a number of advantages to this method as compared to that of NGK or BMZ. One such advantage of the MTY method is that, as we will describe in this note, it is straightforward to include the effect of gravity, though in all cases it gives rise to certain theoretical complications. A thorough treatment of the gravitational interaction is especially important for understanding relativistic corrections to many classes of boson stars, where gravity contributes importantly at leading order.
As the axion field is a Hermitian (real) scalar, there is no symmetry which protects axion stars against number changing interactions. Nevertheless, there is an approximately conserved particle number that renders the configuration extremely long-lived [8,18]. The NREFT presented here renders the calculation of the classical decay rate for emission of high-energy particles extremely straightforward. As a concrete example, we apply this method to estimate the lifetime of dense axion stars.
One way to establish the consistency of an NREFT is to compare a numerical oscillon solution and prediction by EFT; this was established for the MTY method in [8]. Another way to check the consistency is to compare an exact solution and prediction by EFT. There are a few situations where an exact solution is known; one example is spatially homogeneous φ 4 theory, where the exact solution is given by an elliptic function. We will compare this exact solution with the result from the MTY, BMZ, and NGK effective field theory calculations.
The paper is organized as follows. In Section 2, we present the calculation of relativistic corrections in NREFT as an expansion in the small parameters of the theory, building on the MTY formalism of [8]. In Section 3, we describe how to couple this theory to gravity, and compare corrections coming from gravity to the other expansion parameters. We explain how to calculate the imaginary part of the Lagrangian, which gives rise to decay processes even in the classical regime, in Section 4. Finally, we compare various NREFT methods in Section 5, and we conclude in Section 6.
We use natural units throughout, where ħ h = c = 1.

NR EFT without gravity
In this section, we explain the method of calculation in our NREFT, building on our original work [8].
The advantage of this method in the study of oscillons is sixfold: the calculation is easy and straightforward by using Feynman diagrams; the reason for (approximate) stability is clear; it is easy to calculate the background configuration; it is straightforward to include the effect of gravity; the lifetime can be calculated from the imaginary part of the Lagrangian ; and the higher order corrections in the small parameter expansion become small. Each of these points will be illustrated in what follows. We start from the following Lagrangian for a relativistic quantum field theory of a real scalar field φ: where we assume a 2 symmetry for representative simplicity. The case without 2 symmetry has been discussed in Ref. [8], though in that work we did not include the effect of gravity. In this section, we refine the method and clarify the expansion scheme of our NREFT. In Section 3, we will explain how to include the effect of gravity.

Expansion scheme
We first decompose φ into slowly and fast oscillating parts: where the frequency ω ( m φ ) will be determined later. We sometimes expand the fast oscillating parts as where each ψ n is assumed to be a slowly varying field. 3 Note that this is just a Fourier decomposition and each ψ n is the mode that has a positive frequency. The goal is to derive the EFT by assuming ψ 1 is the dominant mode (in the non-relativistic regime) and integrating out the other highly oscillating modes ψ n>1 .
It is instructive to mention the implicit limitation of the mode expansion in Eq. (2.3) here. This expansion of the relativistic component δφ is useful when we compute the classical NREFT, defined by tree-level diagrams. Note however that this is not true if we would like to consider quantum effects because the energy of δφ must be continuous in loop diagrams. For this purpose we have to stop the decomposition of the scalar field at φ = (e −i ωt ψ 1 + h.c.) + δφ. See Ref. [8] for the definition of δφ in this case. The case of quantum decay for axion stars was also considered previously in [19].
Let us also clarify the role of ω here. In Sec. 2.2 we discuss how to get a stationary configuration of ψ 1 . There we will take ω so that it factors out the time dependence of ψ 1 . Nevertheless we sometimes keep both ∂ t ψ 1 and ω, since one may also take ω → m φ and keep ∂ t ψ 1 instead. The latter limit is useful for computing scatterings of free ψ 1 particles as done in Ref. [22]. We emphasize that both are equivalent after solving the equation of motion for ψ 1 . We clarify that we do not need to include a time-dependent term by hand in our EFT, as suggested in Ref. [22]. Such terms are already taken into account in Ref. [8]. We demonstrate the equivalence of each EFT in Sec. 5 by comparing with an exact solution in the limit of ∇ψ 1 → 0. Also, in Sec. 5.3.3, we explicitly show that our EFT contains the time derivative term pointed out in Ref. [22], if one takes ω → m φ and keeps ∂ t ψ 1 .
The EFT method is useful for the study of a clump of oscillating scalar field, like an oscillon or axion star. In the non-relativistic regime, there are three possible small quantities in the equation of motion: For each of them, we can assign pseudo parameters ε x ,V ,t = 1 like Then, we can systematically estimate relativistic modes by using an expansion in ε * . As a result, a term which is proportional to , c ∈ are exponents of the expansion parameters. This fact supports the validity of an expansion in ε * , given that δ * are small. There are certain oscillon solutions for which we find the following relation, In such a case, the resulting effective Lagrangian can be decomposed into powers of a single expansion parameter (ε V for example). In general, the EFT will contain higher derivative terms with respect to time, but we can use the equation of motion to remove them (see Refs. [27,28]). Let us illustrate how to obtain the NREFT here. Our strategy to get the NREFT is straightforward: keep the NR mode, ψ 1 , and integrate out the fast oscillating parts, δφ. The EFT action can be decomposed into kinetic, mass, and interaction parts S N R = S kin + S mass + S int . (2.11) First of all, substituting this decomposition into Eq. (2.1), one may readily find that the kinetic term plus the mass term become: where x = d 4 x is the spacetime volume integration. Note here that there are no cross terms among different n of ψ n . This is because the ψ n s are slowly oscillating by assumption, meaning that their time dependence is weak compared to the exponential prefactor. To make this point clear, let us consider the cross term of ψ † n and ψ n for n = n . There, we expect terms proportional to e i (n −n)ωt ψ † n ψ n . If ψ n and ψ n are time independent, then the integral of such a term over time is proportional to a delta function δ ((n − n )ω), which is nonzero only if n = n. A weak dependence on time effectively shifts this delta function, but in such a way that its argument is never zero, and so such terms integrate to zero in this case as well. Hence, all the terms must appear in a pair of ψ † n and ψ n . Then we move on to the interaction term. To make our discussion concrete, let us consider the following potential as an example: 14) The fact that ψ n s are slowly varying compared to e −i ωt is also essential for rewriting the interaction term by means of our decomposition. Taking this into account, one can decompose the interaction term into two parts: 1 Figure 1: The leading order diagram at ε 2 V for the λ 4 φ 4 interaction is shown. A thick line represents the relativistic mode, whose off-shell propagator is the inverse of the operator in parentheses in Eq. (2.13). with Here we have adopted the following notation ψ −n = ψ † n . The second term, δS , represents the interaction between the NR mode ψ 1 and the relativistic modes contained in δφ. In the second equality, we have used the mode expansion of δφ given in Eq. (2.3) which is useful to derive the NREFT in the classical case. Terms with uncompensated e ±i n ω t prefactors, e.g. e −4 i ω t ψ 4 1 , all integrate to zero, because of the slow time dependence of ψ n (see the discussion surrounding Eq. (2.13)). Now we are in a position to integrate out the relativistic mode δφ. The only diagram that contributes up to ε 2 V is depicted in Fig. 1, namely a 3-to-3 process mediated by one relativistic particle: Note that the propagator for the internal line is off-shell, and can be read off from the expression of Eq. (2.13). Also, note that there exists a cut contribution that represents the production of relativistic particles. We will come back to this point in Sec. 4 later, but for now just recognize that we have an imaginary part in our EFT.
To sum up, the resulting EFT for ψ 1 takes the form of where Γ [ψ 1 ] represents the imaginary part of the action, which is responsible for the lifetime of clumps; this will be discussed in Sec. 4. Note that each term of the resulting EFT must contain the same number of ψ 1 and ψ † 1 . This, again, owes to the slow oscillation of ψ 1 compared to e i ω t : a vertex like ψ n 1 (ψ † 1 ) n e −i ω(n −n)t integrates to zero unless n = n . 4 Consequently, the final result just depends on ψ 1 .
For example, if the original interaction term is given by Eq. (2.14), the effective potential is up to ε 2 V . As there are no corrections of the orders of ε V ε x or ε V ε t , the corrections involving ε t or ε x start from ε 2 V ε t or ε 2 V ε x : where we defined derivative operators as The corrections up to ε 3 x ,V ,t is calculated in appendix A. Before leaving this section, we would like to emphasize that the calculation of the effective potential in this method is extremely straightforward. One can derive the EFT corrections directly using Feynman diagrams, or by considering the ψ n >1 as perturbations on the equation of motion for ψ 1 ; both formulations give identical results. We illustrate the calculation procedure more fully in appendix B, where we compute the effective potential at (ε 3 ).

Stable configuration
If we neglect the imaginary part of the action, there is a global U(1) symmetry and the particle number of the field ψ 1 is conserved. In Ref. [29], Coleman showed that there is a stable configuration in such a complex scalar field theory with a U(1) symmetry under some conditions. The stable solution is known as a Q-ball. Therefore scalar clumps in a real scalar field theory, like oscillon and axion star, can be understood as a projection of a Q-ball by the decomposition of Eq. (2.2). That is, neglecting higher-order modes which give rise to decay, the NR oscillon field is which is the projection on the real axis of the U(1)-symmetric Q-ball [30].
Suppose that the number of particles in a system is given by Q , defined as The most stable and energetically favourable configuration of the scalar field ψ 1 can be calculated by minimizing the energy with its number fixed. Hence we should minimizẽ where ω is a Lagrange multiplier. This can be rewritten as where we set ω = ω . The first term is positive definite, so that it is minimized when ∂ t ψ 1 (x , t ) = 0. In other words, in the stable configuration we define ω so that ∂ t ψ 1 (t , x) = 0.
Assuming spherical symmetry, we can determine the spatial dependence of ψ 1 by the following equation of motion: where r = |x| is the radial coordinate. The boundary condition is ψ 1 = 0 for r → ∞ and ∂ r ψ 1 = 0 at r = 0. The condition for a solution to exist is as follows [29]: where Min ψ 1 [V eff / ψ 1 2 ] represents the minimum of V eff / ψ 1 2 as a function of ψ 1 . The solution of Eq. (2.30) is referred to as a Q-ball in the literature. Given ∂ t ψ 1 = 0, the frequency ω and the particle number Q are related by We can easily show that (see e.g. [31]) which means that the chemical potential inside of the Q-ball is equal to ω. Therefore, if ω < m φ , the Q-ball is stable compared with the free particle state. In addition, the stability condition against small perturbations is simply given by [32][33][34] ω Q We derive this condition, even when gravitational interactions are included, in appendix C. Remembering Eq. (2.2), we can understand the oscillon configuration as a projection of a Q-ball, whose stability is guaranteed by a particle-number conservation. Note that we can determine the configuration of clumps from Eq. (2.30), which can be solved numerically by using the shooting method. The condition for the existence of clumps is also clear through Eq. (2.31). These are some of the advantages of our EFT compared with others used in the literature in the context of bound states of scalars.

Effect of gravity on the scalar field
Here we discuss the inclusion of gravity on the oscillon. The starting point is to promote the spacetime derivatives in the equation of motion to covariant derivatives, Since we are interested in a spherically symmetric configuration, the metric around the oscillon can be written as where A(t , r ) and B (t , r ) are functions determined by the Einstein equation.
The energy-momentum tensor of the (relativistic) scalar field is given by The action of the scalar field is then given by The EFT action for ψ 1 takes the form of where V eff should be calculated as in the previous section. The gravitational corrections come through the factors of A and B in Eq. (3.6). We can write the energy-momentum tensor of the scalar field in the form of where ρ is the energy density, p r is the radial pressure, and p ⊥ is the tangential pressure [12], which in general are understood to be functions of the radial coordinate r as well as time t . The t t -and r r -components of Einstein's equations lead to [35] where G = M P −2 is the Newtonian gravitational constant and M P is the Planck mass. The θ θ -and ϕϕ-components of Einstein's equations are automatically satisfied when the above equations and the equation of motion of the scalar field are satisfied.

NR EFT with gravity corrections
Now we assume that the effect of gravity is small and provide an expansion scheme. At leading order, we write A(r ) = 1 − 2Ψ(r ) and B (r ) = 1 + 2Φ(r ), and define the small parameter We define the expansion in terms of Ψ because, as we will see below, at leading order Φ decouples from the equations and only Ψ is relevant. At higher order, we need to include Ψ and Φ. The parameter controlling the size of gravity corrections is Newton's gravitational constant, so we introduce a pseudo as we did with the other corrections. When we define ρ and p r by Eq. (3.7), they are determined by the scalar field as (3.14) From Eqs. (3.8) and (3.9), the gravitational potentials Φ and Ψ are determined as at the leading order in ε g . 5 At the leading order in ε g ,x ,t , the effective action and Lagrangian is simply given by where ψ ≡ e −i ω t ψ 1 and we impose a boundary condition Ψ(r → ∞) = 0. The equations of motion of this system are given by Note that the time derivative part equivalent to the one obtained from a free field theory, and the conserved charge Q is defined as usual by Eq. (2.27).
The way to obtain a Q -ball (or an axion star) configuration is as follows. For a given ω, we can define three dimensional action can be regarded as a three dimensional action for two coupled scalar fields ψ and Ψ which have canonical kinetic terms, and the potential is given by The existence of the bounce solution is ensured by the following three conditions (see [37] for more details): • The potential V S 3 is non-negative in an open neighborhood which includes the origin ψ, Ψ = 0; • The potential V S 3 becomes negative for some ψ, Ψ; • The effective potential V eff (|ψ|) satisfies lim |ψ|→∞ In general, the solution with the minimal action is ensured to be O (3) symmetric in space even when there are a lot of fields. The first condition is satisfied if and only if m φ − ω > 0. The second condition is automatically satisfied because in the large Ψ limit, the potential becomes negative. The third condition is satisfied as long as the potential is not so pathological. Note that these conditions do not ensure the validity of the perturbative analysis. For example, the obtained solution may have a large Ψ = δ g ∼ (1). We need to check the validity of perturbations separately. The relevant question here is whether the obtained solution is stable against small perturbations or not. Even with the gravitational interaction, stability is ensured provided that which is both a necessary and sufficient condition. We prove this condition in appendix C. Note that ω is still given by Eq. (2.33).

Discussion about expansion parameters
We may roughly estimate the small parameters of the EFT as 6 where R is the size of the oscillon, and we took the self-interaction potential from Eq. (2.14). The total oscillon mass can be estimated as M ∼ m 2 φ φ 2 R 3 , which implies The size of these parameters determines which EFT corrections become important.
We can build intuition about this by looking at stable bound states. Suppose λ 4 < 0 (attractive force). Then a bound state can be supported by a balance of the kinetic pressure against the gravitational and self-interaction terms in the equation of motion (this is the standard picture of a boson star). If δ V δ x ∼ δ g , then kinetic pressure supports the state against gravity. In this regime we have which is the standard mass-radius relation for a boson star with no self-interactions [11,12,38]. On the other hand, if δ g δ V ∼ δ x , then we have which is the mass-radius relation for non-gravitating boson stars [14,15]. These non-interacting (nongravitating) states are stable (unstable) under perturbations. If all three terms are of similar order δ x ∼ δ V ∼ δ g , we have a transition region which includes a critical (maximum) mass at corresponding to a radius of reproducing the standard results [14,17]. This raises another important point about the EFT. When we compute corrections at some order in ε x ,V ,t , this corresponds to corrections of size δ x ,V ,t at this order. If δ g is as large as the other small 6 We will ignore δ t for the purposes of this section. parameters in the problem (as in a stable boson star), then it is not consistent to neglect corrections at higher order in δ g . Another way to say this is that, if gravity contributes in the equation of motion, then post-Newtonian corrections proportional to Ψ 2 in the equations will be as important as those proportional to ∇ 4 φ, λ 2 4 , etc., which are the ones calculated in all formulations of scalar EFT we have been discussing [8,21,23,25]. The application of these EFTs beyond leading order, in a system where gravity contributes in an important way in the equations of motion, is not correct. Therefore, we need to compute corrections coming from the gravity effect ε g up to the same order with ε x if we consider the case where the clump forms due to the gravitational interaction. We leave a full analysis of higherorder gravitational corrections to future work.

Numerical results
Here we show a numerical result for an axion star, where the interaction potential is given by a relativistic periodic potential minus the mass term, In general, such a potential will have contributions at every even power of the field φ, but for simplicity we will truncate at φ 6 . This is the minimal axion potential for which a dense axion star solution exists [26]. The coefficients are where f a is an axion decay constant. In this section, we use the EFT contributions from the φ 6 potential, which we calculate in appendix A, and we take into account ε 2 V as well as ε g corrections, but neglect higher orders in ε x and ε t . In order to perform numerical simulations, it is convenient to rescale variables and fields such that they are dimensionless. We adopt the following rescaling: where we define (3.37) The equations of motion are given bȳ  where the prime denotes the derivative with respect tor . The boundary conditions are now given bȳ

Log[Q]
We numerically solve the equations of motion (3.38) and (3.39) with various initial conditions. Then we calculate the total charge (i.e., the total particle number Q ) and the radius R of the configuration. Note that we define the radius as the one at which 90% of the energy is enclosed. The resulting phase diagram is shown in Fig. 2, where we take g = 10 −4 and 10 −6 . In the plot, we show the rescaled quantitiesR andQ , defined by the rescaling of Eqs. (3.34) and (3.35) For a large radius, the gravity effect dominates and we obtain R ∝ 1/Q , which is consistent with Eq. (3.28). This branch is shown as thick lines in Fig. 2. AsR decreases and reaches to (1) (see Eq. (3.31)), the oscillon enters into the regime of d Q /d ω > 0 and becomes unstable. This is shown as dashed lines in the figure. There is another stable branch in the phase diagram, shown as thin solid lines, where the radius is of the order of (but is still larger than) m −1 φ . The solution in this regime has been well studied in the context of oscillons, where the gravity effect is negligible. In the context of axion clumps, it is sometimes called a dense axion star [24] or axiton [39]. A typical radius of an axion star can be estimated asR ∼ g /δ x from the dimensional analysis. For the dense axion star, δ x is smaller but not much smaller by many orders of magnitude than unity. Therefore we expectR ∼ g for the dense axion star, which is consistent with our numerical results.
If gravity is absent, the condition for the existence of axion star solutions is given by Eq. (2.31). We find that it cannot be satisfied unless m φ − ω is not so suppressed; near the right endpoint of Fig. 2, 7 the expansion parameters are not small and the non-relativistic expansion may not be a good approximation. It is not directly clear how next-to-leading order corrections would affect the endpoint of axion star solutions. In addition, the decay rate in this regime is not suppressed, as we will see in the next section. Our numerical results show that it becomes difficult to find a long-lived dense axion star solution if we take smaller g .
Our results are broadly consistent with the ones in previous literature, as in most works the authors ignore ε x and ε t corrections (as we did in this section) [24,41,42]. A leading-order analysis of the first nontrivial harmonic in a dense axion star was performed in [20], whereas all leading-order relativistic corrections were taken into account in [25]. Here we have taken into account the leading order relativistic corrections at (ε 2 V ) and our calculation can be applicable to relatively large ψ 1 , which is the case for the dense axion star branch with a largeQ . Of course, any leading-order analysis will break down for the very strongly-bound regime of dense axion stars, where m φ −ω is no longer small. In our analysis, stability against small perturbations is also confirmed in both dilute and dense branches, as seen in the lower panel of Fig. 2, which is in very good agreement with the stability analysis of [25]. † 1 · · · · · · ⇢ n Figure 3: The n to n process mediated by one relativistic particle is shown. The off-shell propagator for the internal line is determined by Eq. (2.13).

Classical decay rate
The number conservation is approximate and is violated by relativistic particle production, which leads to a finite lifetime for these localized clumps. Their lifetime can therefore be estimated from the imaginary part of the Lagrangian in the EFT, which is calculated by cutting relativistic propagators as done in Ref. [8]. An analysis of quantum decay for axion stars was performed by [19].
Let us first confirm that the imaginary part Γ of the action breaks the charge conservation. Taking a time derivative of Eq. (2.32), we readily geṫ which clearly shows that Γ breaks the approximate number conservation. Note that there is no time dependence on ψ 1 for the stable configuration if we neglect Γ . Now we are in a position to evaluate how Γ depends on ψ 1 . Since we focus on a classical NREFT, the relevant cutting diagrams which yield Γ must be tree-level. (We will comment on quantum loop diagrams shortly.) Motivated by this observation, we consider an n to n process mediated by one relativistic particle. See also Fig. 3. To be concrete, we consider a potential of the form in the original relativistic theory. After using the EFT expansion of Eq. (2.2), this term gives where δφ represents relativistic modes collectively. Classical decays of this type were first considered in the context of axion stars in [18].
The imaginary part appears when the relativistic propagator hits the pole of where we have taken ∂ t ψ 1 = 0 for the stable configuration, and we use the Feynman boundary condition for the relativistic propagator. 8 It is clear that, if ψ 1 is homogeneous, there is no pole and hence no imaginary part. Thus, the decay rate strongly depends on how ψ 1 localizes in space. This is why we first discuss the profile of ψ 1 by neglecting Γ . This procedure can be justified a posteriori if the decay from Γ is much slower than the typical formation time scale of classical lumps. In other words, if this is the case, one may assume that the scalar field can track the stable solution during the course of its slow decay. Note in addition that this is a fully relativistic phenomenon which cannot be captured by any finite expansion in terms of spatial derivatives of the propagator in Eq. (4.5).
To see the structure more clearly, we move to the Fourier space. Then the imaginary part of the action Γ gives rise to a decay rate Γ n→1 = d 3 x Γ . Then Eq. (4.2) from this particular diagram can be evaluated as which was computed in detail in [8] (though we have chosen a different normalization for the interaction term). We have factorized a typical amplitude of the profile as ψ 1 (x) = ψ 1 j (x), and defined Now it is obvious that the classical decay rate is non-zero if j n (p) has a non-zero value for p = n 2 ω 2 − m 2 φ , which originates from a spatial gradient of the localized lump.
We can calculateQ n→1 numerically by taking a Fourier transformation for the configuration of clump derived from Eq. (2.30). One may assume that the profile is approximated by the Gaussian: (4.12) 8 Strictly speaking, we would like to study the dynamics of ψ 1 and hence it is more appropriate to take the closed-timepath formalism. Then this term looks like where ± denotes the fields residing on the upper/lower contour in the closed-time-path formalism as mentioned in Ref. [8].
The equation of motion for ψ 1 can be derived by taking a derivative with respect to ψ 1 and then ψ (+) [8] for more details. Although the propagator for the relativistic mode here is the retarded one, one can show that the result coincides with Eq. (4.5) for vacuum of the relativistic mode where G ret = G Fyn holds. Here G ret/Fyn represents the retarded/Feynman propagator.
One can clearly see that the rate is suppressed for R ω R m φ 1 which is nothing but the nonrelativistic condition Eq. (2.4). Thus, our treatment is justified a posteriori. For a larger n , it is not easy to produce the relativistic particles because we need more energy to hit the pole. This is the case of the Gaussian profile, but it is possible to use other profiles to calculate the decay rate. Importantly, the decay rate depends on the tail of the momentum distribution, and so a compact function will give an incorrect result. It is not difficult to use the exact numerical solution of Eq. (2.30) for the oscillon profile; this was what was done in [18]. The result is parametrically similar to what we have estimated in Eq. (4.12) using the Gaussian ansatz. 9 One non-trivial example is the case when the profile is rather flat [40]. In such a case, j n (p ) can become zero at a certain radius. As a result, the configuration stays almost all time at the point where the lowest j n , becomes zero since the decay rate is highly suppressed at that point [8]. Note in addition that the decay rate is proportional to exp[− (δ −1 x )] and such effects can not be realized by any order of ε (or δ) expansions. In this sense, this decay process can be regarded as a kind of non perturbative one.
Here, let us estimate a typical time scale due to this classical decay by using a Gaussian profile. In that case we have δ x ≡ 1/(m φ R ) 2 1, where R is the radius parameter in the Gaussian profile. The total number of particles is given by 2Q 4π πR 3 m φ |ψ 1 | 2 . In addition, we define where we used the potential V n of Eq. (4.3). Note that δ x δ V n is expected from the equation of motion. Then, we have the following expression for the typical time scale for decay: If δ x is not so suppressed, this decay process may determine the fate of the configuration. Now we shall consider the axion star, which is numerically calculated in Sec. 3.4. In this case, the typical value of δ x can be estimated from the top panel Fig. 2 because Using δ x ∼ δ t = m φ − ω /m φ , one can also read it from the bottom panel of the figure. In the dense axion star branch, it is easy to see that δ t is of order 0.1−1 and hence by examining Eq. (4.14), we expect that these objects would decay before the present epoch. This was concluded also in the analysis of [18,20]. On the other hand, in the dilute axion branch, δ t is many orders of magnitude smaller. We can see in Figure 2 that the ratio of δ t in the dilute branch compared to that of the dense branch is smaller than g = (f a /M P ) 2 . Since we expect g ∼ 10 −6 for the GUT-scale axion or g ∼ 10 −14 for the (standard) QCD axion, the lifetime of dilute axion star is much longer than the present age of the universe because of the exponential factor in the decay rate. Therefore we conclude that dilute axion stars can survive until the present, if they formed in the early universe. The decay rate can be significantly larger when g is larger, as in models of Fuzzy Dark Matter [43].
Finally, we comment on quantum decay of oscillons. In this work, we have discussed the classical decay of oscillons via relativistic corrections, where the gradient term hits the pole of the propagator in the tree-level diagram. One may calculate loop diagrams to take into account quantum corrections in the EFT. Then there are some corrections to the imaginary term of the EFT from the cutting diagrams, which represent quantum decay of oscillon. If the decay rate is sufficiently small, it can be roughly estimated as the elementary decay rate of scalar field times the number of particles inside the oscillon. However, if the decay rate is sufficiently large, the statistics of daughter particles may become relevant. If the daughter particle is fermionic and obeys Fermi-Dirac statistics, its production rate from the surface of the oscillon has an upper bound by the Pauli exclusion principle [44,45]. If the decay rate is saturated, it is proportional to the number of degrees of freedom for particles that interact with the oscillon.
On the other hand, if the daughter particle is bosonic and obeys Bose-Einstein statistics, its production rate may be enhanced by the Bose-enhancement effect [7,46]. The latter effect may lead to an interesting observable effect for axion star [47]. The classical decay discussed in this paper could in principle have been affected by the Bose-enhancement effect. However, the effect is relevant only if the production rate is sufficiently larger than the escape rate of daughter particles. Therefore this enhancement effect is not important for the classical decay via the self-interaction.

Comparison with other works
In this section we will discuss other methods for calculating relativistic corrections in NREFT.

Method proposed by Namjoo, Guth, and Kaiser
In Ref. [23], NGK proposed an EFT by defining the scalar field with a nonlocal operator, giving rise to an equation of motion with only first-order time derivatives. In this case, the denominator of the propagator is linear in the energy, resulting in propagation which is only forward in time. In other words, in their EFT one integrates out part of the non-relativistic mode with frequency ∼ −m φ . On the other hand, our EFT contains the second-order time derivatives and integrates out only modes with high frequencies > m φ . Though the results agree at the level of T-and S-matrix elements [22,23], our method has the advantage of increased computational efficiency and smaller corrections beyond the leading order. In this section, we consider a polynomial potential with φ 4 and φ 6 order terms, although in [23] only φ 4 was considered. For the purposes of comparison, the relevant φ 6 corrections in the MTY method of the EFT are shown in appendix A.

Canonical transformation
The Hamiltonian of the relativistic theory (A.1) is given by where π is the canonical momentum. The canonical variables are φ and π. The Hamiltonian can be rewritten in terms of new canonical variables ψ and i ψ * by the following canonical transformation: where is defined by This should be included since otherwise the equation of motion of ψ depends on the spatial derivative of ψ * [23]. Note for future reference that the dimension of ψ in the normalization of NGK differs from that of Eq. (2.2) by a factor of m φ . This canonical transformation can be done by the following generating function F [φ, ψ, t ]: Hence the new Hamiltonian is given by where φ and π should be rewritten in terms of ψ and ψ * by using Eqs. (5.2) and (5.3). This is similar to the Hamiltonian in the non-relativistic field theory if we expand by using ∇ 2 /m 2 φ 1. The Lagrangian is therefore given by As a result, the equation of motion looks like a Schrödinger equation in quantum mechanics. Since we have just used the canonical transformation, the commutation relations for ψ and i ψ * are the same as those for φ and π. Up to now, the calculations are exact.

Expansion scheme
Following Ref. [23], we decompose ψ(t , x) as x)e i νm φ t , (5.11) where each of the ψ ν (x, t ) is assumed to be slowly varying. The modes with odd ν do not play any role in the following analysis because of the 2 symmetry in the original Lagrangian.
Here, we note that the zero point energy is shifted by an amount of m φ by the canonical transformation. Hence the ψ 0 represents the non-relativistic field whose kinetic energy is much smaller than its rest mass. Its rest energy m is removed from the Hamiltonian by the constant shift. So one may expect that the non-relativistic Lagrangian for ψ 0 can be constructed by integrating out ψ ν with ν = 0.
However, we should note that the modes ψ ν with ν ≥ 2 represents components of the original scalar field that have negative frequencies. In particular, the mode ψ 2 satisfies iψ 2 +m φ ψ 2 = −m φ ψ 2 , where a factor of m φ comes from the constant shift of the Hamiltonian. This means that the EFT of ψ 0 includes a contribution from integration of the mode ψ 2 with frequency of order −m φ . This leads to a large correction at the first nontrivial order.
To integrate out relativistic modes, we introduce a pseudo parameter ε and expand ψ ν (ν = 0) as We assume that ψ 0 is the zeroth order for ε. Then the equation of motion can be expanded by this small parameter. The result is given by at the first order of ε (= 1). Here and hereafter, we set = 1 for simplicity and focus only on the effective potential. The relativistic modes can be rewritten in terms of ψ 0 by usingψ ν m φ ψ ν : 14) where we define Substituting these into Eq. (5.13), we obtain the effective potential as We found that the effective potentialṼ eff is equivalent to Eq. (A.2) when we take into account the difference of the normalization of ψ. We expect that the form of the equation of motion can be factorized into a combination of ψ 0 and ψ 2 , and then the resulting effective potential is given by the same form asṼ eff and given by Eq. (A.2). This is consistent with the fact that ψ 2 , which is a fast oscillating mode from the standpoint of NGK, is a non-relativistic mode from the standpoint of MTY. 10

Method proposed by Braaten, Mohapatra, and Zhang
Here we explain the method proposed by BMZ in Ref. [21,22]. They first write down the effective Lagrangian (5.20) The effective coupling constants λ 2n are determined by the following procedure. First, they calculate the scattering amplitudes which have only non-relativistic particles in the external lines by using the original relativistic Lagrangian and the effective Lagrangian. Then compare the results and determine λ 2n so that the results are consistent. Note that the vertex factor in the EFT is given by λ 2n /2 n . In addition, there should be an additional factor 2 m for the diagram with m external legs in the EFT to compare the amplitude with the one calculated by the original theory. It is important that the contribution to the amplitude from the non-relativistic internal line in the EFT has a different meaning from the one calculated in the original theory. This is because in the non-relativistic EFT the denominator of the propagator is linear in the energy, resulting in propagation which is only forward in time. The resulting EFT is therefore equivalent to the one calculated in Ref. [23]. In fact, we have checked that the effective potential calculated in Ref. [21] (i.e. Eq. (5.18)) is consistent with the one calculated in Ref. [23] at least for the coefficients of λ 2 4 , λ 4 λ 6 , and λ 2 6 .

Consistency check
Here we check the consistency of EFTs: [8](MTY) and [21][22][23](BMZ-NGK). Since each EFT defines the dominant non-relativistic part in a different way, we should be careful about how to compare different EFTs. One way may be to compare some physical quantities predicted by EFTs which are truncated at the same ε n order. For example, EFT can predict the frequency-dependent total energy of oscillons E (ω), which is an important physical quantity. We expect the difference of the physical quantities to be (ε n+1 ).
Below, we will show such a consistency check of the effective theory at (ε 3 V ) in [8] (MTY) and [21][22][23] (BMZ-NGK). We consider a simple setup where there is only the quartic coupling and no spatial dependence for simplicity. The Lagrangian and equation of motion are then given by where the dot denotes derivative with respect to time. If the amplitude of φ is given by a × m φ with a being some dimensionless number, the frequency of oscillation ω can be written as where K (k 2 ) is elliptic integral of the first kind and given by Note that in the limit λ 4 → 0, we recover ω/m φ = 1.
For this system, we can use each EFT to estimate ω. The procedure is as follows. First, we fix the amplitude of the dominant mode (ψ 1 in MTY and ψ 0 in BMZ-NGK) to some constant value. Then, by using the equations of motion derived from each EFT, we can estimate ω. On the other hand, we can also estimate the amplitude a which is modified by fast oscillating modes. Then, we can use the frequency formula above Eq. (5.23) to estimate ω. If ω derived from an EFT and the one from the frequency formula are the same, the EFT can be regarded as a consistent one.
Below, we estimate ω up to third order in the λ 4 expansion by using formula Eq. (5.23) and EFT in MTY and BMZ-NGK. We will see that both EFTs can reproduce the correct result at least to third order in λ 4 . This means that if we make a prediction of amplitude dependent frequency ω(a ) by using both EFTs at ε 3 , the difference appears at ε 4 . This fact supports the consistency of both EFTs.

MTY
We fix the amplitude of the non-relativistic field to be |ψ 1 |/m φ = 1 in order to set a = 1 + (λ 4 ). This choice is just for representational simplicity. For any value of |ψ 1 |, we can do the same procedure and see consistency.
At the order of (λ 2 4 ), there exist the following fast oscillating modes: As a result, when we evaluate the expression at |ψ 1 /m φ | = 1, the total amplitude is shifted as Then, the frequency formula Eq. (5.23) gives On the other hand, our EFT gives Note here that the term proportional to (ω/m φ −1) comes from the interaction at ε 2 V ε t . This is exactly the same term pointed out in Ref. [22]. We will also clarify this point in Sec. 5.3.3. One may readily solve this equation at |ψ 1 /m φ | = 1 to give 30) which is consistent with the exact result given in Eq. (5.28).

BMZ-NGK
We fix |ψ 0 |/m 3/2 φ = 1/ 2 in order to set a = 1 + (λ). As before, this choice is just for representational simplicity. At (λ 2 4 ), there are ν = −4, −2, 2, 4, 6 modes in ψ (Note that some of these modes are mixed with the non-relativistic mode in MTY sense). Taking their contributions into account, the total amplitude becomes  34) which is consistent with the exact result given in Eq. (5.32). Finally let us comment on the size of corrections. Since MTY and BMZ-NGK are based on different expansion schemes, the sizes of fast oscillation modes are also different. The size of fast oscillation modes are estimated as a ratio of the strength of the source and the size of the propagator. In general, the size of the propagator in MTY is larger than that in BMZ-NGK. This is in part because the size of the propagator in MTY is proportional to n 2 − 1 with n ω being the frequency of fast oscillation modes while that of BMZ-NGK is given by n . As a result, corrections to the effective action in MTY becomes smaller than that in BMZ-NGK. For example, comparing the size of corrections in Eq. (5.30) and Eq. (5.34), we can see that at higher orders, the corrections become increasingly large in the method of NGK compared to MTY. This is one advantage supporting the use of the method of MTY, whose higher order corrections are relatively small.

On interaction with time derivative
In Ref. [22], they argued that there is an error in MTY [8], namely there is a missing term with a time derivative. However, as we have seen in Sec. 5.3, both EFTs give the same result, which suggests their claim is incorrect. In this section, we clarify that our EFT contains the required term. We explicitly show that the interaction with a time derivative arises if one would like to match the EFT with scatterings of NR free particles as done in Ref. [22].
First, we illustrate how to derive the interaction term proportional to (ω − m φ ) in the case of the stationary ψ 1 in Sec. 2.2. Let us start from Eq. (2.20). Thanks to ∂ t ψ 1 = 0, one finds The second term represents the correction at ε 2 V ε t which gives the term proportional to (ω − m φ ) in Eq. (5.29). Now let us confirm that the same term can be reproduced even if we take ω → m φ but keep ∂ t ψ 1 = 0 as done in Ref. [8]. This limit is useful if one would like to compute the scattering amplitude of free NR particles rather than to obtain the stationary solution. One easily gets The second term is nothing but the one pointed out in Ref. [22], which is clearly included in our EFT. We can also see that this expression is equivalent to Eq. (5.35) if one replaces ψ 1 (t , x) with e −i (ω−m φ )t ψ 1 (x).

Conclusion
We have provided a method to construct a classical NREFT in a scalar field theory including the effect of gravity. There are several advantages in our EFT: • the calculation is easy and straightforward by using Feynman diagrams; • the reason of (approximate) stability is clear; • it is easy to calculate the background configuration; • it is straightforward to include the effect of gravity; • the lifetime can be calculated from the imaginary part of the Lagrangian in the NREFT; • the higher order corrections become relatively small.
We have clarified an expansion scheme for gravitational corrections as well as relativistic corrections.
Since the gravitational potential has a radial dependence, we need to solve coupled differential equations to determine the oscillon profile.
Since the number of particles is approximately conserved in the NREFT, the oscillon can be understood as a projection of a Q-ball onto the real axis. Then its (quasi-)stability can be understood by the approximate conservation of the number of particles in the NREFT. If the gradient energy hits the pole of the propagator of the relativistic field, it gives imaginary terms in the Lagrangian of the EFT. This leads to the emission of relativistic particles and hence the decay of the oscillon. We have checked that the lifetime is exponentially suppressed with a large exponent which is inversely proportional to a small expansion parameter.
We have also discussed stability against small perturbations and see that the Q-ball analogy works well. As a result, we have found even with gravitational interaction, a necessary and sufficient condition for the stability against small perturbations is simply given by dQ /dω < 0.
As an example, we have considered an axion-like potential and found oscillon solutions taking into account the gravity effect. Our results are consistent with the ones in the literature, but can be extended to the regime where relativistic corrections become relevant. We have found that a longlived oscillon solution is absent in the limit where the gravity effect vanishes. The lifetime of axion stars can be estimated from these results and we have found that it is much shorter than the present age of the universe for the dense axion star but is much longer for the dilute axion star (in accordance with the previous results of [18]). We have concluded that dilute axion stars survive until the present day for realistic parameters of axion.
We also discuss the consistency of our EFT with the ones proposed by BMZ and NGK in Refs. [21][22][23]. As they have shown in those papers, we have checked that all of the EFTs are consistent, though depending on the application, different approaches have different advantages. In particular, one of the relatively low-frequency modes is integrated out in the EFTs of BMZ and NGK, but is not integrated out in our EFT. Since the low-frequency mode gives a large correction to the NREFT, our method obtains a smaller correction than that of BMZ and NGK. Even if we do not integrate out that relatively low-frequency mode, we can easily calculate the profile of scalar field configuration by solving the equation of motion numerically.
We have pointed out that, in a system where gravity contributes at leading order in the equation of motion, it is not consistent to account for relativistic corrections without taking into account corrections to the gravitational interaction as well. This is the situation in a typical dilute boson star, where gravity is one of the dominant forces determining the bound state configuration. The EFT we have presented here, as well as that of NGK [23] and BMZ [22], does not calculate corrections to the gravitational interaction. Interestingly, these corrections can give rise to new interactions and potentially new decay diagrams mediated by gravity and self-interactions. We leave a full analysis of this topic for a future work.

Diagrammatic derivation.
First we list all the interaction vertices after the decomposition given in Eq. (2.2). The interactions between the NR mode and relativistic mode are given by † and their conjugates. There also exists a self interaction for the relativistic mode, † Then, to get the NREFT, all one has to do is to integrate out the relativistic mode δφ. Throughout this paper, we focus on the classical NREFT, and hence all the diagrams that will be integrated out must be tree diagrams. The leading order term in the coupling expansion is obtained from † where we dropped terms of higher order than ε 3 and also the imaginary part in the second line. The next to leading order term involves three vertices as depicted below † where again we dropped terms of higher order than ε 3 and also the imaginary part. To sum up, the effective potential up to ε 3 is given by Derivation from EoM. Here we present another way to derive the NREFT which is essentially equivalent but still useful. Starting from the equation of motion (EoM), we solve δφ order by order in terms of ε V . In this case, the original EOM is given by We decompose φ as and assume ψ 1 dominates φ. Then, at ε V order, we have and the resulting solution is c.] into EOM, we obtain the EOM at ε 2 V order. In this order, only ψ (ε 2 V ) 3 contributes to effective action at (ε 3 ). We have 3 ) + h.c.] into EOM, we obtain the effective EOM for ψ ≡ e −i ωt ψ 1 Then, we can obtain the effective action at (ε 3 ) order:

C The stability condition against small perturbations
Here, we derive the stability condition against small perturbations when gravitational interactions are included. We fundamentally follow stability discussions done in [33,34]. The Lagrangian of the system can be written as 11 where K denotes a kinetic term and Ψ is a gravitational potential. The dot denotes derivative with respect to time. Suppose that the Q-ball solution is given by e −i ωt σ(r )/ 2 where σ(r ) is a real function. S 3 (ω) is stationarized by the Q-ball solution σ(r )/ 2. In second order in fluctuations, we can find eigenvectors and eigenvalues: We normalize them as follows Note that real and imaginary parts of δψ are separated. We decompose functional space into real and imaginary directions as δψ = δψ R + i δψ I . We know that S 3 (ω) has one negative mode in real direction. Below, we assume S 3 (ω) has only one negative mode, which is ensured if the Q-ball solution is a solution of the equation of motion with the minimal action S 3 . In addition, there are four zero modes for S 3 : three of them are related to spatial shift which we denote δψ S ,k (k = 1, 2, 3) and one of them corresponds to a U (1) shift which we denote δψ θ . δψ S ,k are a real function and δψ θ is a purely imaginary function. We remove such zero modes from δψ space and treat them in a different way.
Since the real direction has one negative mode, the stability is non trivial. We concentrate on eigenvectors with a real δψ a . As we will see, the negative mode can be shifted to a positive mode due to the mixing with the zero mode. We consider fluctuations in the real direction and a zero mode 12 : where θ (t ) denotes some real function, which is a zero mode in S 3 . What is non trivial here is the mixing between the zero mode θ (t ) and the real direction in K . We can solve the equation of motion forθ and the resulting kinetic term becomes We expand the vector (σ(r ), 0) t in terms of eigenvectors: Here δψ a are real functions in eigenmodes. We define matrix B and C as For c a (t ), the Lagrangian can be written as 12 δψ S ,k do not contribute to this analysis (see [33,34]).
Since B a b is a positive matrix, the minimal eigenvalue of the following matrix M a b ≡ λ a δ a b + 4ω 2 I σ t σ, (C. 18) determines the stability. If the minimal eigenvalue is positive, we expect that Q-ball configuration is stable against small perturbations. On the contrary, if the minimal eigenvalue is negative, Q-ball configuration becomes unstable agains small perturbations. We denote the eigenvalues of M by {Λ a }, and require that (Λ 1 < Λ 2 < ...) and also (λ 1 < λ 2 < ...) are different. 13  These relations ensure that there is only one Λ a between λ a and λ a +1 . Note that λ a has only one negative mode which we denote λ 1 . The condition Λ 1 > 0 is now equivalent to As is in [33,34], we can connect G (0) and d Q /d ω and show that the stability condition is given by Below, we derive G (0) = ω Q dQ dω . First, we note that σ(ω) is a solution to the equation of motion derived by S 3 (ω). Taking the derivative of the equations of motion for ψ and Ψ with respect to ω, we have Multiplying a σ a on 1 λ a S 3 (δψ a δΨ a ) t = (δψ a δΨ a ) t , we have Since δψ a and ∂ σ/∂ ω do not contain the zero mode of S 3 , we can multiply the inverse matrix of S 3 in (C.24) and (C.25). Then, we have 13 One may worry about the degeneracy of eigenmodes. However, if we deform the potential infinitesimally, such a degeneracy is broken in general. The physical results are expected to be unaffected by such an infinitesimal deformation. Therefore, we conclude that stability is ensured when Eq. (C.23) is satisfied.