Asymptotically flat, spherical, self-interacting scalar, Dirac and Proca stars

We present a comparative analysis of the self-gravitating solitons arising in the Einstein-Klein-Gordon, Einstein-Dirac and Einstein-Proca models, for the particular case of static, spherically symmetric spacetimes. Differently from the previous study arXiv:1708.05674, the matter fields possess suitable self-interacting terms in the Lagrangians, which allow for the existence of $Q$-ball--type solutions for these models in the flat spacetime limit. In spite of this important difference, our analysis shows that the high degree of universality observed in arXiv:1708.05674 remains, and various spin-independent common patterns are observed.

1 Introduction and motivation 1

.1 General remarks
The (modern) idea of solitons as extended particle-like configurations can be traced back (at least) to Lord Kelvin, around one and half centuries ago, who proposed that atoms are made of vortex knots [2]. The first explicit example of solitons in a relativistic field theory, however, was found by Skyrme [3], (almost) one hundred years later. The latter, dubbed Skyrmions, exist in a model with four (real) scalars subject to a constraint. They were proposed as a field theory realizations of baryons; nonetheless Skyrmions capture the basic properties of a generic soliton.
In the context of this work, solitons are defined as spatially localized, singularity-free solutions of a field theory model, which possess a finite mass and angular momentum 1 . As already suspected in the literature before Skyrme's seminal work 2 many non-linear field theories possess such solutions. The study of solitons has attracted a lot of attention in the last 50 years; they have important applications in various contexts ranging from e.g. the models of condensed matter physics to high energy physics and cosmology. Moreover, solitons are also significant for a nonperturbative quantum description of states. Famous examples of solitons include the sphalerons in the Standard Model [5] and the magnetic monopoles [6] in various extensions thereof.
A review of soliton physics in the non-gravitating case can be found in the textbooks [7,8]; see also the reviews [9,10]. In this work, we shall focus on a special class of solitons in four spacetime dimensions that are time-dependent (at the level of the matter field) but non-dissipative (with a time independent spacetime geometry), and carrying a Noether charge Q associated with a global U (1) symmetry. We shall restrict ourselves to models with a single matter field (except for the fermionic case) and possessing a standard kinetic term. These solitons do not possess topological properties; but the field theory Lagrangian must have suitable self-interactions.
Three different models are considered here: (ungauged) field theories describing particles of spin s = 0, 1/2, 1, with the goal of providing a comparative analysis of their spherically symmetric solitonic solutions. As such, this work is an extension of that in Ref. [1], which dealt with the same matter fields, but without selfinteractions. The latter introduce a novel feature: self-interactions allow the existence of solitonic solutions even on a fixed flat spacetime background, i.e. ignoring gravity. As we shall see, self-interactions lead to a more intricate landscape of (gravitating) solutions, with some new qualitative features as compared to the picture revealed in [1].
Historically, in the non-gravitating case, solitons in a complex scalar field model have been known for more than 50 years, being dubbed Q − balls by Coleman [11] (see also the previous work [12]). These are central for our discussion, sharing most of the basic features with the higher spin generalizations. Within a Dirac field model with a positive quartic self-interaction term, Ivanenko [13], Weyl [14], Heisenberg [15] and Finkelstein et. al. [16,17], have made early attempts to construct particle-like solutions. The first rigorous numerical study of such configurations, however, has been done by Soler in 1976 [18]. The case of a (Abelian) vector field has been ignored until recently; the study of Q-ball-like solutions was initiated by Loginov [19].
In all three cases, a large literature has grown based on these early studies, which includes generalizations in various directions. A particularly interesting case concerns the study of backreaction of these solutions on the spacetime geometry. This is a legitimate question, since the existence of solitons relies on the nonlinearity of the field theory and General Relativity (GR) is intrinsically highly nonlinear. Therefore, gravitational interactions could significantly alter the flat spacetime solitons.
The results in the literature prove that all flat space solitons survive when including gravity effects; however, a more complicated picture emerges as compared to the non-self-interacting case. The study of gravitating generalization of the flat spacetime Q-balls has started with Friedberg et al. [20]. The higher spin fields are less studied; GR generalization of the self-interacting spinors were considered only recently by Dzhunushaliev and Folomeev [21], while the gravitating Q-Proca stars have been studied by Brihaye et al. [22] and Minamitsuji [23].
The main purpose of this work is to review various existing results on this type of solitons, putting them together under a consistent set of notations and conventions. A comparative study is presented, starting with the generic framework and scaling symmetries (Section 2). There and afterwards, the mathematical description of each of the three models is made in parallel to emphasise their similarities. The basic properties of solutions in a flat spacetime background are discussed in Section 3, where we adapt Deser's argument [24] to find virial identities which yield some insight on the existence of such configurations. Gravity is included in Section 4; a curved spacetime generalization of the Deser-type virial identities is also discussed. Extending the models' framework to include GR effects allows for the possibility to add a black hole (BH) horizon at the center of soliton, a possibility which is indeed realized for other well-known solutions, such as Skyrmions, magnetic monopoles and sphalerons (see e.g. the review work [26]). However, as we shall discuss in Section 4 (and Appendix B), the situation is different for the considered (spherically symmetric) solutions, and there are no BH generalizations, regardelss of the presence of self-interactions. We also discuss the physical interpretation of the fermionic solutions. Concluding remarks and some open questions are presented in Section 5.

Conventions
Throughout, we set c = = 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 ∂ µ , ∇ µ andD µ to denote partial, covariant and spinor derivatives, respectively.
The flat spacetime metric in spherical coordinates reads The conventions for scalars are those in Ref. [27]. In the Proca field case, we shall use the notation and conventions in [28,29]. For fermions we shall follow the framework (including the definitions and conventions) in [30], 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, Φ R , Φ I , via Φ = Φ R + 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 cannot be spherically symmetric.
• spin 1: A is a complex 4-potential, with the field strength F = dA. Again, the model can be described in terms of two real vector fields, A = A R + iA I .
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 close to the origin at = 10 −6 for global tolerance 10 −14 , adjusting for shooting parameters (which are some constants which enter the near origin expansion of the solutions) and integrating towards r → ∞. The accuracy of the solutions was also monitored by computing virial relations 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

The action and field equations
We consider Einstein's gravity minimally coupled to spin−s matter fields ψ (s) = {Φ; Ψ; A}, with s = 0, 1/2, 1. The action reads where the kinetic part of the matter Lagrangians reads, respectively, where we have denoted for spin 0, 1/2 and 1, respectively. In all cases, µ corresponds to the mass of the elementary quanta of the field(s). Also, λ, ν are input parameters of the theory, which are not fixed a priori. We shall also denotė Extremizing the action (2.3) leads to the Einstein field equations where the energy-momentum tensor reads, respectively, Proca : In the Proca case, taking the 4-divergence of the field eqs. (2.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 In all cases, the action of the matter fields ψ possesses a global U (1) invariance, under the transformation ψ → e iα ψ, with α constant. This implies the existence of a conserved 4-current, which reads This current is conserved via the field equations, 14) It follows that integrating the timelike component of this 4-current in a spacelike slice Σ yields a conserved quantity -the Noether charge:

The metric and matter fields
The spherically symmetric configurations are most conveniently studied in Schwarzschild-like coordinates, within the following metric Ansatz: The matter field Ansatz which is compatible with a spherically symmetric geometry reads: introducing a single real function φ(r); introducing two real potentials F (r) and H(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 where with κ = ±1 and z(r) a complex function. In what follows, we shall consider the case κ = 1 only (note that qualitatively similar solutions with κ = −1 exist as well). Also, in the above Ansatz it is convenient to define where f (r) and g(r) are two real functions, a choice which simplifies the equations. The reason we choose two independent s = 1/2 fields, with the Ansatz (2.19), (2.20), is the following. For both spinors, the individual energy-momentum tensors are not spherically symmetric 3 such that the full configuration is spherically symmetric, being compatible with the line element (2.16).
Within this framework, the explicit expressions of ψ 2 , as defined in (1.2), are for spin 0, 1/2 and 1, respectively. In all case, w > 0 is the frequency of the matter field.

The explicit equations
The equations for the mass function m(r) read, respectively, The equations for the metric function σ(r) read, respectively, scalar : The equations for the matter fields are An additional supplementary constraint is also present, which is a 2nd order equation for the metric functions m(r) and σ(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 where Given the above framework, the energy density measured by a static observer, ρ = −T t t , is The mass of the flat space solutions is computed as the integral of ρ, while in the self-gravitating case it can be read off from the asymptotic behaviour of the g tt metric potential As for the Noether charge, it reads

Units and scaling symmetries
The matter fields in action (2.3) have the following dimensions: (with L = lenght): Also, the coupling constants which enter the potential are (generically) dimension-full, with: 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 which are not mentioned explicitly remain invariant).
The first one (s0) is very simple In the Dirac case, the corresponding symmetry is more complicated, with In all cases the product mµ and the ratio w/µ are left invariant by the symmetry (s1). Under this transformation, the global quantities behave as The symmetry (s1) is usually employed to work in units of length set by the field mass, 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): The symmetries (s1), (s2) are used in practice to simplify the numerical study of the solitons. First, the symmetry (s2) can be used to setḠ = 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 [1].
However, the probe limit of the solutions becomes less transparent for this choice. An alternative route, employed in this work, is to use (2.38), (2.39) the set the coefficient of the quartic term to unity, 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 scalar and Proca : The ratio of these fundamental mass scales defines the dimensionless coupling constant 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 scalar and Proca : β = νµ 2 λ 2 ; Dirac : β = νµ λ 2 . (2.44) As such, the mass and charge of the non-gravitating solutions is given in units set by µ, λ, while in the gravitating case, to make contact with the previous results without self-interaction, we use units 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 to greatly simplify the numerical study of the solutions. For example, the Einstein equations read with w the scaled frequency. Also, the scaled scalar potential reads For ψ 2 > 0, U is strictly positive for β > 1/4. The case β = 1/4 is special, since the potential becomes U = ψ 2 (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. [31].
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 argument given long ago by Deser [24] (used therein to rule out the existence of finite energy time-independent solutions in Yang-Mills theory) to obtain virial-type identities and a simple relation between the mass of solutions and the trace of the energy-momentum tensor 4 -see also [25].
Working in Cartesian coordinates x a (a = 1, 2, 3), assume the existence of a stationary soliton in some field theory model. Following [24], consider the following (trivial) identity 5 together with its volume integral. The left hand side vanishes from regularity and finite energy requirements 6 . The second term in (3.48) vanishes from energy-momentum conservation (plus staticity) and thus we are left with the virial-type identity 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 When applied to the specific Ansätze in this work, (3.49) results in the following expressions: • scalar field: This relation can be simplified by using the equation for φ in a simpler form 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: which can also be written, via field equations, in the alternative form 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 played by the nonlinear quartic term for the existence of solitons. 4 The same virial identities are found by adapting Derrick's scaling argument [37]. 5 Observe that Deser's argument cannot be extended to a curved spacetime background. 6 Note that in the spin-1/2 case, one considers the total energy-momentum tensor.
• Proca field: After eliminating the kinetic term (F − wG) 2 , via field equations, the above relation takes the suggestive form: 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). But since A 2 may take both positive or negative values, one cannot use the above relation to predict the sign of λ.
It is also interesting to note that the mass of the flat space Proca solitons takes the simple form (3.57)

General remarks
The corresponding equations are found by taking in the corresponding general equations in Section 2, and we shall not display them here. Moreover, the boundary conditions 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 2 acts as an effective tachyonic contribution to the mass term, and thus it can be absorbed into µ 2 , by defining µ 2 − w 2 → µ 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 [32], for a new potential (3.59) Note that although φ formally satisfies the same equation as before, the energy-momentum tensor and the total mass of these solutions are different. Also, 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 [32] M (w = 0) = M (w) − wQ, (3.60) which relates the mass of a static solution (w = 0 and the redefined potential (3.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 F (r) satisfies essentially the same equation as a scalar Q-ball, 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 (3.57) together with the above equations one can prove that the total mass of this configuration is negative 7 In all cases, the sextic self-interaction term is not necessary for the existence of solutions. Since the β = 0 case has some special properties, we shall discuss it separately at the end of this section.

Solutions with a sextic self-interaction term, β > 0
In what follows, to simplify the picture, we shall assume β > 1/4, such that the potential U (ψ 2 ) (as given by (2.6)) is strictly positive (assuming ψ 2 > 0, which, as we shall see, is not necessarily the case).
Since no exact solutions are known, the solitons are constructed numerically. The profile of typical solutions are shown in Figures 1-3 (left panels). One notices that in the Proca case, ψ 2 = A 2 is negative for small enough r, while ψ 2 = iΨ (A) Ψ (A) is positive in the Dirac case. 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 exist only 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 ψ ∼ e − √ µ 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 the 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)). 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.  Fig. 1 for Dirac stars. Note that the single particle condition, Q = 1, is not imposed here.

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 Fig. 1 (right panel, inset). Note that all 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 r−interval, are found for small enough w (although M > 0 always).
In fact, the equation for the scalar field has an interesting form, the frequency parameter being irrelevant. After using the alternative rescaling one can see that no free parameter exists in this case. The scalar field satisfies the kink-like equation The pattern exhibited by the Proca solutions without a sextic term is different, being displayed in Fig.  3. As found in Ref. [22], 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. Since M > Q, one expects all these solutions to be unstable. Also, 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 argument 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 [18] ) seems to be similar to the one found in the β > 0 case, see Figure 2  Let us mention that in each case, one can construct an approximate form of the solutions both at r = 0 and at infinity, compatible with the boundary conditions above.

Virial identities
The reduced action (2.25) allows us to prove that the solutions satisfy the (simple enough) interesting virial identities, which generalize for a curved geometry Deser's relations (3.51), (3.53, (3.55). Following [35,36] (wherein generalizations for a curved geometry of Derrick's argument for flat spacetime [37] were established), assume the existence of a solution described by m(r), σ(r) and ψ(r) (with ψ = {φ; f, g; F, H} the matter fields) with suitable boundary conditions at r = 0 and at infinity. Then each member of the 1-parameter family m Λ (r) ≡ m(Λr), σ Λ (r) ≡ σ(Λr), and ψ Λ (r) ≡ ψ(Λr), assumes the same boundary values at r = 0 and at r = ∞, and the action S Λ ≡ S ef f [m Λ , σ Λ , ψ Λ ] must have a critical point at Λ = 1, i.e. [dS/dΛ] Λ=1 = 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: • Dirac field: • Proca field: This relation can be further simplified after using the Proca field equations to yield 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). Then 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.

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, 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 found in the probe case at the limit(s) of the w-interval. In particular, this regularization implies no mass gap exists for gravitating configurations: M, Q → 0 as w → µ. Also, 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 becomes very diluted and the solution trivializes. At some intermediate frequency, a maximal mass is attained (which may be a global, or only local, maximum, depending on α, β).
The effects of the quartic and sextic self-interacting terms in the potential (2.6) become irrelevant for large enough α. For example, the α = 10 curves displayed in Figures 4-6 are well approximated by the   Fig. 4 for Dirac stars. The single particle condition, Q = 1, is not imposed here.
corresponding Einstein-scalar/Dirac/Proca results with β = ν = 0, the maximal relative difference being of only a few percent (towards the critical value of frequency).
One can also identify also 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 (M, 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 (M, w)-curves appear to end in a critical configuration before describing a spiral. In the Proca case, this feature is discussed in [22].
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 = w c , one finds a secondary set, extending from w = 0 to some maximal value of w < w c .
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 [38] for a scalar field, and in [23] 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

No hair results
It is interesting to inquire if the solutions above may allow the existence of a black hole horizon, inside. Indeed, this is a well known feature found for a variety of other solitons, see e.g. the review work [26]. 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 [39]. Their proof covers also the case of the potential (2.6) considered in this work, and still holds if the hair has the harmonic time-dependence we consider.
• Dirac case No hair results for the Einstein-Dirac case with a massive spinor (no self-interaction ) were proposed in [40,41]. The nonexistence of stationary states for the nonlinear Dirac equation with a quartic self-interaction on the Schwarzschild metric has been proven in [42].
• Proca case A no hair theorem has been proven in [29] for a massive, non-selfinterating Proca field. In appendix B, we generalize it for an arbitrary Proca potential U (A 2 ).

The issue of particle numbers: bosons vs. fermions
In all results displayed above, we have treated equally the bosonic and fermionic fields. However, while the classical treatment of the Dirac equation is mathematically justified 9 , 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 discussed in [1], the one particle condition is imposed as follows. Let us suppose we have a numerical solution with some values for the mass and charge (M (num) , Q (num) ). Then we can use the symmetry (2.34)- (2.36) in order to map it to a 'new' solution with where Q 0 > 0 is arbitrary. This assumption fixes the value of the scaling parameter c, such that the numerical mass of the 'new' solution will be [1] .
Note that the values of w, µ, λ and ν should also be scaled according to (2.34), (2.35).
This transformation is used to normalize the total charge of a fermion to one, Q 0 = 1. With this condition 10 Q = 1 the (M, w)-curves in Fig. 2 (right panels) and Fig. 5 are not sequences of solutions with fixed parameters in the potential and varying Q; instead, they get 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).  The corresponding results are shown in Fig. 7, 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 qualitatively similar to that found in the λ = ν = 0 case [1]. Again, both the total mass, M and the mass of the field µ are bounded and never exceeds, roughly, M P l .

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 [1], the matter fields herein are self-interacting, such that all three models possess (Q-balllike) solutions on the flat spacetime limit. As such, a more complicated landscape of gravitating solutions is found, with some new qualitative features as compared to the picture revealed in [1].
Despite this fact, however, 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 < µ. Also, 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 [43], 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 [46,47]. Indeed, even flat space spinning solitons with spin s = 1/2, 1 are yet 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 [48].
Finally, let us remark that 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 [44], [29]. However, this does not seem to be the case for fermions, regardless of the presence of a self-interacting potential.
Here, Roman and Greek letters are used for tetrad and coordinate indices, respectively. Roman indices are raised and lowered with η ab . It follows that Matrix indices are raised/lowered in the standard way:γ a = η abγ b and γ α = g αβ γ β . One uses the Weyl/chiral representation, in which where I is the 2 × 2 identity, O is the 2 × 2 zero matrix, and σ i are the Pauli matrices Then the matricesγ a are defined aŝ The Dirac four-spinor is written as where ψ + and ψ − are (left-and right-handed) two-spinors, which may be projected out from Ψ with the operators P ± = 1 2 (I ±γ 5 ) whereγ The spinor connection matrices Γ ν are defined, up to an additive multiple of the unit matrix, by the relation where Γ µ νλ is the affine connection. A suitable choice satisfying (A.14) makes use of the spin connection ω α bc , B Self-interacting Proca field: a no hair result As for the previous study [29] 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 = r H > 0, which requires that N (r H ) = 0 , (B.1) while σ(r H ) > 0. Since we are assuming that there are no more exterior horizons, then r > r H =constant are timelike surfaces and N (r H ) > 0. Also, we can choose without loss of generality that σ(r H ) > 0, since the equations of motion are invariant under σ → −σ. It follows that N (r) and σ(r) are strictly positive functions for any r > r H . In establishing the theorem, we shall use the relation (2.12) written in the generic form together with one of the Proca equations 3) 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 Then the function F (r) starts from zero at the horizon and remains strictly positive (or negative) for some r-interval (the case of a negative F (r) can be discussed in a similar way). Now, let us assume F (r) > 0 for r H < r < r 1 . It follows that, in this interval, F (r) is a strictly increasing (and positive) function. Next, we consider the expression (which appears in Eqs. (B.2), (B.3)), P (r) ≡ 1 − 2σ 2 (r)N (r)U w 2 . (B.5) One can see that P (r H ) = 1; actually P becomes negative for large r, since N → 1, σ → 1 as r → ∞, while A 2 → 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 r−interval r H < r < r 2 where P (r) is a strictly positive function. Now, the same reasoning applies also forU (sinceU (r H ) = µ 2 /2 > 0), whileU → 0 asymptotically. This implies the existence of some interval in the vicinity of the horizon whereU , in that interval. Consequently, H(r) < 0, since σ, N are positive everywhere outside the horizon. The last conclusion implies a contradiction: H(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 (B.3), that H(r) > 0. Thus we conclude that F (r) = H(r) = 0 is the only solution compatible with a BH geometry (q.e.d.).