We consider a complex scalar field minimally coupled to gravity with an action given by S = ∫ d 4 − g − R κ + g μ ν ∇ μ Φ * ∇ ν Φ + V | Φ | 2 , (2.1) where R represents the Ricci scalar, g = det (g _μν ), κ = 16π, Φ represents the scalar field, and the star stands for the complex conjugate. V represents the self-interacting potential. The variation of the action 2.1 with respect to the metric tensor g _μν leads to Einstein’s equations. R μ ν − 1 2 R g μ ν = κ T μ ν , (2.2) where the stress–energy tensor is given by T μ ν = 1 2 ∇ μ Φ * ∇ ν Φ + ∇ μ Φ ∇ ν Φ * − 1 2 g μ ν ∇ α Φ * ∇ α Φ + V | Φ | 2 . (2.3)

The conservation of energy applied to the stress–energy tensor (2.3) reduces to the equation □ − d V d | Φ | 2 Φ = 0 (2.4) and its complex conjugate, where □Φ = g ^μν ∇ _μ ∇ _ν Φ.

We are interested in static spherically symmetric spacetimes, so we assume the metric can be written as d s 2 = − α r 2 d t 2 + a r 2 d r 2 + r 2 d Ω 2 , (2.5) where dΩ² = dθ ² + sin² θdφ ² is the metric defined on the two-sphere. In order to construct stationary configurations, we assume a scalar field with a time harmonic dependence of the form Φ r , t = ϕ r e − i ω t . (2.6)

As described in the work of Guerra et al. (2019), axion stars are constructed with the potential V = 2 μ 2 f a 2 B 1 − 1 − 4 B ⁡ sin 2 ϕ f a , (2.7) where B = 0.22 is a numerical factor that depends on the mass of the up and down quarks and μ is identified as the mass of the particles. In the limit f _a ≫ ϕ, the leading term in the power series of the potential leads to a potential with a quartic self-interaction of the form V ϕ ≈ μ 2 ϕ 2 − ℏ μ 2 f a 2 3 B − 1 12 ϕ 4 . (2.8)

From 2.8 we thus expect, the major contributions of the periodic potential appear in the limit of small f _a . With the ansatz for metric 2.5 and form 2.6 of the field, the Einstein-scalar field system 2.2, 2.4, 2.7 in spherical symmetry becomes 1 α d α d r = a − 1 2 r − r κ a μ 2 f a 2 2 B 1 − 1 − 4 B ⁡ sin 2 ϕ 2 f a + r a ϕ 2 ω 2 2 + r 2 d ϕ d r 2 , (2.9) 1 a d a d r = 1 − a r + r κ a μ 2 f a 2 2 B 1 − 1 − 4 B ⁡ sin 2 ϕ 2 f a + r a ϕ 2 ω 2 α 2 + r d ϕ d r 2 , (2.10) d 2 ⁡ ϕ d r 2 = − a ϕ ω 2 α 2 + d ϕ d r 1 2 a d a d r − 2 r − 1 α d α d r + a 4 π μ 2 f a 1 − 4 B ⁡ sin 2 ϕ 2 f a sin ϕ f a . (2.11)

To solve systems 2.9–2.11, one must choose appropriate boundary conditions, which are as follows. In order to guarantee that the spacetime is locally flat at the origin, the condition a (0) = 1 is required. Furthermore, we ask for d ϕ ( 0 ) d r = 0 . Additionally, for solutions that represent an isolated configuration, the scalar field must vanish at infinity. In this limit, asymptotic flatness is reached. These conditions reduce the system to an eigenvalue problem for ω such that for each choice of ϕ _c = ϕ(0), the system has a solution that decays exponentially at infinity. In order to solve the system numerically, it is convenient to use dimensionless quantities defined by r ̃ = μ r , ω ̃ = ω μ , ϕ ̃ = κ 2 ϕ , α ̃ = μ ω α . (2.12)

After substituting this scaling 2.12 in systems 2.9–2.11, the rescaled system is left in units of the mass parameter μ, which fixes the scale. Given a value of the field at the origin ϕ _c as a free parameter, we choose a trial value of the lapse and integrate the system outward from the origin using a fourth-order Runge–Kutta scheme with an adaptive step size. This adaptive scheme allowed us to reach smaller values of f _a . Then, we use a shooting algorithm to find the value of the lapse that corresponds to asymptotically flat solutions. We construct families of axion boson stars for different values of the parameter f _a and different values of the central scalar field. We also focus on solutions with no nodes in the scalar field corresponding to the ground state.

As demonstrated by Guerra et al. (2019), in the limit of large f _a , mini-boson stars are recovered, and this happens for values f a ∼ O ( 10 ) , which, in physical units, corresponds to values of f _a ∼ 1.22 × 10²⁰ GeV. In bosonic stars, the scalar field decays exponentially (since one asks for asymptotically flat spacetime), and consequently, they do not have a well-defined boundary; thus, it is common to describe their effective size in terms of the R ₉₉ radius, which is defined as the radius of the sphere containing 99% of the total mass of the star. In addition, we determine the star’s compactness as C = M R 99 .

The energy density ρ, radial pressure p _r , and tangential pressure p _t are defined in terms of the stress–energy tensor 2.3 as ρ = − T t t = 1 2 ω 2 α 2 ϕ 2 + 1 a 2 ∂ r ϕ 2 + 1 2 V ϕ , (2.13) p r = T r r = 1 2 ω 2 α 2 ϕ 2 + 1 a 2 ∂ r ϕ 2 − 1 2 V ϕ , (2.14) p t = T θ θ = T ϕ ϕ = 1 2 ω 2 α 2 ϕ 2 − 1 a 2 ∂ r ϕ 2 − 1 2 V ϕ . (2.15)

The total mass of the object M is given by the limit r → ∞ of the function M r = 4 π ∫ 0 r r ̃ 2 ⁡ ρ r ̃ d r ̃ . (2.16)

In our numerical calculations, we approach M as the value of this function evaluated at the outer point of the numerical grid. Alternatively, one can also find the total mass of the system by assuming that far away, the metric reduces to the Schwarzschild metric (Misner et al., 1973) so that M = lim r → ∞ r 2 1 − 1 a 2 . (2.17)

This expression converges very rapidly as r grows due to the exponential decay of the scalar field, and nonetheless, we checked that the results obtained with both expressions 2.16, 2.17 agree in the limit of large r.

It is known since the pioneer works of Kaup (1968) and Ruffini and Bonazzola (1969) that the mass of mini-boson stars has a maximum value M Kaup = 0.633 m Pl 2 / μ , where m _Pl represents the Planck mass. For f _a = 10, the system tends to the standard mini-boson star in which there is a local maximum at a critical value of ϕ _c that separates between stable and unstable configurations. This case is recovered for large values of f _a . Figure 2 displays the existence plots of the mass versus the central value of the field ϕ _c , for different configurations with some representative values of the decaying constant: f _a = 10 log₁₀ (f _a ) = −1.5, log₁₀ (f _a ) = −1.7, and log₁₀ (f _a ) = −2.0. The top left panel shown in Figure 2 corresponds to this limiting behavior. However, smaller values of the decay constant f _a lead to the formation of new branches. Furthermore, when the numerical value of f _a decreases, finding solutions becomes more challenging from the numerical point of view, and in this work, we found solutions for values up to log₁₀ f _a = −2.7. As shown in Figure 3, as f _a decreases, the existence diagram becomes more intricate with several local maxima and minima. The existence of these local extreme points indicate new stability branches at higher densities, giving rise to radially stable boson stars. Furthermore, the compactness of such solutions increases when the value of f _a decreases.

For larger f _a , the solutions reduce to the standard mini-boson stars, but in the limit f _a → 0, solutions with the appropriate boundary conditions, representing localized objects, are more difficult to find.

We now move to the description of the stars for a fixed value of f _a . For the sake of simplicity, we show the results for log₁₀ (f _a ) = −1.7. Figure 4 shows the parametric plots of radial and tangential pressures [p _r (r), p _t (r)] as given by (2.14 and 2.15. In the figure, the dashed line corresponds to the identity p _r = p _t . Larger deviations from this line, thus, correspond to larger anisotropy. The plots correspond to configurations marked with dots in Figure 2 and labeled with letters from A–F. Configuration A, with smaller values of ϕ _c , display only a small deviation from isotropy; however, as the value of ϕ _c increases, the stars become more anisotropic. For configuration E, the change of sign in the tangential pressure becomes more evident. However, the anisotropy is not as large as the one presented in other boson stars (Steven, 2012; Alcubierre et al., 2018; Alcubierre et al., 2019). The values of the radius, mass, and compactness for configurations A–F are listed in Table 1. As one moves to the right in the existence plot, with larger values of ϕ _c , axion stars become more compact. In the last column of Table 1, the numerical values of compactness are listed. As one moves to higher values of ϕ _c , the compactness increases. This trend is also present for the smaller values of f _a used in this work. Nonetheless, it seems that the compactness does not reach the maximum compactness limit (Buchdahl, 1966). Figure 5 displays the profiles of radial and tangential pressure as well as the density given in 2.13 of the configuration with the local maximum mass labeled as B and E in Figure 2. For configuration B, the difference in terms of absolute values between p _r and p _t is quite small. For configuration E, the same trend holds, but the tangential pressure becomes negative in the more exterior parts of the star. In both models, the maximum density is attained in the origin as in the usual mini-boson stars. Figure 5 shows the radial profiles of the metric coefficients α and a. In both models, the relation a/α → 1 holds at infinity, which is consistent with a Schwarzschild asymptotic behavior. In the more compact model E, a lower value of the lapse is attained at the origin, which is related with the fact that this configuration is more compact.

In general relativity, the maximum compactness of a self-gravitating, isotropic and spherically symmetric object made of a perfect fluid is M/R = 4/9, where M represents the mass of the object and R represents its radius. However, the above number, known as Buchdahl’s bound, relies strongly on the hypothesis of isotropy (Buchdahl, 1966). Considerable effort has been dedicated to model the properties of anisotropic matter, with the hope of finding physically viable models of compact stars. While anisotropies are generally negligible as compared to the pressure, it has been shown that even small anisotropies in fluid stars may induce significant changes on the mass and compactness of the star (Raposo et al., 2019). To determine the structure of a compact star, a widely followed path is to specify an equation of state and then solve the field equations. Customarily, this is carried out considering hydrodynamical equilibrium; however, in this work, we found that moderate anisotropic configurations exist for the potential (2.7), yielding to highly compact objects. Nonetheless, our results indicate that for the axion periodic potential, it is not possible to surpass Buchdahl’s bound.

It will also be helpful in identifying some general properties of the motion of test particles propagating in the spacetime associated with the previously found stars and, in particular, in determining whether the solutions admit stable circular orbits or light rings. The interest of studying light rings in compact objects has been renewed because of the correspondence between the quasi-normal modes of black holes and light ring oscillations to describe the initial part of the ringdown gravitational-wave signal of black holes (Cardoso et al., 2009; Cardoso and Pani, 2017; Khanna and Price, 2017).

In order to describe the motion of test particles in the spacetime (2.5), let us consider the Lagrangian equation as follows: 2 L = − α 2 t ̇ 2 + a 2 r ̇ 2 + r 2 θ ̇ 2 + r 2 ⁡ sin 2 ⁡ θ φ ̇ 2 , (3.1) where “dot” denotes the derivative with respect to the affine parameter τ. It is possible to consider 2 L = − k so that if k = 1, it corresponds to time-like geodesics and if k = 0, it corresponds to null geodesics.

Since the spacetime is spherically symmetric, we focus on particles in the plane θ = π/2 in 3.1. Associated to the time and angular symmetry, there are two conserved quantities during the motion of the particles, namely, the energy at infinity E and the angular momentum ℓ given by E = ∂ L ∂ t ̇ = − α 2 t ̇ , ℓ = ∂ L ∂ ϕ ̇ = r 2 ϕ ̇ . (3.2)

Normalization of the four velocities u ^μ u _μ = −k can be written in terms of the conserved quantities 3.2 yielding the equation − k = − E 2 α 2 + a 2 r ̇ 2 + ℓ 2 r 2 , (3.3) Equation 3.3 can be written as a 2 α 2 r ̇ 2 = E 2 − α 2 k + ℓ 2 r 2 . (3.4)

Defining a new variable z, through d z d r = a α , equation (3.4) becomes z ̇ 2 = E 2 − V eff , (3.5) where the effective potential in 3.5 is defined as V eff = α 2 k + ℓ 2 r 2 . (3.6)

Circular motion of particles is possible when the conditions r ̇ = 0 and V eff d r = 0 applied to 3.6 are fulfilled. In the time-like case, these conditions completely specify the energy and angular momentum to be E 2 = α 3 α − r α ′ and  ℓ 2 = α ′ r 3 α − r α ′ , (3.7) where the prime denotes the derivative with respect to r. For null geodesics, the existence of circular orbits is given by the condition h(r) = α − rα′ = 0 for which expressions (3.7) are undetermined. Moreover, it has been proven that regular configurations can have two light rings, of which one is stable (Pedro et al., 2017). However, not all configurations considered here admit light rings, and their appearance will depend on the compactness of the star. On the other hand, regarding massive particles, there always exist stable circular orbits. Figure 6 shows the potential (3.6) for null geodesics (with k = 0) for configurations A to F in Figure 2 with log₁₀ (f _a ) = −1.7. For this value of f _a , configurations are not compact enough to have light rings. Since the potential for these configurations is very alike, the motion of null particles around them is quite similar. In the right panel of Figure 6, the function h(r) is displayed, and the change of behavior close the origin for configurations D–F is observed, and there is a slight decay on the value of h(r). It can be seen in the figure that as the configurations become more compact, h decreases in such a way that it will cross the origin twice, giving rise to the existence of a pair of light rings. For configurations with lower values of log₁₀ (f _a ) ≈ − 2.3, this is actually the case, and a pair of light rings appear. Figure 7 displays the effective potential for massive particles (k = 1) with ℓ = 1 for configurations A–F. In all cases, the potential is minimum, indicating the existence of circular orbits.