Axion stars in the infrared limit

Following Ruffini and Bonazzola, we use a quantized boson field to describe condensates of axions forming compact objects. Without substantial modifications, the method can only be applied to axions with decay constant, fa, satisfying δ = (fa/MP)2 ≪ 1, where MP is the Planck mass. Similarly, the applicability of the Ruffini-Bonazzola method to axion stars also requires that the relative binding energy of axions satisfies Δ=1−Ea/ma2≪1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \varDelta =\sqrt{1-{\left({E}_a/{m}_a\right)}^2}\ll 1 $$\end{document}, where Ea and ma are the energy and mass of the axion. The simultaneous expansion of the equations of motion in δ and Δ leads to a simplified set of equations, depending only on the parameter, λ=δ/Δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \lambda =\sqrt{\delta }/\varDelta $$\end{document} in leading order of the expansions. Keeping leading order in Δ is equivalent to the infrared limit, in which only relevant and marginal terms contribute to the equations of motion. The number of axions in the star is uniquely determined by λ. Numerical solutions are found in a wide range of λ. At small λ the mass and radius of the axion star rise linearly with λ. While at larger λ the radius of the star continues to rise, the mass of the star, M , attains a maximum at λmax ≃ 0.58. All stars are unstable for λ > λmax. We discuss the relationship of our results to current observational constraints on dark matter and the phenomenology of Fast Radio Bursts.


Introduction
Scalar fields, which give rise to spin-zero quanta satisfying Bose statistics, are natural in quantum field theories of phenomenological importance. The recently discovered Higgs field, which generates masses for all the other particles, is the prime example. The axion, a yet to be discovered pseudoscalar, postulated to solve the strong-CP problem endemic to QCD, is another well motivated spin-zero particle. It is also generic in string theory for 4-dimensional axion-like degrees of freedom to arise as Kaluza-Klein zero modes from the compactification of anti-symmetric tensor fields defined in 10 space-time dimensions. Quintessence, the almost massless scalar field invoked to drive the late-time acceleration of the universe, is another hypothetical scalar degree of freedom. It would be quite interesting to determine if such spin-zero degrees of freedom could form stable compact structures due to their self-gravitation.
Starting with the seminal works of Kaup [1] and Ruffini and Bonazzola [2], the study of gravitationally bound bosonic degrees of freedom has received a great deal of attention. It has been seen that complex scalar field configurations in the presence of gravity can form stable compact objects, termed boson stars [3][4][5][6][7][8][9][10]. Kaup in his original work solved the Klein-Gordon free field equation in an asymptotically flat spherically symmetric spacetime background and found localized spherically symmetric field configurations which are energy eigenstates. Adopting a slightly different approach, Ruffini and Bonazzola used the expectation value of the energy-momentum tensor operator of a second quantized free hermitian scalar field system, evaluated in an N -particle state where all the particles occupy the lowest energy state, to be the source of gravitational interactions. Here the basis states for second quantization were the wave functions of the Klein-Gordon equation. This approach yielded an equation of motion for the wave function of the N -particle state,

JHEP03(2015)080
which represented a condensate of N bosons in a single state. They found stable localized wave functions for the condensate, indicating the interesting possibility that lumps of selfgravitating scalar field configurations could exist in nature.
In a recent set of publications, [11][12][13][14] analyzed the self-gravitating field theoretical system consisting of a real scalar field with an interaction potential of the form f a 2 m a 2 (1 − cos( φ fa )). They followed the method pioneered by Ruffini and Bonazzola. Such a potential represents the interactions of the axion field at low energies, as derived using the dilute instanton gas approximation, where f a represents the axion decay constant, and m a the axion mass [15][16][17]. In this paper we analyze the same system but present an analytic expansion, which simplifies the equations of motion. We also consider a different range of input parameters. This leads to interesting physical consequences. Without substantial modification, our expansion method can only be applied to axions with f a M P , where M P = 1/ √ 8πG is the Planck mass. In fact, as we will point out later, the field theoretic method of Ruffini and Bonazzola [2] can only be applied in this regime when the relative binding energy of axions, defined by ∆ = √ m a 2 − E a 2 / m a , is small. We expand the Einstein and Klein-Gordon equations in terms of a series in ∆. We find that the contribution of all operators with scaling dimension d > 4 vanish as a power of ∆. The leading order equations depend on the single parameter λ = fa M P ∆ . We find analytic formulas for the mass and the radius of the star as a function of λ, satisfied in the ranges f a / M P λ M P / f a . Parametrizing solutions of the equations of motion by λ, rather than the central density of the star, has the advantage of being able to compare total energies of solutions with equal numbers of axions as shown in figure 2. As expected, the critical temperature of the condensed axion star is very large, T c > 10 10 GeV.

Axion field dynamics in the condensed state
The dynamics of a self-gravitating, free hermitian scalar field was first analyzed in [2]. More recently, the authors of [11,12] used a similar procedure to describe the condensation of interacting hermitian Bose fields. We revisit this analysis by evaluating the expectation value of the axion potential in the N particle condensed state and simplifying the equations of motion using the expansion method described in the previous section. As stated earlier the range of parameters we explore differ from that of [11][12][13][14].

Expectation values
In our analysis, we consider axion dynamics to be described by a hermitian scalar field φ a with potential energy We will see that the leading-order deviation of the space-time metric from flat space is proportional to δ ≡ f 2 a /M 2 P , the effective coupling constant of matter to gravity. 1 If we therefore restrict ourselves to decay constants f a 10 −2 M P , it is sufficient to take gravity into account to leading order in δ 10 −4 .

JHEP03(2015)080
We start with the standard definition of the axion potential (2.1), but replace it by the expectation value of a quantum potential in which the quantum field, expanded in modes having definite radial, polar and azimuthal quantum numbers, is where a † n,l,m is the creation operator for an axion with the appropriate quantum numbers. To a very good approximation every axion in an axion star is in the ground state, so the star is an almost perfect condensate. As we will see later, for QCD axions the critical temperature of the condensate is 10 11 GeV. Consequently, in what follows, we will use the expression and, in appendix A, we calculate the expectation value of the axion potential exactly in the tree approximation. Note that, to justify using the tree approximations one needs to consider the fact that negative energy states and with them loop corrections can be neglected if binding energies are of ∆E a m a . We obtain in the large N limit an effective potential that is different from the classical potential. We thus obtain the equations of motion by taking the expectation value of the Einstein equation and the scalar field equation. For this, we will also need the expectation values (2.5)

Equations of motion
We write the spherically symmetric metric as We consider solutions for (f a / M P ) 2 = δ 1. When δ = 0 then the solution of the Einstein equations is reduced to A = B = 1 (flat metric). When δ is small, we can expand the equations of motion in a power series of δ. We will retain only up to first order contributions in δ. Even for axions of Grand Unified Theories, where δ = O(10 −4 ), second order contributions can be safely neglected.
The classical potential for the radial wave function, R, of the axion field is , but this is replaced by the expectation value of the tree level quantum potential (see appendix A), The expectation value of the Einstein equation and of the nonlinear equation of motion taken between states N − 1| and |N form a closed set of equations, Introducing the dimension-free radial coordinate z = r m, using the rescaled wave function X(z) = 2 √ N R(r) / f a , and defining = E / m, A = 1 + δ a(z) and B = 1 + δ b(z), we obtain in leading order of δ a, b, and X must be regular functions of z at z = 0 and vanish at z = ∞. While the radial wave function X and the metric function b have finite values at z = 0, a must vanish at z = 0. Aside from possible initial conditions, the equations depend on δ and the rescaled energy, .

Expansion in the binding energy of axions
Note that the dimension free coordinate, z = r m, measures the radial distance in units of the Compton wave length of the axion. We will show now that for axion stars with δ 1 and ∆ = √ 1 − 2 2(m − E)/m 1, (2.8) can be reduced to a system of equations depending on and δ through the combination λ = √ δ / ∆ = f a / (M P ∆) only. The parameter λ is not in general small. We will investigate bounds on possible values of λ, consistent with our approximations, later.
To expand our equations in ∆ first we expand the potential in a power series of the radial wave function, as Substituting (2.9) into (2.8) the axion equation of motion takes the form (2.10) It is easy to see that a systematic expansion of (2.10) in powers of ∆ can be performed if we further rescale the dimensionless radial coordinate by introducing x = ∆ z and simultaneously rescale the wave function as X(z) = ∆ Y (x), while keeping the scale of the dimensionless metric components, a and b, unchanged. The powers of ∆ extracted from each term are the "engineering" dimensions of the corresponding operator. All terms of (2.10) are of dimension three, except for the irrelevant, non-renormalizable terms, like X 5 and even higher powers of X, and of the relevant last term on the right hand side JHEP03(2015)080 of (2.10), which is of dimension one. Performing a systematic expansion in ∆ and keeping relevant and marginal terms only is tantamount to taking the infrared limit of the theory.
The leading order equations for the dimensionless field Y (x), depending on the dimensionless coordinate x are As we will see later, λ = √ δ / ∆, which is the only input parameter, determines the mass and radius of the axion star. We will also see that λ is a double valued function of N , the number of axions in the system. N is the natural physical input parameter.
Leading order corrections to (2.11) are of O(δ) and of O(δ λ 2 ), which both should be much less than 1. Then in addition to ∆ = f a / (M P λ) 1 we must have δ λ 2 = (λ f a / M P ) 2 1. These constraints can be combined together to give the range of validity of (2.11) as f a / M P λ M P / f a . Since e.g. for QCD axions f a / M P ∼ 10 −7 , solving the system (2.11) provides a correct solution for a wide range of sizes of axion stars.
We used the shooting method to integrate (2.11) and calculate the function, Y (x). Requiring the regularity of Y , a, and b at x = 0 we are left with two integration constants, which can be chosen as the value of Y (x) and b(x) at the center of the star. The boundary conditions are that Y (x), a(x), and b(x) tend to zero at x → ∞. a(x) and b(x) have newtonian asymptotics of a, b ∼ x −1 at x → ∞. Such boundary conditions are difficult to implement in the numerical calculations. However, notice that (2.11) implies and, consequently a(x) + b(x) tends to zero exponentially when x → ∞ because Y (x) also has similar behavior. Such a boundary condition can be imposed easily in a numerical calculation. We performed the numerical integration of (2.11) for a series of values of λ. As an example, in figure 1 we plot our solution for the wave function, Y (x), as a function of x at λ = 1. For all choices of λ we have considered, the wave function has a similar general shape. The initial values required for attaining the necessary asymptotic behavior are approximately proportional to the inverse of λ.

Analytic approximations to the physical parameters of axion stars
The most important parameters describing axion stars are the total mass, M , the radius inside which 99% of the matter contained in the star is concentrated, R 99 , and the number of axions in the star, N . The mass of the star is given in leading order of the infrared limit (∆ → 0) by where Furthermore, using the standard definition for boson stars, we define x 99 as Then restoring the appropriate scale, the radius of the axion star becomes Combining (3.1) and (3.4) we obtain the relationship between its radius and its mass. . (3.5) Using the function Y (x) found by numerical integration we calculated V (∞) and x 99 for a series of values 0.01 ≤ λ ≤ 10. As they are functionals of Y (x), they are uniquely fixed by the value of λ. The ratio x 99 / V (∞), appearing in (3.5) is larger than 0.1 throughout the range of λ we consider, leading to the relationship where R S is the Schwarzschild radius. This relationship, which implies gravitational stability is satisfied provided δ 1 / 500. As shown by (3.6), our formalism could even be applied to axions with f a ∼ 10 15 − 10 16 GeV, because for axions with decay constant f a ≈ 10 16 GeV, δ ∼ 10 −4 .
At For λ 0.5 the rise of R 99 is steeper, gives an excellent fit. R 99 continues to increase throughout the range of λ considered. The axion stars investigated in [11,12]  We plotted the dependence of the mass, scaled by the factor f a m a / M P , as a function of λ in figure 2. The striking feature of the plot is the maximum of the mass as a function of parameter λ = f a / (M P ∆). Maxima of M , as a function of the central value of the axion wave function have already been found in previous work [3,[6][7][8][9][10][11][12], but the maximum in figure 2 has a physical significance due to the following considerations. 2 Note that the number of axions in the condensate is approximately equal to     The maximum, M max = M (∆ max ) is attained at λ max = √ δ / ∆ max 0.58. Then the maximal number of axions in an axion star is For every N < N max , there are two masses, M (∆ 1 ) and M (∆ 2 ), such that (3.11) is satisfied. Assume ∆ 1 > ∆ 2 , or equivalently, λ 1 < λ max < λ 2 . Then M (∆ 1 ) < M (∆ 2 ), hence the relationship between the masses is illustrated by the dotted line in figure 2. As a result, every star corresponding to the branch of the curve at λ > λ max , which is plotted by a dashed line, is unstable. In fact, using (3.12) it is easy to estimate the mass difference between the two states of the axion star containing N axions. We obtain δM M δM / M 1, but as is shown in table 1, δM , still represents a substantial amount of macroscopic energy. Now, admittedly, there is an energy barrier between states M (∆ 2 ) and M (∆ 1 ). In the present paper we do not attempt to calculate the decay rate. That will be the subject of future research.
In table 1 we provide the mass, mass difference between stars containing the same number of axions (δM ), radius, and density of axion stars for several values of λ between 10 −1 ≤ λ ≤ 10. (Note that this range is well inside the interval f a /M P λ M P /f a , where our approximations are valid.) For the input values of f a and m a we use QCD axions, with parameters satisfying m a = 6 µeV × 10 12 GeV / f a [18,19]. If we pick the

JHEP03(2015)080
value m a = 10 −5 eV, we obtain f a = 6 × 10 11 GeV. Note, however, that in table 1, at fixed m a f a , the values of M , δM and d are increasing functions of f a , while R 99 is independent of the choice of f a .
A final comment concerns the critical temperature of the axion condensate. To see that the axion star is a pure condensate, with very little contribution from excited states one needs to consider its critical temperature. Using the average density, ρ = M / (m a V ) 3 M / (4 π R 3 99 m a ), and the expression for the radius in (3.6) we obtain the following estimate for the critical temperature (neglecting interactions) For a QCD axion, with f a ≈ 6 × 10 11 GeV, this would give T c ∼ 10 11 GeV.

Conclusions
As it has been pointed out in the previous section, to describe compact objects formed from axions one needs to consider that, except possibly in the early universe, the temperature of the object is many orders magnitudes smaller than the critical temperature of boson condensation. In other words, in a good approximation one may assume that all the axions are in the ground state. Therefore, along with Ruffini and Bonazzola [2] we substitute the wave function of the axions by a quantized field. As the resulting complicated nonlinear field theory cannot be solved analytically, one is restricted to use perturbation theory. A similar approach was used by Barranco and Bernal [11,12] in their application to axion stars, but they explore a region in which parameter λ is extremely small, λ = O(10 −6 ). The most important observations of this paper are that the application of the method of [2] to axion stars must be restricted in two different ways: (i) the decay constant of the axions, f a must be much smaller than the Planck mass, or using the notations of this paper δ = f a 2 / M P 2 1; (ii) only weakly bound axions should be considered, i.e. the axion binding energy should satisfy m a − E m a . The reason for requiring small binding energy is a requirement of using the Ruffini-Bonazzola method [2] for an interacting field theory: if the binding energy is of O(m a ) then the perturbation expansion of the axion field theory breaks down. Pairs of positive and negative energy states contribute to the ground state and to the expectation values of physical quantities, as well. To extend calculations beyond the region δ 1 and ∆ = 1 − E 2 / m 2 a 1 requires using non-perturbative field theory, which is beyond the scope of this paper.
Requiring ∆ 1 has allowed us to expand the equations of motion in the scale parameter ∆. If we consider that the scaling dimension of the axion field is unity, the expansion in ∆ is tantamount to taking the infrared limit of the axion field theory. Only marginal terms and a relevant term with coefficient λ 2 = δ / ∆ 2 of the Lagrangian give a contribution in leading order of ∆. Consequently, sixth and higher order terms of the transcendental axion potential do not contribute.
We have calculated the masses, radii and densities of axion stars as a function of the parameter λ = √ δ / ∆ over two orders of magnitude of the variable λ, consistent with the restrictions δ, ∆ 1. We found that the radius of the star increases approximately linearly JHEP03(2015)080 throughout the range of λ we considered. However, while the mass of the star, M (λ), rises linearly at small λ, it reaches a maximum at λ max 0.58 and tends to zero at large λ. The number of axions in the system, N (λ), which is determined uniquely by λ, is the same at two different values, λ 1 and λ 2 , where λ 1 < λ max < λ 2 . We show (see eq. (3.13)) that M (λ 2 ) > M (λ 1 ). This implies that all states on the branch λ > λ max of the M (λ) curve are unstable. For a QCD axion, assuming m a = 10 −5 eV we obtained M max ≡ M (λ max ) 10 19 kg and R max = R(λ max ) 1000 km. 3 For fixed f a m a = 6 × 10 −3 GeV 2 the radius is independent of f a , while the mass is an increasing function of f a . The maximal number of axions in an axion star is N max = 7.3 × 10 58 , but it is a fast increasing function of f a .
In the future we will also investigate rotating axion stars and analyze the maximal mass and the question of stability as a function of angular momentum. We will investigate collisions of axion stars in view of the existence of the maximal mass. Furthermore if an axion star is in isolation, unstable states, with λ > λ max should decay into the corresponding stable state. We will compute rate of transition, and possible signatures of their decay.
The masses of compact axion star solutions found in our work are consistent with the mass bounds derived by Tkachev for condensate formation through gravity [20]. Axion stars, along with free axions may form all or part of dark matter. Various production mechanisms of axions in the early universe are discussed in [21][22][23]. Bounds on the axion decay constant resulting from astrophysical and cosmological observations are discussed in [21] and [24][25][26]. We will investigate how our conclusions change if axion stars are in equilibrium within a cloud of free axions after reaching the maximal mass. We will also investigate the consequence if axions come in a multitude of flavors, as expected in theories derived from compact extra dimensions.
The collisions of axion condensates with neutron stars have been studied in the past [13,14,28]. Recently axion stars have been proposed as progenitors for Fast Radio Bursts (FRBs) [29][30][31][32] via collisions with neutron star atmospheres [33][34][35]. Whatever produces FRBs should be able to generate large amounts of energy, on the order of approximately 10 42 ergs/s, in a fairly tight frequency range around 1.4 GHz over time scales on the order of milliseconds. Assuming that all of the mass of an axion star is converted into radiation during a putative collision with a neutron star, one finds that an axion star mass on the order of 10 −12 M would be required. This mass is indeed compatible with our results listed in table 1. However, as the radius of our boson star is at least an order of magnitude larger that that of a neutron star we feel that the possibility of converting all the axions to radiation in such a collision is highly unlikely. Therefore we plan to perform an accurate estimation of this conversion rate in the near future.
The radius of the axion stars may be estimated from the duration of the bursts and was found in [29][30][31][32] to be on the order of 100 km, which is also compatible with our results in table 1. Moreover, [29][30][31][32] also showed that the frequency of the radiation can be generated by axions of mass approximately 10 −5 eV. We hope to examine this mechanism in greater detail and will report on the results of our investigation elsewhere. 3 Note that Rmax is the radius of the heaviest axion star but not the axion star of largest radius. The radii of unstable axion stars are always larger than Rmax.

JHEP03(2015)080
As mentioned, one result of our calculations is that for a fixed number of axions there are in general two masses corresponding to two possible values of ∆. An intriguing consequence is the possibility of tunneling from the state with a higher mass to the one with a lower mass. We will investigate whether collisions with neutron stars could induce such a tunneling process resulting in a fainter companion burst of approximately 10 29 ergs, where we have taken δM ∼ 10 5 kg for an axion star of mass 10 18 kg and radius 100 km. If detected, such a companion burst could serve to distinguish between axion star progenitors of FRBs and other proposals [36][37][38]. However, considering the faintness of the companion burst, it is unlikely to be observed from events that occur outside our galaxy. The observed frequency of FRB events (about 10 −3 events per galaxy per year) would therefore make this phenomenon even more difficult to observe.
Finally, we would like to emphasize again that there is no fundamental reason to limit calculations to theories with decay constants satisfying f a M P and stars with sufficiently small binding energy. We just state that the method of taking the expectation value of the equations of motion in free particle states, which are defined in flat space is not permissible if the axion decay constant is comparable with the Planck mass and/or the relative binding energy is not much smaller than 1. Using the current formalism to extend the results presented in this paper beyond leading order expansion in δ and ∆ would lead to false results. Such an extension would require using exact solutions of an interacting field theory, possibly in curved space, a calculation, which is beyond the scope of our current analysis. Restricting ourselves to leading order contributions has the added benefit of simplifying results to the extent of being able to obtain a clearer physical interpretation of the properties of axion stars.

JHEP03(2015)080
The Baker-Campbell-Hausdorff lemma implies that, if [X + [X + , X − ] = 0, then cos(X + + X − ) = e − 1 2 [X + ,X − ] e i X + e i X − + c.c. Omitting loop corrections and using free particle states, the right hand side of (A.4) can be readily calculated after expanding the two exponentials into power series of the creation and annihilation operators. Owing to the fact that in the condensate only ground state particles can be annihilated in the expectation values all contributions containing operators of the excited states vanish on the right hand side of (A.4). It follows that where R = R 1,0,0 is the single particle ground state radial wave function. When N → ∞ at fixed R / f a , the sum is dominated by k N . Then the last multiplier on the right hand side of (A.5) tends to N k giving our final result Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.