Prepared for submission to JHEP

Non-topological soliton Bardeen boson stars and its

frozen state
arXiv:2502.14153v1 [gr-qc] 19 Feb 2025

Zhen-Hua Zhao∗,a , Yining Gu a , Shucong Liu b , Long-Xing Huang b , Yong-Qing

Wang ∗,b
a
Department of Applied Physics, Shandong University of Science and Technology, Qingdao 266590,
China

b
School of Physical Science and Technology, Lanzhou University, Lanzhou 730000, China
E-mail: zhaozhh78@[Link], yqwang@[Link], ∗ corresponding
author

Abstract:
We investigate a Bardeen model coupling Einstein gravity with nonlinear electro-
magnetic fields and non-topological soliton complex scalar fields, governed by the mag-
netic charge q̃, the complex scalar field frequency ω̃, and the self-interaction parameter
η̃. Our results reveal that the magnetic charge q̃ exhibits η̃-dependent critical values q̃c ,
beyond which (q̃ > q̃c ) Bardeen boson stars (BBSs) may transition into frozen states
(ω̃ → 0). These frozen states are characterized by a critical horizon whose radius r̃cH
H,RN H,RN H,RN RN denote the inner and outer horizons of
satisfies r̃inner < r̃c < r̃outer , where r̃inner and r̃outer
magnetic Reissner-Nordström (RN) black holes with equivalent mass and magnetic charge.
Notably, the ADM mass of frozen BBSs is independent of η̃. Furthermore, light ring (LR)
solutions exist universally across all tested combinations of q̃ and η̃, with all frozen BBSs
exhibiting LRs whose outer radius r̃outer LR is independent of η̃. Compared to magnetic RN
black holes, frozen BBSs possess a smaller outer LR radius (r̃outer LR < r̃ LR, RN ).
outer
Contents

1 Introduction 1

2 The model 2

3 Numerical scheme 5

4 Numerical results 5
4.1 Solutions 6
4.1.1 Case of q̃ < q̃c 6
4.1.2 Case of q̃ > q̃c 6
4.2 Frozen states and critical horizon 8
4.3 Compactness 14
4.4 Light Rings 14

5 Conclusions 16

6 Acknowledgment 19

1 Introduction

The existence of black holes has been reinforced by the discovery of gravitational waves
produced by the merger of two black holes [1] and the observation of black hole shadows by
the Event Horizon Telescope (EHT) [2–4]. The realisation of these experiments represents
a significant milestone in the process of black hole research. However, there remain a
number of long-standing issues that require further investigation, including the existence of
singularities and singular rings, as well as the so-called Black hole information paradox [5–7].
Penrose’s singularity theorem [8] led to the understanding that the existence of a singularity
in the solution of a black hole in general relativity appeared to be an intrinsic feature of
the theory. Although the weak cosmic censorship conjecture [9] permits singularities of
gravitational collapse to be contained within the event horizon, which effectively shields
them from external observers, this does not alter the fundamental reality of the existence
of singularities.
The existence of a singularity implies that matter can be compressed to a point, which
may be acceptable as a mathematical object, but challenging to accept as a physical entity
[10, 11]. The study of the nature of the singularity has also prompted interest in the study
of singularity-free black holes, known as regular black holes [12–19]. The earliest known
study in this field was proposed by Bardeen in 1968 [12], who suggested a metric for a
black hole without singularities, which is now known as the “Bardeen black holes” [17].

–1–
However, obtaining exact solutions for singularity-free black holes by solving the Einstein
equations remains highly challenging. It was not until 1999 that Ayón-Beato and García
[20] proposed the first exact regular black hole solution using nonlinear electrodynamics
coupled to general relativity. Subsequently, they also derived the exact solution for the
Bardeen black hole [21], interpreting it as a nonlinear magnetic monopole.
The construction of general relativistic objects by introducing electromagnetic fields
and Einstein gravitational coupling can be traced back to the work of Wheeler in the 1950s,
who termed this “gravitational-electromagnetic entity” as geons [22]. In the 1960s, Kaup
[23] suggested massive complex scalar fields as an alternative to the electromagnetic field.
The stable local solutions that resulted from this were known as “Klein-Gordon Geon” and
subsequently as boson stars (BSs).
Unlike black holes, BSs lack event horizons. However, recent studies [24, 25] have shown
that, under specific conditions, BSs can possess an critical horizon that closely resembles the
event horizon of black holes. From the perspective of an observer at infinity, time within the
critical horizon appears to slow down infinitely, causing matter to asymptotically approach
but never cross the horizon within a finite time. BSs in this state are called “frozen stars"
. The concept of frozen stars resulting from gravitational collapse in a vacuum was first
introduced by Oppenheimer and Snyder [26].
The frozen stars studied in Ref. [24] are are constructed within the framework of the
Bardeen model, where Einstein gravity is coupled to a nonlinear electromagnetic field and
a free complex scalar field. In this paper, based on Wang’s work [24], we further explore
scenarios involving self-interactions of the complex scalar field. Specifically, we focus on
the case where the complex scalar field adopts a non-topological configuration. Its coupling
with Einstein gravity can give rise to a non-topological soliton star, a concept first proposed
by Friedberg, Lee, and Pang [27, 28]) and subsequently studied extensively in, for example,
Refs. [29–37].
The organization of this paper is as follows. In Sec. 2, we present the model what we
discussed. In Sec. 3, we outline Numerical Scheme. In Sec. 4, we present the numerical
results and the discussion in four topics. In Sec. 5, we present some conclusions.

2 The model

This section provides a brief introduction to the theoretical framework, which includes the
Einstein nonlinear electrodynamics model [20] coupled to a self-interaction complex scalar
field. We adopt the natural units where ℏ = c = G = 1. The action is given by follows:
√
 
R
Z
S= −gd4 x + L(1) + L(2) , (2.1)
16π
where the Lagrangian densities are defined as
p !5
2
3 2q 2 F
L(1) =− p , (2.2)
2s 1 + 2q 2 F
L(2) = −∇µ Φ∗ ∇µ Φ − U (Φ, Φ∗ ). (2.3)

–2–
Here,

2
U (ΦΦ∗ ) = µ2 ΦΦ∗ 1 − 2η 2 ΦΦ∗ , (2.4)
R is the scalar curvature, Φ is a complex scalar field, Φ∗ is the complex conjugate of Φ, F =
1 µν with the electromagnetic field strength F
4 Fµν F µν = ∂µ Aν − ∂µ Aν , and Aµ represents
the electromagnetic field. The constants q, s, µ, and η are independent parameters, where q
denotes the magnetic charge, µ the scalar field mass parameter, and η the coupling constant
of the scalar fields.
When q = 0 while the complex scalar field remains non-vanishing, the action describes
a non-topological soliton black hole (BS) theory. Conversely, if the complex scalar field
vanishes but q ̸= 0, the model reduces to the Bardeen black hole theory. The line element
of the Bardeen theory in spherical coordinates (t, r, θ, φ) is given by

ds2 = −f (r)dt2 + f (r)−1 dr2 + r2 dθ2 + sin2 θdφ2 ,

(2.5)
where
q3 r2
f (r) = 1 − . (2.6)
s(r2 + q 2 )3/2
For q ≥ qb = 33/4 s/2, the equation f (r) = 0 has a real root, indicating a Bardeen
p

spacetime solution with an event horizon. When q < qb , there is no event horizon, and thus
no black hole solution. The asymptotic behavior of f (r) at infinity is

q3
f (r) = 1 − 2 /r + O(1/r3 ). (2.7)
2s
The Arnowitt-Deser-Misner (ADM) mass M of this model, derived from the second term,
is given by M = q 3 /(2s) [38].
Varying the action (2.1) with respect to gµν , Aµ , and Φ yields the equations of motion:

1 (1) (2)
Rab − gµν R − 8π(Tµν + Tµν )=0, (2.8)
2 !
∂L(1) µν
∇µ F = 0, (2.9)
∂F
□Φ − U̇ (ΦΦ∗ )Φ = 0. (2.10)

Here,
∂U (Φ, Φ∗ )
U̇ (ΦΦ∗ ) = = µ2 1 − 4η 2 (ΦΦ∗ ) + 3η 4 (ΦΦ∗ )2 ,
 
∗
(2.11)
∂(ΦΦ )

(1) ∂L(1)
Tµν =− Fµρ Fν ρ + gµν L(1) , (2.12)
∂F
 
(2) ∗ ∗ 1 λρ ∗ ∗ ∗
Tµν = ∂µ Φ ∂ν Φ + ∂ν Φ ∂µ Φ − gµν g (∂λ Φ ∂ρ Φ + ∂ρ Φ ∂λ Φ) + U (Φ, Φ ) . (2.13)
2
The action (2.1) is invariant under the global U (1) transformation of the complex scalar
field, Φ → eiα Φ, where α is a constant. According to Noether’s theorem, this symmetry
leads to a conserved charge:

–3–
 Z
Q= J 0 dV, (2.14)
Σ

where Σ is a spacelike hypersurface and J 0 is the timelike component of the conserved

current:
J µ = −i (Φ∗ ∂ µ Φ − Φ∂ µ Φ∗ ) . (2.15)

We consider the generic static, spherically symmetric spacetime with the line element

dr2
ds2 = −N (r)σ 2 (r)dt2 + + r2 dθ2 + sin2 θdφ2 ,


(2.16)
N (r)

where N (r) = 1 − 2m(r)

r . The functions N (r), m(r), and σ(r) depend only on the radial
coordinate r. The ansatzes for the electromagnetic field [20] and the scalar field [27, 28] are
given by
θ φ ϕ(r)
Fµν = 2δ[µ δν] q sin(θ), Φ = √ e−iωt , (2.17)
2
where ϕ(r) is a real scalar field and ω is the frequency of Φ.
The Noether charge obtained from Eq. (2.14) is written as
∞
ωϕ2
Z
Q = 8π r2 dr . (2.18)
0 N σ

Substituting the ansatzes (2.16) and (2.17) into the field equations (2.8), we obtain the
following equations:

2 N ′ σ′
   2 
′′ ′ ω ϕ
ϕ + + + ϕ + 2
− U̇ (ϕ) = 0, (2.19)
r N σ Nσ N
ω 2 rϕ2 N 12πr 1
N ′ + 4π 2
+ 4πN rϕ′2 + + r8πU (ϕ) + 5/2
− = 0, (2.20)
Nσ r s (1 + r2 /q 2 ) r
′ 2 2
 
σ ω ϕ
− 4πr ϕ′2 + 2 2 = 0, (2.21)
σ N σ

where
µ2 2
2
U (ϕ) = ϕ 1 − η 2 ϕ2 ,
2
∂U (ϕ)
= µ2 1 − 4η 2 ϕ2 + 3η 4 ϕ4 ,


U̇ (ϕ) = 2
∂(ϕ )
2m′ (r) 2m(r)
N′ = − + ,
r r2
and the prime denotes differentiation with respect to r.
To solve these ordinary differential equations, we impose the following boundary con-
ditions to ensure regularity at the origin and asymptotic flatness at infinity:

m(0) = 0, σ(∞) = 1, m(∞) = M, σ(0) = σ0 . (2.22)

–4–
where M is the ADM mass, and σ0 is a constant determined by solving the differential
equations. The scalar field boundary conditions are:

dϕ(r)
ϕ(∞) = 0, = 0. (2.23)
dr r=0

3 Numerical scheme

To facilitate numerical computations, we apply the following scaling transformations:

r̃ q̃ √ 4πs̃ m̃ ϕ̃
r→ , q → , η → η̃ 4π, ω → ω̃µ, s → 2 , m → , ϕ → √ , (3.1)
µ µ µ µ 4π
M̃ Q̃
M→ , and Q → 2 , (3.2)
µ µ

where
∞
ω̃ ϕ̃2
Z
Q̃ = 2 r̃2 dr̃ . (3.3)
0 N σ
After applying these transformations, the equations (2.19)–(2.21) simplify by setting:

µ → 1 and 4π → 1. (3.4)

This demonstrates that the equations are scale-invariant with respect to the mass parameter
µ.
To further simplify the numerical implementation, we introduce a compactified radial
coordinate:
r̃
x= . (3.5)
1 + r̃
This transformation maps the domain [0, ∞) to [0, 1], with x = 0 corresponding to the
origin and x = 1 to spatial infinity.
The open-source software Fenicsx (version 0.8) [39–42] is employed to solve the differ-
ential equations (2.19), (2.20), and (2.21) subject to the boundary conditions (2.22) and
(2.23). The computational domain [0, 1] is discretized using at least 3000 grid points, en-
suring a relative numerical error of less than 10−5 . In the subsequent analysis, we focus on
the influence of the parameters q̃, η̃, and ω̃, while setting s̃ = 0.2 without loss of generality.

4 Numerical results

In this section, we present our numerical results. When η̃ = 0, the model reduces to the one
studied in Ref. [24]. It is known that when the value of q̃ exceeds a critical value, denoted
as q̃c , a solution with frequency ω̃ → 0 exists. In our numerical calculations, the smallest
value of ω̃ is set to 0.0001. However, in principle, solutions with even smaller values of ω̃
can be obtained by reducing the step sizes of x and ω̃.

–5–
 η̃ 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
q̃c 0.59979 0.60015 0.60317 0.6028 0.6021 0.60165 0.60136 0.60115 0.60098

Table 1. The relation between η̃ and q̃c . There are different values of q̃c for different values of η̃.

The model under consideration includes a self-interaction term for the complex scalar
field. Our calculations reveal that the critical value q̃c varies with the self-coupling con-
stant η̃. The values of q̃c for η̃ = 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 are computed and
summarized in Table 1. From the table, we observe that the variation in η̃ has a subtle but
non-negligible effect on q̃c , with the maximum difference between q̃c values not exceeding
0.004. Although this difference is small, its implications cannot be overlooked and will be
further discussed in Subsection 4.1.2.

4.1 Solutions
In general, if the values of q̃ and ω̃ are fixed, the solutions of the equations (2.19), (2.20),
and (2.21) can cover a wide range of values of η̃. For example, η̃ can take values in the range
[0, 100]. We present some representative solutions in Fig. 1, with q̃ = 0.3 and ω̃ = 0.9.

4.1.1 Case of q̃ < q̃c

In Fig. 2(a), we illustrate the relationship between the minimum value of ω̃, denoted by
ω̃min , and η̃. As η̃ increases, ω̃min generally decreases. Additionally, as ω̃ decreases, the
problem becomes increasingly stiff and more challenging to solve. To demonstrate this,
the solutions for η̃ = 5 are presented in Fig. 2(b), Fig. 2(c), and Fig. 2(d). Notably, the
slope of the solutions for ϕ̃(x), m̃(x), and σ(x) around x = 0.95 increases significantly as ω̃
decreases.
When η̃ is relatively small, for instance, when η̃ = 1, the M̃ versus ω̃ curves maintain
a classical spiral shape, as depicted in Fig. 3(a) with q̃ = 0.0, 0.3, 0.5. As η̃ increases,
such as when η̃ = 2, the curve exhibits a bimodal structure, as shown in Fig. 3(b) with
q̃ = 0.0, 0.3. Furthermore, there is a significant increase in the maximum value of M̃ as
η̃ increases, as illustrated in Fig. 3(c). As both q̃ and η̃ increase, the M̃ versus ω̃ curves
may adopt a very different shape from the classical one, for example, as shown in Fig. 3(d),
where cusps appear at the left end of the curves, and the curves almost fold into each other.
Additionally, the lower branch of the curve intersects with the upper branch when q̃ = 0.5
and η̃ = 4.
The Q̃ versus ω̃ curves exhibit similar characteristics to the M̃ versus ω̃ curves, as
shown in Fig. 4(a). Notably, the upper and lower branches of the curve do not cross when
q̃ = 0.5 and η̃ = 4, as illustrated in Fig. 4(b).

4.1.2 Case of q̃ > q̃c

In this case, the shapes of the M̃ versus ω̃ and Q̃ versus ω̃ curves differ significantly from
those in the q̃ < q̃c case, with the most notable difference being the absence of the inflection
point at the left end of the curves. For example, in Fig. 5 and Fig. 6, we present diagrams

–6–
 0.25
η̃ =0.0 0.5 η̃ =0.0
η̃ =1.0 η̃ =1.0
0.20
η̃ =10.0 0.4 η̃ =10.0
0.15 η̃ =100.0 η̃ =100.0
0.3

m̃(x)
φ̃(x)

0.10 0.2

0.05 0.1

0.00 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x

(a) (b)

1.00

0.95
η̃ =0.0
η̃ =1.0
σ(x)

0.90
η̃ =10.0
η̃ =100.0
0.85

0.80
0.0 0.2 0.4 0.6 0.8 1.0
x

(c)

Figure 1. Solutions of the fields functions with q̃ = 0.3, ω̃ = 0.9 and different values of η̃.

of the M̃ versus ω̃ and Q̃ versus ω̃ curves, respectively, for q̃ = 0.63, 0.65, 0.67, 0.7. It is
evident that for all values of η̃, there are no inflection points at the left end of the curves.
Furthermore, the M̃ versus ω̃ curves converge to the same point at both the left and right
endpoints, with the values of ω̃ corresponding to the left endpoints being 0.0001 for all
curves.
As previously mentioned, it is theoretically possible to obtain solutions with smaller
values of ω̃ by decreasing the step sizes of x and ω̃. However, the stiffness of the problem
increases as ω̃ decreases. As shown in Fig. 7, with q̃ = 0.61, η̃ = 3, and ω̃ = 0.1, 0.01, 0.0001,
the slopes of the solutions for the functions ϕ̃ and σ increase dramatically with decreasing
ω̃ around the point x = 0.35. This indicates that the smaller the value of ω̃, the more
difficult the equations are to solve. Therefore, in this paper, we limit the minimum value
of ω̃ to 0.0001.
In Table 1, there is a very slight difference in q̃c for different values of η̃, but it cannot
be ignored. For example, the value of q̃ = 0.6 is the critical value for η̃ = 0, but not for
η̃ = 1 and η̃ = 2. For η̃ = 0, the lowest frequency value is 0.0001, whereas for η̃ = 1 and
η̃ = 2, the lowest frequency values are 0.01586 and 0.01104, respectively. Consequently, as

–7–
 0.8
q̃ =0.0
0.7 q̃ =0.3 0.20
q̃ =0.5
0.6 0.15
ω̃min

0.5

φ̃(x)
0.10
0.4
ω̃ =0.4
0.3 0.05
ω̃ =0.2
0.2 ω̃ =0.14021
0.00
0 1 2 3 4 0.0 0.2 0.4 0.6 0.8 1.0
η̃ x

(a) (b)

ω̃ =0.4 1.0
5
ω̃ =0.2
4 ω̃ =0.14021 0.8
m̃(x)

3
σ(x)

0.6
2
ω̃ =0.4
1 0.4
ω̃ =0.2
ω̃ =0.14021
0 0.2
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x

(c) (d)

Figure 2. (a)ω̃min versus η̃. (b)-(d) Solutions of the field functions with η̃ = 5 and q̃ = 0.3.

shown in Fig. 8, the M̃ versus ω̃ curves corresponding to η̃ = 1 and 2 exhibit an inflection

point at the left end.
Another consequence is the effect on the metrics. To illustrate this effect, we present the
images of the metric components g rr and −gtt in Fig. 9. From these images, it is evident
that the minima of g rr and −gtt at η̃ = 0 are approximately three orders of magnitude
smaller than the minima of g rr and −gtt at η̃ = 1 and 2. Moreover, the minima of g rr and
−gtt at η̃ = 0 can be further reduced relative to the case of η̃ = 1 and 2 by increasing the
number of coordinate intervals and decreasing the step size of ω̃.

4.2 Frozen states and critical horizon

In the case of q̃ > q̃c , as ω̃ → 0, the value of −g00 also tends to zero in the central region
inside the BBSs. As illustrated in Fig. 10(a), within this region, the value of −g00 can be
reduced by an order of magnitude, reaching 10−7 . The boundary of this region, denoted
by xc , is determined by the minimum value of g rr , referred to as g rr min . Notably, g rr min is
also reduced by an order of magnitude, reaching 10−7 , as shown in Fig. 10(b).

–8–
 0.5 0.7 q̃=0.0
η̃=2.0 q̃=0.3
0.6
0.4 q̃=0.5
0.5
η̃=1.0
M̃

M̃
0.3 0.4

q̃=0.0 0.3
0.2
q̃=0.3 0.2
q̃=0.5
0.1 0.1
0.4 0.6 0.8 1.0 0.4 0.6 0.8 1.0
ω̃ ω̃

(a) (b)

η̃=0.0 2.5
3 η̃=3.0
q̃=0.3 η̃=1.0 q̃=0.5 η̃=4.0
η̃=2.0
2.0
η̃=3.0
2
η̃=4.0
M̃

1.5
M̃

1 1.0

0.5
0
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8
ω̃ ω̃

(c) (d)

Figure 3. M̃ versus ω̃ with q̃ < q̃c .

12
q̃=0.3, η̃=0.0 q̃=0.5, η̃=3.0
10 q̃=0.3, η̃=1.0 q̃=0.5, η̃=4.0
6
q̃=0.3, η̃=2.0
8
q̃=0.3, η̃=3.0
q̃=0.3, η̃=4.0 4
Q̃

6
Q̃

4
2
2

0 0
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8
ω̃ ω̃

(a) (b)

Figure 4. Q̃ versus ω̃ with q̃ < q̃c .

–9–
 0.75
η̃=0.0 0.74 η̃=0.0
η̃=1.0 η̃=1.0
η̃=2.0 0.72 η̃=2.0
0.70 η̃=3.0 η̃=3.0
η̃=4.0 0.70 η̃=4.0
M̃

M̃
0.65 0.68

0.66
0.60
q̃=0.61 0.64 q̃=0.63
0.62
0.0 0.2 0.4 0.6 0.0 0.2 0.4 0.6
ω̃ ω̃

(a) (b)

η̃=0.0 0.885 η̃=0.0

0.76 η̃=1.0 η̃=1.0
η̃=2.0 0.880 η̃=2.0
η̃=3.0 η̃=3.0
0.74 η̃=4.0 0.875 η̃=4.0
M̃

M̃

0.870
0.72
0.865
0.70
q̃=0.65 0.860 q̃=0.7

0.0 0.2 0.4 0.6 0.0 0.1 0.2 0.3

ω̃ ω̃

(c) (d)

Figure 5. M̃ versus ω̃ and Q̃ versus ω̃ with q̃ > q̃c .

Fig. 10(c) demonstrates that within the specified range where x < xc , −g00 is less than
g rr
min . This phenomenon is universal, as evidenced by Fig. 10(d). Consequently, when the
minimum value of g rr reaches the order of 10−7 , the sphere defined by the radial coordinate
x = xc is termed the “critical horizon” [25, 43, 44].
This implies that, from the perspective of an observer at infinity, time nearly halts
inside the critical horizon, and matter approaching the surface of the star appears to stop
near the critical horizon. Thus, the BBSs are identified as frozen at this moment.
To gain an intuitive understanding of the critical horizon radius, we compare it with
the event horizon radius of the magnetic Reissner-Nordström (RN) black hole. The line
element of the magnetic RN black hole is given by

ds2 = −∆(r̃)dt2 + ∆(r̃)−1 dr̃2 + r̃2 dθ2 + sin2 θdφ2 ,

(4.1)

where
2M̃ q̃ 2
∆(r̃) = 1 − + 2.
r̃ r̃

– 10 –
 0.25
η̃=0.0 0.125 η̃=0.0
q̃=0.61 q̃=0.63
η̃=1.0 η̃=1.0
0.20
η̃=2.0 0.100 η̃=2.0
0.15 η̃=3.0 η̃=3.0
η̃=4.0 0.075 η̃=4.0
Q̃

Q̃
0.10
0.050

0.05 0.025

0.00 0.000
0.0 0.2 0.4 0.6 0.0 0.2 0.4 0.6
ω̃ ω̃

(a) (b)

0.08
η̃=0.0 0.025
q̃=0.65 η̃=0.0
q̃=0.7
η̃=1.0 η̃=1.0
0.06 0.020
η̃=2.0 η̃=2.0
η̃=3.0 η̃=3.0
0.015
0.04 η̃=4.0 η̃=4.0
Q̃

Q̃

0.010
0.02
0.005

0.00 0.000
0.0 0.2 0.4 0.6 0.0 0.1 0.2 0.3
ω̃ ω̃

(c) (d)

Figure 6. M̃ versus ω̃ and Q̃ versus ω̃ with q̃ > q̃c .

1.0
ω̃ =0.1
0.5 ω̃ =0.01
0.8 0.05
ω̃ =0.0001
0.4
0.6
0.00
σ(x)
φ̃(x)

0.3 0.05
0.350 0.375
0.4
0.2
ω̃ =0.1
0.00 0.2
0.1 ω̃ =0.01
0.375 0.400
ω̃ =0.0001
0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x

(a) (b)

Figure 7. Solutions of ϕ̃ and σ with q̃ = 0.61 and η̃ = 3.

– 11 –
 η̃=0.0
0.64 η̃=1.0
η̃=2.0
0.62

0.60

M̃
0.649
0.58

0.56 0.648

0.54 0.000 0.025

0.0 0.2 0.4 0.6 0.8
ω̃

Figure 8. M̃ versus ω̃ with q̃ = 0.6. The minimum frequencies corresponding to the curves with
η̃ = 0, η̃ = 1.0, and η̃ = 2.0 are ω̃ = 0.0001, ω̃ = 0.01586, and ω̃ = 0.01104, respectively.

100
η̃ =0.0
−1
10 η̃ =1.0
×10−6 10−2 η̃ =2.0
1.0
−gtt (x)

10−3 10−4
g rr (x)

0.5

0.3400 0.3425 10−6

10−5 η̃ =0.0
η̃ =1.0
10−8
η̃ =2.0
10−7
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x

(a) (b)

Figure 9. Solutions of g rr (x) and −gtt (x) with q̃ = 0.6. The frequencies corresponding to the
curves with η̃ = 0, η̃ = 1.0, and η̃ = 2.0 are ω̃ = 0.0001, ω̃ = 0.01586, and ω̃ = 0.01104, respectively.

The event horizon radii of the magnetic RN black hole can be determined from the condition
∆ = 0:
 s 
q
q̃ 2
H,RN
r̃inner/outer = M̃ ∓ M̃ 2 − q̃ 2 = M̃ 1 ∓ 1 − ,
M̃ 2

H,RN H,RN
where r̃inner and r̃outer denote the inner and outer event horizons of the magnetic Reissner-
Nordström (RN) black hole, respectively.
Our calculation results are summarized in Table 2, where the critical horizon radius
of a BBS is denoted by r̃cH . From the table, we identify two key trends: (1) the ADM
mass of the frozen BBSs increases monotonically with q̃; and (2) r̃c consistently satisfies
H,RN H,RN H,RN H
the inequality r̃inner < r̃cH < r̃outer . Additionally, similar to r̃outer , r̃c exhibits a monotonic
increase with q̃.

– 12 –
 10−1 10−1
q̃=0.63 q̃=0.63 g rr min =8.74 × 10−7
−3
10
−gtt (x)

10−3

g rr (x)
10−5 ω̃ =0.1 ω̃ =0.1
ω̃ =0.01 10−5 ω̃ =0.01
10−7 ω̃ =0.001 ω̃ =0.001
ω̃ =0.0001 ω̃ =0.0001
10−7
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x

(a) (b)

10−1 10−1
q̃=0.63
−3
10 10−3 q̃ =0.63, g rr
q̃ =0.63, −gtt
q̃ =0.65, g rr
10−5 10−5 q̃ =0.65, −gtt
g rr q̃ =0.67, g rr
−7 q̃ =0.67, −gtt
10 −gtt 10−7
q̃ =0.7, g rr
−gtt < g rr min q̃ =0.7, −gtt

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x

(c) (d)

Figure 10. Solutions of −g00 and g rr , the vertical axes of are represented on a logarithmic scale,
and η̃ = 1.

H,RN H,RN
M̃ r̃cH r̃outer r̃inner
q̃ = 0.61 0.67 0.59 0.95 0.39
q̃ = 0.63 0.72 0.69 1.06 0.37
q̃ = 0.65 0.77 0.77 1.17 0.36
q̃ = 0.67 0.81 0.84 1.27 0.35
q̃ = 0.70 0.86 0.98 1.43 0.34

Table 2. Radii of the critical horizons of frozen BBSs and the event horizons of magnetic RN
block holes. All numerical values are rounded to two decimal places.

– 13 –
4.3 Compactness
Each M̃ versus ω̃ curve exhibits a maximum value of M̃ , denoted by M̃max . The compact-
ness discussed in this section is specifically evaluated at this maximum mass value.
The definition of compactness adopted here is presented in Ref. [45]. The radius X is
determined by solving the equation:
m̃(X)
= 0.99. (4.2)
M̃max
Additional definitions of compactness can be found in Ref. [35]. Using Eq. (3.5), the corre-
sponding physical radius R̃ can be derived as
X
R̃ = , (4.3)
1−X
and the compactness is then calculated by

M̃max
C= . (4.4)
R̃
Simultaneously, the frequency corresponding to M̃max , denoted by ω̃ ∗ , can also be identified.
Fig. 11 illustrates the variations of M̃max and ω̃ ∗ as functions of η̃. For the case of
q̃ < q̃c , as η̃ increases, M̃max initially decreases and then increases, with the maximum
difference in M̃max exceeding 2.0. When η̃ > 1, the value of ω̃ ∗ decreases monotonically
with increasing η̃. On the other hand, for q̃ > q̃c , M̃max remains nearly constant as η̃
increases, and the maximum difference in M̃max is less than 0.1. Additionally, except for
the case of q̃ = 0.7 where ω̃ ∗ remains constant, ω̃ ∗ increases monotonically with η̃ beyond
a certain value of η̃.
Fig. 12 depicts the compactness C and the radius R̃ of BBSs as functions of η̃. For q̃ <
q̃c , the compactness C decreases initially and subsequently increases as η̃ increases, whereas
the radius R̃ initially exhibits fluctuations and then undergoes a significant increase overall.
In contrast, for q̃ > q̃c , both C and R̃ demonstrate negligible variation with increasing η̃,
with the compactness values remaining consistently around C ≈ 0.1.

4.4 Light Rings

The geodesic equation for photons in a gravitational field is given by

gµν ẋµ ẋν = 0, (4.5)

where the dot denotes a derivative with respect to an affine parameter along the geodesic
and xµ = (t, r, θ, ϕ).
Due to the static spherical symmetry, we assume photon orbits lie in the equatorial
plane θ = π/2, implying θ̇ = 0. The metric (2.16) possesses two Killing vectors ∂t and
∂φ , corresponding to conserved quantities: the energy E = −g00 ṫ and angular momentum
L = r̃2 φ̇. Equation (4.5) then reduces to
1 ˙2
r̃ + V (r̃) = 0, (4.6)
2

– 14 –
 3.5 q̃ =0.0 q̃ =0.0
q̃ =0.3 0.8 q̃ =0.3
3.0 q̃ =0.5 q̃ =0.5
q̃ =0.63 q̃ =0.63
2.5 0.6 q̃ =0.65
q̃ =0.65
M̃max

q̃ =0.67 q̃ =0.67
2.0

ω̃ ∗
q̃ =0.7 0.4 q̃ =0.7

1.5

1.0 0.2

0.5 0.0
0 1 2 3 4 0 1 2 3 4
η̃ η̃

(a) (b)

Figure 11. (a) M̃max versus η̃ with s̃ = 0.2. (b) ω̃ ∗ versus η̃ with s̃ = 0.2. The dots in the figure
correspond to the calculated values of M̃max and ω̃ ∗ .

12
0.30 q̃ =0.0 q̃ =0.65
q̃ =0.0
q̃ =0.3 q̃ =0.67
q̃ =0.3
q̃ =0.5 10 q̃ =0.5 q̃ =0.7
0.25 q̃ =0.63
q̃ =0.63
q̃ =0.65
0.20 q̃ =0.67 8
R̃
C

q̃ =0.7
0.15
6
0.10
4
0.05
0 1 2 3 4 0 1 2 3 4
η̃ η̃

(a) (b)

Figure 12. (a) Compactness C versus η̃. (b) Radius R̃ versus η̃. The dots in the figure correspond
to the calculated values of C and R̃.

with the effective potential

1 L2 −g00
 
1
V (r̃) = 2
− 2 , (4.7)
2 −g00 g11 r̃ b
where b = L/E is the impact parameter. This formulation describes a classical particle of
unit mass moving in the potential V (r̃) with zero energy [38].
Light rings (LRs) correspond to circular photon orbits satisfying V (r̃) = 0 and V ′ (r̃) =
0, yielding:
−g00
 
1
− 2 = 0, (4.8)
r̃2 b
 g ′
00
− 2 = 0. (4.9)
r̃

– 15 –
We define the reduced effective potential [43, 44, 46]:
g00
Veff (r̃) = − (4.10)
r̃2
to characterize LR positions. Similarly, the radial stability condition V ′′ (r) > 0 can be
′′ (r̃) > 0. The stability of the orbit along the θ-direction has been rigorously
replaced by Veff
proved by Cunha, Berti and Herdeiro [47].
Our calculation reveals LR solutions exist for all parameter combinations:

q̃ ∈ {0.0, 0.1, 0.3, 0.5, 0.6, 0.61, 0.63, 0.65, 0.67, 0.7},

and
η̃ ∈ {0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 3.5, 4.0}.
Fig. 13 displays representative effective potentials for η̃ = 1 with varying q̃ and ω̃.
LR solutions exist only within specific frequency ranges. Fig. 14 illustrates this by shad-
ing LR regions on M̃ -ω̃ curves. For small q̃ and η̃, LRs primarily occupy counterclockwise-
end regions (subfigure 14(a)). As these parameters increase, LR regions expand significantly
(subfigures 14(b)-14(d)).
Comparing frozen BBSs with magnetic RN black holes, the RN effective potential is:

RN ∆(r̃)
Veff (r̃) = − , (4.11)
r̃2
yielding two LRs at:
 
1
q
LR,RN
r̃inner = 3M̃ − 9M̃ 2 − 8q̃ 2 , (4.12)
2
 
1
q
LR,RN 2 2
r̃outer = 3M̃ + 9M̃ − 8q̃ . (4.13)
2
Boson stars also exhibit two light rings (LRs) [31], with the inner and outer radii
LR
denoted by r̃inner LR , respectively. Fig. 15 plots the radii of the LRs of BBSs
and r̃outer
as functions of ω̃ for q̃ = 0.61, 0.63, 0.67, 0.7. As illustrated in Fig. 15, as ω̃ → 0, the
LR
outer radius r̃outer LR does not. This
converges to a fixed value, while the inner radius r̃inner
convergence of the outer radius to a fixed point implies that the position of the outer LRs
of frozen BBSs is independent of η̃.
Fig. 16 illustrates the relationship between the light ring radii of frozen BBSs and RNs
black holes for η̃ = 0.0, 1.0, 2.0, 4.0. For the tested values of q̃, the outer light ring radius
LR , is consistently smaller than that of RN black holes, r̃ LR,RN . Conversely,
of BBSs, r̃outer outer
LR , is larger than that of RN black holes, r̃ LR,RN ,
the inner light ring radius of BBSs, r̃inner inner
LR and r̃ LR increase with
for all tested values of q̃ except q̃ = 0.61. Additionally, both r̃inner outer
increasing q̃.

5 Conclusions

In this work, we have investigated a Bardeen model that couples Einstein’s gravity with
nonlinear electromagnetic fields and non-topological soliton complex scalar fields. When

– 16 –
 3.5 0.20
q =0.0, ω =0.80816 q =0.63, ω =0.355
3.0
q =0.3, ω =0.83976 q =0.67, ω =0.006
2.5 q =0.5, ω =0.46141 0.15 q =0.7, ω =0.00031
2.0
Veff (r)

Veff (r)
0.10
1.5

1.0
0.05
0.5

0.0 0.00
0.0 0.1 0.2 0.3 0.4 0 1 2 3 4
r r

(a) (b)

Figure 13. Effective potentials for photon circular orbits with η̃ = 1 and varying q̃, ω̃.

0.6 q=0.0, η=3.0

q=0.0, η=4.0
3
q=0.3, η=3.0
0.5 q=0.3, η=4.0

3.0
0.4 2
M
M

0.3
2.5
q=0.0, η=0.5 1
q=0.0, η=1.0 0.200 0.225
0.2
q=0.3, η=0.5
q=0.3, η=1.0
0.1 0
0.8 0.9 1.0 0.2 0.4 0.6 0.8 1.0
ω ω

(a) (b)

0.7 2.5 q=0.5, η=3.0

q=0.5, η=4.0
q=0.63, η=3.0
0.6 2.0 q=0.63, η=4.0

1.5
M

0.5 q=0.5, η=0.0

q=0.5, η=0.5
q=0.5, η=1.0 1.0
0.4 q=0.63, η=0.0
q=0.63, η=0.5
0.5
q=0.63, η=1.0
0.3
0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8
ω ω

(c) (d)

Figure 14. M̃ -ω̃ relations with LR regions (gray) for various q̃, η̃.

– 17 –
 1.6 1.6
1.4 LR
η =0.0, r̃inner 1.4 LR
η =0.0, r̃inner
η =0.0, LR
r̃outer LR
η =0.0, r̃outer
1.2
η =2.0, LR
r̃inner 1.2 LR
η =2.0, r̃inner
r̃

r̃
1.0 η =2.0, LR
r̃outer LR
η =2.0, r̃outer
LR
1.0 LR
η =4.0, r̃inner η =4.0, r̃inner
0.8 LR LR
η =4.0, r̃outer 0.8 η =4.0, r̃outer
0.6 q̃=0.61 q̃=0.63
0.6
0.4
0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4
ω ω

(a) (b)

1.8 2.0

1.6 LR
1.8 LR
η =0.0, r̃inner η =0.0, r̃inner
LR
η =0.0, r̃outer LR
η =0.0, r̃outer
1.4 1.6
LR
η =2.0, r̃inner LR
η =2.0, r̃inner
r̃

r̃

1.2
LR
η =2.0, r̃outer 1.4 LR
η =2.0, r̃outer q̃=0.7
η =4.0, LR
r̃inner LR
η =4.0, r̃inner
1.0 LR 1.2 LR
η =4.0, r̃outer η =4.0, r̃outer

0.8
q̃=0.67 1.0

0.8
0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3
ω ω

(c) (d)

LR LR
Figure 15. Inner (r̃inner ) and outer (r̃outer ) LR radii versus ω̃ for BBSs.

the magnetic charge q̃ = 0 and the complex scalar field self-interaction parameter η̃ = 0,
the model reduces to the mini-boson stars model, serving as a baseline for our study. To
explore the effects of magnetic charge and scalar field self-interaction, we systematically
increase the values of q̃ and η̃. We pay particular attention to parameter values leading to
frozen star solutions.
Our findings reveal that the magnetic charge q̃ plays a crucial role in determining the
frozen states of Bardeen boson stars (BBSs), with different critical values for varying η̃.
Notably, the ADM mass M̃ of frozen BBSs is independent of η̃, and these frozen solutions
exhibit a critical horizon whose radius lies between the inner and outer horizon radii of a
magnetic Reissner-Nordström (RN) black hole.
For fixed q̃ and η̃, we have determined the maximum ADM mass of BBSs from the
M̃ -ω̃ curves and computed the corresponding compactness. Notably, for q̃ > q̃c , both the
compactness and the maximum M̃ remain nearly constant with respect to η̃.
We also investigate the existence of light rings for BBSs. Our calculations demonstrate
that light ring solutions exist for all tested combinations of q̃ and η̃, and all frozen stars

– 18 –
 2.0 2.0

LR
r̃inner LR
r̃inner
1.5 LR 1.5 LR
r̃outer r̃outer
r̃

r̃
LR,RN LR,RN
r̃inner r̃inner

1.0 η̃=0.0 LR,RN

r̃outer 1.0 η̃=1.0 LR,RN
r̃outer

0.5 0.5

0.62 0.64 0.66 0.68 0.70 0.62 0.64 0.66 0.68 0.70
q̃ q̃

(a) (b)

2.0 2.0

LR
r̃inner LR
r̃inner
1.5 LR
r̃outer
1.5 LR
r̃outer
r̃

r̃

LR,RN LR,RN
r̃inner r̃inner

1.0 η̃=2.0 LR,RN

r̃outer
1.0 η̃=4.0 LR,RN
r̃outer

0.5 0.5
0.62 0.64 0.66 0.68 0.70 0.62 0.64 0.66 0.68 0.70
q̃ q̃

(c) (d)

Figure 16. Comparison of LR radii between frozen BBSs and RN black holes versus q̃.

possess light rings. The outer light ring radius of frozen stars is independent of η̃. Compared
to magnetic RN black holes, frozen BBSs exhibit a smaller outer light ring radius.
Furthermore, we have examined the existence of light rings (LRs) for BBSs. Our results
demonstrate that LR solutions exist for all tested combinations of q̃ and η̃, and all frozen
BBSs possess LRs. The outer LR radius of frozen stars is independent of η̃, and compared
to magnetic RN black holes, frozen BBSs exhibit a smaller outer LR radius.

6 Acknowledgment

This work is supported by the National Natural Science Foundation of China (Grant No.
12275110 and No. 12247101) and the National Key Research and Development Program of
China (Grant No. 2022YFC2204101 and 2020YFC2201503).

References

[1] LIGO Scientific, Virgo collaboration, B. Abbott, R. Abbott, T. Abbott et al.,

– 19 –
 Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116
(2016) 061102 [1602.03837].
[2] Event Horizon Telescope collaboration, K. Akiyama, A. Alberdi, W. Alef et al., First
M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole,
Astrophys. J. Lett. 875 (2019) L1 [1906.11238].
[3] Event Horizon Telescope collaboration, K. Akiyama, A. Alberdi, W. Alef et al., First
M87 Event Horizon Telescope Results. V. Physical Origin of the Asymmetric Ring,
Astrophys. J. Lett. 875 (2019) L5 [1906.11242].
[4] Event Horizon Telescope collaboration, K. Akiyama, A. Alberdi, W. Alef et al., First
Sagittarius A* Event Horizon Telescope Results. I. The Shadow of the Supermassive Black
Hole in the Center of the Milky Way, Astrophys. J. Lett. 930 (2022) L12 [2311.08680].
[5] S. W. Hawking, Black hole explosions?, Nature 248 (1974) 30.
[6] R. M. Wald, On particle creation by black holes, Commun. Math. Phys. 45 (1975) 9.
[7] S. W. Hawking, Breakdown of predictability in gravitational collapse, Phys. Rev. D 14 (1976)
2460.
[8] R. Penrose, Gravitational Collapse and Space-Time Singularities, Phys. Rev. Lett. 14 (1965)
57.
[9] R. Penrose, Gravitational collapse: The role of general relativity, Riv. Nuovo Cim. 1 (1969)
252.
[10] A. Einstein, On a Stationary System With Spherical Symmetry Consisting of Many
Gravitating Masses, Annals Math. 40 (1939) 922 [1968902].
[11] R. P. Kerr, Do Black Holes have Singularities?, arXiv:2312.00841 [gr-qc] (2023)
[2312.00841].
[12] J. Bardeen, Non-singular general relativistic gravitational collapse, in Proceedings of the 5th
International Conference on Gravitation and the Theory of Relativity, (Tiflis, U.S.S.R.),
Tbilisi University Press, 1968.
[13] K. A. Bronnikov, V. N. Melnikov, G. N. Shikin and K. P. Staniukovich, Scalar,
electromagnetic, and gravitational fields interaction: Particlelike solutions, Annals Phys. 118
(1979) 84.
[14] A. Borde, Open and closed universes, initial singularities, and inflation, Phys. Rev. D 50
(1994) 3692 [gr-qc/9403049].
[15] C. Barrabès and V. P. Frolov, How many new worlds are inside a black hole?, Phys. Rev. D
53 (1996) 3215 [hep-th/9511136].
[16] M. Mars, M. M. Martín-Prats and J. M. M. Senovilla, Models of regular Schwarzschild black
holes satisfying weak energy conditions, Class. Quant. Grav. 13 (1996) L51.
[17] A. Borde, Regular black holes and topology change, Phys. Rev. D 55 (1997) 7615
[gr-qc/9612057].
[18] A. Cabo and E. Ayon, About Black Holes Without Trapping Interior, Int. J. Mod. Phys. A
14 (1999) 2013 [gr-qc/9704073].
[19] S. A. Hayward, Formation and evaporation of non-singular black holes, Phys. Rev. Lett. 96
(2006) 031103 [gr-qc/0506126].

– 20 –
[20] E. Ayón-Beato and A. García, Regular Black Hole in General Relativity Coupled to
Nonlinear Electrodynamics, Phys. Rev. Lett. 80 (1998) 5056 [gr-qc/9911046].
[21] E. Ayón-Beato and A. García, The Bardeen model as a nonlinear magnetic monopole, Phys.
Lett. B 493 (2000) 149 [gr-qc/0009077].
[22] J. A. Wheeler, Geons, Phys. Rev. 97 (1955) 511.
[23] D. J. Kaup, Klein-Gordon Geon, Phys. Rev. 172 (1968) 1331.
[24] X. E. Wang, From Bardeen-boson stars to black holes without event horizon,
arXiv:2305.19057 [gr-qc] (2023) [2305.19057].
[25] Y. Yue and Y.-Q. Wang, Frozen Hayward-boson stars, arXiv:2312.07224 [gr-qc] (2023)
[2312.07224].
[26] J. R. Oppenheimer and H. Snyder, On Continued Gravitational Contraction, Phys. Rev. 56
(1939) 455.
[27] T. D. Lee, Soliton stars and the critical masses of black holes, Phys. Rev. D 35 (1987) 3637.
[28] R. Friedberg, T. D. Lee and Y. Pang, Scalar soliton stars and black holes, Phys. Rev. D 35
(1987) 3658.
[29] B. W. Lynn, Q-stars, Nucl. Phys. B 321 (1989) 465.
[30] P. Jetzer, Boson stars, Phys. Rept. 220 (1992) 163.
[31] C. F. B. Macedo, P. Pani, V. Cardoso and L. C. B. Crispino, Astrophysical signatures of
boson stars: Quasinormal modes and inspiral resonances, Phys. Rev. D 88 (2013) 064046
[1307.4812].
[32] Y. Brihaye, A. Cisterna, B. Hartmann and G. Luchini, From topological to nontopological
solitons: Kinks, domain walls, and $Q$-balls in a scalar field model with a nontrivial vacuum
manifold, Phys. Rev. D 92 (2015) 124061 [1511.02757].
[33] C. F. B. Macedo, Compact Objects in General Relativity and Beyond, Ph.D. thesis, Para U.,
2015.
[34] M. Bošković and E. Barausse, Soliton boson stars, Q-balls and the causal Buchdahl bound,
JCAP 02 (2022) 032 [2111.03870].
[35] L. G. Collodel and D. D. Doneva, Solitonic boson stars: Numerical solutions beyond the
thin-wall approximation, Phys. Rev. D 106 (2022) 084057 [2203.08203].
[36] N. Siemonsen, Nonlinear Treatment of a Black Hole Mimicker Ringdown, Phys. Rev. Lett.
133 (2024) 031401 [2404.14536].
[37] T. Ogawa and H. Ishihara, Gravastars as nontopological solitons, Phys. Rev. D 110 (2024)
124003 [2409.07818].
[38] S. M. Carroll, Spacetime and Geometry: An Introduction to General Relativity. Pearson
Education, Harlow, pearson new international edition ed., 2014.
[39] I. A. Baratta, J. P. Dean, J. S. Dokken, M. Habera, J. S. Hale, C. N. Richardson et al.,
DOLFINx: The next generation FEniCS problem solving environment, Dec., 2023.
10.5281/zenodo.10447666.
[40] M. W. Scroggs, J. S. Dokken, C. N. Richardson and G. N. Wells, Construction of arbitrary
order finite element degree-of-freedom maps on polygonal and polyhedral cell meshes, ACM
Transactions on Mathematical Software 48 (2022) 18:1.

– 21 –
[41] M. W. Scroggs, I. A. Baratta, C. N. Richardson and G. N. Wells, Basix: A runtime finite
element basis evaluation library, Journal of Open Source Software 7 (2022) 3982.
[42] M. S. Alnaes, A. Logg, K. B. Ølgaard, M. E. Rognes and G. N. Wells, Unified Form
Language: A domain-specific language for weak formulations of partial differential equations,
ACM Transactions on Mathematical Software 40 (2014) .
[43] L.-X. Huang, S.-X. Sun and Y.-Q. Wang, Frozen Bardeen-Dirac stars and light ball,
2312.07400.
[44] J.-R. Chen and Y.-Q. Wang, Hayward spacetime with axion scalar field, arXiv:2407.17278
[hep-th] (2024) [2407.17278].
[45] P. A. Cano, L. Machet and C. Myin, Boson stars with nonlinear sigma models, Phys. Rev. D
109 (2024) 044043 [2311.03433].
[46] X.-X. Zeng, C.-Y. Yang, Y.-X. Huang, K.-J. He, G.-P. Li and S. Guo, Optical Images of Mini
Boson Stars in Palatini $f(R)$ Gravity, arXiv:2501.13764 [gr-qc] (2025) [2501.13764].
[47] P. V. P. Cunha, E. Berti and C. A. R. Herdeiro, Light ring stability in ultra-compact objects,
Phys. Rev. Lett. 119 (2017) 251102 [1708.04211].

– 22 –