# Soliton boson stars, Q-balls and the causal Buchdahl bound

Mateja Bošković  SISSA, Via Bonomea 265, 34136 Trieste, Italy and INFN Sezione
di Trieste  IFPU - Institute for Fundamental Physics of the Universe, Via
Beirut 2, 34014 Trieste, Italy  Enrico Barausse  SISSA, Via Bonomea 265, 34136
Trieste, Italy and INFN Sezione di Trieste  IFPU - Institute for Fundamental
Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy

###### Abstract

Self-gravitating non-topological solitons whose potential admits multiple
vacua are promising candidates for exotic compact objects. Such objects can
arise in several extensions of the Standard Model and could be produced in the
early Universe. In this work, we focus on objects made from complex scalars
(gravitating Q-balls/soliton boson stars), deriving analytic solutions in
spherical symmetry and comparing them with fully numerical ones. In the high-
compactness limit we find that these objects present an effectively linear
equation of state, thus saturating the Buchdahl limit with the causality
constraint. Far from that limit, these objects behave either as flat space-
time Q-balls or (in the low-compactness limit) as mini boson stars stabilized
by quantum pressure. We establish the robustness of this picture by analyzing
a variety of potentials (including cosine, quartic and sextic ones).

###### Contents

  1. I Introduction
  2. II Q-balls: a review
    1. II.1 Simplest Q-ball potential: analytic description
  3. III Soliton boson stars
    1. III.1 Scaling arguments
    2. III.2 Structure equations in GR
    3. III.3 Representative configurations
    4. III.4 Origin of maximum compactness
    5. III.5 Analytic construction
      1. III.5.1 Exterior zone
      2. III.5.2 Interior zone
      3. III.5.3 Boundary zone
      4. III.5.4 Energy balance
      5. III.5.5 Finale: Semi-analytic solution
  4. IV Parameter space of Soliton boson stars
    1. IV.1 Planck scale regime
    2. IV.2 Low compactness stable branch
    3. IV.3 Cosine potential
    4. IV.4 General sextic potential
    5. IV.5 Non-polynomial quartic potential
  5. V Conclusion
  6. A Numerical solutions of Einstein-Klein-Gordon system
  7. B Definitions of boson star radius
  8. C Analytic construction of Soliton Boson stars: technical details
    1. C.1 $u_{B}$ metric coefficient on the boundary
    2. C.2 Details of the energy balance calculation
    3. C.3 On the errors of estimating the radius

##  I Introduction

Gravitational wave (GW) Taylor and Weisberg (1982); Abbott  _et al._ (2016,
2017, 2019, 2021a, 2021b) and electromagnetic (EM) Giacconi  _et al._ (1962);
Bowyer  _et al._ (1965); Reid  _et al._ (1999); Ghez  _et al._ (2000); Akiyama
_et al._ (2019a, b) observations provide a unique opportunity to understand
the formation, dynamics and environment of compact objects in the Universe, as
well as to test gravity in the hitherto unexplored strong-field regime.
Besides known compact objects, such as white dwarfs, neutron stars and black
holes (BHs), more exotic compact objects (ECOs) Giudice  _et al._ (2016);
Cardoso and Pani (2019) may also exist as a result of extensions of the
Standard Model, and might represent a fraction of the dark matter (DM) or be
relics from processes in the early Universe Bertone  _et al._ (2020). In
particular, ECOs with compactness $\mathcal{C}\equiv GM/(Rc^{2})>1/3$ ($M$ and
$R$ being respectively the mass and radius of the body), also referred to as
“ultra-compact objects”, may be hard to distinguish from BHs if they present a
light ring Cardoso and Pani (2019). However, toy models for ECOs (such as
wormholes Visser (1995), gravastars Mazur and Mottola (2001), anisotropic
fluid stars Raposo  _et al._ (2019) etc) have point at a rich range of
possible observational features distinguishing these objects from BHs, e.g.
non-vanishing Love numbers Cardoso  _et al._ (2017), distinct post-merger
phase for binary systems Toubiana  _et al._ (2021), GW echoes Barausse  _et
al._ (2014); Cardoso  _et al._ (2016a, b); Cardoso and Pani (2019) etc.

Many ECO toy models, however, are not realistic candidates, as they contradict
some of the following reasonable requirements: stability (on relevant
astrophysical/cosmological scales); existence of possible production
mechanisms and astrophysical/cosmological formation channels; consistency with
known and tested physics; and embedding in a plausible beyond-Standard Model
theory. A more promising way to construct realistic ECO models is to consider
solitons, i.e. localized, finite-energy and stable solutions of the equations
of motion of a field theory. Derrick’s theorem (e.g. Năstase (2019)) is a
powerful constraint on the existence of these solutions, as it prohibits
static non-trivial solutions for a set of real scalar fields in $D\geq 2$,
where $D$ is the number of space dimensions. In order to evade Derrick’s
theorem, one must either consider topologically non-trivial configurations
(topological solitons) or consider a theory with conserved charge (non-
topological solitons) Lee and Pang (1992); Năstase (2019). In this work, we
will consider the latter type of solitons111Examples of topological solitons
in the strong gravity context are cosmic strings Helfer  _et al._ (2019) and
sine-Gordon stars Franzin  _et al._ (2018)..

The simplest examples of non-topological solitons are boson stars, i.e.
configurations made from complex scalars $\Phi$ minimally coupled to gravity
and admitting a $U(1)$ symmetry. Boson stars are regular horizonless solutions
in general relativity. These objects have been thoroughly studied as potential
ECOs Jetzer (1992); Schunck and Mielke (2003); Liebling and Palenzuela (2017);
Cardoso and Pani (2019); Visinelli (2021). Among the different models of boson
stars, the most investigated are mini-boson stars Kaup (1968) (MBSs; which
only have a mass term in the potential), self-interacting boson stars Colpi
_et al._ (1986) (SIBSs; whose potential also includes a repulsive
$\lambda\Phi^{4}$ term) and solitonic boson stars Friedberg  _et al._ (1987a)
(SBSs), which present the potential

$\displaystyle
V=\mu^{2}|\Phi|^{2}\Big{(}1-2\frac{|\Phi|^{2}}{\sigma_{0}^{2}}\Big{)}^{2}\,,$
(1)

where $\sigma_{0}$ is the degenerate vacuum and $\mu$ is a scalar mass222More
complicated examples of boson stars are either made from multiple Bezares and
Palenzuela (2018) or higher-spin fields Brito  _et al._ (2016a), have a bigger
Herdeiro  _et al._ (2019) or gauged symmetry group Liebling and Palenzuela
(2017), or include both bosonic and fermionic components Henriques  _et al._
(1989)..

These models describe three distinct physical mechanisms of stabilizing the
configurations (see Sec. III.1). In MBSs/SIBSs, that mechanism is provided by
the equilibrium between gravity and respectively quantum pressure/repulsive
pressure. For SBSs, which have a bubble-like structure in the most compact
part of the parameter space, we will see that accumulation of the energy near
the surface gives rise to a surface tension between different vacua, which
allows for developing massive and highly compact configurations. Because of
this, SBSs are the only model among the aforementioned three that exists also
in the $G\to 0$ $(M_{\rm Pl}\to\infty)$ limit. (In that limit, SBSs are also
referred to as Q-balls Coleman (1985)). The three models mentioned above
should be considered as an illustration of the different mechanisms that can
bind boson stars together. Models more motivated from the particle physics
point of view would follow from renormalizable Freitas  _et al._ (2021) or
effective field theories for the scalar Krippendorf  _et al._ (2018); Barranco
_et al._ (2021). However, we expect that the macroscopic behaviour of the
resulting boson stars should still fall, at least qualitatively, in one of the
three classes mentioned above, as we will argue here in the case of SBSs.

Boson stars generally satisfy the soundness checks that we mentioned above.
Non-rotating boson stars have at least one stable branch in the low-
compactness (MBS) limit Seidel and Suen (1990); Liebling and Palenzuela
(2017). Although rotation can lead to instabilities Di Giovanni  _et al._
(2020), it has been recently shown that sufficiently strong self-interactions
can sustain it Siemonsen and East (2021); Dmitriev  _et al._ (2021). As these
objects are constructed from covariant Lorentz-invariant actions, they are
consistent with causality from the start (something often not clear with ad
hoc ECO models), and the modelling of their dynamics is in principle
accessible. Boson stars can form from scalar collapse, with the excess field
‘evaporating’ (gravitational cooling) Seidel and Suen (1994). Once a
population of non-rotating boson stars is formed, it can reproduce in binaries
through the channel $\rm{BS}+\rm{BS}\to\rm{BS}$ Palenzuela  _et al._ (2017).
In this work we will focus on SBSs, as they stand out as the most compact
ones, with $\mathcal{C}_{\rm max}\approx 0.33$ Cardoso  _et al._ (2016b) or
even $\mathcal{C}_{\rm max}\approx 0.348$ Friedberg  _et al._ (1987a).

SBSs are the simplest representatives of non-topological solitons that exist
in the $M_{\rm Pl}\to\infty$ limit and admit a false or degenerate vacuum in
one of the bosonic degrees of freedom333The study of these objects was
initiated by T. D. Lee and collaborators in the 1980s Lee (1987a); Friedberg
_et al._ (1987b, a); Lee and Pang (1987); Lee (1987b); Lee and Pang (1992)..
Thus forth, we will focus on this specific model, in order to get an insight
on the general features of the whole class of objects and potentially on those
of more realistic and complicated models that we plan to address in the future
Lee and Pang (1987); Bahcall  _et al._ (1989); Hong  _et al._ (2020); Gross
_et al._ (2021); Bošković  _et al._ .

Besides simplicity, there are also at least two additional reasons to be
interested in this kind of objects. First, Q-balls/SBSs can be formed in the
early Universe, through scalar fragmentation after inflation Cotner  _et al._
(2019); Amin and Mocz (2019), thermal phase transitions Krylov  _et al._
(2013) or solitosynthesis Postma (2002); Croon  _et al._ (2020). Therefore,
the possible observation of SBSs or similar ECOs can shed light on beyond-
Standard Model physics, from baryogenesis to dark matter Rubakov and Gorbunov
(2017). The charge of these objects could be protected by an approximate, low-
energy global symmetry Krippendorf  _et al._ (2018), and depending on the
specific model a fraction of the compact-object population of the universe
could be in this form Troitsky (2016). Note that even if SBSs are low-
compactness at formation, their subsequent interaction and merger could lead
to the formation of more compact configurations.

SBSs could also be considered as a proxy for similar self-gravitating
structures made from (real) bosons unprotected by a symmetry. Among the latter
are axion stars Helfer  _et al._ (2017), but also oscillatons Nazari  _et al._
(2021), moduli stars Krippendorf  _et al._ (2018); Muia  _et al._ (2019) etc.
These objects are not strictly speaking solitons, and both the energy-momentum
tensor and metric are time dependent Seidel and Suen (1991); Helfer  _et al._
(2017); Bošković  _et al._ (2018). However, owing to the large number of
particles, their decay is exponentially suppressed Kasuya  _et al._ (2003);
Hertzberg (2010); Zhang  _et al._ (2020). Notwithstanding their pseudo-
solitonic nature, these objects have macroscopic properties qualitatively
similar to the corresponding boson star models. For example, the $M$-$R$ curve
of scalar stars with a mass-term potential (oscillaton) differs only by (at
most) a few percent from the corresponding MBS model, and only in the most
compact branch Brito  _et al._ (2016b). Axion stars and similar objects are
expected to form in a wider variety of settings - through cosmological
evolution of axion dark matter Marsh (2016); Hui (2021); Schive  _et al._
(2014); Levkov  _et al._ (2018); Widdicombe  _et al._ (2018); Arvanitaki  _et
al._ (2020), as inflation relics Amin  _et al._ (2014); Niemeyer and Easther
(2020); Arvanitaki  _et al._ (2020) or through other early Universe processes
similar to those giving rise to complex configurations Krippendorf  _et al._
(2018); Cotner  _et al._ (2019); Amin and Mocz (2019). Similar objects can
also arise in certain scalar-tensor theories Barranco  _et al._ (2021).
Although an exact parallel between boson star models and pseudo-solitonic
configurations warrants a more in-depth study, in particular when a
false/degenerate vacuum exists, the static nature of boson star spacetimes
allows for a technically easier preliminary analysis of these objects.

Besides the rich literature on Q-balls (starting from the seminal paper by
Coleman Coleman (1985); see also Nugaev and Shkerin (2020) for a recent
review), the structure of relativistic SBSs has been investigated for the
simplest potential (1) Friedberg  _et al._ (1987a); Kesden  _et al._ (2005);
Macedo  _et al._ (2013); Cardoso  _et al._ (2016b); Siemonsen and East (2021),
for more general sextic potentials Kleihaus  _et al._ (2005, 2012); Siemonsen
and East (2021) and for related potentials, such as the cosine one Guerra  _et
al._ (2019); Siemonsen and East (2021). Quasi-normal modes and geodesics
around SBSs have also been studied Macedo  _et al._ (2013), their Love numbers
have been calculated Cardoso  _et al._ (2017); Sennett  _et al._ (2017), and
they have been simulated in binary systems Cardoso  _et al._ (2016b); Bezares
_et al._ (2017); Palenzuela  _et al._ (2017); Bezares and Palenzuela (2018);
Helfer  _et al._ (2021); Bezares  _et al._ (2022).

The relation between Q-balls and SBSs has been at least partially discussed
previously Multamaki and Vilja (2002); Tamaki and Sakai (2010, 2011a, 2011b)
and the most dramatic effects of gravity are non-perturbative in
$\sigma_{0}/M_{\rm Pl}$. In this regime, as $M_{\rm Pl}<\infty$, Q-ball limit
is not realized. However, there is no major qualitative difference between
Q-balls and SBSs in the part of the parameter space where the effects of
gravity are not important and perturbative. For clarity, in the rest of the
paper we will refer to the $M_{\rm Pl}\to\infty$ limit as Q-balls and to the
gravitating case as SBSs.

Notwithstanding all these studies, to the best of our knowledge some basic
aspects of SBSs and non-topological solitons that admit degenerate vacua in
general have not been appropriately addressed in the literature. In this work
we aim to fill this gap by constructing analytic descriptions of SBSs for the
simplest potential (1) in the sub-Planck limit $\sigma_{0}/M_{\rm Pl}\ll 1$,
by understanding the physics that sets the maximal compactness of these
objects, and by illuminating the non-trivial connection with the Q-ball limit.
We will establish the robustness and the limits of this picture by exploring
the parameter space in $\sigma_{0}/M_{\rm Pl}$, and also by considering
cosine, sextic and quartic potentials. We will show that diverse choices of
the potential do not correspond to dramatically distinct macroscopic
behaviours of SBSs, and that the structure of SBSs mostly depends on the
distance between the central field and the scalar false vacuum, and on whether
the false vacuum is sufficiently deep.

This paper is organized as follows. In Sec. II we will provide a brief review
of Q-balls, upon which we build up description of SBS structure and properties
in Sec. III. In Sec. IV we will explore the parameter space of these objects
and the effect of the scalar potential on their structure. The technical
details of our numerical methods are presented in App. A, while in App. B we
discuss various definitions of the SBS radii and in App. C we provide some
additional details on the analytic construction of SBSs. Throughout this
paper, we will employ a metric signature $-+++$ and natural units $c=\hbar=1$,
with $M^{2}_{\rm Pl}=1/8\pi G$ and $m^{2}_{\rm Pl}=1/G$. Note that in
geometric units $(G=c=1)$ $M^{2}_{\rm Pl}/\hbar=1/8\pi$ and $[\hbar]=[M]^{2}$.

##  II Q-balls: a review

The Lagrangian of a scalar field with $U(1)$ symmetry in flat space-time is
given by

$\mathcal{L}_{\Phi}=-\partial_{\mu}\Phi^{\dagger}\partial^{\mu}\Phi-V(|\Phi|)\,.$
(2)

The scalar field can be decomposed in Fourier modes as

$\Phi=\phi(r)e^{-i\omega t}\,,$ (3)

with time translations corresponding to a change in phase. Because of the
$U(1)$ symmetry, and as long as only one Fourier mode is excited, this
superficial time dependence does not propagate to any observable quantity,
such as thermodynamics parameters: the density $\rho=-T^{t}_{t}$ , the radial
pressure $P_{\rm rad}=T^{r}_{r}$ etc., and Derrick’s theorem is circumvented.

Varying the action (2) with respect to $\phi$, one obtains the Klein-Gordon
equation

$\displaystyle\phi^{\prime\prime}+\frac{2}{r}\phi^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{dU_{\omega}}{d\phi}\,,$ (4) $\displaystyle U_{\omega}$
$\displaystyle=$ $\displaystyle\frac{1}{2}(\omega^{2}\phi^{2}-V(\phi))\,.$ (5)

Note that this equation can be interpreted as the equation of motion of a
Newtonian particle under a friction term and an effective potential
$U_{\omega}$. The energy of the Newtonian particle is

$\displaystyle\mathcal{E}=\frac{1}{2}(\phi^{\prime})^{2}+U_{\omega}\,,$ (6)

and it is conserved only in $D=1$ dimensions, when the friction is absent.

In order for the field energy (i.e. the integrated energy density) $E=\int
dV\rho$ to be finite, one must have $\phi(r\to\infty)\to 0$, which implies
that the Q-ball is a localized object. Let us assume that at infinity the
scalar is free, with leading mass term in the potential
$V=\frac{1}{2}\mu^{2}\phi^{2}+\mathcal{O}(\phi^{4})$. As
$U_{\omega}\sim(\omega^{2}-\mu^{2})\phi^{2}$ , $\phi\to 0$, we require
$\omega<\mu$ for the scalar to converge to the vacuum state at infinity. The
leading order behaviour is of the form
$\phi\sim\exp{(-\sqrt{\mu^{2}-\omega^{2}}r)}$, which implies zero “velocity”
$\phi^{\prime}$ at infinity and $\mathcal{E}(\infty)=0$. The asymptotic
Yukawa-like behaviour implies that the Q-ball radius is ill defined. We follow
the common practice (e.g. Macedo  _et al._ (2013)) and arbitrarily define the
radius as that enclosing 99$\%$ of the total mass $R_{99}$. (Note however that
different conventions and definitions can be found in the literature, see App.
B).

At the initial point (in the particle perspective), in order for the friction
to be overcome one must have $\phi^{\prime}(r\to 0)\propto
r^{1+\epsilon}\,,\epsilon>0$ i.e. $\phi^{\prime}(0)=0$. As $U_{\omega}(0)=0$,
the particle needs to be released at some point
$\phi_{c}|U_{\omega}(\phi_{c})>0\,,U_{\omega}^{\prime}(\phi_{c})>0$ in order
to overcome the friction and arrive at infinity with finite energy. The last
condition implies that $U_{\omega}$ must admit an additional hill, or more
formally there must be a non-trivial minimum of $V/\phi^{2}$ at some point
$\phi_{0}\neq 0$ and $\omega\geq\omega_{0}\equiv\text{min}[V/\phi^{2}]$
Coleman (1985). The simplest potential is thus of the form

$\displaystyle
V=\mu^{2}|\Phi|^{2}\Big{(}1-2\frac{|\Phi|^{2}}{\sigma_{0}^{2}}\Big{)}^{2}\,,$
(7)

with $\phi_{0}=\sigma_{0}/\sqrt{2}$. A theory with this potential is non-
renormalizable, and thus must interpreted in the context of effective or
thermal field theories444If one considers non-Abelian groups, instead,
renormalizable potentials can also admit stable Q-balls Safian  _et al._
(1988)..

A sketch of the corresponding $U_{\omega}$ is shown in Fig. 1. In the limit
$\omega=\omega_{0}$, the fictitious particle is at the second hill
$\phi_{c}\to\phi_{0}$, and the scalar would need an infinite amount of time to
travel to the first hill (trivial vacuum). This scenario, in the Q-ball
perspective, corresponds to an infinitely large configuration. The impact of
the friction in this case can be neglected, as it is proportional to $1/R\to
0$. For a small but non-zero value of $\omega$ , the traveling time/Q-ball
radius becomes finite, but still large (thin-wall regime). In the opposite
case when $\omega\to\mu$, the impact of friction is pronounced and the
transition to the origin is smoother (thick-wall regime). The particle picture
also demonstrates how one can construct these solutions numerically - through
a shooting algorithm (see Sec. III.2 and App. A for the numerical
formulation). If one releases the particle on the left of $\phi_{c}$ (e.g.
somewhere around the valley for $\omega/\mu=0.5$), it will not have enough
energy to reach the first hill. However, releasing the particle sufficiently
near the second hill will let it reach the trivial vacuum with an excess of
energy. In this way, one can iteratively find the true solution.

From the invariance under internal $U(1)$ symmetry, the Noether current and
charge can be found as

$\displaystyle j^{\mu}$ $\displaystyle=$
$\displaystyle-i(\Phi^{\dagger}\partial^{\mu}\Phi-\Phi\partial^{\mu}\Phi^{\dagger})\,,$
(8) $\displaystyle Q$ $\displaystyle\equiv$ $\displaystyle-\int
dVj^{t}=4\pi\omega\int dr\,r^{2}\phi^{2}\,.$ (9)

In order for the Q-ball not to decay to constituent scalars, the energy per
unit charge should be lower than the particle mass:

$\displaystyle E<\mu Q\,.$ (10)

In Coleman (1985), Coleman proved the stability of Q-balls, under the above
restrictions on the form of the potential in classical field theory (necessary
condition) and the requirement (10) . Other aspects of stability are discussed
in Paccetti Correia and Schmidt (2001); Sakai and Sasaki (2008).

To obtain a rough understanding of the Q-ball properties, let us consider the
potential (1). When the central field is near the degenerate vacuum
$\sigma_{0}$ (thin-wall regime), the main contribution to the mass of the
object comes from the bulk $\sigma^{2}_{0}R^{3}$ and the surface tension
$\delta^{-1}\sigma^{2}_{0}R^{2}$, where $\delta\sim\mu^{-1}$ is the size of
the potential wall. In equilibrium one has [with
$(\Phi\sim\sigma_{0}\exp{i\omega t})$]

$R\sim\frac{\mu}{\omega^{2}}\,.$ (11)

This rough argument provides the leading order behaviour of Q-balls in the
thin-wall limit.

Figure 1: (Left) Effective particle potential $U_{\omega}$ for various values
of the field frequency $\omega$. The dots correspond to initial field values
for physical configurations and all curves represent the potential (1), except
for the purple one which represents the quadratic potential
$V=\mu^{2}|\Phi|^{2}$ (that doesn’t allow for solitonic solutions in flat
space-time). (Right) Radial field profile (numerical results) for several
physical configurations. Note that the $\omega/\mu=0$ case corresponds to the
infinitely large Q-ball with scalar $\phi=\sigma_{0}/\sqrt{2}$.

###  II.1 Simplest Q-ball potential: analytic description

We will now focus on the simplest potential with a global $U(1)$ symmetry and
which admits stable Q-balls (1). For this potential, we can use the following
scaling:

$\mathsf{r}=\mu
r\,,w=\frac{\omega}{\mu}\,,\varphi=\phi/(\sigma_{0}/\sqrt{2})\,,$ (12)

which makes the equations of motion dimensionless.

In Heeck  _et al._ (2021a), analytic profiles for Q-balls have been
constructed (building on Paccetti Correia and Schmidt (2001); Tsumagari  _et
al._ (2008)) by matching solutions in three regimes – interior (perturbative
solution around $\varphi_{c}$), boundary (expansion around the radius) and
exterior (asymptotic). Here, we only mention these results, as we will
partially review them (together with their curved space-time generalization)
in the next section:

$\displaystyle\varphi_{<}$ $\displaystyle\approx$
$\displaystyle\varphi_{+}\Big{(}1-c_{<}\frac{\sinh(\alpha\mathsf{r})}{\mathsf{r}}\Big{)}\,,\,$
(13) $\displaystyle\varphi_{B}$ $\displaystyle\approx$
$\displaystyle\varphi_{+}\frac{1}{\sqrt{1+2e^{2(\mathsf{r}-\mathsf{R^{\ast}})}}}\,,$
(14) $\displaystyle\varphi_{>}$ $\displaystyle\approx$
$\displaystyle\frac{\varphi_{+}c_{>}}{\mathsf{r}}e^{-\sqrt{1-w^{2}}\mathsf{r}}\,,$
(15)

where

$\displaystyle\alpha^{2}=\frac{4}{3}\big{(}1+3w^{2}+2\sqrt{1+3w^{2}}\big{)}\,,$
$\displaystyle\varphi^{2}_{+}=\frac{1}{3}(2+\sqrt{1+3w^{2}})\,.$ (16)

Here, $\varphi_{+}$ corresponds to the maximum of $U_{\omega}$. For $w\ll 1$,
the maximum of $U_{\omega}$ lies close to the degenerate vacuum
$\varphi\approx 1$; $\mathsf{R^{\ast}}=\mu R^{\ast}$ is the inflection point
of the field $\varphi^{\prime\prime}(\mathsf{R^{\ast}})=0$, taken as the
estimate of the size of the Q-ball. We confirm the observation from Heeck  _et
al._ (2021a) that the function (14) does a good job at describing the
numerical profile of the Q-ball, even outside its range of validity, in
contrast with the interior and asymptotic approximants, which only work
sufficiently close/far to the centre, respectively.

By matching the fields (13)-(15) and their derivatives at
$\\{\mathsf{r_{<}},\mathsf{r_{>}}\\}$, one finds Heeck  _et al._ (2021a):

$\displaystyle\,c_{<}\approx\mathsf{R}^{\ast}e^{-2\mathsf{R}^{\ast}}\,,\,c_{>}\approx\sqrt{2}\mathsf{R^{\ast}}e^{\mathsf{R^{\ast}}}\,.$
(17)

However, this matching is not sufficient to close the system. One needs also
the energy balance condition

$\mathcal{E}(\infty)-\mathcal{E}(0)=-\int^{\infty}_{0}dr\big{[}\frac{2}{r}(\phi^{\prime})^{2}\big{]}\,.$
(18)

In summary, we have five equations in five unknowns
$\\{c_{<},c_{>},\mathsf{R^{\ast}},\mathsf{r_{<}},\mathsf{r_{>}}\\}$ as
functions of $\omega$, or equivalently
$\varphi_{c}(\omega)=\varphi_{+}(1-c_{<})$. From Eq. (17) we see that in the
thin-wall regime (large Q-balls), the deviation of the field from the central
value in the interior zone is exponentially suppressed. In the outer zone, the
field itself is exponentially suppressed, and thus the integral on the right
hand side of Eq. (18) is dominated by the boundary zone.

From invariance under time reversal, we expect $R^{\ast}$ to be a Laurent
polynomial in even powers of $\omega$. From Eq. (18) and introducing
$z=\mathsf{r}-\mathsf{R^{\ast}}$, we then have

$\displaystyle\left(e^{-2\mathsf{R^{\ast}}}\mathsf{R^{\ast}}-1\right)^{2}(2+\sqrt{1+3w^{2}})\times$
$\displaystyle\Big{[}w^{2}\mathsf{R^{\ast}}-\mathsf{R^{\ast}}\Big{(}1-\frac{1}{3}\big{(}e^{-2\mathsf{R^{\ast}}}\mathsf{R^{\ast}}-1\big{)}^{2}(\sqrt{3w^{2}+1}+2)\Big{)}^{2}\Big{]}$
$\displaystyle=4\varphi_{+}\int^{+\infty}_{-\infty}\frac{dz}{1+\frac{z}{\mathsf{R^{\ast}}}}\frac{4e^{4z}}{(1+2e^{2z})^{3}}\,.$
(19)

If we take the $\mathsf{R^{\ast}}\to\infty,\omega\to 0$ limits, the integral
on the right hand side is convergent. In order for the left hand side (i.e.
the term in square brackets) to be finite, we must have

$R\omega^{2}=1+\sum^{\infty}_{n=0}c_{2n}\omega^{2n+2}\,,$ (20)

which confirms the preliminary expectation given by formula (11), i.e.
$R^{\ast}\propto\omega^{-2},\omega\to 0$. By expanding the denominator for
$z\ll\mathsf{R^{\ast}}$ and ignoring exponentially suppressed terms, we find

$c_{0}=\frac{1}{4}(2\log 2-1)\,,\,c_{2}=\frac{1}{48}\left(4\pi^{2}-27\right)$
(21)

for the leading order behaviour.

In order to connect $R^{\ast}$ to the mass-based radius $R\equiv R_{99}$, we
must numerically solve $4\pi\int^{R}_{0}dr\,r^{2}\rho(r)=0.99M$ for $R$. Here,
we construct an analytic approximate solution by considering the density
support of the Q-ball. By construction of the profile,
$\mathsf{R}-\mathsf{R^{\ast}}\equiv\lambda>0$ and we will assume that
$\lambda\equiv\mathsf{R}-\mathsf{R^{\ast}}$ corresponds to the tail of the
boundary zone555The most compact regime corresponds to $\omega\to 0$ when
$\lambda\to 0+$. All other, physically reasonable, profiles would present a
smoother decay of the scalar i.e. fatter tails. Thus, they would reach the
inflection point well before most of $99\%$ of the energy density has been
accumulated. We have checked this explicitly, for all of the models considered
in this work, both with and without gravity.. The width of the boundary zone
can be found by taking the derivative of the field (14) and computing the
standard deviation of that symmetric function, which we generalize by changing
$2z\to\delta z$ in the exponent, having in mind the inclusion of gravity in
Sec. III. This yields

$\lambda=3\frac{1}{2}\sqrt{-\int^{+\infty}_{-\infty}z^{2}d\Big{(}\frac{1}{\sqrt{1+2e^{\delta
z}}}\Big{)}}=\frac{3}{2}\frac{2.66}{\delta}\,.$ (22)

The prefactor arises from the fact that
$\int^{+3\delta}_{-3\delta}d(1+2e^{\delta z})^{-1/2}\approx
0.99\int^{+\infty}_{-\infty}d(1+2e^{\delta z})^{-1/2}$. For Q-balls (in
Minkowski space) $\delta=2$ and

$\mathsf{R}\approx\mathsf{R^{\ast}}+\lambda\,,\lambda=2.66.$ (23)

From the solution, other macroscopic parameters can also be calculated, such
as the energy (mass) $E=M$ and the Noether charge $Q$ Heeck  _et al._ (2021a)
:

$\displaystyle\bar{Q}$ $\displaystyle\equiv$
$\displaystyle\frac{\mu^{2}Q}{m^{2}_{\rm Pl}}$ (24) $\displaystyle\approx$
$\displaystyle\frac{\Lambda^{2}}{6}(\mathsf{R^{\ast}})^{3}w\Big{[}1-\frac{3\ln{2}}{2\mathsf{R^{\ast}}}+\mathcal{O}{\big{(}(\mathsf{R^{\ast}})^{-2}\big{)}}\Big{]}\,,$
$\displaystyle\bar{M}$ $\displaystyle\equiv$ $\displaystyle\frac{\mu
M}{m^{2}_{\rm Pl}}$ (25) $\displaystyle\approx$ $\displaystyle
w\bar{Q}+\frac{\Lambda^{2}}{24}(\mathsf{R^{\ast}})^{2}\Big{[}1+\frac{1-\ln{2}}{\mathsf{R^{\ast}}}+\mathcal{O}{\big{(}(\mathsf{R^{\ast}})^{-2}\big{)}}\Big{]},$

where we have introduced

$\Lambda=\frac{\sigma_{0}}{M_{\rm Pl}}\,,$ (26)

anticipating the connection with the curved space-time case in the next
section. Note that the previous relations imply the following scaling of the
compactness in the thin-wall regime

$\mathcal{C}\propto\Lambda^{2}\frac{1}{\omega^{2}}\,.$ (27)

In Fig. 2 we show the $\mathcal{C}-M$ curve from both the numerical and
analytic calculations. From the plots, we also see that the analytic
approximation for $R$ agrees well with the numerically calculated value. The
strong gravity regime is expected to be attained for values of
$\mu\,,M\,,\sigma_{0}$ such that
$\mathcal{C}\,\mbox{{\scriptsize$\stackrel{{\scriptstyle>}}{{\sim}}$}}\,0.1$:

$\displaystyle\Big{(}\frac{\mathcal{{C}}}{0.1}\Big{)}$ $\displaystyle\approx
1.6\times\Big{(}\frac{\mu}{10^{-8}\rm{eV}}\Big{)}^{1/2}\Big{(}\frac{M}{1M_{\odot}}\Big{)}^{1/2}\Big{(}\frac{\sigma_{0}}{10^{17}\rm{GeV}}\Big{)}^{-1}$
$\displaystyle\approx
1.6\times\Big{(}\frac{\mu}{10^{-14}\rm{eV}}\Big{)}^{1/2}\Big{(}\frac{M}{10^{6}M_{\odot}}\Big{)}^{1/2}\Big{(}\frac{\sigma_{0}}{10^{17}\rm{GeV}}\Big{)}^{-1}\,.$

The threshold frequency between the stable and the unstable branch can be
found from Eqs. (10), (24) and (25), and is given by Heeck  _et al._ (2021a)

$w<w_{\rm{Qb-s}}\approx 0.82\,.$ (29)

Figure 2: Mass-compactness relations for Q-balls. Dots correspond to numerical
configurations, while the solid line is found from the analytic expressions
(23), (25).

##  III Soliton boson stars

Turning gravity on introduces a new scale $\Lambda$ in the problem, with both
perturbative and non-perturbative effects. In the low-compactness limit, SBSs
effectively “see” a quadratic potential and are thus stabilized by quantum
pressure (Sections III.1 and IV.2). In the high-compactness limit, Eq. (27)
shows that for any given $\Lambda$ there exists a sufficiently low $\omega$
such that at some point the Schwarzschild compactness will be reached, and the
scalar will collapse to a BH, simply because of the hoop conjecture Thorne
(1972); Cardoso and Pani (2019). Thus, from the one stable and the one
unstable branch in the flat space-time limit, we expect two stable and two
unstable branches for any value of $\Lambda\ll 1$. This has already been
established numerically and using catastrophe theory arguments Tamaki and
Sakai (2011a); Kleihaus  _et al._ (2012). In the following, we will confirm
previous results and supplement them with a physical interpretation and also
an analytic model, focusing on the compact stable branch in the perturbative
regime $\Lambda\ll 1$.

###  III.1 Scaling arguments

The structure of macroscopic objects is determined by the physics that
stabilizes them. In this section, partially inspired by Lee (1987b); Burrows
and Ostriker (2014); Freivogel  _et al._ (2020), we will provide rough scaling
arguments to demonstrate what kind of configuration properties are to be
expected.

A polytropic equation of state $P\sim\kappa\rho_{\rm m}^{\gamma}\,$,
$[\kappa]=[M]^{-n}$ provides pressure support counteracting the attractive
gravitational force. In equilibrium, one has

$\displaystyle M^{2-\gamma}\sim\kappa m^{2}_{\rm Pl}R^{4-3\gamma}\,.$ (30)

If these configurations can reach the Schwarzshild compactness scale $R_{\rm
Sch}\sim M/m^{2}_{\rm Pl}$, where strong-gravity effects are important, the
maximum Chandrasekhar mass is given by

$\displaystyle M_{\rm Ch}^{2\gamma-2}\sim\kappa m^{6\gamma-6}_{\rm Pl}\,.$
(31)

SIBSs are pressure-supported by the repulsive interaction
$P\sim\lambda|\Phi|^{4}$ for sufficiently large values of $\lambda>0$. As the
matter density is given by $\rho_{\rm m}\sim\mu^{2}|\Phi|^{2}$ (with $\mu$ the
scalar mass), the equation of state is $P\sim\frac{\lambda}{\mu^{4}}\rho_{\rm
m}^{2}$. SIBSs thus have the same mass-radius scaling as fermionic pressure-
supported objects [scaling (30)], and the maximum mass is given by the
Chandrasekhar limit

$\displaystyle M_{\rm Ch}\sim\sqrt{\lambda}\frac{m^{3}_{\rm Pl}}{\mu^{2}}\,.$
(32)

Pressure-supported nature of the SIBS is also the reason why the maximal
attainable compactness for this model $\mathcal{C}_{\rm max}[\rm{SIBS}]\approx
0.16$ Amaro-Seoane  _et al._ (2010) is very close to the neutron star value
$\mathcal{C}_{\rm max}[\rm{NS}]\approx 0.19$ Chen and Piekarewicz (2015).

The stability of MBSs stems from the interplay of the kinetic term in the
Lagrangian (“quantum pressure”) and gravity. Microscopically, quantum pressure
originates from Heisenberg’s uncertainty principle. Let $R$ be the
characteristic size of a MBS and $v_{\text{vir}}\sim\sqrt{M/(m^{2}_{\rm
Pl}R)}$ its virialized velocity. From the uncertainty principle, one has

$\displaystyle\mu M\sim\frac{m^{2}_{\rm Pl}}{\mu R}.$ (33)

The maximum mass, when $R\sim R_{\rm Sch}$, corresponds to the Kaup limit Kaup
(1968)

$\displaystyle M_{\rm Kp}\sim\frac{m^{2}_{\rm Pl}}{\mu}.$ (34)

Note that the difference between the Kaup and Chandrasekhar scaling originates
from the fact that the quantum pressure is not a polytropic radial pressure,
but rather an anisotropic non-local stress Chavanis and Harko (2012).

In contrast to MBSs and SIBSs, SBSs are stabilized already in flat space-time,
as discussed in Sec. II. We can estimate the maximum mass by equating
expression (11) with the Schwartzshild scale, which gives the Lee limit
Friedberg  _et al._ (1987a); Lee (1987b)

$\displaystyle M_{\rm Lee}\sim\frac{m^{4}_{\rm Pl}}{\sigma^{2}_{0}\mu}.$ (35)

###  III.2 Structure equations in GR

The action (2) can be generalized to curved space-time via a minimal coupling
to gravity:

$S=\int d^{4}x\sqrt{-g}(M^{2}_{\rm Pl}R+\mathcal{L}_{\Phi})\,.$ (36)

We will consider spherically symmetric and static space-times of the form

$\displaystyle ds^{2}=-e^{v}dt^{2}+e^{u}dr^{2}+r^{2}d\Omega^{2}\,.$ (37)

The radial profile of the metric coefficients follows from the Einstein field
equations, which together, with the Klein-Gordon equation for the scalar,
describe the full structure of the object:

$\displaystyle\frac{1}{r^{2}}\left(r\,e^{-u}\right)^{\prime}-\frac{1}{r^{2}}$
$\displaystyle=$ $\displaystyle-\frac{1}{M^{2}_{\rm Pl}}\rho\,,$ (38)
$\displaystyle
e^{-u}\left(\frac{v^{\prime}}{r}+\frac{1}{r^{2}}\right)-\frac{1}{r^{2}}$
$\displaystyle=$ $\displaystyle\frac{1}{M^{2}_{\rm Pl}}P_{\text{rad}}\,,$ (39)
$\displaystyle\phi^{\prime\prime}+\left(\frac{2}{r}+\frac{v^{\prime}-u^{\prime}}{2}\right)\phi^{\prime}$
$\displaystyle=$ $\displaystyle
e^{u}\left({\frac{dV}{d|\Phi|^{2}}}-\omega^{2}e^{-v}\right)\phi,$ (40)

with the density $\rho$, the radial pressure $P_{\rm rad}$ and the tangential
pressure $P_{\rm tan}$ defined by

$\displaystyle\rho=e^{-v}\omega^{2}\phi^{2}+e^{-u}(\phi^{\prime})^{2}+V\,,$
(41) $\displaystyle P_{\rm
rad}=e^{-v}\omega^{2}\phi^{2}+e^{-u}(\phi^{\prime})^{2}-V\,,$ (42)
$\displaystyle P_{\rm
tan}=e^{-v}\omega^{2}\phi^{2}-e^{-u}(\phi^{\prime})^{2}-V\,.$ (43)

We will now focus on the simplest potential (1) and adopt the rescaling

$\mathsf{r}=\mu r\,,\bar{m}(\mathsf{r})=\frac{\mu m(r)}{m^{2}_{\rm
Pl}}\,,w=\frac{\omega}{\mu}\,,\varphi=\phi/(\sigma_{0}/\sqrt{2})$ (44)

[c.f. relations (12)]. With this parameterization, the structure equations
become dimensionless. We further assume that ${}^{\prime}\equiv d/d\mathsf{r}$
unless we state otherwise.

Note that a shift $v\to\tilde{v}=v-v(0)$ corresponds to a rescaling of the
time coordinate and thus to a redefinition of the scalar frequency
$\tilde{\omega}=\exp{(-v(0)/2)}\omega$. Fixing $\varphi(0)\equiv\varphi_{c}$
and $\varphi(\infty)=0$, let us use this gauge freedom to also set
$\tilde{v}(0)=0$. This specifies666We are here focusing on the ground state
(nodeless) solutions. It has been previously established that the excited
solutions in various boson star models decay to the ground state Balakrishna
_et al._ (1998); Collodel  _et al._ (2017). However, there are recent
indications that sufficiently strong self-interactions can stabilise excited
states over long time scales at least for SIBSs Sanchis-Gual  _et al._ (2021).
a boundary value problem, with eigenvalue $\tilde{\omega}$, which we determine
with a shooting method. Once that the configuration is calculated in this
gauge, we can rescale the time coordinate at will, so that $v(\infty)\equiv
v_{\infty}=0$. More details on the numerical method are given in the App. A.

The Arnowitt-Deser-Misner (ADM) mass Carroll (2004) can be extracted from the
asymptotic behaviour of the metric functions:

$\exp{[-u(\mathsf{r})]}=1-\frac{2\bar{M}}{\mathsf{r}}+\mathcal{O}\Big{(}\frac{1}{\mathsf{r}^{2}}\Big{)}\,.$
(45)

From Eq. (39), one can also obtain the integral representation

$\bar{M}=4\pi\int
d\mathsf{r}\,\mathsf{r}^{2}\frac{\rho(\mathsf{r}^{\prime})}{\mu^{2}\sigma_{0}^{2}}\,.$
(46)

Finally, the Noether charge is given by

$\displaystyle Q$ $\displaystyle\equiv$ $\displaystyle-\int
dVj^{t}=4\pi\omega\int drd\theta\,\sqrt{-g}e^{-v}\phi^{2}\,.$ (47)

###  III.3 Representative configurations

In the rest of this Section, we will set $\Lambda=0.186$ in order to
illustrate our results. In the gauge $\tilde{v}(0)=0$, the fictitious particle
initially does not feel the impact of gravity, and consequently, in the thin-
wall regime, we should expect that the central field is close to
$\varphi_{+}$.

Like in the flat space-time case, in curved space-time numerical results
produce two types of configurations - thin-wall ones for small $w$, and thick-
wall ones for larger $w$. Unlike in the Q-ball case, the thin-wall regime of
SBSs has two sub-classes, one perturbatively close to the flat space-time
case, and a non-perturbative one close to the maximum compactness
configuration. In Fig. 3 we show the $w-\varphi_{c}$, $\tilde{w}-\varphi_{c}$
and $\mathcal{C}-\varphi_{c}$ curves for our choice of $\Lambda$.

Various authors have used catastrophe theory arguments to assess the (local)
stability of Q-balls and boson stars Sakai and Sasaki (2008); Tamaki and Sakai
(2010, 2011a); Kleihaus  _et al._ (2012). We will not review these arguments
here, but let us mention that they show that stability can be checked by
analyzing the position of the turning points, for fixed $\Lambda$, in the
$\bar{M}-\varphi_{c}$ curve, up to a final (unstable) “spiral” that occurs for
$\varphi_{c}\,\mbox{{\scriptsize$\stackrel{{\scriptstyle>}}{{\sim}}$}}\,1$ in
the gravitating case. Thus, one can start from the stable MBS limit and track
the turning points as $\varphi_{c}$ increases. Already from Fig. 3, we see two
turning points that mark the boundaries of a middle stable branch, with two
unstable branches on the left and on the right of it (a low compactness stable
branch is not shown on this plot, see Fig. 14). Note however that the turning
point in the $\bar{M}-\varphi_{c}$ parameter space does not map exactly to the
same point in the $\mathcal{C}-\varphi_{c}$ parameter space. We have
established that the most massive configurations, which separate the compact
stable branch from the unstable one, have compactnesses slightly bellow the
maximum one. For example, for our representative case, the maximum mass
configuration has $\mathcal{C}=0.336$, while the maximally compact one has
$\mathcal{C}=0.346$.

We will now focus on the three representative configurations for each of the
three aforementioned (sub-)classes: the thick-wall regime (I); the thin-wall
regime perturbatively close to the Q-ball limit (II); and the thin-wall regime
close to the non-perturbative $\mathcal{C}_{\rm max}$ cut-off (III) for our
$\Lambda=0.186$. Configuration III is also the maximally massive one for our
$\Lambda$. The parameters of these configurations are given in Table 1. The
field and density profiles are shown in Figs. 5 and 6, while the metric
coefficients are displayed in Fig. 7.

| $\phi(0)/(\sigma_{0}/\sqrt{2})$ | $\tilde{\omega}/\mu$ | $\omega/\mu$ | $\tilde{v}(\infty)$ | $\mu R$ | $\mu M/m^{2}_{\rm Pl}$ | $\mathcal{C}$
---|---|---|---|---|---|---|---
I | 1.05 | 0.488 | 0.454 | 0.147 | 7.243 | 0.188 | 0.026
II | 1.014 | 0.241 | 0.146 | 1.009 | 29.821 | 5.719 | 0.192
III | 1.032 | 0.380 | 0.088 | 2.924 | 41.182 | 13.84 | 0.336
Table 1: Parameters of the representative configurations of SBSs: thick-wall
regime (I); thin-wall regime perturbatively close to the Q-ball limit (II);
strong gravity, thin-wall non-perturbative branch (III). Configuration III has
also the maximal value of mass for $\Lambda=0.186$ and thus is at the border
between the compact stable and the unstable branch. Figure 3: (Left) Frequency
vs central field $\varphi_{c}\equiv\phi(0)/(\sigma_{0}/\sqrt{2})$ for
$\Lambda=0.186$, in two gauges $\omega\,,\tilde{\omega}$ compared to the
Q-ball prediction appropriate in the $\tilde{v}(0)=0$ gauge (dark red line).
Note that two almost identical segments of the thin-wall regime in the
$\tilde{v}(0)=0$ gauge, denoted by the pink rectangles and the black stars,
correspond to separate segments in the $v(\infty)=0$ gauge, denoted by the
blue rectangles and the green stars, respectively. (Right) Compactness vs
central field for $\Lambda=0.186$ (cyan circles) and the Q-ball result (dark
red line). The three representative configurations (denoted by the orange club
for I, purple triangle for II and black square for III) are indicated.

###  III.4 Origin of maximum compactness

The maximum compactness of an object is usually discussed in the context of
Buchdahl’s theorem Buchdahl (1959). The latter states, in its most general
form, that GR self-gravitating objects with energy density monotonically
decreasing outwards $\rho^{\prime}\leq 0$ and positive (or vanishing) pressure
anisotropy (i.e. radial pressure larger than or equal to the tangential one)
$P_{\rm rad}\geq P_{\rm tan}$ must have $C\leq 4/9$ Buchdahl (1959); Urbano
and Veermäe (2019). In the case of fluid stars, isotropic “normal matter”
(i.e. one satisfying the weak energy condition $\rho\geq 0\,\land\,\rho+P\geq
0$ and the micro-stability condition $P\geq 0\,\land\,dP/d\rho\geq 0$) is
consistent with these assumptions. However, Buchdahl’s limit is not a robust
concept, as it can be easily evaded by violating the theorem’s assumptions,
producing even more compact objects. A discussion of the parameterized bounds
on $\mathcal{C}$ can be found in Guven and O’Murchadha (1999), and various toy
models that exceed the Buchdahl bound have been constructed in the literature
Cardoso and Pani (2019).

Buchdahl’s theorem can be additionally strengthened by requiring that the
equation of state be consistent with causality. In Urbano and Veermäe (2019)
(see also Alho  _et al._ (2021)), it has been shown that a useful toy model
for understanding the maximally compact and causal configurations is given by
objects described by a linear equation of state (LinEoS)

$\displaystyle\rho=\rho_{c}+\frac{P}{c^{2}_{s}}\,,$ (48)

where $c^{2}_{s}\equiv\partial P/\partial\rho>0$ (the speed of sound) and
$\rho_{c}$ are constant. Hence, allowing for $c^{2}_{s}>1$ accounts for
violations of causality. The limit $c^{2}_{s}\to\infty$ corresponds to
constant density stars, whose maximum compactness $\mathcal{C}_{\rm B}=4/9$
follows from Buchdahl’s theorem. For a given $c_{s}^{2}$, the LinEoS describes
the stiffest possible matter, and consequently it yields the most compact
configurations. Maximally compact and causal configurations consistent with
the assumptions of Buchdahl’s theorem cannot surpass
$\mathcal{C}_{\rm{B+C}}=0.354$ when $c^{2}_{s}=1$. The dependence of the
compactness of LinEoS configurations on $c^{2}_{s}$ is approximately given (to
within a $3.6\%$ error) by the fitting formula Urbano and Veermäe (2019)

$\frac{4}{9\mathcal{C}_{\rm{LinEoS}}}-1\approx\frac{0.77+0.51c^{2}_{s}}{c^{2}_{s}(4.18+c^{2}_{s})}\,.$
(49)

SBSs in the thin-wall regime are a physical example of objects with a
LinEoS777This is, as far as we know, original insight. A comparison between
SBSs and constant density stars $(c^{2}_{s}\to\infty)$ was discussed in Macedo
_et al._ (2013)., because in the bulk of the star $\varphi\approx 1$ and hence
$\varphi^{\prime}\approx V\approx 0$, making in turn $P_{\rm rad}\approx\rho$.
This argument implies that the maximal compactness of SBS is $\mathcal{C}_{\rm
max}\,\mbox{{\scriptsize$\stackrel{{\scriptstyle<}}{{\sim}}$}}\,\mathcal{C}_{\rm
B+C}$, which is consistent with our numerical results and the values reported
in the previous work Friedberg  _et al._ (1987a); Kesden  _et al._ (2005);
Cardoso  _et al._ (2016b); Visinelli (2021). Now we address two apparent
loopholes in the previous argument. First, SBSs do not have monotonic energy
density profiles in their surface regions (Fig. 6, right), which violates the
assumptions of Buchdahl’s theorem Urbano and Veermäe (2019). The region where
this violation occurs, however, is parametrically smaller than the size of the
SBS bulk in the thin-wall regime, and thus it does not affect
$\mathcal{C}_{\rm{B+C}}$ significantly. Secondly, although SBSs have
anisotropic pressure, the radial pressure is larger than the tangential one
[Eqns. (42), (43)], which does not allow for violating the Buchdahl
compactness bound (unlike the opposite case in which the tangential pressure
is larger than the radial one Raposo  _et al._ (2019); Alho  _et al._ (2021)).

In more detail, we can estimate how well the LinEoS describes the matter
inside the SBS in the thick wall regime. Assuming the Q-ball results (c.f.
Section III.5 for a justification) we will ignore exponentially suppressed
scalar derivatives [Eq. (17)] and approximate Eqns. (41), (42) to obtain

$\displaystyle\frac{\rho}{\frac{\mu^{2}\sigma^{2}_{0}}{2}}\approx
w_{+}^{2}\varphi^{2}_{c}+V(\varphi_{c})\,,\frac{P_{\rm
rad}}{\frac{\mu^{2}\sigma^{2}_{0}}{2}}\approx
w_{+}^{2}\varphi^{2}_{c}-V(\varphi_{c})\,,$

from where it follows that

$\displaystyle(c^{2}_{s})_{a}\approx\frac{\varphi_{c}^{2}\left(3\varphi_{c}^{2}-2\right)}{6\varphi_{c}^{4}-6\varphi_{c}^{2}+1}\,,$
(50)

with $w_{+}$ being the inverse of Eq. (16) and
$\varphi_{+}\approx\varphi_{c}$. For $\varphi_{c}=1.032$ (configuration III)
from Eq. (50) we find $(c_{s})_{a}=0.95$. Using this value in Eq. (49) we get
$\mathcal{C}_{\rm LinEoS}=0.350$, close to $\mathcal{C}_{\rm III}=0.336$.
Taking the limit $\varphi_{c}\to 1$, Eq. (50) implies $(c^{2}_{s})_{a}\to 1$
and hence $\mathcal{C}_{\rm max}\to\mathcal{C}_{\rm B+C}$. Note however that
the exact limit $\varphi_{c}=1$ in the thin-wall regime is attainable only in
the Minkowski space-time ($\Lambda=0$). This is a singular limit as the
absence of gravity implies $\mathcal{C}\to\infty$ when $\varphi\to 1$, as
elaborated in Section II.1.

In order to scrutinize previous analysis further, we have calculated (for our
representative $\Lambda=0.186$) numerically the average

$\langle c^{2}_{s}\rangle=\frac{1}{R_{<}}\int^{R_{<}}_{0}dr\frac{\partial
P_{\rm rad}}{\partial\rho}$ (51)

over several configurations [with $R_{<}$ the boundary of the bulk, c.f.
expression (68)], and we have compared the compactness predicted by Eq. (49)
with the actual one. The results are shown in Fig. 4 , and clearly indicate
good agreement in the thin wall regime (i.e. at high compactnesses). We have
also presented results for the position dependence of $c^{2}_{s}$ in Fig. 5 ,
for two specific configurations (I and III). Note that negative and seemingly
non-causal values of $c^{2}_{s}$ appear in the boundary zone, signaling a
breakdown of the hydrodynamic description. This breakdown also occurs in flat
space-time Nugaev and Shkerin (2020).

Finally, in the thick-wall regime, commensurability between the bulk and the
boundary does not allow for an effective LinEoS description (Fig. 5). As a
result, in this regime $\mathcal{C}\ll\mathcal{C}_{\rm{B+C}}$.

Figure 4: Relative difference between the numerically determined compactness
$\mathcal{C}$ and the prediction from the LinEoS $\mathcal{C}_{\rm LinEoS}$
(49) [using the numerically found average speed of sound (51)]. Figure 5:
Scalar field $\varphi\equiv\phi/(\sigma_{0}/\sqrt{2})$ and the speed of sound
$c^{2}_{s}$ radial profile for the two representative configurations: I
(orange) and III (black). The blue dot-dashed line represents the estimate of
the speed of sound (for III) $(c^{2}_{s})_{a}$ from Eq. (50).

###  III.5 Analytic construction

As already mentioned, the structure of SBSs can also be interpreted in light
of the dynamics of a fictitious Newtonian particle in a time dependent
potential, with $r$ playing the role of time888This interpretation was
mentioned in Tamaki and Sakai (2011b) at the qualitative level, but without
exploring its implementation and consequences. A similar approach was also
applied to (Minkowski) Q-balls with gauged $U(1)$ symmetry while this work was
well into preparation Heeck  _et al._ (2021b). , as can be seen from the
Klein-Gordon equation

$\displaystyle\varphi^{\prime\prime}+\Big{(}\frac{2}{\mathsf{r}}-\frac{W^{\prime}}{W}\Big{)}\varphi^{\prime}=\Big{[}\mathsf{m}^{2}(1-4\varphi^{2}+3\varphi^{4})-W^{2}\Big{]}\varphi\,,$
$\displaystyle\mu W=\omega e^{(u-v)/2}\,,\mu\mathsf{m}=\mu e^{u/2}\,,$ (52)

Equations (38) and (39) can be formulated as dynamical equations for “time-
dependent” frequency $\mu W$ and scalar mass $\mu\mathsf{m}$. Their dynamics
depends on $\varphi$ through $\rho$ and $P_{\rm rad}$. However, in certain
regimes we can approximately evaluate $W,\mathsf{m}$ without knowing the full
behaviour of $\varphi$.

If we work in the gauge where $\tilde{v}(0)=0$, the Q-ball does not “feel” the
gravitational field initially (i.e. at the centre). Like in the flat space-
time case, for a given $\varphi_{c}$ we need to find $W(0)=w$ such that the
“SBS-particle” rolls over the time-dependent potential and reaches $\varphi=0$
in an infinite time ($\mathsf{r}\to\infty$) with $\mathcal{E}=0$. Note that
$W$ is a gauge invariant object (i.e. it is left unaffected by a change of the
time coordinate).

The particle analogy sheds light, already at the qualitative level, on some of
the properties of SBSs: as Q-balls have a unique thin-wall regime, for SBSs
this regime will be described by a curve in the $\tilde{w}-\varphi_{c}$ space,
very close to the Q-ball result. However, for some SBS configurations, the
effective particle will exhibit the thin-wall regime perturbatively close to
the Q-ball one (like for configuration II), while others will be in the non-
perturbative part of the parameter space (like in the case of configuration
III). Thus, the thin wall regime in the $\tilde{w}-\varphi_{c}$ representation
consists of two almost degenerate curves (denoted by black stars and pink
rectangles in Fig. 3). For example, there are configurations where the
difference in $\tilde{w}$ can be as low as $\sim 10^{-51}$, with the
corresponding difference in $w$ $\sim 10^{-6}$. These curves split into two
separate curves in the $w-\varphi_{c}$ parameter space, one close to the
Q-ball line, while the other is the horizontal asymptote, represented by green
stars and blue rectangles on Fig. 3, respectively.

In the following subsections, we will give an analytic description of the
scalar in the strong-field regime in three zones (analogous to the flat-
spacetimes ones introduced in Sec. II), starting from the simplest one.

####  III.5.1 Exterior zone

The simplest regime is the asymptotic one, where the space-time is to a good
approximation Schwarzschild, owing to the fast decay of the scalar:

$\displaystyle u_{>}=-\log{\Big{(}1-\frac{2\bar{M}}{\mathsf{r}}\Big{)}}\,,$
(53)
$\displaystyle\tilde{v}_{>}=\tilde{v}_{\infty}+\log{\Big{(}1-\frac{2\bar{M}}{\mathsf{r}}\Big{)}}\,.$
(54)

The evolution of the scalar is found by solving the Klein-Gordon equation (40)
with these $u,v$, expanding in powers of $1/r$:

$\displaystyle[\mathcal{O}_{\rm L}+\mathcal{O}_{\rm
NL}]\varphi(\mathsf{r})=0\,,$ (55) $\displaystyle\mathcal{O}_{\rm
L}=\frac{d^{2}}{d\mathsf{r}^{2}}+\frac{2}{\mathsf{r}}\frac{d}{d\mathsf{r}}-\Big{(}1-e^{-\tilde{v}_{\infty}}\tilde{w}^{2}+\frac{2\bar{M}}{\mathsf{r}}\big{(}1-2e^{-\tilde{v}_{\infty}}\tilde{w}^{2}\big{)}\Big{)}\,,$
$\displaystyle\mathcal{O}_{\rm
NL}=\frac{\bar{M}}{\mathsf{r}}\Big{(}6\varphi^{5}-8\varphi^{3}\Big{)}+3\varphi^{5}-4\varphi^{3}\,.$

Note that the non-linear terms can be treated as a perturbation. While it
appears that there is no simple analytic estimate of these corrections, the
leading term gives a good description of the asymptotic behaviour, because the
scalar is suppressed (more than exponentially).

Requiring $\varphi\to 0$ as $r\to\infty$, we obtain the solution in terms of
hypergeometric functions, or equivalently in terms of the Wittaker function
Kling and Rajaraman (2017)

$\displaystyle\varphi_{\infty}=\frac{\mathcal{A}}{2\mathsf{r}\sqrt{1-e^{-\tilde{v}_{\infty}}\tilde{w}^{2}}}\times$
$\displaystyle\,W\Big{[}\bar{M}\frac{2e^{-\tilde{v}_{\infty}}\tilde{w}^{2}-1}{\sqrt{1-e^{-\tilde{v}_{\infty}}\tilde{w}^{2}}},-\frac{1}{2},2\sqrt{1-e^{-\tilde{v}_{\infty}}\tilde{w}^{2}}\Big{]}\,.$
(56)

This asymptotic solution is parameterized by $\tilde{v}_{\infty}$, $\bar{M}$,
$\tilde{w}$ and the normalization amplitude $\mathcal{A}$. The linear equation
$\mathcal{O}_{\rm L}\varphi(\mathsf{r})=0$ is invariant under the rescaling
$\varphi\to\mathsf{C}\varphi$, but the matching condition breaks the
invariance and selects $\mathsf{C}=\mathcal{A}$. Expanding the Whittaker
function in $1/\mathsf{r}$, the leading behaviour gives

$\displaystyle\varphi_{\infty}\simeq\frac{\mathcal{A}}{\mathsf{r}^{1+\beta_{>}}}e^{-\alpha_{>}\mathsf{r}}\,,$
(57)
$\displaystyle\alpha_{>}=\sqrt{1-e^{-\tilde{v}_{\infty}}\tilde{w}^{2}}\,,$
(58)
$\displaystyle\beta_{>}=\frac{\bar{M}}{\alpha_{>}}(1-2e^{-\tilde{v}_{\infty}}\tilde{w}^{2})\,.$
(59)

Note that for $\Lambda\to 0$, $\bar{M}\to 0$ and we recover the Q-ball result
(15). We show a comparison between the numerical and the asymptotic field
behaviour in Fig. 6.

Figure 6: (Left) Numerical results vs asymptotic approximation for the ratio
between the field radial profile and its derivative $\phi/\phi^{\prime}$ for
configurations I, II. (Right) Numerical results vs boundary zone analytic
approximation for the field energy density radial profile
$\rho/(\frac{\mu^{2}\sigma^{2}_{0}}{2})$ for configurations I, II, III.
Numerically determined parameters are used as input.

####  III.5.2 Interior zone

In the gauge $\tilde{v}(0)=0$, the flat space-time result (17) should be valid
sufficiently close to the origin, and the scalar field derivative is
suppressed by $\tilde{w}$. Thus, we will calculate the metric coefficients in
the interior perturbatively in $\tilde{w}$, and approximate
$\varphi^{\prime}\approx 0$ and $V_{<}\approx V(\varphi_{c})$. This
description provides a good approximation for configurations similar to II in
the strong-field and the thin-wall regime and close to the Q-balls limit,
where $\tilde{w}$ is small and controls the size of the star.

The perturbative expansion in $\tilde{w}$ must be appropriately resummed (or
performed from the start in a suitable form) so that it can be matched with
the exterior. Taking

$\displaystyle\tilde{u}_{<}=\log[\mathsf{\tilde{u}}_{0}(r)+\mathsf{\tilde{u}}_{2}(r)\tilde{w}^{2}+\mathsf{\tilde{u}}_{4}(r)\tilde{w}^{4}+\mathcal{O}(\tilde{w}^{6})]\,,$
(60)
$\displaystyle\tilde{v}_{<}=\log[\mathsf{\tilde{v}}_{0}(r)+\mathsf{\tilde{v}}_{2}(r)\tilde{w}^{2}+\mathsf{\tilde{v}}_{4}(r)\tilde{w}^{4}+\mathcal{O}(\tilde{w}^{6})]\,,$
(61)

from Eqns. (38) and (39) one gets

$\displaystyle\mathsf{\tilde{u}}_{0}=1\,,\,\,\mathsf{\tilde{u}}_{2}=\frac{1}{6}\Lambda^{2}r^{2}\,,\mathsf{\tilde{u}}_{4}=\frac{1}{360}\left(45\Lambda^{2}r^{2}-2\Lambda^{4}r^{4}\right)\,$
(62)
$\displaystyle\mathsf{\tilde{v}}_{0}=1\,,\mathsf{\tilde{v}}_{2}=\frac{1}{3}\Lambda^{2}r^{2}\,,\mathsf{\tilde{v}}_{4}=\frac{1}{120}\left(15\Lambda^{2}r^{2}+4\Lambda^{2}r^{4}\right)\,.$
(63)

In Fig. 7 we show the numerical results and the interior and asymptotic
approximations for the metric coefficients.

Figure 7: (Left) Numerical vs. approximate analytic (interior, boundary,
asymptotic) results for the $u$ metric coefficient radial profile for II.
(Right) Numerical vs. approximate analytic (interior, asymptotic) results for
the $v$ radial profile metric coefficient for II. Numerically determined
parameters are used as input.

####  III.5.3 Boundary zone

In the transition region, let us expand around $(\mathsf{R^{\ast}})^{-1}$ as
in Heeck  _et al._ (2021a). Thus, we neglect the friction $1/R^{\ast}$ term in
the Klein-Gordon equation (40). We can estimate the contribution from the
time-dependent frequency as

$\frac{W^{\prime}}{W}\sim\frac{1}{1-\frac{2\bar{M}}{\mathsf{R^{\ast}}}}\frac{2\bar{M}}{(\mathsf{R^{\ast}})^{2}}\sim\frac{2}{\mathsf{R^{\ast}}(\mathcal{C}^{-1}-1)}\,.$
(64)

As $\mathcal{C}_{\rm{max}}\sim 0.35$ (cf. Sec. III.4), we find that this term
is of the order of $1/R^{\ast}$ and can be neglected. Now, neglecting friction
and assuming in the first iteration $W\approx W_{\ast}\equiv
W(\mathsf{R^{\ast}})\,,\mathsf{m}\approx\mathsf{m}_{\ast}\equiv\mathsf{m}(\mathsf{R^{\ast}})$
(as we are interested only in a tiny strip around of the thin wall), we can
use the conservation of energy and the fact that $\mathcal{E}(\infty)=0$ to
reduce the equations of motion to Heeck  _et al._ (2021a)

$\varphi^{\prime}=\pm\varphi\sqrt{\Big{[}\mathsf{m}^{2}_{\ast}(1-\varphi^{2})^{2}-W^{2}_{\ast}\Big{]}}\,.$
(65)

In the thin-wall regime, the expectation $W^{2}_{\ast}\ll 1$ leads to

$\varphi_{B}=\frac{1}{\sqrt{1+2\exp{[2\mathsf{m}_{\ast}(\mathsf{r}-\mathsf{R^{\ast}})]}}}\,,$
(66)

where we have specified the integration constant by requiring
$\varphi^{\prime\prime}(\mathsf{R^{\ast}})=0$. Note that the second
integration constant is determined by the value of the energy and the
requirement that the function is monotonously decreasing. Furthermore, Heeck
_et al._ (2021a) includes an ad hoc prefactor $\varphi_{+}$ in $\varphi_{B}$
for Q-balls [formula (14)], because it slightly improves the analytic
description in the thick-wall limit. (Note that in the thin-wall regime
$\varphi_{+}\approx 1$ in any case.) We have not included this prefactor in
the SBS context, as it tends to worsen the model.

From the analysis around Eq. (22), we find the estimate of the width of the
density support in the boundary zone to be

$\displaystyle\lambda=\frac{3}{4}\frac{2.66}{\mathsf{m}_{\ast}}\,,$ (67)
$\displaystyle\mathsf{R_{>/<}}=\mathsf{R^{\ast}}\pm\lambda$ (68)

In the flat space-time limit $\mathsf{m}_{\ast}=1$ and we recover Eq. (14).

A plot of the field density is given in Fig. 6, where the metric coefficients
are taken from the interior zone perturbative series (60), (61). Note that the
boundary-zone solution works (somewhat surprisingly) even far from its region
of a priori validity, down to the transition region (like in the flat space-
time case). The worst agreement occurs for the configuration I, as expected
since this configuration is in the thick-wall regime.

Having a preliminary understanding of the field behaviour in the boundary
zone, as well as of the radial profile of the metric coefficients in the
internal and external zones, we can understand the junction conditions for the
metric coefficients. From expression (66):

$\displaystyle\bar{\rho}_{\rm
kin}=e^{-v}\tilde{w}^{2}\varphi^{2}\approx\frac{e^{-v_{\ast}}\tilde{w}^{2}}{2e^{2\mathsf{m}_{\ast}z}+1}\,$
(69) $\displaystyle\bar{\rho}_{\rm
st}=e^{-u}(\varphi^{\prime})^{2}\approx\frac{4e^{4\mathsf{m}_{\ast}z}}{\left(2e^{2\mathsf{m}_{\ast}z}+1\right)^{3}}\,$
(70) $\displaystyle\bar{\rho}_{\rm
pot}=V\approx\frac{4e^{4\mathsf{m}_{\ast}z}}{\left(2e^{2\mathsf{m}_{\ast}z}+1\right)^{3}}\,,$
(71)

with $z=\mathsf{r}-\mathsf{R^{\ast}}$ and
$\bar{\rho}_{i}=\Big{(}\frac{\mu^{2}\sigma^{2}_{0}}{2}\Big{)}^{-1}\rho_{i}$.
In the thin wall limit one has $\bar{\rho}_{\rm kin}\to H(-z)$, where $H$ is
the Heaviside step function, while $\bar{\rho}_{\rm st}\to\delta(z)$. From
Eqns. (38), (39) and the limiting cases, we can expect that $u$ will have a
step-like behaviour at $R^{\ast}$, while $v$ will present a smoother
transition. This type of behaviour was noted already in Friedberg  _et al._
(1987a). Solving Eq. (38) we find a complicated expression that we report in
App. C.1. For the configuration II we show $u_{B}$ in Fig. 7 (Left).

####  III.5.4 Energy balance

Like in the flat space-time case [Eq. (18)], we can determine the inflection
point in the case of SBSs by using energy balance arguments. We now need to
account for the “time” dependence of the potential parameters as

$\frac{dU_{\omega}}{dr}=\frac{\partial
U_{\omega}}{\partial\phi}\phi^{\prime}+\frac{\partial U_{\omega}}{\partial
r}\,.$ (72)

One then finds

$\displaystyle\frac{2}{\mu^{2}\sigma^{2}_{0}}\mathcal{E}(0)=$ (73)
$\displaystyle\int^{\infty}_{0}d\mathsf{r}\big{[}\frac{2}{\mathsf{r}}(\varphi^{\prime})^{2}-\frac{W^{\prime}}{W}(\varphi^{\prime})^{2}+\mathsf{m}^{\prime}\mathsf{m}\varphi^{2}(1-\varphi^{2})-W^{\prime}W\varphi^{2}\big{]}\,.$

Note that in the last expression, only the first term is present in the
Minkowski limit (18), because $\mu W\to\omega$, $\mu\mathsf{m}\to\mu$ and
$\omega,\mu$ do not run in “time”. Like the first one, the other terms are
also dominant in the particular zones, as can be inferred from their form. The
second and third terms are important in the boundary zone where the field
interpolates between $\varphi_{c}$ and the exponential tail, while the fourth
term receives important contributions both from the interior and the boundary
zone. The leading order behaviour of all these terms is provided separately by

$\displaystyle A_{\mathcal{E}}$ $\displaystyle\equiv$
$\displaystyle\int^{\mathsf{R_{>}}}_{\mathsf{R_{<}}}d\mathsf{r}\frac{2}{\mathsf{r}}(\varphi^{\prime})^{2}\approx\frac{\mathsf{m}_{\ast}}{2\mathsf{R^{\ast}}}\,,$
(74) $\displaystyle B_{\mathcal{E}}$ $\displaystyle\equiv$
$\displaystyle-\int^{\mathsf{R_{>}}}_{\mathsf{R_{<}}}d\mathsf{r}\frac{W^{\prime}}{W}(\varphi^{\prime})^{2}\approx-\frac{\mathsf{m}_{\ast}}{\mathsf{R^{\ast}}}\Big{[}\frac{1-\mathsf{m}^{2}_{\ast}}{4}+\frac{1}{80}\mathsf{m}^{2}_{\ast}\Lambda^{2}\mathsf{R^{\ast}}^{2}\Big{]}\,,$
(75) $\displaystyle C_{\mathcal{E}}$ $\displaystyle\equiv$
$\displaystyle\int^{\mathsf{R_{>}}}_{\mathsf{R_{<}}}d\mathsf{r}\mathsf{m}^{\prime}\mathsf{m}\varphi^{2}(1-\varphi^{2})\approx\frac{\mathsf{m}_{\ast}}{\mathsf{R^{\ast}}}\Big{[}\frac{1-\mathsf{m}_{\ast}^{2}}{8}+\mathsf{m}_{\ast}^{2}\Lambda^{2}\mathsf{R^{\ast}}^{2}\Big{(}\frac{1}{80}+\frac{1}{48}w^{2}e^{-v_{\ast}}\Big{)}\Big{]}\,,$
(76) $\displaystyle D_{\mathcal{E}\,<}$ $\displaystyle\equiv$
$\displaystyle-\int^{\mathsf{R_{<}}}_{0}d\mathsf{r}W^{\prime}W\varphi^{2}\approx\frac{w^{2}\varphi_{c}^{2}}{2}\big{[}1-e^{u_{<}(\mathsf{R_{<}})-v_{<}(\mathsf{R_{<}})}\big{]}\,,$
(77) $\displaystyle D_{\mathcal{E}\,B}$ $\displaystyle\equiv$
$\displaystyle-\int^{\mathsf{R_{>}}}_{\mathsf{R_{<}}}d\mathsf{r}W^{\prime}W\varphi^{2}\approx-2\frac{\mathsf{m}_{\ast}}{\mathsf{R^{\ast}}}w^{2}e^{-v_{\ast}}\Big{[}\frac{1}{2}\log\Big{(}\frac{3}{2}\Big{)}(1-\mathsf{m}^{2}_{\ast})+\frac{7}{648}\mathsf{m}^{2}_{\ast}\Lambda^{2}\mathsf{R^{\ast}}^{2}\Big{]}\,.$
(78)

The left hand side of Eq. (73) is therefore given as in Eq. (II.1), i.e.,
neglecting exponentially suppressed terms, by

$\displaystyle\bar{\mathcal{E}}(0)\approx\frac{\left(\sqrt{3w^{2}+1}+2\right)}{27}\left(3w^{2}+\sqrt{3w^{2}+1}-1\right).$
(79)

From the numerical results of Fig. 3, it is clear that the effect of (strong)
gravity is to introduce a new branch. As a result, for a subset of
$\varphi_{c}\sim 1$ we have two different stable configurations [and possibly
one more (un)stable one] for the same $\varphi_{c}$. For the class of
configurations that contains I and II, we expect the Q-ball result
$\mathsf{R^{\ast}}\sim 1/\tilde{w}^{2}$. For the most compact configurations,
from Eqns. (53), (60) and ignoring the jumping conditions, one naively gets

$\displaystyle\mathsf{R^{\ast}}\sim\frac{\sqrt{12\mathcal{C}_{\rm
max}}}{\Lambda\tilde{w}}\sim\frac{2}{\Lambda\tilde{w}}\,.$ (80)

where we assumed in the last relation
$\mathcal{C}_{\max}\approx\mathcal{C}_{\rm B+C}$. Although this expectation is
too simplistic, it suggests the useful variable

$\displaystyle T\equiv\mathsf{R^{\ast}}\Lambda\tilde{w}\,.$ (81)

Because of the aforementioned degeneracy, $T(\tilde{w})$ is not a single-
valued function. Instead, it is more useful to look for
$\tilde{w}=\tilde{w}(T)$. The minimum of this curve, $\tilde{w}_{\cup}$,
separates the Q-ball-like branch from the non-perturbative strong-gravity one,
and provides the lowest $\tilde{w}$ for a given $\Lambda$. In flat space-time,
one has $\tilde{w}_{\cup}=0$ $(\mathsf{R^{\ast}}\to\infty)$, so we expect
$\tilde{w}_{\cup}\sim\mathcal{O}(\Lambda)$.

Approximating the complicated algebraic expressions in $\mathsf{R^{\ast}}$ we
find (see App. C.2 for the details)

$\displaystyle\tilde{w}\approx\Lambda\frac{\sqrt{-\frac{T^{4}}{5}+6T^{2}+36}\left(T^{4}+10T^{2}+30\right)}{T\left(-T^{4}+30T^{2}+180\right)}\,.$
(82)

The minimum of this function corresponds to $\tilde{w}_{\cup}/\Lambda\approx
1.1$, while for $T\to 0$ we recover the Q-ball asymptotics
$\tilde{w}\sim\Lambda/T$. From there the approximate behaviour of $w$ follows:

$\displaystyle\frac{w}{\tilde{w}}\approx$
$\displaystyle\Big{(}1+\frac{5T^{2}\left(3\tilde{w}^{2}+4\right)-(2/3)T^{4}}{(2\sqrt{30})^{2}}\Big{)}^{-1/2}\times\,$
(83)
$\displaystyle\Big{(}1+\frac{5T^{2}\left(3\tilde{w}^{2}+8\right)+4T^{4}}{(2\sqrt{30})^{2}}\Big{)}^{-1/2}\,,$

which we show in Fig. 8. Our rough approximations have therefore given us a
simple analytic description that predicts the horizontal branching off of the
$w-\varphi$ curve for any $\Lambda\ll 1$ (corresponding to the non-
perturbative effect of strong-gravity), with $\mathcal{C}$ getting close to
LinEoS limit discussed in Sec. III.4.

####  III.5.5 Finale: Semi-analytic solution

Instead of ignoring sub-leading terms, which led us to Eq. (83), we can also
solve the full master (algebraic) Eqns. (74)-(79) numerically, using
expression (82) as a guess for each $T$. Now, we use $u_{B}$ when calculating
$\mathsf{m}_{\ast}$ in Eqns. (74)-(78), except in $u_{B}$ itself [expression
(123)], where we iteratively take $\mathsf{m}_{\ast}\approx$
$\exp{(u_{<}(\mathsf{R}^{\ast})/2)}$ . In contrast to the costly numerical
solution of the boundary-value Einstein-Klein-Gordon system (App. A), the
semi-analytic approach provides a solution in a few seconds on a laptop
computer.

This procedure leads to a very good description of the configurations in the
stable compact branch, as shown in Fig. 8. The semi-analytic results are in
excellent agreement with the numerical ones for $\bar{M}$, $v_{\infty}$, and
consequently for $w$. For both $\mathsf{R^{\ast}}$ and $\mathsf{R}$, the
agreement is almost perfect in the Q-ball-like strong-gravity branch (i.e. for
configurations similar to II), while in the non-perturbative strong-gravity
branch (similar to III), there is a systematic deviation, up to a few percent
relative error. This error increases as $\varphi$ increases, but note that
this occurs mostly in the unstable branch, for which accurate approximations
are not of crucial importance. We elaborate more on the reason for the
systematic error in App. C.3.

Once the parameters are determined in the described way, one can use the
expansions from the previous Subsections to reconstruct both the scalar and
the gravitational field throughout the space-time, as well as the
thermodynamic functions: density, pressure(s), speed of sound etc.

Figure 8: (Left) Approximate analytic, semi-analytic (pink diamonds) and
numerical calculations (cyan circles) of the $w-\varphi$ behaviour for SBSs.
The purple dashed line shows both branches from expression (83), while the
dark red line represents the the Q-ball limit [expression (20), using
expansion (60), (61)]. (Right) Semi-analytic vs. numerical calculation of the
$M-R$ curves. Both plots correspond to the benchmark scenario $\Lambda=0.186$.
Three representative configurations from Table 1 are also shown.

##  IV Parameter space of Soliton boson stars

In this Section, we will move away from the benchmark scenario of Sec. III,
where we only considered the compact stable branch of SBSs with the simple
potential (1) for $\Lambda\ll 1$. We will now consider the Planck limit
$\Lambda\sim 1$ in Sec. IV.1 and the low-compactness stable branch (for
generic $\Lambda$) in Sec. IV.2. We will then adopt potentials with multiple
degenerate vacua in Sec. IV.3, ones with a false vacuum instead of a
degenerate one in Section IV.4, and we will test the robustness of our
conclusions in Sec. IV.5, by considering a non-polynomial effective potential.

###  IV.1 Planck scale regime

To understand the qualitative impact of a large “control parameter” $\Lambda$,
let us assume that the Q-ball description is valid up to
$\mathcal{C}_{\rm{B+C}}$ i.e. $\mathcal{C}_{\rm{max}}\sim$
$\frac{\Lambda^{2}}{16\pi}\varphi^{2}_{c}\mathsf{R}^{2}w^{2}$
$\propto\Lambda^{2}/w_{\rm{min}}^{2}$ $<\mathcal{C}_{\rm{B+C}}$. As $\Lambda$
increases, in order for $\mathcal{C}_{\rm{max}}$ to asymptote to
$C_{\rm{B+C}}$ in the thin wall regime, $w_{\rm{min}}$ has to increase.
However, as $w_{\rm{min}}$ increases, the thin wall regime is superseded by
the thick wall one, and when $w_{\rm{min}}\simeq w_{\rm{Qb-s}}$ [cf. Eq. (29)]
the stable branch inherited from the flat space-time limit disappears.

We have presented numerical results for $\tilde{w}\,,w$ as functions of
$\Lambda$ in Fig. 9. Note that the dip in the $w(\varphi_{c})$ curve increases
with $\Lambda$. This corresponds to the growth of the height of the horizontal
asymptote of the $w(\varphi_{c})$ and hence minimal possible $w$ as argued
above. Consequently the Planck scale regime leads to less compact
configurations (Fig. 11). As strong deviations from the Q-ball description
occur in the compact stable branch at $\Lambda\simeq 0.4$, we cannot make
quantitative predictions for when the compact stable branch will disappear,
but numerically we find that this occurs at $\Lambda\simeq 1.1$. The mass-
radius relation is shown in Fig. 10, while Fig. 11 represents how the maximum
achievable compactness changes999It is interesting to note that scalar stars
in Horndeski’s theory have been recently constructed and present very similar
a behaviour for $\mathcal{C}$ (c.f. Fig. 10 in Barranco  _et al._ (2021)).
However, we are not aware of a simple mapping between SBSs in GR and this kind
of configurations. from $C_{\rm B+C}$ to $C_{\rm MBS}\approx 0.11$ over
$\Lambda\sim 0.1-1$ range.

Figure 9: (Left) Frequency vs central scalar field value in $\tilde{v}(0)=0$
gauge and (Right) in the $v(\infty)=0$ gauge, for SBSs with different values
of $\Lambda$, along with the (gauge-independent) analytical result for Q-balls
(16).  Figure 10: Mass-radius relation for SBSs with different values of
$\Lambda$. The dots represent numerical results, while the dashed black line
represents the analytic result for NBS (88). The shapes and colours of dots
are the same as in Fig. 9. Figure 11: Maximum compactness of SBSs for given
$\Lambda$. We have also indicated the Buchdahl bound $\mathcal{C}_{\rm
B}=0.44$, the causal Buchdahl bound $\mathcal{C}_{\rm B+C}=0.354$, the
condition for the photon sphere $\mathcal{C}=0.33$ and the maximal compactness
for MBSs $\mathcal{C}=0.11$. In the ultra-compact domain (cyan region), the
numerical results are supplemented with semi-analytical ones. Note that for
$\Lambda=0$: $\mathcal{C}_{\rm{max}}\to\infty$ (Q-ball limit). The orange
region represents the MBS part of the parameter space and extends to
$\Lambda\to\infty$.

###  IV.2 Low compactness stable branch

The low compactness stable branch of SBSs is supported by quantum pressure,
and as a result in this limit SBSs behave as MBSs. The exact MBS limit of SBSs
is $\Lambda\to\infty$, and the appropriate field parametrization is

$\varphi_{\rm MBS}=\phi/M_{\rm Pl}\,,$ (84)

as there are no scales on which the structure equations depend, except for the
Planck scale. These configurations were thoroughly studied in the original
paper by Kaup Kaup (1968) and are reviewed in Liebling and Palenzuela (2017);
Schunck and Mielke (2003). The numerical procedure to obtain MBSs is analogous
to the SBS case (App. A).

In the SBS setting, at $\Lambda\simeq 1.1$ the compact stable branch of SBS
vanishes, and only the MBS one is left. For small $\Lambda$, the low
compactness stable branch is replaced by the unstable Q-ball branch in the
weak field regime. As the control parameter increases from $\Lambda\simeq
0.7$, more and more configurations in the low compactness stable branch
develop higher compactness, up to the MBS limit of $C_{\rm MBS}\approx 0.11$.
The low compactness branch presents different behaviour than the compact one,
where the Planck scale regime leads to lower compactness. Note that the
highest compactness is achieved in the unstable branch, while in the stable
one $C^{(s)}_{\rm MBS}\approx 0.08$.

In the weak-field approximation of MBSs, the Einstein-Klein-Gordon system
reduces to the Schrödinger-Poisson system (Newtonian boson stars) Seidel and
Suen (1990); Guzman and Urena-Lopez (2003); Liebling and Palenzuela (2017);
Kling and Rajaraman (2017); Bošković  _et al._ (2018); Annulli  _et al._
(2020)

$\displaystyle e\phi=-\frac{1}{2\mu}\nabla^{2}\phi+\mu\phi\Omega\,,$ (85)
$\displaystyle\nabla^{2}\Omega=\frac{1}{2}M^{-2}_{\rm Pl}\mu^{2}\phi^{2}\,,$
(86)

where $e^{v}\approx 1+2\Omega$, $\omega\approx\mu+e$ and $\Omega\ll 1$,
$e\ll\mu$. Like in the general discussion, the scalar mass can be factored
out, and in addition system admits a scaling symmetry

$\displaystyle\phi\to k^{2}\phi\,,e\to k^{2}e\,,r\to r/k\,,\Omega\to
k^{2}\Omega\,.$ (87)

This allows for a universal description of these objects, and for performing
the numerical integration only once (for the analytic solution see Kling and
Rajaraman (2017, 2018)).

Fixing the scale with $-k^{2}=e/2$, we can find several useful relations
between the macroscopic parameters, which will be compared with the
relativistic numerical results, e.g.

$\displaystyle\bar{M}=\frac{\beta Z}{\bar{R}}\,,$ (88)
$\displaystyle\bar{M}=\sqrt{2\beta^{2}\Big{(}1-w\Big{)}}\,,$ (89)
$\displaystyle\varphi_{c}=\frac{2s_{0}}{\Lambda}\frac{1}{\sqrt{w}}\Big{(}1-w\Big{)}\,,$
(90)

where $s_{0}=1.022$ and $\beta=1.753$ Kling and Rajaraman (2017) and we find
the scale-invariant radius (that encloses 99% of the BS mass) to be
$Z=5.6741$.

In Fig. 10, we see that the Newtonian boson star scaling gives a good
description of SBS configurations in the $\varphi_{c}\to 0$ limit.

###  IV.3 Cosine potential

In axionic physics, a cosine potential $V\sim\cos{(a/f_{a})}$, where $a$ is
the axion field and $f_{a}$ is a decay constant, is often considered. This
potential arises from non-perturbative effects that generate small masses for
the (initially massless) Goldstone boson associated with the spontaneous
breaking of the Peccei-Quinn symmetry Marsh (2016); Hui  _et al._ (2017);
Helfer  _et al._ (2017). In the strong gravity context, an axion potential can
produce non-trivial effects on the stability of axion stars Helfer  _et al._
(2017).

Inspired by axion star solutions with the axion potential, some authors have
considered boson star models with similar potentials Guerra  _et al._ (2019);
Delgado  _et al._ (2020); Siemonsen and East (2021). Note that in the absence
of beyond standard model physics that could motivate such potentials for
complex scalars with $U(1)$ symmetry, one should consider these models only as
proxies to understand (pseudo-real) axion stars (and only if different minima
are physically sensible). As these potentials develop multiple minima and
having in mind the Taylor expansion of the $\cos$ function, our previous
discussion would imply that “axion boson stars” would periodically replicate
SBSs for the different values of $\Lambda$ corresponding to different minima,
up to the Planckian threshold.

For concreteness we will consider a specific form of the potential, from
Guerra  _et al._ (2019):

$V=\frac{2\mu^{2}f^{2}_{a}}{B}\Big{[}1-\sqrt{1-4B\sin^{2}\Big{(}\frac{|\Phi|}{2f_{\rm
a}}\Big{)}}\Big{]}\,,$ (91)

where $B$ is a model dependent constant [taken to be $B\approx 0.22$ in Guerra
_et al._ (2019)]. The minimum of potential (91) occurs at

$\phi_{\rm min}=f_{a}2n\pi\,,n\in\mathbb{N}\,.$ (92)

This gives an $n$-dependent $\Lambda$ scale

$\Lambda_{n}=\frac{f_{a}}{m_{\rm Pl}}2n\pi\sqrt{16\pi}\,,n\in\mathbb{N}\,.$
(93)

In Fig. 12, we have compared numerical results from Guerra  _et al._ (2019)
with a set of SBSs specified by the potential (1) and the control parameter
(93). It is clear that the periodic features for the “axion boson stars” occur
for the appropriate field values given by (92). The $\bar{M}-\phi(0)$ plot for
the first minimum is in excellent quantitative agreement with the SBS results.
The agreement, however, is not perfect as $\cos$ can only locally be
approximated with the sextic polynomial. As we progress in $n$, the agreement
worsens. This should come as no surprise, because in the analogue perspective
an axion boson star “particle” has to go through an effective potential that
has several peaks and troughs, unlike in the SBS case. Finally, a sufficiently
large $n$ field is in the Planck scale regime, and the final unstable branch
is reached, like in the SBS case (Sec. IV.1).

---
Figure 12: Mass-central field dependence for the cosine potential (91). The
black line represents numerical results from Guerra  _et al._ (2019), while
dots correspond to a set of SBSs with $\Lambda_{n}$ from Eq. (93) and
$n=1,2,3,4$. The left panel represents cosine potential boson stars with
$f_{a}/m_{\rm Pl}=10^{-2.2}$, while the right one assumes $f_{a}/m_{\rm
Pl}=10^{-2.3}$, with the other parameters set to the values used in Guerra
_et al._ (2019). Note the different rescaling of the scalar with respect to
the rest of this work, in accordance with the conventions of Guerra  _et al._
(2019).

The cosine potential example is illustrative also for the following reason:
the sextic potential is non-renormalizable and thus not valid for arbitrary
field values (as one needs to be within the limit of validity of the effective
field theory). Close to this limit, higher-dimensional operators become
relevant. The cosine potential illustrates that such terms do not modify
qualitatively the macroscopic behaviour, as long as degenerate vacua (and even
false ones, as we will argue next) are present.

###  IV.4 General sextic potential

The case of degenerate vacua is somewhat special, while a more generic
scenario would allow for a non-zero false vacuum:

$\displaystyle V=\mu^{2}|\Phi|^{2}-\beta|\Phi|^{4}+\xi|\Phi|^{6}\,,$ (94)
$\displaystyle\beta>0\,,\xi>0\,.$

There are two useful reparametrizations of this potential. The first one, used
in Heeck  _et al._ (2021a), parametrizes the potential as a deviation from the
degenerate vacuum case:

$\displaystyle
V_{6}=\phi^{2}_{0}\Big{[}(\mu^{2}-\omega^{2}_{0})\varphi^{2}(1-\varphi^{2})^{2}+\omega^{2}_{0}\varphi^{2}\Big{]}\,,$
(95)
$\displaystyle\varphi=\frac{|\Phi|}{\phi_{0}}\,,\,\phi_{0}=\sqrt{\frac{\beta}{2\xi}}\,,\,w_{0}=\frac{\omega_{0}}{\mu}=\sqrt{1-\frac{\beta^{2}}{4\xi\mu^{2}}}$

The parameter choices $\beta=4\mu^{2}/\sigma_{0}^{2}$,
$\xi=4\mu^{2}/\sigma_{0}^{4}$ imply $\omega^{2}_{0}=0$ and
$\phi_{0}=\sigma_{0}/\sqrt{2}$ and reproduce the potential (1) i.e. the
benchmark scenario of this work, while $\omega_{0}$ parametrizes deviation
from the simplest potential (1). Another useful approach is to relate
$\omega_{0}$ to the ratio between the potential barrier $\phi_{\rm B}:$
$dV/d\phi|_{\rm B}=0$ , $d^{2}V/d\phi^{2}|_{\rm B}<0$ and the non-trivial
minimum (false vacuum) $\phi_{\rm F}$:

$\displaystyle V_{6}$ $\displaystyle=$
$\displaystyle\frac{\mu^{3}}{6\sqrt{3\xi}}\frac{\varphi_{x}^{2}}{x^{3}}\Big{(}6x^{2}+\left(-3x^{2}-3\right)\varphi_{x}^{2}+2\varphi_{x}^{4}\Big{)}\,,$
(96) $\displaystyle\varphi_{x}$ $\displaystyle\equiv$
$\displaystyle\frac{\phi}{\phi_{\rm B}}\,,\phi_{\rm
F/B}=\frac{1}{\sqrt{3\xi}}\sqrt{\beta\pm\sqrt{\beta^{2}-3\xi\mu^{2}}}\,,$
$\displaystyle x$ $\displaystyle\equiv$ $\displaystyle\frac{\phi_{\rm
B}}{\phi_{\rm
F}}=\sqrt{\frac{2-\sqrt{1-3\mathsf{w}_{0}^{2}}}{2+\sqrt{1-3\mathsf{w}_{0}^{2}}}}\,,\mathsf{w}_{0}=\frac{\omega_{0}}{\sqrt{\mu^{2}-\omega_{0}^{2}}}\,.$

In this parameterization, the limits $x=1/\sqrt{3}$ ($\mathsf{w}_{0}=0$) and
$x=1$ ($\mathsf{w}_{0}=1/\sqrt{3}$) interpolate between the degenerate vacua
case and the scenario where the potential develops an exact stationary
inflection point. Allowing for $x<1/\sqrt{3}$ (and hence imaginary $w_{0}$)
makes $\varphi_{\rm F}$ a true vacuum, but we will not consider that scenario
in this work101010See Sakai and Sasaki (2008) for the flat space-time case and
Tamaki and Sakai (2011a) for the gravitating case.. In principle, Coleman’s
(necessary) stability criterion (c.f. Sec. II) $w_{0}\leq w<1$ allows even for
the cases where the second minimum disappears $\mathsf{w}_{0}>1/\sqrt{3}$
$(w_{0}>1/2)$ and the parametrization given in Eq. (96) is not applicable.
However, such configurations are in the deep thick wall regime, as we will
argue below. Thus, defining the class of non-topological solitons that we have
examined in this work by the presence of the false/degenerate vacuum in the
potential is not completely rigorous. However, not only does that describe the
largest part of the parameter space, but also the region where the
phenomenology of these objects differs most significantly from “normal
matter”.

We can define (as in Heeck  _et al._ (2021a))

$\varkappa^{2}=\frac{w^{2}-w^{2}_{0}}{1-w_{0}^{2}}\,,$ (97)

so that the Minkowski Klein-Gordon equation has the same form as in Sec. II.1
if one substitutes $w\to\varkappa$ and
$\mathsf{r}\to\sqrt{\mu^{2}-\omega_{0}^{2}}r$. One can thus use both analytic
and numerical results for the scalar profile of Q-balls with the simplest
potential. The macroscopic properties $R,M,Q...$, however, depend explicitly
on $w_{0}$ (see Heeck  _et al._ (2021a) for the relevant expressions). As
$w_{0}$ increases from $0$ to $1$, the length scale of the boundary
$\propto(1-w^{2}_{0})^{-1/2}$ increases even if $\varkappa\ll 1$. Thus, as the
false vacuum departs from $0$, the thick-wall regime increasingly dominates
the Q-ball behaviour.

In curved spacetime, the Klein-Gordon equation is not invariant under the
above reparametrization111111This is expected from the equivalence principle:
in Minkowski space it is enough to know the energy difference between the two
vacua, while in GR information about both vacua is needed.. Instead, the
system (38) - (40) can be formulated as:

$\displaystyle\frac{1}{\mathsf{r}^{2}}\left(\mathsf{r}\,e^{-u}\right)^{\prime}-\frac{1}{\mathsf{r}^{2}}=-\frac{\Lambda^{2}}{2}\times$
(98)
$\displaystyle\Big{[}e^{-v}\varkappa^{2}\varphi^{2}+e^{-u}(\varphi^{\prime})^{2}+\varphi^{2}(1-\varphi^{2})^{2}+\mathsf{w}_{0}^{2}\varphi^{2}(1+e^{-v})\Big{]}\,,$
$\displaystyle
e^{-u}\left(\frac{v^{\prime}}{\mathsf{r}}+\frac{1}{\mathsf{r}^{2}}\right)-\frac{1}{\mathsf{r}^{2}}=\frac{\Lambda^{2}}{2}\times$
(99)
$\displaystyle\Big{[}e^{-v}\varkappa^{2}\varphi^{2}+e^{-u}(\varphi^{\prime})^{2}-\varphi^{2}(1-\varphi^{2})^{2}-\mathsf{w}_{0}^{2}\varphi^{2}(1-e^{-v})\Big{]}\,,$
$\displaystyle\varphi^{\prime\prime}+\left(\frac{2}{\mathsf{r}}+\frac{v^{\prime}-u^{\prime}}{2}\right)\varphi^{\prime}=$
(100) $\displaystyle
e^{u}\left[(1-4\varphi^{2}+3\varphi^{4})-\varkappa^{2}e^{-v}+\mathsf{w}_{0}^{2}(1-e^{-v})\right]\varphi\,,$

where

$\mathsf{r}=\sqrt{1-w_{0}^{2}}\,\mu
r\,,\bar{m}(\mathsf{r})=\frac{\sqrt{1-w_{0}^{2}}\mu m(r)}{m^{2}_{\rm
Pl}}\,,\Lambda=\frac{\sqrt{2}\phi_{0}}{M_{\rm Pl}}$

and the other conventions from Eqs. (95), (96) and (97) apply, while the prime
′ denotes spatial derivatives with respect to $\mathsf{r}$. In the Minkowski
limit, one has $u\to 0,v\to 0$, and the explicit dependence on
$\mathsf{w}_{0}^{2}$ disappears in the Klein-Gordon equation.

Self-gravitating configurations with the potential (94) have been considered
for particular values of the coefficients in Kleihaus  _et al._ (2005); Tamaki
and Sakai (2011a); Kleihaus  _et al._ (2012); Siemonsen and East (2021). The
parameterization outlined here allows us to perform a systematic exploration
of the parameter space, by varying $w_{0}$ from $0$ to $1$ in discrete steps
using the same numerical approach as in the rest of this work (App. A). In
Fig. 13 (left) we show how the compactness decreases (for fixed $\Lambda$)
from the most compact configurations $w_{0}=0$ to $w_{0}\to 1$. This is
represented in the $w-\varphi$ parameter space by the growth of the height of
the horizontal asymptote, or by that of the tipping point of the two branches
in the $\tilde{\omega}-\varphi$ parameter space, as demonstrated in Fig. 13
(right).

The interpretation of these results is straightforward: increasing the height
of the false vacuum implies thicker walls and hence larger minimal
frequency/smaller maximal mass and compactness. The picture outlined in Sec.
III.5, where the analogue particle in the $\tilde{v}(0)=0$ gauge does not
initially feel the presence of gravity, is valid also for the general sextic
potential. Henceforth, the curve $\tilde{w}-\varphi$ in the thin-wall regime
is given by Eq. (16) (with $w\to\varkappa$):

$\tilde{w}^{2}=(1-4\varphi^{2}_{c}+3\varphi^{4}_{c})(1-w^{2}_{0})+w^{2}_{0}\,.$
(101)

The compactness dependence on $w_{0}$ can also be understood in terms of the
LinEoS. Using an arguments analogous to those that lead to Eq. (50) one gets,
for the general sextic potential, the following estimate of the speed of sound

$\displaystyle(c^{2}_{s})_{a}=\frac{\varphi_{c}^{2}\left(w_{0}^{2}\left(3\varphi_{c}^{2}-4\right)+6\varphi_{c}^{2}-4\right)}{2+4\left(w_{0}^{2}-3\right)\varphi_{c}^{2}-3\left(w_{0}^{2}-4\right)\varphi_{c}^{4}}\,.$
(102)

This equation predicts that even for $\varphi_{c}=1$, the speed of sound will
be subluminal when $w_{0}\neq 0$. For example, the compactness of the maximum
mass configuration with $w_{0}=0.24$ is numerically found to be
$\mathcal{C}=0.321$, while Eqns. (102) and (49) predict a similar value for
the corresponding $\varphi_{c}=1.083$: $\mathcal{C}_{\rm{LinEoS}}=0.326$.

In agreement with the discussion in Sec. IV.2, for
$\Lambda\,\mbox{{\scriptsize$\stackrel{{\scriptstyle>}}{{\sim}}$}}\,1$ the
configurations exhibit a MBS-like behaviour, irrespective of the value of
$w_{0}$.

Figure 13: Compactness (left) and scalar frequency (right) of SBSs for the
general sextic potential (94), in the $v(\infty)=0$ gauge, as a function of
the central field. Various values of the parameter $\omega_{0}$ that describes
the deviation from the degenerate vacuum case are considered, and we set
$\Lambda=0.186$, comparing also with the vanila case $\omega_{0}/\mu=0$.

###  IV.5 Non-polynomial quartic potential

As one last departure from the benchmark scenario of this work, we will now
consider a real scalar field $\phi$ with the renormalizable potential

$\displaystyle V_{4}(\phi)=\mu^{2}\phi^{2}-g\phi^{3}+\lambda\phi^{4}\,,$ (103)
$\displaystyle g>0\,,\lambda>0\,.$

Q-balls with this effective potential can form in the presence of other fields
Kusenko (1997a); Postma (2002); Bošković  _et al._ , or we can consider this
scenario as a proxy for a pseudo-soliton composed of real scalars. Formally
(and in line with the rest of this work) we will take this model to originate
from the non-polynomial potential of a $U(1)$ complex scalar

$V_{4}(|\Phi|)=\mu^{2}|\Phi|^{2}-g(|\Phi|^{2})^{3/2}+\lambda|\Phi|^{4}\,,$
(104)

with the ansatz (3) [$\Phi=\phi(r)e^{-i\omega t}$] giving a real scalar $\phi$
in Eq. (103).

Like for the generic sextic potential of Sec. IV.4, the potential of Eq. (103)
admits false and degenerate vacua, and can be parameterized as a deviation
from the degenerate case

$\displaystyle
V_{4}=\phi^{2}_{0}\Big{[}(\mu^{2}-\omega^{2}_{0})\varphi^{2}(1-\varphi)^{2}+\omega^{2}_{0}\varphi^{2}\Big{]}\,,$
(105)
$\displaystyle\varphi=\frac{|\Phi|}{\phi_{0}}\,,\,\phi_{0}=\frac{g}{2\lambda}\,,\,w_{0}=\frac{\omega_{0}}{\mu}=\sqrt{1-\frac{g^{2}}{4\lambda\mu^{2}}}\,,$

or parameterizing the two vacua

$\displaystyle V_{4}$ $\displaystyle=$
$\displaystyle\frac{\lambda\phi^{4}_{\rm
F}}{3}\varphi_{x}^{2}\left(6x-4(x+1)\varphi_{x}+3\varphi_{x}^{2}\right)\,,$
(106) $\displaystyle\varphi_{x}$ $\displaystyle\equiv$
$\displaystyle\frac{\phi}{\phi_{\rm B}}\,,\phi_{\rm
F/B}=\frac{3g\pm\sqrt{9g^{2}-32\lambda\mu^{2}}}{8\lambda}\,,\frac{1}{\mu^{2}\phi^{2}_{\rm
F}}\frac{\lambda\phi^{4}_{\rm F}}{3}=\frac{1}{6x}\,,$ $\displaystyle x$
$\displaystyle\equiv$ $\displaystyle\frac{\phi_{\rm B}}{\phi_{\rm
F}}=\frac{3-\sqrt{1-8\mathsf{w}_{0}^{2}}}{3+\sqrt{1-8\mathsf{w}_{0}^{2}}}\,,\mathsf{w}_{0}=\frac{\omega_{0}}{\sqrt{\mu^{2}-\omega_{0}^{2}}}\,.$

In contrast to the $V_{6}$ case in flat spacetime, all configurations are
(classically) stable provided that121212Quantum effects can influence
stability in part of the parameter space for small Q-balls Postma (2002). Note
that stable solutions (under small perturbations) can exist also for
$w^{2}_{0}<0$ Paccetti Correia and Schmidt (2001); Sakai and Sasaki (2008).
$w^{2}_{0}\geq 0$ Kusenko (1997b); Paccetti Correia and Schmidt (2001); Sakai
and Sasaki (2008). Otherwise, the discussion is similar to the generic $V_{6}$
case: for $1\gg w^{2}\sim w^{2}_{0}$ Q-balls are in the thin wall regime,
while for $w\sim 1$ they are in the thick wall one.

The gravitating case for this potential was considered in Tamaki and Sakai
(2010), but only for one value of (in our parametrization) $w_{0}\approx 0.4$,
which is in the intermediate thick wall regime. Expecting a similar
phenomenology with respect to $w_{0}\neq 0$ as in Sec. IV.4, we have focused
only on the $w_{0}=0$ case (degenerate vacua), noting that the above
parametrization allows for a straightforward systematic exploration with
respect to the height of the false vacuum.

In the thin-wall gravitating limit, the relations among the macroscopic
parameters are very similar to those of SBSs with $V_{6}$ [Eq. (1)], as can be
seen in Fig. 14 for the mass and Fig. 15 for the compactnesses. This should
come as no surprise as the presence of the non-trivial vacuum renders the
equation of state linear in the thin-wall regime and the same arguments from
Sec. III.4 apply. For example, the estimate of the speed of sound, in analogy
with the Eq. (50), is

$\displaystyle(c^{2}_{s})_{a}\approx\frac{\varphi_{c}(4\varphi_{c}-3)}{12\varphi_{c}^{2}-15\varphi_{c}+4}\,.$
(107)

From the above and the LinEoS results (49) we can predict $\mathcal{C}_{\rm
LinEoS}=0.337$ for the compactness of the maximaum mass configuration with
$\Lambda=0.186$ ($\varphi_{c}=1.053$), whose true value is
$\mathcal{C}=0.331$.

The presence of only one turning point in the $\bar{M}-\varphi_{c}$ diagram
[Fig. 14] indicates that only one stable branch is inherited from flat
spacetime, in contrast with the $V_{6}$ case. This branch is, as in the
$V_{6}$ case, succeeded by the compact unstable branch, as elaborated using
catastrophe theory arguments in Tamaki and Sakai (2010). The fact that both
SBSs with a quartic potential and MBSs have only one stable and one unstable
branch, noted in Tamaki and Sakai (2010), is in fact accidental, as the MBS
stable zone originates from the quantum pressure, while the stable zone of
SBSs with the $V_{4}$ potential is inherited from the corresponding Q-balls.
This is the reason why the $\bar{M}-\varphi_{c}$ relation in the low
compactness stable branch for the sextic potential is described by the NBS
scaling (90) in Fig. 14, while the quartic is not. This difference can be also
illustrated by the $w-\varphi$ diagram in Fig. 16: for the quartic potential
one can observe a much sharper decline of the $w(\varphi_{c})$ curve from
$w=1$ than with the $V_{6}$ potential. Convergence between the two potentials
occurs, manifest in both $\bar{M}-\varphi_{c}$ and $w-\varphi_{c}$
representations (Figs. 14, 16), in the compact limit, where gravity does not
discriminate between the highest powers of the scalar potential. In agreement
with Sec. III.5, the $\tilde{w}(\varphi_{c})$ curve in the thin wall regime
matches well the Q-ball result:

$\tilde{w}=\sqrt{1-3\varphi_{c}+2\varphi_{c}^{2}}\,.$ (108)

Figure 14: Mass vs central field of SBSs in the case of the quartic potential
(103), compared with the benchmark case (1) for the same control parameter
$\Lambda=0.186$. We also show the NBS scaling (90) (black, dashed). Figure 15:
Compactness vs central field of SBSs in the case of the quartic potential
(103), compared with the benchmark case (1) for the same control parameter
$\Lambda=0.186$. The causal Buchdahl bound $\mathcal{C}_{\rm B+C}$ is
indicated as a gray and dashed line. Figure 16: SBSs scalar frequency as a
function of the central field in the case of the quartic potential (103),
compared with the benchmark case (94) for $\Lambda=0.186$. Q-ball (with the
quartic potential) results, both numerical and analytical, are indicated with
the purple triangles and the purple line, respectively.

##  V Conclusion

In this work we have provided a comprehensive analysis of the self-gravitating
solitonic objects made of complex scalars that obey potentials with
false/degenerate vacua. We have built on previous studies by reinterpreting
and improving them, and also by providing novel results. In more detail, we
find that in the thin wall regime, because of the presence of the non-trivial
vacuum, these objects can achieve high compactness, saturating the Buchdahl
limit with the causality constraint $\mathcal{C}_{\rm B+C}=0.354$ (Sec.
III.4). These values are the highest that have been found so far for motivated
ECO models. This results also provides a new perspective on SBSs - objects
with an incredibly stiff equation of state, where the information on local
disturbances is transmitted at almost the speed of light.

We have established the robustness of this picture by checking various
potentials considered in the literature - general sextic, quartic and cosine
potentials (Sec. IV). Although in this work we have stressed the general
features of SBSs, particular models considered in the literature present some
differences. For the ease of navigating amongst different models, we summarise
in Table 2 the potentials that we have considered in this work. In addition,
the bubble-like (in the thin-wall regime) behaviour of the scalar field allows
for an analytic description of these configurations, which we have presented
for the simplest case of degenerate vacua (Sec. III.5), although we expect
that this description can be extended (in a straightforward manner) to the
other potentials considered in Sec. IV. This analytic solution can be used as
an approximation to the numerical one, and it has helped us obtain analytic
control of SBS solutions for arbitrarily small values of $\sigma_{0}/M_{\rm
Pl}$.

Boson star model | Potential | $M_{\rm Pl}\to\infty$ | $|\Phi|\ll M_{\rm Pl}$ | $\mathcal{C}_{\rm max}$
---|---|---|---|---
Mini (MBS) | $\mu^{2}|\Phi|^{2}$ | / | Newtonian (NBS) | 0.11
Self-interacting (SIBS) | $\mu^{2}|\Phi|^{2}+|\lambda||\Phi|^{4}$ | / | NBS: self-interacting | 0.16
Soliton (SBS): simplest | $\mu^{2}|\Phi|^{2}\Big{(}1-2\frac{|\Phi|^{2}}{\sigma_{0}^{2}}\Big{)}^{2}$ | Q-ball: simplest | NBS/MBS | $\,\mbox{{\scriptsize$\stackrel{{\scriptstyle<}}{{\sim}}$}}\,0.354$
SBS: sextic | $\mu^{2}|\Phi|^{2}-|\beta||\Phi|^{4}+|\xi||\Phi|^{6}$ | Q-ball: sextic | NBS/MBS | $\,\mbox{{\scriptsize$\stackrel{{\scriptstyle<}}{{\sim}}$}}\,0.354-0.06w_{0}^{2}$
SBS: cosine [“axion BS”] | $|\alpha|\Big{[}1-\sqrt{1-|\beta|\sin^{2}\Big{(}\frac{|\Phi|}{\sigma_{0}}\Big{)}}\Big{]}$ | Q-ball: cosine | NBS/MBS | $\,\mbox{{\scriptsize$\stackrel{{\scriptstyle<}}{{\sim}}$}}\,0.354$
SBS: quartic | $\mu^{2}|\Phi|^{2}-|g|(|\Phi|^{2})^{3/2}+|\lambda||\Phi|^{4}$ | Q-ball: quartic | Q-ball: quartic | $\,\mbox{{\scriptsize$\stackrel{{\scriptstyle<}}{{\sim}}$}}\,0.354-0.2w^{2}_{0}$
Table 2: Models considered in this work, their potentials and main properties
(including the non-gravitating limit $M_{\rm Pl}\to\infty$, the low
compactness limit $|\Phi|\ll M_{\rm Pl}$ and the maximal compactness
$\mathcal{C}_{\rm max}$). For the quartic SBS we have quoted $\mathcal{C}_{\rm
max}$ for the general potential (105) and both general quartic and sextic
potential expressions are valid for small $w_{0}$. Value for $\mathcal{C}_{\rm
max}[\rm{SIBS}]$ is taken from Amaro-Seoane  _et al._ (2010).

In the low compactness limit, the configurations are stabilized either by
gravity (thus behaving as MBS or NBS, cf. Sec. IV.2), or by self-interactions
(like in the case of the quartic potential, cf. Sec. IV.5). For field values
close to the Planck scale, the compact stable branch shrinks, and only the
low-compactness one, described by the MBS model, survives (cf. Sec. IV.1).

One follow up to this project is to generalize our formalism to study scalar-
fermion solitonic configurations, which could form through different channels
and possibly develop other interesting phenomenology Lee (1987b); Bahcall  _et
al._ (1989); Hong  _et al._ (2020); Gross  _et al._ (2021). This work is
currently in progress and we will present it elsewhere Bošković  _et al._ .
There are several other topics that are a natural continuation of this work
and could be of interest: investigating the applicability of our results for
non-Abelian global or gauged gravitating Q-balls or higher-spin fields;
comparing properties of (gravitating) Q-balls to configurations made of real
bosons , motivated by axion or vector DM Hui  _et al._ (2017); Hui (2021),
with a false/degenerate vacuum in the potential; obtaining analytic control
and the parameter space behaviour of the SBS quasi-normal modes Macedo  _et
al._ (2013) and tidal Love numbers Cardoso  _et al._ (2017); Sennett  _et al._
(2017); probing the I-Love-Q and other universal relations Yagi and Yunes
(2013, 2017); describing GW echoes Urbano and Veermäe (2019) etc.

From the astro-particle physics perspective, it is also important to
understand what class of Q-ball production mechanisms can lead to a
significant fraction of compact objects, and how can the presence or absence
of SBS signatures in present and future GW detectors help us constrain their
formation mechanism and possibly gain information on the early universe. On
the phenomenological side, it is therefore imperative to systematically
explore the behaviour of SBSs in binaries and their GW signatures Palenzuela
_et al._ (2017), also for low-$\mathcal{C}$ configurations Chia and Edwards
(2020). Work addressing the last topic is reported elsewhere Bezares  _et al._
(2022).

###### Acknowledgements.

We acknowledge financial support provided under the European Union’s H2020 ERC
Consolidator Grant “GRavity from Astrophysical to Microscopic Scales” grant
agreement no. GRAMS-815673. This work was supported by the EU Horizon 2020
Research and Innovation Programme under the Marie Sklodowska-Curie Grant
Agreement No. 101007855. We would also like to thank Alfredo Urbano for an
initial collaboration on these topics, Miguel Bezares for useful conversations
on boson stars and for checking our numerical solutions with an independent
evolution code and Paolo Pani and Steve Liebling for comments on the
manuscript.

##  Appendix A Numerical solutions of Einstein-Klein-Gordon system

In the following, we present the technical details regarding the well-
posedness and numerical solution of the Einstein-Klein-Gordon system given by
Eqs. (38), (39) and (40), which we reproduce here for clarity:

$\displaystyle\frac{1}{r^{2}}\left(r\,e^{-u}\right)^{\prime}-\frac{1}{r^{2}}$
$\displaystyle=$ $\displaystyle-\frac{1}{M^{2}_{\rm Pl}}\rho\,,$ (109)
$\displaystyle
e^{-u}\left(\frac{v^{\prime}}{r}+\frac{1}{r^{2}}\right)-\frac{1}{r^{2}}$
$\displaystyle=$ $\displaystyle\frac{1}{M^{2}_{\rm Pl}}P_{\text{rad}}\,,$
(110)
$\displaystyle\phi^{\prime\prime}+\left(\frac{2}{r}+\frac{v^{\prime}-u^{\prime}}{2}\right)\phi^{\prime}$
$\displaystyle=$ $\displaystyle
e^{u}\left({\frac{dV}{d|\Phi|^{2}}}-\omega^{2}e^{-v}\right)\phi.$ (111)

First, let us notice that the structure equations have poles at $r=0$, and we
thus have to impose regularity there. In addition, we require solutions to be
asymptotically flat. The resulting boundary eigenvalue problem uniquely
determines the eigenvalue $\omega$. In fact, on the one hand, a local
expansion of the fields around $r=0$ up to $\mathcal{O}(\epsilon^{4})$ yields

$\displaystyle u(\epsilon)$ $\displaystyle\approx$ $\displaystyle
0+\frac{1}{6}\Lambda^{2}\epsilon^{2}\varphi_{c}^{2}\left(1+e^{-v_{c}}w^{2}+\varphi_{c}^{4}-2\varphi_{c}^{2}\right)\,,$
(112) $\displaystyle v(\epsilon)$ $\displaystyle\approx$ $\displaystyle
v_{c}-\frac{1}{6}\Lambda^{2}\epsilon^{2}\varphi_{c}^{2}\left(1-2e^{-v_{c}}w^{2}+\varphi_{c}^{4}-2\varphi_{c}^{2}\right)\,,$
(113) $\displaystyle\varphi(\epsilon)$ $\displaystyle\approx$
$\displaystyle\varphi_{c}+\frac{1}{6}\epsilon^{2}\varphi_{c}\left(1-e^{-v_{c}}w^{2}+3\varphi_{c}^{4}-4\varphi_{c}^{2}\right)\,,$
(114)

where we can set $\tilde{v}_{c}=0$ and $w\to\tilde{w}$ with a rescaling of the
time coordinate. On the other hand, the leading order asymptotic behaviour, as
discussed in Sec. III.5.1, is given by Eqs. (53), (54), (57), which we
reproduce here for the reader’s convenience:

$\displaystyle u_{>}$ $\displaystyle=$
$\displaystyle-\log{\Big{(}1-\frac{2\bar{M}}{\mathsf{r}}\Big{)}}\,,$ (115)
$\displaystyle\tilde{v}_{>}$ $\displaystyle=$
$\displaystyle\tilde{v}_{\infty}+\log{\Big{(}1-\frac{2\bar{M}}{\mathsf{r}}\Big{)}}\,$
(116) $\displaystyle\varphi_{\infty}$ $\displaystyle\simeq$
$\displaystyle\frac{\mathcal{A}}{\mathsf{r}^{1+\beta_{>}}}e^{-\alpha_{>}\mathsf{r}}\,,$
(117) $\displaystyle\alpha_{>}$ $\displaystyle=$
$\displaystyle\sqrt{1-e^{-\tilde{v}_{\infty}}\tilde{w}^{2}}\,,$ (118)
$\displaystyle\beta_{>}$ $\displaystyle=$
$\displaystyle\frac{\bar{M}}{\alpha_{>}}(1-2e^{-\tilde{v}_{\infty}}\tilde{w}^{2})\,,$
(119)

One can then match the numerical solution from the interior, obtained with the
initial conditions determined by (112)-(114) at some small but non-zero radius
$\epsilon$, with the asymptotic expansion at the infinity, at a finite but
sufficiently large matching radius (“direct shooting”). This can be done by
solving the four junction conditions

$\displaystyle\Delta u|_{\bullet}=0\,,\Delta
v|_{\bullet}=0\,,\Delta\varphi|_{\bullet}=0\,,\Delta\varphi^{\prime}|_{\bullet}=0\,,$
(120)

where $\Delta x\equiv x_{>}-x_{<}$ and $\mathsf{r}_{\bullet}$ is the matching
radius, in the unknowns $\varphi_{c}$, $\tilde{w}$, $\mathcal{A}$ and
$\tilde{v}_{\infty}$. Numerical integrations were performed using
Mathematica’s Inc.  default stiff solver. The stiffness of the system in the
thin wall regime requires an extraordinary amount of fine tuning for the
eigenvalues. Examples of the precision levels needed in order to produce
compact configurations with a shooting method are given in Macedo  _et al._
(2013).

Note that this procedure is applicable to SBSs, MBSs, to the potentials
considered in Sec. IV, and by taking $v\to 0$, $u\to 0$ also to Q-balls. (The
asymptotics of Q-balls are discussed in Sec. II.1). Note that the local
expansion (112)-(114) is given for the SBS model and is different for the
other potentials. We have validated our results by successfully reproducing
the MBS and SBS configurations of Macedo  _et al._ (2013) and by verifying
that the static configurations that we find do not change when used as initial
data in the evolution code of Bezares  _et al._ (2022). The radius of our
solutions is found by inverting $\bar{M}(\mathsf{R})=0.99\bar{M}$ (using
bisection) and the Noether charge is calculated numerically by integrating Eq.
(47). An alternative method for numerical calculation is outlined in Kleihaus
_et al._ (2005); Siemonsen and East (2021).

##  Appendix B Definitions of boson star radius

Figure 17: Relative difference between $\mathcal{C}$ calculated with the
benchmark definition of the radius enclosing 99% of the mass and
$\mathcal{C}_{x}=M/R_{x}$, where $R_{x}=\\{R_{95},R^{\ast},\langle
r\rangle,\sqrt{\langle r^{2}\rangle}\\}$.

As the radius of boson stars is not well defined, various definitions have
been used in the literature Schunck and Mielke (2003). Besides the definition
used in this work (corresponding to 99% enclosed mass), the radius $R_{95}$
enclosing 95% of the total mass has also been used Barranco  _et al._ (2021),
as well as the following moments of the density distribution:

$\displaystyle\langle r\rangle=\frac{1}{M}\int^{\infty}_{0}r\rho(r)dV\,,$
(121) $\displaystyle\sqrt{\langle
r^{2}\rangle}=\sqrt{\frac{1}{M}\int^{\infty}_{0}r^{2}\rho(r)dV}\,.$ (122)

Some authors have also considered radii enclosing a given fraction (95% or
99%) of the total Noether charge $Q$ Palenzuela  _et al._ (2017), or the
moments of $j^{t}$ Kleihaus  _et al._ (2012). We will not discuss these
definitions, because in this paper we are interested in the behaviour of
Q-balls/SBSs as compact objects, for which energy density based radii are more
relevant. Finally, the inflection point $R^{\ast}$ was also taken as a Q-ball
radius in Heeck  _et al._ (2021a, b).

In Fig. 17, we show the difference between our benchmark definition of
compactness $\mathcal{C}=M/R$ and the compactness calculated with
(respectively) $R_{95}$ , $R^{\ast}$ , $\langle r\rangle$ ,$\sqrt{\langle
r^{2}\rangle}$ for $\Lambda=0.186$. The cutoff for the numerical integrals in
Eqs. (121) and (122) was taken to be the domain of the integration. As can be
seen, while $R_{95}$ and $R^{\ast}$ produce relative differences in the value
of the compactness in the compact stable branch of respectively $\sim
10^{-1}\,-\,10^{-2}$ and $\sim 3\cdot 10^{-1}\,-\,2\cdot 10^{-2}$ for
$\mathcal{C}\,\mbox{{\scriptsize$\stackrel{{\scriptstyle>}}{{\sim}}$}}\,0.05$,
$\langle r\rangle$ and $\langle r^{2}\rangle$ yield $\sim 0.2\,-\,0.7$
relative differences. This difference can be understood in the following way:
in the flat space-time and $\omega\to 0$ limit, the scalar profile approaches
a step function and $R^{\ast}$ approaches the hard surface radius. Neglecting
the surface tension and the potential in that limit, one finds $\langle
r\rangle\approx(3/4)R^{\ast}$ and $\sqrt{\langle
r^{2}\rangle}\approx\sqrt{(3/5)}R^{\ast}$.

##  Appendix C Analytic construction of Soliton Boson stars: technical details

Here we provide some additional information on the analytic construction of
SBS.

###  C.1 $u_{B}$ metric coefficient on the boundary

On the basis of the discussion from Sec. III.5.3 we find the $\ln(g_{tt})$
metric coefficient jump:

$\displaystyle
u_{B}=-2\mathsf{m}_{\ast}\mathsf{R}_{\ast}+v_{<}(\mathsf{R^{\ast}})-\log{\Big{[}\frac{\Lambda^{2}}{r}(e^{-2\mathsf{m}_{\ast}\mathsf{R}_{\ast}+v_{<}(\mathsf{R^{\ast}})}\\{\chi_{0}+\chi_{1}r+\chi_{2}r^{2}+\chi_{2}r^{3}\\}+c_{\rm
int})\Big{]}}\,,$ (123)
$\displaystyle\chi_{0}=\frac{1}{24}\left(-\frac{3e^{-v_{<}(\mathsf{R^{\ast}})}w^{2}\text{Li}_{3}\left(-2e^{2\mathsf{m}_{\ast}(r-\mathsf{R}_{\ast})}\right)}{\mathsf{m}_{\ast}^{3}}+\frac{3\text{Li}_{2}\left(-2e^{2\mathsf{m}_{\ast}(r-\mathsf{R}_{\ast})}\right)}{\mathsf{m}_{\ast}^{3}}+\frac{3\log\left(2e^{2\mathsf{m}_{\ast}(r-\mathsf{R}_{\ast})}+1\right)}{\mathsf{m}_{\ast}^{3}}\right)\,,$
$\displaystyle\chi_{1}=\frac{1}{24}\left(\frac{24}{\Lambda^{2}}+\frac{6e^{-v_{<}(\mathsf{R^{\ast}})}w^{2}\text{Li}_{2}\left(-2e^{2\mathsf{m}_{\ast}(r-\mathsf{R}_{\ast})}\right)}{\mathsf{m}_{\ast}^{2}}+\frac{6}{\mathsf{m}_{\ast}^{2}\left(2e^{2\mathsf{m}_{\ast}(r-\mathsf{R}_{\ast})}+1\right)}+\frac{6\log\left(2e^{2\mathsf{m}_{\ast}(r-\mathsf{R}_{\ast})}+1\right)}{\mathsf{m}_{\ast}^{2}}-\frac{6}{\mathsf{m}_{\ast}^{2}}\right)\,,$
$\displaystyle\chi_{2}=\frac{1}{24}\left(\frac{6e^{-v_{<}(\mathsf{R^{\ast}})}w^{2}\log\left(2e^{2\mathsf{m}_{\ast}(r-\mathsf{R}_{\ast})}+1\right)}{\mathsf{m}_{\ast}}-\frac{6}{\mathsf{m}_{\ast}\left(2e^{2\mathsf{m}_{\ast}(r-\mathsf{R}_{\ast})}+1\right)^{2}}+\frac{12}{2\mathsf{m}_{\ast}e^{2\mathsf{m}_{\ast}(r-\mathsf{R}_{\ast})}+\mathsf{m}_{\ast}}-\frac{6}{\mathsf{m}_{\ast}}\right)\,,$
$\displaystyle\chi_{3}=-\frac{1}{6}e^{-v_{<}(\mathsf{R^{\ast}})}w^{2}\,.$

The integration constant $c_{\rm int}$ is determined by matching the previous
solution to $u_{<}$ at $\mathsf{R_{<}}$ and $\text{Li}_{n}$ is the
polylogarithm function.

###  C.2 Details of the energy balance calculation

Here we provide details on obtaining analytic approximations in Sec. III.5.4 .
We will first approximate the complicated algebraic expressions (74) - (78).
The sum $B_{\mathcal{E}}+C_{\mathcal{E}}+D_{\mathcal{E}B}$ is subleading, as
the dominant contribution, proportional to
$\mathsf{m}^{2}_{\ast}\Lambda^{2}\mathsf{R^{\ast}}^{2}$, cancels out. In the
Q-ball limit, one also has $\mathsf{m}^{2}_{\ast}\sim 1$, which suppresses the
term proportional to $(1-\mathsf{m}^{2}_{\ast})$. In the more compact branch,
we have numerically established that
$(1-\mathsf{m}^{2}_{\ast})\,\mbox{{\scriptsize$\stackrel{{\scriptstyle<}}{{\sim}}$}}\,2.5$;
however, in that regime the terms $\sim 1/\mathsf{R^{\ast}}$ are subdominant
relative to the volume originating terms in $D_{\mathcal{E}\,<}$, because
$(\mathsf{R^{\ast}}\Lambda\tilde{w})^{2}\sim 4$, while
$1/\mathsf{R^{\ast}}\sim\Lambda\tilde{w}/2\ll 1$ . Finally, the contributions
with $w^{2}e^{-v_{\ast}}$ are suppressed, both due to $v_{\ast}>0$ and to the
numerical pre-factors.

Leaving only the terms $A_{\mathcal{E}}+D_{\mathcal{E}\,<}$ on the right hand
side of our master equation Eq. (73), approximating
$\mathsf{m}_{\ast}\approx\exp{(u_{<}(\mathsf{R^{\ast}})/2)}$ and expanding in
$\Lambda$, we the find Eq. (82) which we reproduce here

$\displaystyle\tilde{w}\approx\Lambda\frac{\sqrt{-\frac{T^{4}}{5}+6T^{2}+36}\left(T^{4}+10T^{2}+30\right)}{T\left(-T^{4}+30T^{2}+180\right)}\,.$
(124)

Inverting this equation, we find complicated expressions for the two branches,
which can be approximated as

$T\approx
T(\tilde{w}_{\cup}/\Lambda)\pm\sqrt{3.1\,(\tilde{w}-\tilde{w}_{\cup})/\Lambda}\,,\tilde{w}\sim\tilde{w}_{\cup}$
(125)

where $T(\tilde{w}_{\cup}/\Lambda)\approx 1.6$. Using Eq. (125), the simple
expectation for $v_{\infty}\approx
v_{<}(\mathsf{R^{\ast}})+u_{<}(\mathsf{R^{\ast}})$ allows us to find an
approximate behaviour for $w$, given in Eq. (83).

###  C.3 On the errors of estimating the radius

In Section III.5.5, the semi-analytic calculation develops a few percent
systematic error near the maximum mass. There are two reasons for this. First,
the expansion of $u$ and $v$ in the interior [expressions (60)-(61)] works
well only for small $\tilde{w}$, i.e. for configurations similar to $\rm{II}$.
In the thick-wall regime, the perturbative expansion outlined in Section
III.5.2 does not work well by construction. Although the frequency $\tilde{w}$
is higher for the configurations similar to III, the field derivative is still
exponentially suppressed and physically these configurations have even thinner
walls than the ones similar to II. Thus, one can expect that some different
perturbative scheme could improve the analytic description of the metric
coefficients for the configurations similar to III.

The second reason for the systematic error is that Eq. (22), which describes
the deviation between $R$ and $R^{\ast}$, does the best job for the
configurations similar to II. For configurations similar to III, we have
established numerically that the choice

$\lambda_{\rm III}=\frac{4}{2}\frac{2.66}{\delta}\,$ (126)

is the most appropriate one, while for configurations similar to I:

$\lambda_{\rm I}=\frac{2}{2}\frac{2.66}{\delta}\,.$ (127)

This issue can be circumvented by performing the (cheap) numerical inversion
$4\pi\int^{R}_{0}dr\,r^{2}\rho(r)=0.99M$ around the trial value determined by
the semi-analytic algorithm, with a piecewise analytic approximation for
$\rho$ (obtained from the field and the metric coefficients approximants).
Finally, the limits of the analytic description of the thin-wall Q-ball-like
region in the parameter space (i.e. configurations similar to I) are discussed
in Heeck  _et al._ (2021a).

Regarding the compact unstable branch: the numerical calculations indicate
additional step-like features in the bulk of the scalar profile. This feature,
not accounted for in our analytical description, is probably the origin of the
increasing error in this branch.

## References

  * Taylor and Weisberg (1982) J. H. Taylor and J. M. Weisberg, Astrophys. J. 253, 908 (1982).
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 100, 104036 (2019), arXiv:1903.04467 [gr-qc] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, KAGRA, VIRGO), Astrophys. J. Lett. 915, L5 (2021b), arXiv:2106.15163 [astro-ph.HE] .
  * Giacconi  _et al._ (1962) R. Giacconi, H. Gursky, F. R. Paolini, and B. B. Rossi, Phys. Rev. Lett. 9, 439 (1962).
  * Bowyer  _et al._ (1965) S. Bowyer, E. T. Byram, T. A. Chubb, and H. Friedman, Science 147, 394 (1965).
  * Reid  _et al._ (1999) M. J. Reid, A. C. S. Readhead, R. C. Vermeulen, and R. N. Treuhaft, Astrophys. J. 524, 816 (1999), arXiv:astro-ph/9905075 .
  * Ghez  _et al._ (2000) A. Ghez, M. Morris, E. E. Becklin, T. Kremenek, and A. Tanner, Nature 407, 349 (2000), arXiv:astro-ph/0009339 .
  * Akiyama  _et al._ (2019a) K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L1 (2019a), arXiv:1906.11238 [astro-ph.GA] .
  * Akiyama  _et al._ (2019b) K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L6 (2019b), arXiv:1906.11243 [astro-ph.GA] .
  * Giudice  _et al._ (2016) G. F. Giudice, M. McCullough, and A. Urbano, JCAP 1610, 001 (2016), arXiv:1605.01209 [hep-ph] .
  * Cardoso and Pani (2019) V. Cardoso and P. Pani, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Bertone  _et al._ (2020) G. Bertone _et al._ , SciPost Phys. Core 3, 007 (2020), arXiv:1907.10610 [astro-ph.CO] .
  * Visser (1995) M. Visser,  _Lorentzian wormholes: From Einstein to Hawking_ (1995).
  * Mazur and Mottola (2001) P. O. Mazur and E. Mottola, (2001), arXiv:gr-qc/0109035 [gr-qc] .
  * Raposo  _et al._ (2019) G. Raposo, P. Pani, M. Bezares, C. Palenzuela, and V. Cardoso, Phys. Rev. D99, 104072 (2019), arXiv:1811.07917 [gr-qc] .
  * Cardoso  _et al._ (2017) V. Cardoso, E. Franzin, A. Maselli, P. Pani, and G. Raposo, Phys. Rev. D95, 084014 (2017), [Addendum: Phys. Rev.D95,no.8,089901(2017)], arXiv:1701.01116 [gr-qc] .
  * Toubiana  _et al._ (2021) A. Toubiana, S. Babak, E. Barausse, and L. Lehner, Phys. Rev. D 103, 064042 (2021), arXiv:2011.12122 [gr-qc] .
  * Barausse  _et al._ (2014) E. Barausse, V. Cardoso, and P. Pani, Phys. Rev. D89, 104059 (2014), arXiv:1404.7149 [gr-qc] .
  * Cardoso  _et al._ (2016a) V. Cardoso, E. Franzin, and P. Pani, Phys. Rev. Lett. 116, 171101 (2016a), [Erratum: Phys. Rev. Lett.117,no.8,089902(2016)], arXiv:1602.07309 [gr-qc] .
  * Cardoso  _et al._ (2016b) V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela, and P. Pani, Phys. Rev. D94, 084031 (2016b), arXiv:1608.08637 [gr-qc] .
  * Năstase (2019) H. Năstase,  _Classical Field Theory_ (Cambridge University Press, 2019).
  * Lee and Pang (1992) T. Lee and Y. Pang, Phys.Rept. 221, 251 (1992).
  * Helfer  _et al._ (2019) T. Helfer, J. C. Aurrekoetxea, and E. A. Lim, Phys. Rev. D99, 104028 (2019), arXiv:1808.06678 [gr-qc] .
  * Franzin  _et al._ (2018) E. Franzin, M. Cadoni, and M. Tuveri, Phys. Rev. D97, 124018 (2018), arXiv:1805.08976 [gr-qc] .
  * Jetzer (1992) P. Jetzer, Phys.Rept. 220, 163 (1992).
  * Schunck and Mielke (2003) F. Schunck and E. Mielke, Class.Quant.Grav. 20, R301 (2003), arXiv:0801.0307 [astro-ph] .
  * Liebling and Palenzuela (2017) S. L. Liebling and C. Palenzuela, Living Rev. Rel. 20, 5 (2017), arXiv:1202.5809 [gr-qc] .
  * Visinelli (2021) L. Visinelli, (2021), arXiv:2109.05481 [gr-qc] .
  * Kaup (1968) D. J. Kaup, Phys. Rev. 172, 1331 (1968).
  * Colpi  _et al._ (1986) M. Colpi, S. Shapiro, and I. Wasserman, Phys.Rev.Lett. 57, 2485 (1986).
  * Friedberg  _et al._ (1987a) R. Friedberg, T. D. Lee, and Y. Pang, Phys. Rev. D35, 3658 (1987a).
  * Bezares and Palenzuela (2018) M. Bezares and C. Palenzuela, Class. Quant. Grav. 35, 234002 (2018), arXiv:1808.10732 [gr-qc] .
  * Brito  _et al._ (2016a) R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, Phys. Lett. B752, 291 (2016a), arXiv:1508.05395 [gr-qc] .
  * Herdeiro  _et al._ (2019) C. Herdeiro, I. Perapechka, E. Radu, and Ya. Shnir, JHEP 02, 111 (2019), arXiv:1811.11799 [gr-qc] .
  * Henriques  _et al._ (1989) A. Henriques, A. R. Liddle, and R. Moorhouse, Phys.Lett. B233, 99 (1989).
  * Coleman (1985) S. R. Coleman, Nucl. Phys. B262, 263 (1985), [Erratum: Nucl. Phys.B269,744(1986)].
  * Freitas  _et al._ (2021) F. F. Freitas, C. A. R. Herdeiro, A. P. Morais, A. Onofre, R. Pasechnik, E. Radu, N. Sanchis-Gual, and R. Santos, (2021), arXiv:2107.09493 [hep-ph] .
  * Krippendorf  _et al._ (2018) S. Krippendorf, F. Muia, and F. Quevedo, JHEP 08, 070 (2018), arXiv:1806.04690 [hep-th] .
  * Barranco  _et al._ (2021) J. Barranco, J. Chagoya, A. Diez-Tejedor, G. Niz, and A. A. Roque,  (2021), 10.1088/1475-7516/2021/10/022, arXiv:2108.01679 [gr-qc] .
  * Seidel and Suen (1990) E. Seidel and W.-M. Suen, Phys. Rev. D42, 384 (1990).
  * Di Giovanni  _et al._ (2020) F. Di Giovanni, N. Sanchis-Gual, P. Cerdá-Durán, M. Zilhão, C. Herdeiro, J. A. Font, and E. Radu, Phys. Rev. D 102, 124009 (2020), arXiv:2010.05845 [gr-qc] .
  * Siemonsen and East (2021) N. Siemonsen and W. E. East, Phys. Rev. D 103, 044022 (2021), arXiv:2011.08247 [gr-qc] .
  * Dmitriev  _et al._ (2021) A. S. Dmitriev, D. G. Levkov, A. G. Panin, E. K. Pushnaya, and I. I. Tkachev, Phys. Rev. D 104, 023504 (2021), arXiv:2104.00962 [gr-qc] .
  * Seidel and Suen (1994) E. Seidel and W.-M. Suen, Phys. Rev. Lett. 72, 2516 (1994), arXiv:gr-qc/9309015 [gr-qc] .
  * Palenzuela  _et al._ (2017) C. Palenzuela, P. Pani, M. Bezares, V. Cardoso, L. Lehner, and S. Liebling, Phys. Rev. D96, 104058 (2017), arXiv:1710.09432 [gr-qc] .
  * Lee (1987a) T. D. Lee, Phys. Rev. D35, 3637 (1987a).
  * Friedberg  _et al._ (1987b) R. Friedberg, T. D. Lee, and Y. Pang, Phys. Rev. D35, 3640 (1987b).
  * Lee and Pang (1987) T. D. Lee and Y. Pang, Phys. Rev. D35, 3678 (1987).
  * Lee (1987b) T. D. Lee, Comments Nucl. Part. Phys. 17, 225 (1987b).
  * Bahcall  _et al._ (1989) S. Bahcall, B. W. Lynn, and S. B. Selipsky, Nucl. Phys. B 325, 606 (1989).
  * Hong  _et al._ (2020) J.-P. Hong, S. Jung, and K.-P. Xie, Phys. Rev. D 102, 075028 (2020), arXiv:2008.04430 [hep-ph] .
  * Gross  _et al._ (2021) C. Gross, G. Landini, A. Strumia, and D. Teresi, JHEP 09, 033 (2021), arXiv:2105.02840 [hep-ph] .
  * (56) M. Bošković, E. Barausse, and A. Urbano, Work in progress.
  * Cotner  _et al._ (2019) E. Cotner, A. Kusenko, M. Sasaki, and V. Takhistov, JCAP 1910, 077 (2019), arXiv:1907.10613 [astro-ph.CO] .
  * Amin and Mocz (2019) M. A. Amin and P. Mocz, Phys. Rev. D100, 063507 (2019), arXiv:1902.07261 [astro-ph.CO] .
  * Krylov  _et al._ (2013) E. Krylov, A. Levin, and V. Rubakov, Phys. Rev. D 87, 083528 (2013), arXiv:1301.0354 [hep-ph] .
  * Postma (2002) M. Postma, Phys. Rev. D 65, 085035 (2002), arXiv:hep-ph/0110199 .
  * Croon  _et al._ (2020) D. Croon, A. Kusenko, A. Mazumdar, and G. White, Phys. Rev. D 101, 085010 (2020), arXiv:1910.09562 [hep-ph] .
  * Rubakov and Gorbunov (2017) V. A. Rubakov and D. S. Gorbunov,  _Introduction to the Theory of the Early Universe: Hot big bang theory_ (World Scientific, Singapore, 2017).
  * Troitsky (2016) S. Troitsky, JCAP 11, 027 (2016), arXiv:1510.07132 [hep-ph] .
  * Helfer  _et al._ (2017) T. Helfer, D. J. E. Marsh, K. Clough, M. Fairbairn, E. A. Lim, and R. Becerril, JCAP 1703, 055 (2017), arXiv:1609.04724 [astro-ph.CO] .
  * Nazari  _et al._ (2021) Z. Nazari, M. Cicoli, K. Clough, and F. Muia, JCAP 05, 027 (2021), arXiv:2010.05933 [gr-qc] .
  * Muia  _et al._ (2019) F. Muia, M. Cicoli, K. Clough, F. Pedro, F. Quevedo, and G. P. Vacca, JCAP 1907, 044 (2019), arXiv:1906.09346 [gr-qc] .
  * Seidel and Suen (1991) E. Seidel and W. M. Suen, Phys. Rev. Lett. 66, 1659 (1991).
  * Bošković  _et al._ (2018) M. Bošković, F. Duque, M. C. Ferreira, F. S. Miguel, and V. Cardoso, Phys. Rev. D98, 024037 (2018), arXiv:1806.07331 [gr-qc] .
  * Kasuya  _et al._ (2003) S. Kasuya, M. Kawasaki, and F. Takahashi, Phys. Lett. B 559, 99 (2003), arXiv:hep-ph/0209358 .
  * Hertzberg (2010) M. P. Hertzberg, Phys. Rev. D82, 045022 (2010), arXiv:1003.3459 [hep-th] .
  * Zhang  _et al._ (2020) H.-Y. Zhang, M. A. Amin, E. J. Copeland, P. M. Saffin, and K. D. Lozanov, JCAP 07, 055 (2020), arXiv:2004.01202 [hep-th] .
  * Brito  _et al._ (2016b) R. Brito, V. Cardoso, C. F. B. Macedo, H. Okawa, and C. Palenzuela, Phys. Rev. D93, 044045 (2016b), arXiv:1512.00466 [astro-ph.SR] .
  * Marsh (2016) D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Hui (2021) L. Hui, (2021), arXiv:2101.11735 [astro-ph.CO] .
  * Schive  _et al._ (2014) H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Phys. 10, 496 (2014), arXiv:1406.6586 [astro-ph.GA] .
  * Levkov  _et al._ (2018) D. G. Levkov, A. G. Panin, and I. I. Tkachev, Phys. Rev. Lett. 121, 151301 (2018), arXiv:1804.05857 [astro-ph.CO] .
  * Widdicombe  _et al._ (2018) J. Y. Widdicombe, T. Helfer, D. J. E. Marsh, and E. A. Lim, JCAP 1810, 005 (2018), arXiv:1806.09367 [astro-ph.CO] .
  * Arvanitaki  _et al._ (2020) A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, Phys. Rev. D101, 083014 (2020), arXiv:1909.11665 [astro-ph.CO] .
  * Amin  _et al._ (2014) M. A. Amin, M. P. Hertzberg, D. I. Kaiser, and J. Karouby, Int. J. Mod. Phys. D24, 1530003 (2014), arXiv:1410.3808 [hep-ph] .
  * Niemeyer and Easther (2020) J. C. Niemeyer and R. Easther, JCAP 2007, 030 (2020), arXiv:1911.01661 [astro-ph.CO] .
  * Nugaev and Shkerin (2020) E. Ya. Nugaev and A. V. Shkerin, J. Exp. Theor. Phys. 130, 301 (2020), arXiv:1905.05146 [hep-th] .
  * Kesden  _et al._ (2005) M. Kesden, J. Gair, and M. Kamionkowski, Phys.Rev. D71, 044015 (2005), arXiv:astro-ph/0411478 [astro-ph] .
  * Macedo  _et al._ (2013) C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Phys. Rev. D88, 064046 (2013), arXiv:1307.4812 [gr-qc] .
  * Kleihaus  _et al._ (2005) B. Kleihaus, J. Kunz, and M. List, Phys. Rev. D72, 064002 (2005), arXiv:gr-qc/0505143 [gr-qc] .
  * Kleihaus  _et al._ (2012) B. Kleihaus, J. Kunz, and S. Schneider, Phys. Rev. D85, 024045 (2012), arXiv:1109.5858 [gr-qc] .
  * Guerra  _et al._ (2019) D. Guerra, C. F. B. Macedo, and P. Pani, JCAP 1909, 061 (2019), [Erratum: JCAP2006,no.06,E01(2020)], arXiv:1909.05515 [gr-qc] .
  * Sennett  _et al._ (2017) N. Sennett, T. Hinderer, J. Steinhoff, A. Buonanno, and S. Ossokine, Phys. Rev. D 96, 024002 (2017), arXiv:1704.08651 [gr-qc] .
  * Bezares  _et al._ (2017) M. Bezares, C. Palenzuela, and C. Bona, Phys. Rev. D 95, 124005 (2017), arXiv:1705.01071 [gr-qc] .
  * Helfer  _et al._ (2021) T. Helfer, U. Sperhake, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, (2021), arXiv:2108.11995 [gr-qc] .
  * Bezares  _et al._ (2022) M. Bezares, M. Bošković, S. Liebling, C. Palenzuela, P. Pani, and E. Barausse, (2022), arXiv:2201.06113 [gr-qc] .
  * Multamaki and Vilja (2002) T. Multamaki and I. Vilja, Phys. Lett. B 542, 137 (2002), arXiv:hep-ph/0205302 .
  * Tamaki and Sakai (2010) T. Tamaki and N. Sakai, Phys. Rev. D81, 124041 (2010), arXiv:1105.1498 [gr-qc] .
  * Tamaki and Sakai (2011a) T. Tamaki and N. Sakai, Phys. Rev. D83, 044027 (2011a), arXiv:1105.2932 [gr-qc] .
  * Tamaki and Sakai (2011b) T. Tamaki and N. Sakai, Phys. Rev. D84, 044054 (2011b), arXiv:1108.3902 [gr-qc] .
  * Safian  _et al._ (1988) A. M. Safian, S. R. Coleman, and M. Axenides, Nucl. Phys. B297, 498 (1988).
  * Paccetti Correia and Schmidt (2001) F. Paccetti Correia and M. G. Schmidt, Eur. Phys. J. C21, 181 (2001), arXiv:hep-th/0103189 [hep-th] .
  * Sakai and Sasaki (2008) N. Sakai and M. Sasaki, Prog. Theor. Phys. 119, 929 (2008), arXiv:0712.1450 [hep-ph] .
  * Heeck  _et al._ (2021a) J. Heeck, A. Rajaraman, R. Riley, and C. B. Verhaaren, Phys. Rev. D103, 045008 (2021a), arXiv:2009.08462 [hep-th] .
  * Tsumagari  _et al._ (2008) M. I. Tsumagari, E. J. Copeland, and P. M. Saffin, Phys. Rev. D 78, 065021 (2008), arXiv:0805.3233 [hep-th] .
  * Thorne (1972) K. S. Thorne, in  _Magic Without Magic_ (1972) pp. 231–258.
  * Burrows and Ostriker (2014) A. Burrows and J. P. Ostriker, Proc. Nat. Acad. Sci. 111, 2409 (2014), arXiv:1401.1814 [astro-ph.SR] .
  * Freivogel  _et al._ (2020) B. Freivogel, T. Gasenzer, A. Hebecker, and S. Leonhardt, SciPost Phys. 8, 058 (2020), arXiv:1912.09485 [hep-th] .
  * Amaro-Seoane  _et al._ (2010) P. Amaro-Seoane, J. Barranco, A. Bernal, and L. Rezzolla, JCAP 1011, 002 (2010), arXiv:1009.0019 [astro-ph.CO] .
  * Chen and Piekarewicz (2015) W.-C. Chen and J. Piekarewicz, Phys. Rev. Lett. 115, 161101 (2015), arXiv:1505.07436 [nucl-th] .
  * Chavanis and Harko (2012) P.-H. Chavanis and T. Harko, Phys. Rev. D86, 064011 (2012), arXiv:1108.3986 [astro-ph.SR] .
  * Balakrishna  _et al._ (1998) J. Balakrishna, E. Seidel, and W.-M. Suen, Phys. Rev. D 58, 104004 (1998), arXiv:gr-qc/9712064 .
  * Collodel  _et al._ (2017) L. G. Collodel, B. Kleihaus, and J. Kunz, Phys. Rev. D 96, 084066 (2017), arXiv:1708.02057 [gr-qc] .
  * Sanchis-Gual  _et al._ (2021) N. Sanchis-Gual, C. Herdeiro, and E. Radu, (2021), arXiv:2110.03000 [gr-qc] .
  * Carroll (2004) S. Carroll,  _Spacetime and Geometry: An Introduction to General Relativity_ (Addison-Wesley, San Francisco, 2004).
  * Buchdahl (1959) H. A. Buchdahl, Phys. Rev. 116, 1027 (1959).
  * Urbano and Veermäe (2019) A. Urbano and H. Veermäe, JCAP 1904, 011 (2019), arXiv:1810.07137 [gr-qc] .
  * Guven and O’Murchadha (1999) J. Guven and N. O’Murchadha, Phys. Rev. D60, 084020 (1999), arXiv:gr-qc/9903067 [gr-qc] .