*1.2. Conventions*

Throughout, we set *c* = *h*¯ = 1 and adopt a metric signature +2. Greek letters *α*, *β*, *γ* ... are used for coordinate indices, whereas latin letters *a*, *b*, *c*, ... are used for tetrad basis indices. Symmetrization and antisymmetrization of indices is denoted with round and square brackets, () and [], respectively. We use *∂μ*, <sup>∇</sup>*μ*, and *D*ˆ *μ* to denote partial, covariant, and spinor derivatives, respectively.

The flat spacetime metric in spherical coordinates reads

$$ds^2 = -dt^2 + dr^2 + r^2(d\theta^2 + \sin^2\theta d\phi^2),\tag{1}$$

with *x*0 = *t*, *x*1 = *r*, *x*2 = *θ*, and *x*3 = *ϕ*.

The conventions for scalars are those presented in Ref. [28]. In the Proca field case, we shall use the notation and conventions in [29,30]. For fermions, we shall follow the framework (including the definitions and conventions) in [31], as reviewed in Appendix A.

For the matter content, we shall consider three different fields *ψ*, with the specific form:

• **spin 0:** 

Φis a complex scalar field, which is equivalent to a model with two real scalar fields, <sup>Φ</sup>*R*, <sup>Φ</sup>*I*, via Φ= Φ*<sup>R</sup>* + *i*Φ*I*.

 • **spin 1/2:**

> Ψ(*A*) are massive spinors, with four complex components, the index (*A*) corresponding to the number of copies of the Lagrangian. For a spherically symmetric configuration, one should consider (at least) two spinors (*A* = 1, 2), with equal mass *μ*. A model with a single spinor necessarily possesses a nonzero angular momentum density and it cannot be spherically symmetric.

• **spin 1:**

> A is a complex 4-potential, with the field strength F = *d*A. Again, the model can be described in terms of two real vector fields, A = A*<sup>R</sup>* + *i*A*<sup>I</sup>*.

Some relations ge<sup>t</sup> a simpler form by defining

$$\psi^2 \equiv \{ \Phi^\* \Phi; \ \mathcal{A}\_a \vec{\mathcal{A}}^a; \ i \overleftarrow{\Psi}^{(A)} \Psi^{(A)} \} \,, \tag{2}$$

for spin 0, 1/2 ,and 1, respectively, while we denote *ψ*4 ≡ (*ψ*<sup>2</sup>)<sup>2</sup> and *ψ*6 ≡ (*ψ*<sup>2</sup>)3.

The numerical construction of the solutions reported in this work is standard. In our approach, we use a Runge–Kutta ordinary differential equation solver. For each model, we evaluate the initial conditions that are close to the origin at = 10−<sup>6</sup> for global tolerance <sup>10</sup>−14, adjusting for shooting parameters (which are some constants that enter the near origin expansion of the solutions) and integrating towards *r* → ∞. The accuracy of the solutions was also monitored by computing virial relations that were satisfied by these systems, as discussed below. For a given set of input parameters, the solutions form a discrete set labelled by the number of nodes, *n*, of (some of) the matter function(s). Only fundamental states, which have the minimal number of radial nodes, are considered in this work.

#### **2. The General Framework**

#### *2.1. The Action and Field Equations*

We consider Einstein's gravity that is minimally coupled to spin−*<sup>s</sup>* matter fields *ψ*(*s*) = {<sup>Φ</sup>; Ψ; A}, with *s* = 0, 1/2, 1. The action reads

$$S = \frac{1}{4\pi} \int d^4 \mathbf{x} \sqrt{-\mathbf{g}} \left[ \frac{1}{4\mathbf{G}} \mathbf{R} - L\_{(s)} \right],\tag{3}$$

where

$$L\_{\left(s\right)} = L\_{\left(s\right)}^{\left(0\right)} + \mathcal{U}\_{\left(s\right)}^{\left(\text{int}\right)};\tag{4}$$

the kinetic part of the matter Lagrangians reads, respectively,

$$\begin{aligned} \text{scalar}: \quad &L^{(0)}\_{(0)} = \mathfrak{g}^{a\beta} \Phi^{\*}\_{,a} \Phi\_{,\beta}, \\ \text{Dirac}: \quad &L^{(A)(0)}\_{(1/2)} = i \left[ \frac{1}{2} \left( \{ \mathcal{D} \overline{\Psi}^{(A)} \} \mathcal{V}^{(A)} - \overline{\Psi}^{(A)} \mathcal{D} \Psi^{(A)} \right) \right], \\ \text{Proca}: \quad &L^{(0)}\_{(1)} = \frac{1}{4} \mathcal{F}\_{a\beta} \mathcal{F}^{a\beta}. \end{aligned} \tag{5}$$

*U*(int) (*s*) in (4) is a potential term, which includes self-interactions. In what follows, we shall consider the simplest form of *U*(int) (*s*) , which allows for solitonic solutions in a flat spacetime background, with polynomial terms. As such, apart from the mass term, *U*(int) (*s*) contains quartic and sextic terms only. Afterwards, all three cases can be treated in a unitary way by defining

$$\mathcal{U}\_{\left(\ast\right)}^{\left(\text{int}\right)} = \mathcal{M}\psi^2 - \lambda\psi^4 + \nu\psi^6 \equiv \mathcal{U}\ ,\tag{6}$$

where we have denoted

$$\mathcal{M} \equiv \left\{ \mu^2; \; \mu; \; \frac{1}{2}\mu^2 \right\} \; , \tag{7}$$

for spin 0, 1/2 and 1, respectively. In all cases, *μ* corresponds to the mass of the elementary quanta of the field(s). Additionally, *λ*, *ν* are the input parameters of the theory, which are not fixed *a priori*.

We shall also denote

$$
\delta \mathcal{U} \equiv \frac{\partial \mathcal{U}}{\partial \psi^2} = \mathcal{M} - 2\lambda \psi^2 + 3\nu \psi^4, \qquad \mathcal{U} \equiv \frac{\partial \mathcal{U}}{\partial \psi^2} = -2\lambda + 6\nu \psi^2. \tag{8}
$$

Extremizing the action (3) leads to the Einstein field equations

$$R\_{a\beta} - \frac{1}{2} R\_{\mathcal{S}a\beta} = 2G \ T\_{a\beta} \, , \tag{9}$$

where the energy–momentum tensor reads, respectively,

$$\begin{aligned} \text{scalar}: & \ T\_{\sharp\beta} = \Phi\_{\mathcal{A}}^{\*}\Phi\_{\beta} + \Phi\_{,\beta}^{\*}\Phi\_{,\mathcal{A}} - \mathcal{g}\_{a\beta}L\_{(0)},\\ \text{Dirac}: & T\_{a\beta} = \sum\_{A}T\_{a\beta}^{(A)}, \text{ with } T\_{a\beta}^{(A)} = -\frac{i}{2}\left[\overline{\mathbf{Y}}^{(A)}\gamma\_{(a}\mathcal{D}\_{\beta)}\mathbf{Y}^{(A)} - \left\{\mathcal{D}\_{(a}\overline{\mathbf{Y}}^{(A)}\right\}\gamma\_{\beta}\right]\mathbf{Y}^{(A)} - \mathcal{g}\_{a\beta}L\_{(1/2)}\right],\\ \text{Proca}: & T\_{a\beta} = \frac{1}{2}(\mathcal{F}\_{a\nu}\mathcal{F}\_{\beta\gamma} + \mathcal{F}\_{a\nu}\mathcal{F}\_{\beta\gamma})\mathcal{g}^{\nu\gamma} + \mathcal{U}(\mathcal{A}\_{a}\mathcal{A}\_{\beta} + \mathcal{A}\_{a}\mathcal{A}\_{\beta}) - \mathcal{g}\_{a\beta}L\_{(1)},\end{aligned} \tag{10}$$

and to the matter field equations,

$$\begin{aligned} \text{scalar}: \ (\nabla^2 - \mathcal{U})\Phi &= 0, \\ \text{Dirac}: \ (\mathcal{D} - \mathcal{U})\Psi^{(A)} &= 0, \\ \text{Proca}: \ \frac{1}{2}\nabla\_a \mathcal{F}^{a\beta} - \mathcal{U}\mathcal{A}^{\beta} &= 0. \end{aligned} \tag{11}$$

In the Proca case, taking the four-divergence of the field Equation (11) results in a generalization of the Lorenz condition, which is a dynamical requirement, rather than a gauge choice. This condition takes a particularly simple form in the Ricci-flat case, with

$$(\nabla\_{\mathbf{a}}\mathcal{A}^{\mathbf{a}})\mathcal{U} + \mathcal{U}\nabla\_{\mathbf{a}}(\mathcal{A}\_{\boldsymbol{\beta}}\mathcal{A}^{\boldsymbol{\beta}}) = \mathbf{0} \,. \tag{12}$$

In all cases, the action of the matter fields *ψ* possesses a global *U*(1) invariance, under the transformation *ψ* → *<sup>e</sup><sup>i</sup>αψ*, with *α* constant. This implies the existence of a conserved four-current, which reads

$$\begin{aligned} \text{scalar}: \ j^{\mathbb{A}} &= -i(\Phi^{\*} \partial^{\mathbb{A}} \Phi - \Phi \partial^{\mathbb{A}} \Phi^{\*}), \\ \text{Dirac}: \ j^{\mathbb{A}} &= \Psi \gamma^{\mathbb{A}} \Psi, \\ \text{Proca}: \ j^{\mathbb{A}} &= \frac{i}{2} \left[ \bar{\mathcal{F}}^{\alpha \beta} \mathcal{A}\_{\beta} - \mathcal{F}^{\alpha \beta} \mathcal{A}\_{\beta} \right]. \end{aligned} \tag{13}$$

This current is conserved via the field equations,

$$f\_{;a}^{a} = 0.\tag{14}$$

It follows that integrating the timelike component of this four-current in a spacelike slice Σ yields a conserved quantity—the Noether charge:

$$Q = \int\_{\Sigma} j^t \,. \tag{15}$$

#### *2.2. The Ansätze and Explicit Equations*

#### 2.2.1. The Metric and Matter Fields

The spherically symmetric configurations are most conveniently studied in Schwarzschild-like coordinates, within the following metric Ansatz:

$$ds^2 = -N(r)\sigma^2(r)dt^2 + \frac{dr^2}{N(r)} + r^2(d\theta^2 + \sin^2\theta d\phi^2), \quad \text{with } N(r) \equiv 1 - \frac{2m(r)}{r}.\tag{16}$$

The matter field Ansatz, which is compatible with a spherically symmetric geometry, reads:

$$\text{scalar}: \ \Phi = \phi(r)e^{-i\omega t},\tag{17}$$

introducing a single real function *φ*(*r*);

$$\text{Proca}: \mathcal{A} = \left[ F(r)dt + iH(r)dr \right]e^{-i\text{rot}},\tag{18}$$

introducing two real potentials *<sup>F</sup>*(*r*) and *<sup>H</sup>*(*r*).

The case of a Dirac field is more involved. For a spherically symmetric configurations, we have to consider two Dirac fields, *A* = 1, 2, with

$$\text{Dirac}: \, \Psi^{(A)} = \mathbf{e}^{-\text{iwt}} \mathcal{R}^{(A)}(r) \otimes \Theta^{(A)}(\theta, \phi) \,, \tag{19}$$

where

$$\mathcal{R}^{(1)} = -i\mathcal{R}^{(2)} = \begin{bmatrix} z(r) \\ -i\mathbb{Z}(r) \end{bmatrix}, \quad \mathcal{O}^{(1)} = \begin{bmatrix} -\kappa\sin\frac{\theta}{2} \\ \cos\frac{\theta}{2} \end{bmatrix} \mathbf{e}^{i\frac{\theta}{2}}, \quad \mathcal{O}^{(2)} = \begin{bmatrix} \kappa\cos\frac{\theta}{2} \\ \sin\frac{\theta}{2} \end{bmatrix} \mathbf{e}^{-i\frac{\theta}{2}},$$

with *κ* = ±1 and *z*(*r*) a complex function. In what follows, we shall only consider the case *κ* = 1 (note that qualitatively similar solutions with *κ* = −1 also exist). Additionally, in the above Ansatz, it is convenient to define

$$z(r) = e^{i\pi/4} f(r) - e^{-i\pi/4} g(r) \,, \tag{20}$$

where *f*(*r*) and *g*(*r*) are two real functions, a choice that simplifies the equations. *Symmetry* **2020**, *12*, 2032

The reason that choose two independent *s* = 1/2 fields, with the Ansatz (19) and (20), is the following. For both spinors, the *individual* energy-momentum tensors are not spherically symmetric (this feature is present regardless of the self-interaction potential in the Lagrangian), since *T<sup>t</sup>*(*A*) *ϕ* = 0 (while the other nonzero components *T<sup>r</sup>*(*A*) *r* , *T<sup>θ</sup>*(*A*) *θ* = *Tϕ*(*A*) *ϕ* and *T<sup>t</sup>*(*A*) *t* only depend on *r*). However, *T<sup>t</sup>*(1) *ϕ* + *T<sup>t</sup>*(2) *ϕ* = 0, such that the full configuration is spherically symmetric, being compatible with the line element (16).

Within this framework, the explicit expressions of *ψ*2, as defined in (2), are

$$\psi^2 = \left\{ \phi^2, \ 4(g^2 - f^2), \ H^2N - \frac{F^2}{N\sigma^2} \right\} \tag{21}$$

for spin 0, 1/2 and 1, respectively. In all cases, *w* > 0 is the frequency of the matter field.

#### 2.2.2. The Explicit Equations

The equations for the mass function *m*(*r*) read, respectively,

$$\begin{aligned} \text{scalar}: \; m' &= Gr^2 \left( N\phi'^2 + \frac{w^2 \phi^2}{N\sigma^2} + \mathcal{U} \right), \\ \text{Dirac}: \; m' &= 2Gr^2 \left( 4\sqrt{N}(\mathcal{g}f' - f\mathcal{g}') + \frac{8fg}{r} + \mathcal{U} \right), \\ \text{Proca}: \; m' &= Gr^2 \left[ \frac{(F' - wH)^2}{2\sigma^2} + (\mu^2 - 6\lambda \mathcal{A}^2 + 10\nu \mathcal{A}^4) \frac{F^2}{2N\sigma^2} + \frac{\mathcal{U}}{\mathcal{A}^2} N\mathcal{G}^2 \right]. \end{aligned} \tag{22}$$

The equations for the metric function *σ*(*r*) read, respectively,

$$\begin{aligned} \text{scalar}: \quad \frac{\sigma'}{\sigma} &= 2Gr\left(\phi'^2 + \frac{w^2 \phi^2}{N^2 \sigma^2}\right), \\ \text{Dirac}: \quad \frac{\sigma'}{\sigma} &= 8G\frac{r}{\sqrt{N}} \left(\xi f' - f\xi' + \frac{w(f^2 + g^2)}{N\sigma}\right), \\ \text{Proca}: \quad \frac{\sigma'}{\sigma} &= \frac{2Gr}{N} \left(H^2 N + \frac{F^2}{N\sigma^2}\right) \dot{\mathcal{U}}. \end{aligned} \tag{23}$$

The equations for the matter fields are

$$\begin{aligned} \text{scalar}: \quad & \phi'' + \left(\frac{2}{r} + \frac{N'}{N} + \frac{\sigma'}{\sigma}\right) \phi' + \frac{w^2}{N^2 \sigma^2} \Phi - \mathcal{U} \frac{\Phi}{N} = 0, \\\ \text{Dirac} - \mathbf{f}: \quad & f' + \left(\frac{N'}{4N} + \frac{\sigma'}{2\sigma} + \frac{1}{r\sqrt{N}} + \frac{1}{r}\right) f - \frac{w\varrho}{N\sigma} + \frac{\mathcal{S}}{\sqrt{N}} \mathcal{U} = 0, \\\ \text{Dirac} - \mathbf{g}: \ g' + \left(\frac{N'}{4N} + \frac{\sigma'}{2\sigma} - \frac{1}{r\sqrt{N}} + \frac{1}{r}\right) f + \frac{wf}{N\sigma} + \frac{f}{\sqrt{N}} \mathcal{U} = 0, \\\ \text{Proca} - \mathbf{F}: \qquad & F' - wH + \frac{2N\sigma^2 H}{w} \mathcal{U} = 0, \\\ \text{Proca} - \mathbf{H}: \quad \frac{d}{dr} \left\{ \frac{r^2 \left[wH - F'\right]}{\sigma} \right\} + \frac{2r^2 F}{N\sigma} \mathcal{U} = 0 \ . \end{aligned} \tag{24}$$

An additional supplementary constraint is also present, which is a second order equation for the metric functions *m*(*r*) and *<sup>σ</sup>*(*r*), with first order derivatives of matter fields. However, this equation is a differential consequence of the above field equations; it is used to check the accuracy of the numerical results.

Let us remark that the above equations can be derived from the following effective action

$$\mathcal{S}\_{\rm eff} = \int d\mathbf{r} \mathcal{L}\_{\rm eff} \, \prime \qquad \text{with} \quad \mathcal{L}\_{\rm eff} = \frac{1}{\overline{\mathcal{G}}} \sigma m' - \mathcal{L}\_{\rm (s)} \, \prime \tag{25}$$

where

$$\begin{split} \mathcal{L}\_{(0)} &= r^2 \sigma \left( N \mathfrak{g}'^2 - \frac{w^2 \mathfrak{g}^2}{N \sigma^2} + \mathcal{U} \right), \\ \mathcal{L}\_{(1/2)} &= 8r^2 \sigma \left( \sqrt{N} (\mathfrak{g}f' - f \mathfrak{g}') - \frac{w}{\sigma \sqrt{N}} (f^2 + \mathfrak{g}^2) + \frac{2f\mathfrak{g}}{r} + \frac{1}{4} \mathcal{U} \right), \\ \mathcal{L}\_{(1)} &= r^2 \sigma \left( -\frac{(F' - wH)^2}{2\sigma^2} + \mathcal{U} \right). \end{split} \tag{26}$$

Given the above framework, the energy density measured by a static observer, *ρ* = − *Tt t* , is

$$\begin{aligned} \text{Scalar}: \,\rho &= N\phi'^2 + \frac{w^2}{N\sigma^2}\phi^2 + \mathcal{U}, \\ \text{Dirac}: \,\rho &= 8\left( (\mathcal{g}f' - f\mathcal{g}')\sqrt{N} + \frac{2f\mathcal{g}}{r} + \frac{1}{4}\mathcal{U} \right). \\ \text{Proca}: \,\rho &= \frac{(F'-wH)^2}{2\sigma^2} + \frac{\mathcal{U}}{\mathcal{A}^2}NH^2 + \frac{F^2}{N\sigma^2}\left(\frac{1}{2}\mu^2 - 3\lambda\mathcal{A}^2 + 5\nu\mathcal{A}^4\right). \end{aligned} \tag{27}$$

The mass of the flat space solutions is computed as the integral of *ρ*,

$$M = \int dr r^2 \rho\_{\prime} \tag{28}$$

while, in the self-gravitating case, it can be read off from the asymptotic behaviour of the *gtt* metric potential

$$\log\_{tt} = -N\sigma^2 = -1 + \frac{2MG}{r} + \dots \tag{29}$$

As for the Noether charge, it reads

$$\begin{aligned} \text{scalar}: \ Q &= 2w \int\_0^\infty dr \, r^2 \frac{\Phi^2}{Nr'}\\ \text{Dirac}: \ Q &= 4 \int\_0^\infty dr \, r^2 \frac{(f^2 + g^2)}{\sqrt{N}}, \\ \text{Proca}: \ Q &= 2 \int\_0^\infty dr \, r^2 \frac{(wH - F')N}{\sigma}. \end{aligned} \tag{30}$$

#### *2.3. Units and Scaling Symmetries*

The matter fields in action (3) have the following dimensions: (with *L* = *lenght*):

$$\begin{aligned} \text{scalar}: \ [\Phi] &= \frac{1}{L'}\\ \text{Dirac}: \ [\Psi] &= \frac{1}{L^{3/2}},\\ \text{Proca}: \ [\mathcal{A}\_{\theta}] &= \frac{1}{L}. \end{aligned} \tag{31}$$

Additionally, the coupling constants that enter the potential are (generically) dimension-full, with:

$$\begin{aligned} \text{scalar}: \ [\mu] &= \frac{1}{L}, \ [\lambda] = L^0, \ [\nu] = L^2, \\ \text{Dirac}: \ [\mu] &= \frac{1}{L}, \ [\lambda] = L^2, \ [\nu] = L^5, \\ \text{Proca}: \ [\mu] &= \frac{1}{L}, \ [\lambda] = L^0, \ [\nu] = L^2. \end{aligned} \tag{32}$$

Turning now to scaling symmetries, we notice the existence of *three* different transformations, which, in all cases, leave invariant the equations of motion (in the relations below, the functions or constants that are not mentioned explicitly remain invariant).

The first one (s0) is very simple

$$(s0): \quad \sigma \to c\sigma, \ w \to cw,\tag{33}$$

with *c* some arbitrary positive constant. However, this invariance is fixed when imposing the asymptotic flatness condition *σ* → 1 as *r* → ∞.

More importantly, the equations of motion are invariant under a scaling of the radial coordinate, together with other functions and parameters of the model. In the scalar and Proca cases, this transformation reads:

$$(s1):\quad \text{scalar and Proca}:\ r = c\overline{r},\ \ w = \frac{1}{c}\overline{\upsilon},\ \ \mu = \frac{1}{c}\overline{\mu},\ \ \lambda = \frac{\lambda}{c^2},\ \ \upsilon = \frac{\upsilon}{c^4},\ \ m = c\overline{m}.\tag{34}$$

In the Dirac case, the corresponding symmetry is more complicated, with

$$(s1): \quad \text{Dirac}: \; r = c\mathfrak{r}, \; w = \frac{1}{c}\mathfrak{d}, \; \mu = \frac{1}{c}\mathfrak{d}, \; \nu = c\mathfrak{v}; \; m = c\mathfrak{n}, \; f = \frac{f}{\sqrt{c}}, \; g = \frac{g}{\sqrt{c}}.\tag{35}$$

In all cases, the product *mμ* and the ratio *w*/*μ* are left invariant by the symmetry (*s*<sup>1</sup>). Under this transformation, the global quantities behave as

$$
\mathcal{M} = c\bar{\mathcal{M}},\ \mathcal{Q} = c^2 \bar{\mathcal{Q}}.\tag{36}
$$

The symmetry (*s*<sup>1</sup>) is usually employed in order to work in units of length set by the field mass,

$$
\mu = 1, \text{ i.e. } c = \frac{1}{\mu}. \tag{37}
$$

Finally, the equations are also invariant under a suitable scaling of the matter field(s) functions, together with some coupling constants (while the radial coordinate or the field frequency are not affected):

$$\mathbf{r}(s\mathbf{2}):\ \ \lambda = \frac{\lambda}{c^{\mathbf{2}}},\ \ \nu = \frac{\nu}{c^{\mathbf{2}}},\ \ G = \frac{\mathbf{G}}{c^{\mathbf{2}}},\tag{38}$$

together with

$$\text{s(s2)}: \quad \text{scalar}: \text{ } \phi = \text{\"c\"} \text{c}, \quad \text{Dirac}: f = f \text{\"c\"}, \text{ } g = \text{\"g\text{c}}, \quad \text{Proca}: F = F \text{\"c\"}, \text{ } H = H \text{\"c\"}, \tag{39}$$

while the global quantities transform as

$$M = c^2 \ddot{M}, \ Q = c^2 Q. \tag{40}$$

The symmetries (*s*<sup>1</sup>) and (*s*<sup>2</sup>) are used, in practice, in order to simplify the numerical study of the solitons. First, the symmetry (*s*<sup>2</sup>) can be used to set *G*¯ = 1; i.e., essentially one absorbs Newton's constant in the expression of the matter field(s). This is the usual approach for the models without self-interaction, see, e.g., the discussion in [12].

However, the probe limit of the solutions becomes less transparent for this choice. An alternative route, employed in this work, is to use (38) and (39) to set the coefficient of the quartic term to unity,

$$
\bar{\lambda} = 1, \text{ i.e. } \mathfrak{c} = \frac{1}{\sqrt{\lambda}}.\tag{41}
$$

It follows that two mass scales naturally emerge, one set by gravity *M*Pl = 1/√*G* and the other one, *M*0, set by the field(s) coupling constants, with

$$\text{scalar and Proca}:\ M\_0 = \frac{\mu}{\sqrt{\lambda}};\quad \text{Dirac}:\ M\_0 = \frac{1}{\sqrt{\lambda}}.\tag{42}$$

The ratio of these fundamental mass scales defines the dimensionless coupling constant

$$
\kappa = \frac{M\_0}{M\_{\rm Pl}} \,\,\,\tag{43}
$$

which is relevant in the physics of the solutions.

Another dimensionless input parameter is the scaled constant for the sextic term in the potential *U*, with

$$\text{scalar and Proca}: \ \beta = \frac{\nu \mu^2}{\lambda^2}; \quad \text{Dirac}: \ \beta = \frac{\nu \mu}{\lambda^2}. \tag{44}$$

As such, the mass and charge of the non-gravitating solutions is given in units set by *μ*, *λ*, while, in the gravitating case, in order to make contact with the previous results without self-interaction, we use units that are set by *G* and *μ*.

*A priori*, the range of *α* is unbounded, 0 *α* < ∞. The limit *α* → 0 corresponds to *G* → 0, i.e., the probe limit—one solves the matter field(s) equations on a fixed geometry, which should be a solution of the vacuum Einstein equations. The limit *α* → ∞ corresponds to *λ* → 0, i.e., no quartic, or sextic term (when using *β*, if *ν* is finite) in the action. Thus, solutions of the Einstein-matter field equations without self-interaction are approached in the second limit.

This choice of units with *μ* = *λ* = 1 (after employing the above scaling symmetries) has the advantage of greatly simplifying the numerical study of the solutions. For example, the Einstein equations read

$$R\_{a\beta} - \frac{1}{2}g\_{a\beta} = 2a^2 T\_{a\beta\prime} \tag{45}$$

the only input parameters being

$$\{\mathfrak{a},\ \beta \text{ and } w\},\tag{46}$$

with *w* the scaled frequency. Additionally, the scaled scalar potential reads

$$
\mathcal{U}I = \psi^2 - \psi^4 + \beta \psi^6,\tag{47}
$$

For *ψ*2 > 0, *U* is strictly positive for *β* > 1/4. The case *β* = 1/4 is special, since the potential becomes *U* = *ψ*<sup>2</sup>(1 − *ψ*2/2)2, and, thus, possesses three degenerate minima, at *ψ* = {0, ±√2}. A discussion of these aspects (for a spin-zero field) can be found in Ref. [32].

#### **3. The Probe Limit: Flat Spacetime Solutions**

#### *3.1. Deser's Argument and Virial-Type Identities*

Before performing a numerical study of the solutions, it is useful to derive virial-type identities. For a flat spacetime background, we can adapt a simple argumen<sup>t</sup> given long ago by Deser [26] (used therein to rule out the existence of finite energy time-independent solutions in Yang–Mills theory) in order to obtain virial-type identities and a simple relation between the mass of solutions and the trace of the energy-momentum tensor (the same virial identities are found by adapting Derrick's scaling argumen<sup>t</sup> [33].)—see also [34].

Working in Cartesian coordinates *x<sup>a</sup>* (*a* = 1, 2, 3), assume the existence of a stationary soliton in some field theory model. Following [26], consider the following (trivial) identity observe that Deser's argumen<sup>t</sup> cannot be extended to a curved spacetime background)

$$\frac{\partial}{\partial \mathbf{x}^a} \left( \mathbf{x}^b T^a\_b \right) = T^a\_a + \mathbf{x}^b \frac{\partial T^a\_b}{\partial \mathbf{x}^a} \,, \tag{48}$$

together with its volume integral. The left hand side vanishes from regularity and finite energy requirements (note that, in the spin-1/2 case, one considers the total energy-momentum tensor). The second term in (48) vanishes from energy-momentum conservation (plus staticity) and, thus, we are left with the virial-type identity

$$\int d^3 \mathbf{x} \, T\_a^a = 0. \tag{49}$$

It follows that the total mass-energy of a static soliton in *d* = 3 + 1 dimensions is determined by the trace of the energy-momentum tensor

$$M = -\int d^3 \mathbf{x} \,\, T\_t^t = -\int d^3 \mathbf{x} \,\, T\_\mu^\mu \,. \tag{50}$$

When applied to the specific Ansätze in this work, (49) results in the following expressions:

• scalar field:

$$\int\_0^\infty dr \, r^2 \left[ \frac{1}{3} \phi'^2 + (\mu^2 - w^2) \phi^2 + \nu \phi^6 \right] = \lambda \int\_0^\infty dr \, r^2 \phi^4. \tag{51}$$

This relation can be simplified by using the equation for *φ* in a simpler form

$$\int\_0^\infty dr \, r^2 (\mu^2 - w^2) \phi^2 = \frac{\lambda}{2} \int\_0^\infty dr \, r^2 \phi^4. \tag{52}$$

This clearly shows that the (flat space) Q-ball solutions are supported by the quartic self-interacting term, with *λ* > 0 (i.e., the sextic term is not relevant at this level).

• Dirac field:

$$\int\_0^\infty dr \, r^2 \left( \mathcal{g}f' - f\mathcal{g}' + \frac{2f\mathcal{g}}{r} + \frac{3}{8}lI - \frac{3}{2}w(f^2 + \mathcal{g}^2) \right) = 0 \tag{53}$$

which can also be written, via field equations, in the alternative form

$$\int\_0^\infty dr \, r^2 \left(\mu \psi^2 + \lambda \psi^4\right) = \int\_0^\infty dr \, r^2 \left(4w(f^2 + g^2) + 3\nu \psi^6\right). \tag{54}$$

One concludes that, for *ν* = 0 and *λ* > 0, the solutions are supported by the harmonic time dependence. However, the above relation does not clarify the role that is played by the nonlinear quartic term for the existence of solitons.

• Proca field:

$$\int\_0^\infty dr \, r^2 \left(\frac{1}{2} (F' - wG)^2 - 3lI + 2G^2 \frac{dII}{d\mathcal{A}^2}\right) = 0 \,\,. \tag{55}$$

After eliminating the kinetic term (*F* − *wG*)2, via field equations, the above relation takes the suggestive form:

$$\int\_0^\infty dr \, r^2 \left[ \mu^2 (3F^2 + (\frac{\mu^2}{w^2} - 1)G^2 + \frac{16\lambda^2 G^2}{w^2} \mathcal{A}^4 + \frac{6\nu}{w^2} \mathcal{A}^4 \left( w^2 (F^2 + G^2) + 2G^2 (\mu^2 + 3\nu \mathcal{A}^2) \right) \right] = 0$$

$$2\lambda \int\_0^\infty dr \, r^2 \mathcal{A}^2 \left( 3F^2 + G^2 + \frac{4G^2}{w^2} (\mu^2 + 6\nu \mathcal{A}^2) \right) . \tag{56}$$

From the bound state condition *μ*2 *w*2, it is clear that the existence of finite mass solutions requires a quartic term, *λ* = 0 (*ν* being irelevant). However, since A<sup>2</sup> may take both positive or negative values, one cannot use the above relation in order to predict the sign of *λ*.

It is also interesting to note that the mass of the flat space Proca solitons takes the simple form

$$M = \int\_0^\infty dr \, r^2 \left(\mu^2 \mathcal{A}^2 - 2\nu \mathcal{A}^6\right). \tag{57}$$

#### *3.2. Numerical Results*

#### 3.2.1. General Remarks

The corresponding equations are found by taking

$$N = \sigma = 1\,,\tag{58}$$

in the corresponding general equations presented in Section 2, and we shall not display them here. Moreover, the boundary conditions that are satisfied by the functions *ψ* at *r* = 0, ∞ are similar to those in the gravitating case, as given in the next Section.

The case of a scalar field is special for a flat spacetime metric. The frequency parameter *w* is not relevant, since *w*<sup>2</sup> acts as an effective *tachyonic* contribution to the mass term and, thus, it can be absorbed into *μ*2, by defining *μ*2 − *w*<sup>2</sup> → *μ*2. After this redefinition, the scalar field is static, Φ = *φ*(*r*) and, thus, the Noether charge vanishes. Therefore, all *Q*-ball solutions in a flat spacetime background can be interpreted as *static* scalar solitons, in a model with a shifted scalar field mass term [35], for a new potential

$$
\mathcal{U} = \mathcal{U}\_{(w)} - w^2 \phi^2. \tag{59}
$$

Note that, although *φ* formally satisfies the same equation as before, the energy-momentum tensor and the total mass of these solutions are different. Additionally, the virial identities imply that the redefined potential *U* is necessarily negative for some range of *φ*, which is realised by the solutions, *U* < 0. Moreover, one can prove the following relation [35]

$$M(w=0) = M(w) - wQ\_\prime \tag{60}$$

which relates the mass of a static solution (*w* = 0 and the redefined potential (59)) to the mass and Noether charge of solutions with a given *w* (the other parameters in the potential are kept fixed).

No similar relation seems to exist for higher spin fields. However, it is interesting to notice the existence in this case of a curious static, purely electric solution, i.e., with *H* = 0, *w* = 0. The electric potential *<sup>F</sup>*(*r*) essentially satisfies the same equation as a scalar *Q*-ball,

$$(F'' + \frac{2F'}{r} - (\mu^2 + 4\lambda F^2 + 6\nu F^4)F = 0\,\text{.}\tag{61}$$

The existence of solutions with proper asymptotics requires *λ* < 0, being constructed within the same scheme as in the generic case.

To the best of our knowledge, this case has not been discussed in the literature. However, they possess some unphysical features. For example, by using the virial identity (57) together with the above equations, one can prove that the total mass of this configuration is negative (moreover, we have verified that *M* remains negative when including gravity effects)

$$M = -\frac{1}{3} \int\_0^\infty dr r^2 F'^2 < 0. \tag{62}$$

In all cases, the sextic self-interaction term is not necessary for the existence of solutions. Because the *β* = 0 case has some special properties, we shall discuss it separately at the end of this section.

#### 3.2.2. Solutions with a Sextic Self-Interaction Term, *β* > 0

In what follows, in order to simplify the picture, we shall assume *β* > 1/4, such that the potential *U*(*ψ*<sup>2</sup>) (as given by (6)) is strictly positive (assuming *ψ*2 > 0, which, as we shall see, is not necessarily the case).

Because no exact solutions are known, the solitons are constructed numerically. Figures 1–3 (left panels) show the profile of typical solutions. One notices that, in the Proca case, *ψ*2 = A<sup>2</sup> is negative for small enough *r*, while *ψ*2 = *i*Ψ(*A*)Ψ(*A*) is positive in the Dirac case.

**Figure 1.** (**Left**): the radial profile of a typical non-gravitating scalar soliton. (**Right**): the mass and Noether charge are shown *vs*. the scalar field frequency for the fundamental family of non-gravitating scalar solitons.

The numerical results indicate that the *s* = 0, 1/2, 1 flat spacetime solitons follow the same pattern, which can be summarized, as follows. First, in all cases, the solutions only exist in a certain frequency range, *w*min < *w* < *w*max = *μ*, with *w*min determined by *β*. Second, the solutions with *w* < *μ* decay exponentially in the far field, with no radiation, as *ψ* ∼ *<sup>e</sup>*−√*μ*2−*w*2*r*. Finally, the most interesting qualitative feature is perhaps the existence in all cases of a mass gap. That is, at a critical value of the frequency, both the mass and charge of the solutions assume their minimal (nonzero) value, from where they monotonically increase towards both limiting values of the frequency (see Figures 1–3 (right panels)). When considering the mass of the solutions as a function of the charge, there are thus

two branches, merging at the minimal charge/mass. One expects that a subset of solutions with *M* smaller than the mass of *Q* free particles (bosons or fermions) may be stable, which is possible along the lower frequency branch.

**Figure 2.** Same as Figure 1 for Dirac stars. Note that the single particle condition, *Q* = 1, is not imposed here.

**Figure 3.** Same as Figure 1 for Proca stars.

3.2.3. Solutions without a Sextic Self-Interaction Term, *β* = 0

This case is likely physically less relevant (since the potential *U* is not positive definite); however, it possesses some interesting properties, which depend, to some extent, on the spin of the field.

Starting with the scalar case, some numerical results are displayed in Figure 1 (right panel, inset). Note that all of these solutions are unstable, since *M* > *Q*, and, thus, they have excess energy. We remark that solutions with a negative energy density region, *ρ* < 0, found on a small *<sup>r</sup>*−interval, are found for small enough *w* (although *M* > 0 always).

In fact, the equation for the scalar field has an interesting form, with the frequency parameter being irrelevant. After using the alternative rescaling

$$r \to r / \sqrt{\mu^2 - w^2}, \quad \phi \to \phi \sqrt{\frac{\mu^2 - w^2}{\lambda}},\tag{63}$$

one can see that no free parameter exists in this case. The scalar field satisfies the kink-like equation

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d\phi}{dr}\right) = \phi - 2\phi^3,\tag{64}$$

which does not seem to possess an exact solution (Equation (64) is discussed by many authors, being relevant for the issue of false vacuum decay [36]—e.g., an existence proof can be found in Ref. [37]). Additionally, one can show that the following relation holds

$$M(w; \mu, \lambda) = \frac{M(w = 0; \mu, \lambda)}{\sqrt{\mu^2 - w^2}},\tag{65}$$

with *<sup>M</sup>*(*w* = 0; 1, 1) 1.503.

The pattern that is exhibited by the Proca solutions without a sextic term is different, as displayed in Figure 3. As found in Ref. [24], the limit *w* = 0 is not approached in this case. Instead, a minimal value of *w* is reached, with a backbending towards a critical solution possessing finite charges. Because *M* > *Q*, one expects all of these solutions to be unstable. Additionally, the mass is still positive, *M* > 0, and we did not notice the existence of negative energy densities (although we could not find an analytical argumen<sup>t</sup> to show that *ρ* > 0).

Finally, the pattern of Dirac solitons without a sextic self-interaction (which was the original case in the pioneering work by Soler [20] ) seems to be similar to the one that is found in the *β* > 0 case, see Figure 2 (right panel). In particular, both *M* and *Q* still diverge at the limits of the *w*-interval.

#### **4. Including the Gravity Effects**

#### *4.1. The Boundary Conditions*

The boundary conditions satisfied by the metric functions at the origin are

$$
\sigma(0) = 0, \ \sigma(0) = \sigma\_{0\prime} \tag{66}
$$

while at infinity, one imposes

$$m(\infty) = M, \ \sigma(\infty) = 1,\tag{67}$$

with *<sup>σ</sup>*(0), *M* numbers fixed by numerics. The matter functions should vanish as *r* → ∞

$$\phi(\infty) = f(\infty) = g(\infty) = F(\infty) = H(\infty) = 0,\tag{68}$$

while the boundary conditions at the origin are

$$\begin{aligned} \text{scalar}: \left. \frac{d\phi(r)}{dr} \right|\_{r=0} &= 0, \\ \text{Dirac}: \left. f(0) = 0, \left. \frac{d\lg(r)}{dr} \right|\_{r=0} = 0, \\ \text{Proca}: \left. \frac{dF(r)}{dr} \right|\_{r=0} &= 0, \left. H(0) = 0. \right. \end{aligned} \tag{69}$$

Let us mention that, in each case, one can construct an approximate form of the solutions both at *r* = 0 and at infinity, which is compatible with the boundary conditions above.

#### *4.2. Virial Identities*

The reduced action (25) allows for us to prove that the solutions satisfy the (simple enough) interesting virial identities, which generalize for a curved geometry Deser's relations (51), (53) and (55).

Following [38,39] (wherein generalizations for a curved geometry of Derrick's argumen<sup>t</sup> for flat spacetime [33] were established), assume the existence of a solution that is described by *<sup>m</sup>*(*r*), *σ*(*r*) and *ψ*(*r*) (with *ψ* = {*φ*; *f* , *g*; *F*, *H*} the matter fields) with suitable boundary conditions at *r* = 0 and

at infinity. Subsequently, each member of the one-parameter family *<sup>m</sup>*Λ(*r*) ≡ *<sup>m</sup>*(<sup>Λ</sup>*r*), *<sup>σ</sup>*Λ(*r*) ≡ *<sup>σ</sup>*(<sup>Λ</sup>*r*), and *ψ*Λ(*r*) ≡ *ψ*(<sup>Λ</sup>*r*), assumes the same boundary values at *r* = 0 and at *r* = <sup>∞</sup>, and the action *S*Λ ≡ *Seff* [*<sup>m</sup>*Λ, *σ*Λ, *ψ*Λ] must have a critical point at Λ = 1, i.e., [*dS*/*d*Λ]<sup>Λ</sup>=<sup>1</sup> = 0. Thus, the putative solution must satisfy the following virial relations (here, it is useful to work with unscaled variables, i.e., before considering the transformations (s2), (s3)):

• scalar field:

$$\int\_0^\infty dr \, r^2 \sigma \left[ \phi'^2 + \frac{w^2 (1 - 4N) \phi^2}{N^2 \sigma^2} + 3lI \right] = 0 \; ; \tag{70}$$

• Dirac field:

$$\int\_0^\infty dr \, r^2 \sigma \left[ \sqrt{H} (\mathcal{g} f' - f \mathcal{g}') \left( 2 + \frac{m}{rN} \right) + \frac{4fg}{r} - \frac{w(f^2 + g^2)}{\sigma \sqrt{N}} \left( 3 - \frac{m}{rN} \right) - 3U \right] = 0 \; ; \tag{71}$$

• Proca field:

$$\begin{split} &\int\_{0}^{\infty} dr \, r^{2} \sigma \left[ \mu^{2} \left( \mathcal{A}^{2} \left( 4 - \frac{1}{N} \right) + 2H^{2} (1 - N) \right) + 2\lambda \mathcal{A}^{2} \left( \mathcal{A}^{2} \left( 5 - \frac{2}{N} \right) + 4H^{2} (1 - N) \right) \right. \\ &\left. + 6\nu \mathcal{A}^{4} \left( \mathcal{A}^{2} \left( 2 - \frac{1}{N} \right) + 2H^{2} (1 - N) + 4G^{2} (1 - N) \right) - \frac{(wH - F') (3wH - F')}{\sigma^{2}} \right] = 0. \end{split} \tag{72}$$

This relation can be further simplified after using the Proca field equations to yield

$$\int\_0^\infty dr \, r^2 \sigma \left[ \mu^2 \left( 4 - \frac{1}{N} \right) \mathcal{A}^2 + 2\lambda \left( 5 - \frac{2}{N} \right) \mathcal{A}^4 + 6\nu \left( 2 - \frac{1}{N} \right) \mathcal{A}^6 \right] = \tag{73}$$
 
$$\left[ 2 \int\_0^\infty dr \, r^2 \sigma H^2 (\mu^2 + 4\lambda \mathcal{A}^2 + 6\nu \mathcal{A}^4) \left[ 2N - 1 + \frac{N^2 \sigma^2}{2w^2} (\mu^2 + 4\lambda \mathcal{A}^2 + 6\nu \mathcal{A}^4) \right] \right]$$

These expressions are not transparent, with the effects of gravity and self-interaction being mixed. As such, their main use is to test the accuracy of the numerical results. However, the situation changes once we set *λ* = *ν* = 0 (i.e., no self-interaction). Subsequently, one can see that the solutions are supported by an harmonic time dependence of the matter fields, *w* = 0, except for the Dirac case, where we could not prove a similar result.

#### *4.3. General Features*

The flat space solitons can be generalized to curved spacetime. The presence of higher order self-interacting terms in the potential is not crucial for the existence of gravitating solutions. However, they affect some of their quantitative features.

A common pattern emerges again, with the basic generic properties of the gravitating solutions being summarized, as follows. First, the solutions are topologically trivial, with 0 *r* < ∞. They possess no horizon; the size of the *S*2-sector of the metric shrinks to zero as *r* → 0. At infinity, a Minkowski background is approached. Second, perhaps the most interesting feature is that gravity regularizes the divergencies of the mass and charge that are found in the probe case at the limit(s) of the *w*-interval. In particular, this regularization implies that no mass gap exists for gravitating configurations: *M*, *Q* → 0 as *w* → *μ*. Additionally, assuming the existence of a sextic self-intertion term, with *β* > 1/4, in all cases, the family of solutions describes a continuous curve in a mass *M* (or charge *Q*), *vs*. frequency, *w*, diagram. This curve starts from *M* = *Q* = 0 for *w* = *μ*, in which limit the fields become very diluted and the solution trivializes. At some intermediate frequency, a maximal mass is attained (that may be a global, or only local, maximum, depending on *α*, *β*).

The effects of the quartic and sextic self-interacting terms in the potential (6) become irrelevant for large enough *α*. For example, the *α* = 10 curves that are displayed in Figures 4–6 are well approximated by the corresponding Einstein-scalar/Dirac/Proca results, with *β* = *ν* = 0, the maximal relative difference being of only a few percent (towards the critical value of frequency).

**Figure 4.** ADM mass and Noether charge of the gravitating scalar boson stars *vs*. the scalar field frequency for families of solutions with three different values of the coupling constant *α*. The solutions in the right panel do not possess a sextic self-interacting term.

**Figure 5.** The same as Figure 4 for Dirac stars. The single particle condition, *Q* = 1, is not imposed here.

**Figure 6.** Same as Figure 4 for Proca stars.

One can also identify some specific features, as follows. The scalar solutions with *β* > 1/4 (i.e., a positive potential, *U* > 0) and the Proca starts with large *α* describe a spiral in a (*<sup>M</sup>*, *w*)-diagram. As in similar cases, likely, these spirals approach, at their centre, a critical singular solution. The Dirac solutions with *β* = 0 also describe a spiral. For general Proca solutions with *β* = 0 and Dirac solutions with *β* = 0, the (*<sup>M</sup>*, *w*)-curves appear to end in a critical configuration before describing a spiral. In the Proca case, this feature is discussed in [24].

The scalar field solutions with *β* = 0 possess a more complicated pattern. For small *α*, they possess a static limit. Moreover, two disconnected branches of solutions exist for some intermediate range of *α* e.g., *α* = 0.5 in Figure 4 (right panel). In addition to the familiar spiral starting at *w*/*μ* = 1 and ending for some critical nonzero *w* = *wc*, one finds a secondary set, which extends from *w* = 0 to some maximal value of *w* < *wc*.

Finally, let us mention the existence of another interesting possibility, with quartic interaction only and *λ* < 0, in which case no flat spacetime solitons are found. The corresponding solutions were discussed in [40] for a scalar field, and in [25] for a Proca field. No similar study exists so far for a fermionic field. In the bosonic case, perhaps the most interesting feature is that their maximal mass is proportional with |*λ*|, and, thus, can increase dramatically.

#### **5. Other Aspects**

#### *5.1. No Hair Results*

It is interesting to inquire whether the solutions above may allow the existence of a black hole horizon, inside. Indeed, this is a well known feature that is found for a variety of other solitons, see e.g., the review work [27]. However, this is not the case for the (relatively) simpler solutions in this work.

The situation can be summarized, as follows:

• **Scalar case**

> Peña and Sudarsky established a no-scalar-hair theorem, ruling out a class of spherically symmetric BHs with scalar hair [41]. Their proof also covers the case of the potential (6) considered in this work, and it still holds if the hair has the harmonic time-dependence that we consider.

• **Dirac case**

> No hair results for the Einstein–Dirac case with a massive spinor (no self-interaction ) were proposed in [42,43]. The nonexistence of stationary states for the nonlinear Dirac equation with a quartic self-interaction on the Schwarzschild metric has been proven in [44].

• **Proca case**

> A no hair theorem has been proven in [30] for a massive, non-selfinterating Proca field. In Appendix B, we generalize it for an arbitrary Proca potential *<sup>U</sup>*(A<sup>2</sup>).

#### *5.2. The Issue of Particle Numbers: Bosons vs. Fermions*

In all of the results displayed above, we have equally treated the bosonic and fermionic fields. However, while the classical treatment of the Dirac equation is mathematically justified (a discussion on classical spinors and their possible physical justification can be found, e.g., in Ref. [45]), physically, its fermionic nature should be imposed at the level of the occupation number: at most, a single, particle, in accordance to Pauli's exclusion principle.

Similarly to the non-self-interacting case that is discussed in [12], the one particle condition is imposed, as follows. Let us suppose that we have a numerical solution with some values for the mass and charge (*M*(num), *Q*(num)). Subsequently, we can use the symmetry (34)–(36) in order to map it to a 'new' solution with

$$Q\_{\text{(num)}} = Q\_{0\prime} \tag{74}$$

where *Q*0 > 0 is arbitrary. This assumption fixes the value of the scaling parameter *c*,

$$c = \sqrt{\frac{Q\_{\text{(num)}}}{Q\_0}},\tag{75}$$

such that the *numerical* mass of the 'new' solution will be [12]

$$
\dot{M}\_{\text{(num)}} = M\_{\text{(num)}} \sqrt{\frac{Q\_0}{Q\_{\text{(num)}}}} \,. \tag{76}
$$

Note that the values of *w*, *μ*, *λ*, and *ν* should also be scaled according to (34) and (35).

This transformation is used to normalize the total charge of a fermion to one, *Q*0 = 1. With this condition *Q* = 1, the (*<sup>M</sup>*, *w*)-curves in Figure 2 (right panels) and Figure 5 are not sequences of solutions with fixed parameters in the potential and varying *Q*; instead, they ge<sup>t</sup> mapped into sequences with fixed *Q* and varying *μ*, *ν*. Consequently, one is discussing a sequence of solution of different models (since *μ*, *ν* are input parameters in the action).

Figure 7 shows the corresponding results, for both the probe limit case and including gravity effects. An interesting feature here is that the mass of *Q* = 1 non-gravitating configuration can still take arbitrary large values. However, as expected, gravity effects lead to a picture that is qualitatively similar to that found in the *λ* = *ν* = 0 case [12]. Again, both the total mass, *M* and the mass of the field *μ* are bounded and never exceed, roughly, *MPl*.

**Figure 7.** (**Left panel**) Soliton mass *vs*. the mass of the elementary quanta of the field, for non-gravitating solutions of the Dirac equations with three different values of *β*. (**Right panel**) The same for the gravitating solutions with several values of the coupling constant *α* and *β* = 0. The single particle condition, *Q* = 1, is imposed here.

#### **6. Further Remarks. Conclusions**

The main purpose of this work was to provide a comparative analysis of three different types of solitonic solutions of GR coupled with matter fields of spin *s* = 0, 1/2, 1. A unified framework has been proposed, analysing these cases side-by-side under a consistent set of notations and conventions. Differently from the previous work [12], the matter fields herein are self-interacting, such that all three models possess ( *Q*-ball-like) solutions on the flat spacetime limit. As such, a more complicated landscape of gravitating solutions is found, with some new qualitative features when compared to the picture that is revealed in [12].

However, despite this fact, our study shows that there is, again, a certain universality in the properties of the solitons, being to some extent independent of their (fundamental) spin. As with linear matter fields, the basic ingredients are again: (**i**) a harmonic time dependence of the matter fields; (**ii**) complex field(s)/multiplets such that the energy–momentum is still real; and, (**iii**) the existence of of mass term as a trapping mechanism, creating bound states with *w* < *μ*. Additionally, if one requires the presence of a well defined flat space limit, then (**iv**) the fields should possess at least a quartic self-interaction term.

As an avenue for future work, it would be interesting to go beyond the case of spherical symmetry and consider a comparative study of axially symmetric, spinning solutions. In the non-self-interacting case, this was the subject of the recent work [46], where a common pattern was again revealed to exist. The situation in the presence of self-interactions is less studied; only the scalar field case has, so far, been discussed in the literature [47,48]. Indeed, even flat space spinning solitons with spin *s* = 1/2, 1 are ye<sup>t</sup> unreported in the literature. Moreover, even in the static case, new families of non-spherically symmetric solitons should exist, generalizing, for a self-interacting potential (and possibly for a higher spin), the *s* = 0 multipolar boson stars recently reported in [49].

Finally, let us remark that, for a bosonic field, *s* = 0, 1, the no-hair theorems (as reviewed in the previous Section) can be circumvented for an horizon that rotates synchronously with the field, leading to BHs with scalar or Proca hair [30,50]. However, this does not seem to be the case for fermions, regardless of the presence of a self-interacting potential.

**Author Contributions:** The two authors contributed equally to the conceptualization, methodology, software, validation, formal analysis, investigation, resources, data curation, writing–original draft preparation, writing–review and editing, visualization, supervision, project administration, funding acquisition. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported by the Center for Research and Development in Mathematics and Applications (CIDMA) through the Portuguese Foundation for Science and Technology (FCT - Fundacao para a Ciência e a Tecnologia), references UIDB/04106/2020 and UIDP/04106/2020 and by national funds (OE), through FCT, I.P., in the scope of the framework contract foreseen in the numbers 4, 5 and 6 of the article 23, of the Decree-Law 57/2016, of August 29, changed by Law 57/2017, of July 19. We acknowledge support from the projects PTDC/FIS-OUT/28407/2017, CERN/FIS-PAR/0027/2019 and PTDC/FIS-AST/3041/2020. This work has further been supported by the European Union's Horizon 2020 research and innovation (RISE) programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740. The authors would like to acknowledge networking support by the COST Action CA16104. Computations were performed at the Blafis cluster, in Aveiro University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. The Dirac Field: Conventions**

Since this case is more complicated, we shall include here the basic relations.

Following the framework in [31], we consider a general four dimensional metric *gαβ*, and introduce a tetrad of vectors

$$\mathcal{e}\_a^a = \{e\_0^a, e\_1^a, e\_2^a, e\_3^a\},\tag{A1}$$

which we take to be an orthonormal basis, i.e.,

$$
\epsilon\_{a\beta}\epsilon^a\_a\epsilon^\beta\_b = \eta\_{a\beta\prime} \tag{A2}
$$

where

$$
\eta\_{ab} = \text{diag}(-1, 1, 1, 1). \tag{A3}
$$

Here, Roman and Greek letters are used for tetrad and coordinate indices, respectively. Roman indices are raised and lowered with *ηab*. It follows that

$$
\epsilon^a\_a = \eta^{ab} \varrho\_{a\beta} \epsilon^\beta\_b \quad \text{and} \quad \lg\_{a\beta} = \eta\_{ab} \epsilon^a\_a \epsilon^b\_\beta. \tag{A4}
$$

The next step is to define two sets of 4 × 4 matrices *γ<sup>α</sup>* and *γ*ˆ *a* satisfying the anticommutation relations

$$\{\gamma^a, \gamma^\beta\} = 2g^{a\beta}I\_{4\prime} \quad \{\gamma^a, \gamma^b\} = 2\eta^{ab}I\_{4\prime} \tag{A5}$$

(where, as usual, {*<sup>A</sup>*, *B*} = *AB* + *BA*). Note that the former set *γ<sup>α</sup>* are functions of spacetime position, whereas the latter *γ*ˆ *a* have constant components. The two sets may be related with any orthonormal tetrad,

$$
\gamma^a = \mathfrak{e}\_a^a \hat{\gamma}^a. \tag{A6}
$$

Matrix indices are raised/lowered in the standard way: *γ*ˆ *a* = *ηabγ*<sup>ˆ</sup> *b* and *γα* = *<sup>g</sup>αβγβ*. OneusestheWeyl/chiralrepresentation,inwhich

 

$$
\tilde{\gamma}^0 = \begin{pmatrix} O & I \\ I & O \end{pmatrix}, \quad \tilde{\gamma}^i = \begin{pmatrix} O & \sigma\_i \\ -\sigma\_i & O \end{pmatrix}, \qquad i = 1, 2, 3,\tag{A7}
$$

where *I* is the 2 × 2 identity, *O* is the 2 × 2 zero matrix, and *σi* are the Pauli matrices

$$
\sigma\_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma\_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \sigma\_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. \tag{A8}
$$

Then the matrices *γ*ˆ *a* are defined as

$$
\gamma^1 = \mathbf{i}\,\gamma^3,\ \gamma^2 = \mathbf{i}\,\tilde{\gamma}^1,\ \hat{\gamma}^3 = \mathbf{i}\,\tilde{\gamma}^2,\ \hat{\gamma}^0 = \mathbf{i}\,\tilde{\gamma}^0.\tag{A9}
$$

The Dirac four-spinor is written as

$$\Psi = \begin{pmatrix} \psi\_-\\ \psi\_+ \end{pmatrix},\tag{A10}$$

where *ψ*+ and *ψ*− are (left- and right-handed) two-spinors, which may be projected out from Ψ with the operators *P*± = 12(*I* ± *γ*ˆ5) where

$$
\hat{\gamma}^5 = i\hat{\gamma}^0 \hat{\gamma}^1 \hat{\gamma}^2 \hat{\gamma}^3 = \begin{pmatrix} -I & O \\ O & I \end{pmatrix}.\tag{A11}
$$

We define also the Dirac conjugate

$$
\overline{\Psi} \equiv \Psi^\dagger a,\tag{A12}
$$

with Ψ† denoting the usual Hermitian conjugate and *α* = −*γ*<sup>ˆ</sup> 0.

The spinor covariant derivative *D* ˆ *ν* which enters the Dirac Equation (11) is

$$
\hat{D}\_V = \partial\_V - \Gamma\_V. \tag{A13}
$$

The spinor connection matrices Γ*ν* are defined, up to an additive multiple of the unit matrix, by the relation

$$
\partial\_\nu \gamma^\mu + \Gamma^\mu\_{\,\nu\lambda} \gamma^\lambda - \Gamma\_\nu \gamma^\mu + \gamma^\mu \Gamma\_\nu = 0 \,,\tag{A14}
$$

where Γ*μνλ* is the affine connection. A suitable choice satisfying (A14) makes use of the spin connection *ωα bc*,

$$
\Gamma\_a = -\frac{1}{4} \omega\_{a\,bc} \hat{\gamma}^b \hat{\gamma}^c,\tag{A15}
$$

the spin-connection *ωμab* being defined as

$$
\omega\_{\mu}{}^{a}{}\_{b} = \epsilon^{a}\_{\nu} \varepsilon^{\lambda}\_{b} \Gamma^{\nu}{}\_{\mu \lambda} - \varepsilon^{\lambda}\_{b} \partial\_{\mu} \varepsilon^{a}\_{\lambda}. \tag{A16}
$$

Also, the covariant derivative of the conjugate spinor is

$$
\hat{D}\_{\mu}\overline{\Psi} = \partial\_{\mu}\overline{\Psi} + \overline{\Psi}\Gamma\_{\mu}.\tag{A17}
$$

#### **Appendix B. Self-Interacting Proca Field: A No Hair Result**

As for the previous study [30] for a non-self-interacting field, the theorem is established by contradiction. Let us assume the existence of a regular BH solution of the Einstein-Proca equations. The general Ansatz and the field equations derived in Section 2 apply also to this case. Differently from the globally regular case, the geometry would possess a non-extremal horizon at, say, *r* = *rH* > 0, which requires that

*<sup>N</sup>*(*rH*) = 0 , (A18)

while *σ*(*rH*) > 0. Since we are assuming that there are no more exterior horizons, then *r* > *rH* =constant are timelike surfaces and *<sup>N</sup>* (*rH*) > 0. Also, we can choose without loss of generality that *σ*(*rH*) > 0, since the equations of motion are invariant under *σ* → <sup>−</sup>*σ*. It follows that *<sup>N</sup>*(*r*) and *σ*(*r*) are strictly positive functions for any *r* > *rH*.

In establishing the theorem, we shall use the relation (12) written in the generic form

$$\frac{d}{dr}\left(r^{2}N\sigma H\dot{\mathcal{U}}\right) = -\frac{wr^{2}F\dot{\mathcal{U}}}{N\sigma},\tag{A19}$$

together with one of the Proca equations

$$F' = wH\left(1 - \frac{2N\sigma^2\bar{\mathcal{U}}}{w^2}\right). \tag{A20}$$

The regularity of the horizon implies that the energy density of the Proca field is finite there and also the norm of the Proca-potential A. One can easily see that this condition implies

$$F(r\_H) = 0 \,. \tag{A21}$$

Then the function *<sup>F</sup>*(*r*) starts from zero at the horizon and remains strictly positive (or negative) for some *r*-interval (the case of a negative *<sup>F</sup>* (*r*) can be discussed in a similar way). Now, let us assume *<sup>F</sup>* (*r*) > 0 for *rH* < *r* < *r*1. It follows that, in this interval, *<sup>F</sup>*(*r*) is a strictly increasing (and positive) function.

Next, we consider the expression (which appears in Equations (A19) and (A20)),

$$P(r) \equiv 1 - \frac{2\sigma^2(r)N(r)l\dot{l}}{w^2} \,. \tag{A22}$$

One can see that *<sup>P</sup>*(*rH*) = 1; actually *P* becomes negative for large *r*, since *N* → 1, *σ* → 1 as *r* → <sup>∞</sup>, while A<sup>2</sup> → 0 and *μ* > *w* (which is a bound state condition necessary for an exponential decay of the Proca field at infinity). But the important point is the existence of an *<sup>r</sup>*−interval *rH* < *r* < *r*2 where *<sup>P</sup>*(*r*) is a strictly positive function. Now, the same reasoning applies also for *U* ˙ (since *U* ˙ (*rH*) = *μ*2/2 > 0), while *U* ˙ → 0 asymptotically. This implies the existence of some interval in the vicinity of the horizon where *U* ˙ , *F* and *P* are all positive.

At this point, let us consider an arbitrary value of *r* in this interval. Then we observe that (A19) implies

$$r^2 \sigma(r) N(r) H(r) \dot{\mathcal{U}} = -w \int\_{r\_H}^r dx \frac{\mathbf{x}^2}{\sigma(\mathbf{x}) N(\mathbf{x})} F(\mathbf{x}) \dot{\mathcal{U}}(\mathbf{x}) < 0 \tag{A23}$$

in that interval. Consequently, *<sup>H</sup>*(*r*) < 0, since *σ*, *N* are positive everywhere outside the horizon.

The last conclusion implies a contradiction: *<sup>H</sup>*(*r*) < 0 is not compatible with *F* (*r*) > 0, in that interval. In fact, *F* (*r*) > 0 together with *P* > 0 and *w* > 0, from (A20), that *<sup>H</sup>*(*r*) > 0. Thus we conclude that *<sup>F</sup>*(*r*) = *<sup>H</sup>*(*r*) = 0 is the only solution compatible with a BH geometry (*q*.*e*.*d*.).
