Dynamical formation and stability of fermion-boson stars
Fabrizio Di Giovanni, Saeed Fakhry, Nicolas Sanchis-Gual, Juan Carlos Degollado, José A. Font
DDynamical formation and stability of fermion-boson stars
Fabrizio Di Giovanni, Saeed Fakhry,
2, 1
Nicolas Sanchis-Gual, Juan Carlos Degollado, and Jos´e A. Font
1, 5 Departamento de Astronom´ıa y Astrof´ısica, Universitat de Val`encia,Dr. Moliner 50, 46100, Burjassot (Val`encia), Spain Department of Physics, Shahid Beheshti University, G. C., Evin, Tehran 19839, Iran Centro de Astrof´ısica e Gravita¸c˜ao - CENTRA,Departamento de F´ısica, Instituto Superior T´ecnico - IST,Universidade de Lisboa - UL, Avenida Rovisco Pais 1, 1049-001, Portugal Instituto de Ciencias F´ısicas, Universidad Nacional Aut´onoma de M´exico,Apdo. Postal 48-3, 62251, Cuernavaca, Morelos, M´exico Observatori Astron`omic, Universitat de Val`encia,C/ Catedr´atico Jos´e Beltr´an 2, 46980, Paterna (Val`encia), Spain
Gravitationally bound structures composed by fermions and scalar particles known as fermion-boson stars are regular and static configurations obtained by solving the coupled Einstein-Klein-Gordon-Euler (EKGE) system. In this work, we discuss one possible scenario through which thesefermion-boson stars may form by solving numerically the EKGE system under the simplifying as-sumption of spherical symmetry. Our initial configurations assume an already existing neutron starsurrounded by an accreting cloud of a massive and complex scalar field. The results of our simula-tions show that once part of the initial scalar field is expelled via gravitational cooling the systemgradually oscillates around an equilibrium configuration that is asymptotically consistent with astatic solution of the system. The formation of fermion-boson stars for large positive values of thecoupling constant in the self-interaction term of the scalar-field potential reveal the presence of anode in the scalar field. This suggests that a fermionic core may help stabilize configurations withnodes in the bosonic sector, as happens for purely boson stars in which the ground state and thefirst excited state coexist.
PACS numbers: 95.30.Sf, 04.70.Bw, 04.40.Nr, 04.25.dg
I. INTRODUCTION
Identifying the relevance scalar fields may have forastrophysics and cosmology, in particular as potentialcomponents of the dark matter content of the universe,has long received considerable attention [1–4]. Differentscalar fields have been considered, namely the dilaton instring theories [5, 6], the Higgs boson in the standardmodel of particle physics [7, 8], the inflaton in studiesof the early universe [9, 10], or the axion as a possiblecomponent of cold dark matter [11–14].It has been argued that ultralight bosons form local-ized and coherently oscillating configurations very similarto Bose-Einstein condensates [15, 16]. When the mass ofthe bosonic particle is around 10 − eV [17, 18] thesecondensates provide an alternative to the standard ap-proach to explain large-scale structure formation throughdark-matter seeds. For heavier bosons, the bound con-figurations are smaller and may have the typical size andmass of a sellar compact object such as a neutron star.These objects are generically known as boson stars [19].Boson stars are gravitationally bound configurationsof scalar particles. Since the seminal works of Kaup [20]and Ruffini and Bonnazola [21] their description has beengeneralised in several ways including self-interaction [22],charge [23], rotation [24, 25], oscillating soliton stars [26],stars with more than a single scalar field [27, 28], andeven vector fields (in which case the bosonic star is knownas a Proca star [29]). Reviews on the subject can befound in references [30, 31]. If such bosonic configurations could form from someprimordial gas, it is natural to assume that other par-ticles, such as fermions, could also be present duringthe condensation. Therefore, it would seem theoreti-cally possible that objects made out of a mixture of bothbosons and fermions might also form. Even if the originalconfigurations were mainly composed by either bosonsor fermions, they could be susceptible to further cap-ture fermions and bosons through accretion giving riseto mixed configurations. It is thus a theoretically in-teresting question to investigate the properties of thesemacroscopic composites of fermions and bosons, referredin the literature as fermion-boson stars [32–36] and todiscuss possible means by which they might form. Thisis the focus of this paper. Here we propose a dynam-ical scenario in which a fermionic star (modelled as apolytropic star for simplicity) accretes part of the scalarfield, while part of it is radiated to infinity, and a mixedfermion-boson star forms.The gravitational condensation of a primordial gasand the subsequent radiation of part of the bosonic fieldhas been dubbed gravitational cooling and has been ad-dressed in [37] for purely scalar fields and in [38] for vec-tor fields. Using numerical-relativity simulations thosestudies have shown the dynamical formation of bosonstars and Proca stars, respectively, under the assump-tion of spherical symmetry. In order to be astrophysicallyrelevant, a gravitationally bound system that forms dy-namically must be stable for times much longer than itscharacteristic dynamical timescale. The stability prop- a r X i v : . [ g r- q c ] J un erties of boson stars have been considered in [39–46].In Ref. [47] Seidel and Suen discussed the dynamicalevolution of perturbed boson stars finding, in particu-lar, that unstable stars migrate to the stability region ofstatic configurations which suggests the formation of bo-son stars under generic initial conditions. Further studieson the formation of boson stars were performed in [37]in general relativity and in [48, 49] in the Newtonianregime. These studies concluded that self-gravitating,scalar-field stellar systems settle down into equilibriumconfigurations. We note that this conclusion does notonly apply to the scalar case but it is also valid for thevector counterparts of boson stars, i.e. Proca stars, ashas recently been reported in [38].The purpose of this work is twofold: on the one handwe aim to describe the dynamical formation of fermion-boson stars; on the other hand, we will analyse the stabil-ity properties of those configurations considering a strongself-interaction term in the Klein-Gordon potential of thebosonic part. For this study, and for the sake of sim-plicity, we shall focus on fermion-boson stars assumingspherical symmetry. The starting point of our analy-sis assumes a preexisting neutron star (described with apolytropic equation of state) surrounded by a cloud ofscalar field. Different initial configurations are evolved intime using numerical-relativity simulations. We find thatthe fermionic star is able to capture part of the scalar fieldand the new system evolves toward an almost static con-figuration giving rise to a stable fermion-boson star. Inaddition to show that the dynamical formation of mixedstars is possible we also obtain the corresponding equilib-rium configurations for fermion-boson stars with differentvalues of the self-interaction potential and we study theirstability properties under spherical perturbations.This paper is organized as follows: in Section II we in-troduce the matter model we employ to describe fermion-boson stars and set up the basic equations. Section IIIaddresses the initial data for the dynamical formation ofthe mixed stars and the initial static configurations con-sidering a self-interaction term in the bosonic sector. Thenumerical framework for our simulations is described inSection IV while in Section V the results of the evolu-tions are presented. Finally, our conclusions and finalremarks are reported in Section VI. Our units are suchthat the relevant fundamental constants are equal to one( G = c = (cid:126) = 1). II. SETUP
In this study we consider that bosonic and fermionicmatter only interact through gravity. Therefore, ourmodel is described by a total stress-energy tensor whichis the sum of two contributions, one from a perfect fluidand one from a complex scalar field: T µν = T fluid µν + T φµν , (1) where T fluid µν = [ ρ (1 + (cid:15) ) + P ] u µ u ν + P g µν , (2) T φµν = − g µν ∂ α ¯ φ∂ α φ − V ( φ )+ 12 ( ∂ µ ¯ φ∂ ν φ + ∂ µ φ∂ ν ¯ φ ) . (3)The perfect fluid is described by its pressure P , itsrest-mass density ρ , and its internal energy (cid:15) , while u µ isthe fluid 4-velocity. We consider a quartic self-interactionpotential for the scalar field φV ( φ ) = 12 µ ¯ φφ + 14 λ ( ¯ φφ ) , (4)where µ is the mass of the bosonic particle and λ is theself-interaction parameter; the bar symbol in the last twoequations denotes complex conjugation. The equationsof motion are given by the conservation laws of the stress-energy tensor and the baryonic particles ∇ µ T µν fluid = 0 , (5) ∇ µ ( ρu µ ) = 0 , (6)for the fermionic matter, and by the Klein-Gordon equa-tion ∇ µ ∇ µ φ = µ φ + λ | φ | φ (7)for the complex scalar field, together with the Einsteinequation G µν = 8 πT µν governing the spacetime dynam-ics. Differential operator ∇ µ is the covariant derivativewith respect to the 4-metric g µν . The set of equations (5)-(6) is closed by an equation of state (EoS) for the fluid.We consider both the polytropic EoS and the ideal-gasEoS, P = Kρ Γ = (Γ − ρ(cid:15) . (8)The polytropic EoS is employed to build the equilibriuminitial data while the Γ-law is used for the evolutions as itwould allow to take into account eventual shock-heating(thermal) effects. All equilibrium models we consider areconstructed using K = 100 and Γ = 2. In the nextsubsections we specify our choice for the metric and therelevant equations for both the construction of the staticmodels and the evolution. A. Basic equations for the equilibriumconfigurations
Our formalism for the construction of equilibrium con-figurations of fermion-bosn stars relies on the choice ofa spherically symmetric metric in Schwarzschild coordi-nates ds = − α ( r ) dt + ˜ a ( r ) dr + r ( dθ + sin θ dϕ ) , (9)written in terms of two geometrical functions ˜ a ( r ) and α ( r ). We set a harmonic time dependence ansatz forthe complex scalar field φ ( t, r ) = φ ( r ) e − iωt where ω is its eigenfrequency, and we consider the quartic self-interaction potential for the field given by Eq. (4). Wereplace the self-interaction parameter λ by the dimen-sionless variable Λ, defined asΛ = M λ πµ , (10)in which M Pl = (cid:113) (cid:126) cG indicates the Planck mass (which isone in our units). In the following we consider a scaled ra-dial coordinate r → rµ (together with M → M µ , t → tµ , ω → ω/µ ). Assuming a static fluid, u µ = ( − /α, , , d ˜ adr = ˜ a (cid:18) − ˜ a r + 4 πr (cid:20)(cid:18) ω α + µ + λ φ (cid:19) ˜ a φ +Ψ + 2˜ a ρ (1 + (cid:15) ) (cid:21) (cid:19) , (11) dαdr = α (cid:18) ˜ a − r + 4 πr (cid:20)(cid:18) ω α − µ − λ φ (cid:19) ˜ a φ +Ψ + 2˜ a P (cid:21) (cid:19) , (12) dφdr = Ψ , (13) d Ψ dr = − (cid:18) a − πr ˜ a ( µ φ + λ φ + ρ (1 + (cid:15) ) − P ) (cid:19) Ψ r − (cid:18) ω α − µ − λφ (cid:19) ˜ a φ , (14) dPdr = − [ ρ (1 + (cid:15) ) + P ] α (cid:48) α , (15)where the prime indicates the derivative with respectto r . The system is closed by the EoS (8). To solvethese equations it is necessary to apply certain initial andboundary conditions that are consistent with the geome-try and physical behavior of the mixed stars. In SectionIII B we will introduce these conditions. B. Basic equations for the evolution
For the numerical evolutions we consider a sphericallysymmetric metric in isotropic coordinates ds = − α (ˆ r ) dt + ψ (ˆ r ) γ ij ( dx i + β i dt )( dx j + β j dt ) , (16) where α is the lapse function and β i is the shift vector.The spatial 3-dimensional metric components are γ ij dx i dx j = a (ˆ r ) d ˆ r + b (ˆ r )ˆ r ( dθ + sin θ dϕ ) . (17)We note that a and ˜ a should not be confused as they aredifferent functions; a (ˆ r ) and b (ˆ r ) are the metric functionsfor the isotropic metric, ˆ r denotes the isotropic radial co-ordinate (see section V B for details) and ψ ≡ e χ is theconformal factor. To simplify the notation we will substi-tute ˆ r → r in the following, keeping in mind that all equa-tions and definitions refer nonetheless to the isotropic ra-dial coordinate.Our choice of evolution equations for the spacetimevariables follows Brown’s covariant form [50, 51] of theBaumgarte-Shapiro-Shibata-Nakamura (BSSN) formula-tion of Einstein’s equations [52–54]. The evolved quan-tities used in this work are the spatial metric γ ij , theconformal factor χ , the trace of the extrinsic curvature K , its traceless part A a = A rr , A b = A θθ = A ϕϕ , and theradial component of the so-called conformal connectionfunctions ∆ r (see [53, 54] for definitions).We will not report here explicitly the full system ofevolution equations as it can be found e.g. in Ref. [55].We remind the reader that the equations involve mattersource terms arising from suitable projections of the totalstress-energy tensor T µν , namely the energy density E ,the momentum density j i measured by a normal observer n µ , and the spatial projection of the energy-momentumtensor S ij . These quantities read as E = n µ n ν T µν , (18) j i = − γ µi n ν T µν , (19) S ij = γ µi γ νj T µν . (20)In our setup these quantities are obtained by adding upthe contributions of both the fluid and the scalar field.The explicit expressions we use are listed at the end ofthis section.The gauge conditions we employ in our simulations arethe so-called “non-advective 1+log” gauge condition forthe lapse function α and a variation of the Gamma-drivercondition for the shift vector β r . Further details regard-ing the BSSN evolution equations, gauge conditions, andthe formalism for the hydrodynamic equations can befound in [55].Following our previous work [56] we use two auxiliaryvariables Π = 1 α ( ∂ t − β r ∂ r ) φ, (21)Ψ = ∂ r φ, (22)to cast the Klein-Gordon equation (7) as a first-ordersystem of evolution equations: ∂ t φ = β r ∂ r φ + α Π , (23) ∂ t Π = β r ∂ r Π + αae χ (cid:20) ∂ r Ψ + Ψ (cid:18) r − ∂ r a a + ∂ r bb + 2 ∂ r χ (cid:19)(cid:21) + Ψ ae χ + αK Π − α ( µ + λφ ¯ φ ) φ, (24) ∂ t Ψ = β r ∂ r Ψ + Ψ ∂ r β r + ∂ r ( α Π) . (25)Finally, the system of equations is closed by two con-straint equations, namely the Hamiltonian constraint andthe momentum constraint, which read as H = R − ( A a + 2 A b ) + 23 K − π E = 0 , (26) M r = ∂ r A a − ∂ r K + 6 A a ∂ r χ +( A a − A b )( 2 r + ∂ r bb ) − πj r = 0 , (27)where R is the Ricci scalar.The bosonic contribution to the matter source termsare E φ = 12 (cid:18) ¯Π Π + ¯ΨΨ e χ a (cid:19) + 12 µ ¯ φφ + 14 λ ( ¯ φφ ) (28) j φr = −
12 ( ¯ΠΨ + ¯ΨΠ) , (29) S φa = 12 (cid:18) ¯Π Π + ¯ΨΨ e χ a (cid:19) − µ ¯ φφ − λ ( ¯ φφ ) (30) S φb = 12 (cid:18) ¯Π Π − ¯ΨΨ e χ a (cid:19) − µ ¯ φφ − λ ( ¯ φφ ) (31)where S a = S rr and S b = S θθ = S ϕϕ . Correspondingly, thefermionic contribution to those source terms read E fluid = [ ρ (1 + (cid:15) ) + P ] W − P, (32) j fluid r = e χ a [ ρ (1 + (cid:15) ) + P ] W v r , (33) S fluid a = e χ a [ ρ (1 + (cid:15) ) + P ] W v r + P, (34) S fluid b = P, (35)where W = αu t is the Lorentz factor and v r is the radialcomponent of the fluid 3-velocity. III. INITIAL DATA
As mentioned in the introduction we consider two dif-ferent physical situations in this paper, namely the dy-namical formation of a fermion-boson star and the sta-bility properties of different equilibrium models of suchstars. In the following we discuss the corresponding ini-tial data for either situation.
A. Dynamical formation
To study the dynamical formation of a mixed star webegin with a stable fermionic star (FS) model surroundedby a dilute cloud of bosonic particles. This cloud ac-cretes on to the FS under the gravitational pull of thelatter. Suitable initial data describing this system are se-cured after solving the Hamiltonian constraint (26) andthe momentum constraint (27). To do so we assume aharmonic time dependence for the scalar field and choosea Gaussian radial distribution for the cloud, yielding φ ( r, t ) = A e − r σ e − iωt , (36)where parameters A and σ are the amplitude and thewidth of the Gaussian profile, respectively, and ω is theinitial frequency of the field.To solve the constraints we initially consider the space-time of an isolated spherically symmetric FS by solv-ing the Tolman-Oppenheimer-Volkoff equation. Next, weadd to this solution the dilute cloud of bosonic matter de-scribed by (36). The time symmetry condition, K ij = 0,and the conformally flat condition, a = b = 1, yield thefollowing initial values for a set of spacetime variables β r = 0 ,K = 0 ,A a = A b = 0 , ∆ r = 0 . (37)while the values of the conformal factor ψ and of the lapsefunction α are inferred directly from the FS spacetime.Starting with these initial conditions we solve numeri-cally the Hamiltonian constraint (26) using the proce-dure described in [56]. This yields an updated value ofthe conformal factor ψ and of the γ rr metric component.Due to the harmonic time dependence of the scalarfield, it follows that j φr defined by (29) is zero. Thismeans that the scalar field does not contribute to the mo-mentum constraint equation (27). Therefore, consider-ing (37) the momentum constraint is analytically solved. B. Equilibrium configurations
In Section II A we introduced the basic equations toconstruct the static models of mixed stars. To solve theset of equations, namely equations (11)-(15) and the EoS(8), we need to construct suitable initial data which arecompatible with the physical and geometrical conditionsof the stellar configurations. The system of ODEs be-comes an eigenvalue problem for the frequency ω , whichis a function of two parameters, the central value of thescalar field, φ c , and of the fermionic density, ρ c . Wemake use of the two-parameter shooting method to findthe solution for ω . Once this is found and the centralvalues of all variables are available, we use a 4th-orderRunge-Kutta method to solve the ODEs and reconstructthe radial profiles of the solution.We require the condition of regularity at the origin tobe satisfied for the metric functions. At the outer bound-ary we employ the values provided by the Schwarzschildsolution at the outer radius, which do not depart muchfrom the values of a flat metric, together with a vanish-ing scalar field value. Hence, the boundary conditions forsolving the set of ODEs can be defined as follows˜ a (0) = 1 , φ (0) = φ c ,α (0) = 1 , lim r →∞ α ( r ) = lim r →∞ a ( r ) , Ψ(0) = 0 , lim r →∞ φ ( r ) = 0 ,ρ (0) = ρ c , P (0) = Kρ Γ c , lim r →∞ P ( r ) = 0 . (38)Once the solution is found, one can define the totalgravitational mass based on the value of the metric coef-ficients at infinity M T = lim r −→∞ r (cid:18) − a (cid:19) , (39)which coincides with the Anowitt-Desser-Misner (ADM)mass at infinity. As the Klein-Gordon Lagrangian for acomplex scalar field exhibits invariance under global U(1)transformations φ → φ e iδ , Noether’s theorem predictsthe existence of a conserved charge which can be associ-ated with the number of bosonic particles N B ; moreover,the conservation of the baryonic number provides a def-inition of the number of fermionic particles N F . Thesetwo quantities can be evaluated by integrating their vol-ume density as follows N B = 4 π (cid:90) ˜ aωφ r α dr, N F = 4 π (cid:90) ˜ aρr dr. (40)These quantities will be used to determine the conserva-tion of the number of particles, both bosons and fermionsduring the numerical evolutions. Finally, we evaluatethe radius of the bosonic (fermionic) contribution to themixed star, R B ( R F ), as the radius of the sphere contain-ing 99% of the corresponding particles.As mentioned before, the construction of the static so-lutions for the fermion-boson stars depends on two pa-rameters, namely the central fluid density ρ c and the cen-tral value of the scalar field φ c . We can therefore expressthe mass of the system (39) as a function of these twoparameters M T ( ρ c , φ c ) as we depict in figure 1 for threedifferent values of Λ. In the case of non-rotating bosonstars the parameter space is 1-dimensional and stabil-ity theorems [57] indicate that for each value of Λ thereexists a critical mass such that dM T /dφ c = 0. These crit-ical points indicate the transitions between the stabilityand the instability region of the parameter space. Analo-gous transitions in stability occur in fermionic stars (seee.g. [58, 59]). In the case of fermion-boson stars, as theparameter space is 2-dimensional, the analysis is more φ c ρ c M T =0.6M T =0.8M T =1.0M T =1.2M T =1.4M T =1.5M T =1.55M T =1.6M c =1.637M c =0.248 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 φ c ρ c M T =0.5M T =0.6M T =0.613M T =0.63M T =1.0M T =1.2M T =1.4M T =1.6M c =1.637M c =0.633 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 φ c ρ c M T =0.6M T =0.8M T =1.05M T =1.2M T =1.3M T =1.4M T =1.5M T =1.6M c =1.637M c =1.336 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 FIG. 1: Equilibrium configurations of fermion-boson stars forΛ = −
30 (top), Λ = 0 (middle), and Λ = 30 (bottom).Dashed lines correspond to models with the same total mass M T . The black solid line depicts the boundary between sta-ble and unstable models, and the solid yellow line for the caseΛ = −
30 indicates the maximum value of φ c that assures thenon-negativity of the scalar field potential V ( φ ) in the entirespatial domain. involved. Following [32] we define the critical points asthe values of the pair ( ρ c , φ c ) such that the conditions ∂N B ∂ρ c (cid:12)(cid:12)(cid:12) M =constant = ∂N F ∂ρ c (cid:12)(cid:12)(cid:12) M =constant = 0 ,∂N B ∂φ c (cid:12)(cid:12)(cid:12) M =constant = ∂N F ∂φ c (cid:12)(cid:12)(cid:12) M =constant = 0 , (41)are satisfied. In Fig. 1 we show several curves of constantmass in the parameter space (dashed colored lines). Foreach point of the curves we evaluate the number of bosons N B and fermions N F . If we start from a purely FS config-uration (a point on the horizontal axis in Fig. 1) and wemove along a curve of fixed mass changing the values of φ c and ρ c , the number of bosons increases and the num-ber of fermions deacreases up to a critical point in theparameter space where a maximum is found for N B anda minimum for N F . If we start from a pure boson star (apoint on the vertical axis), the behavior is the opposite,with N B decreasing up to a minimum and N F increas-ing up to a maximum. For each value of the mass, thesecritical points signal the boundary between the stabilityand instability regions. The black solid line in Fig. 1 rep-resents these boundaries in the parameter space for thevalues of Λ = {− , , } . This construction follows thesame approach laid out in [33, 36].As FS do not depend on Λ their threshold mass isconstant for all values of Λ and equal to M c = 1 . M c = 0 . . . −
30, 0, and 30, respectively. For fermion-bosonstars, one can observe that for the same point in theparameter space with fixed values of φ c and ρ c , the totalmass decreases (increases) for positive (negative) valuesof Λ, with respect to the Λ = 0 case.We point out that considering negative values of λ raises the issue that the scalar potential V ( | φ | ) = µ | φ | + λ ( | φ | ) is not bounded from below andcan become negative, breaking the weak-energy condi-tion (see e.g. the discussion in [60]). For Λ = −
30 weevaluate the maximum central value of φ that ensuresthe non-negativity of the scalar field potential, yielding φ c = 0 . IV. NUMERICAL FRAMEWORK
The numerical evolutions of the Einstein-Klein-Gordon-Euler system are performed with the numerical-relativity code originally developed by [55] and up-graded to take into account the complex scalar-field equa-tions in [61]. This computational infrastructure hasbeen extensively used by our group in studies of fun-damental bosonic fields in strong-gravity spacetimes (seee.g. [38, 45, 56, 62, 63]).The time update of the evolution equations is evalu-ated using a Partially Implicit Runge-Kutta method de-veloped by [64, 65]. In this scheme the operators in theright-hand-side of the BSSN evolution equations are di-vided into operators which are evaluated explicitly, andoperators carrying geometrical singularities which areevolved implicitly using the updated values of the firstones. This allows to handle potential numerical instabil- ities arising from 1 /r terms in the equations. While theconstruction of the equilibrium configurations employsSchwarzschild coordinates and an equally spaced lineargrid, the dynamical evolutions make use of isotropic coor-dinates and a logarithmic grid. More precisely, the com-putational domain of the simulations is covered with anisotropic grid which is composed by two different patches,a geometrical progression up to a certain radius and anhyperbolic cosine in the exterior part. This allows toplace the outer boundary sufficiently far from the ori-gin and prevent the effects of reflections. Further de-tails about the computational grid can be found in [62].The minimum resolution we employ in our simulationsis ∆ r = 0 . r min = ∆ r/ r max = 6000.The time step is given by ∆ t = 0 . r in order to obtainlong-term stable simulations. We add 4th-order Kreiss-Oliger numerical dissipation terms to the evolution equa-tions to damp out spurious, high-frequency numericalnoise. All advection terms (such as β r ∂ r f ) are treatedwith an upwind scheme. At the outer boundary we im-pose radiative boundary conditions. V. RESULTSA. Dynamical formation of fermion-boson stars
As described in section III A we start with an initialconfiguration describing a bosonic cloud of matter sur-rounding an already formed FS, and we study the ac-cretion of the bosonic matter on to the FS. The bosoniccloud loses part of its energy through gravitational cool-ing and plunges towards the center of the FS. Intuitively,this process can lead to two possible outcomes: eitherto the formation of a fermion-boson star or, if the massof the entire system is above a certain threshold, to theformation of a Schwarzschild black hole.During the evolutions we compute useful physicalquantities in order to keep track of the formation processand to evaluate the features of the final object. Those willbe used below to compare with some of our static models.We define the bosonic and fermionic energy contained inspheres of different radii r ∗ as E fluid r ∗ = 4 π (cid:90) r ∗ E fluid √ γdr, (42) E φr ∗ = 4 π (cid:90) r ∗ E φ √ γdr, (43)where √ γ = ψ √ abr is the spatial volume element forthe metric (16). Note that we will refer to E fluid /φr max whenreferring to the total energy in the computational grid.Other useful quantities we evaluate along the numericalevolution are the number of bosonic and fermionic par-ticles within spheres of radii r ∗ , computed by means of TABLE I: Initial models leading to stable fermion-boson stars. The vertical lines divide the initial parameters (left), from thephysical quantities evaluated at the end of the time evolution (center) and from the physical quantities of the correspondingequilibrium configuration (right). Note that as model MS5 forms an excited state, we cannot compare it with a nodeless staticsolution. Columns on the left indicate the central rest-mass density ρ c , the self-interaction parameter Λ, and the amplitude ofthe scalar field profile A . Columns at the center indicate the scalar field frequencies ω n , the fermionic energy E fluid30 containedin a sphere of radius r = 30, the bosonic energy E φ , and the ratio between number of bosons and fermions N B /N F . Columnson the right side indicate the frequency ω , the fermionic energy E fluid , the bosonic energy E φ and the ratio between number ofbosons and fermions N B /N F of the corresponding equilibrium configuration.Model ρ c Λ A ω ω E fluid30 E φ N B /N F ω E fluid E φ N B /N F MS1 1 . × − . × − . × −
30 4 . × − . × − -30 4 . × − . × − . × − . × −
30 4 . × − ρ c , the self-interaction parameter Λ, and the amplitude of the scalar fieldprofile A . Model ρ c Λ A MS6 3 . × − -30 3 . × − MS7 3 . × − . × − MS8 3 . × −
30 3 . × − the following integrals N Br ∗ = 4 π (cid:90) r ∗ g ν J ν α √ γdr, (44) N Fr ∗ = 4 π (cid:90) r ∗ ρ √ γdr, (45)where J ν = i ( ¯ φ ∂ ν φ − φ ∂ ν ¯ φ ) is the conserved current as-sociated with the transformation of the U(1) group. Wealso extract the scalar-field frequency ω by performing aFast Fourier transform (FFT) of the real/imaginary partof the scalar field φ . The time window for the FFT ischosen at a sufficiently late time of the evolutions, oncethe bosonic cloud has already accreted on to the FS andthe final object oscillates around an equilibrium configu-ration.For our study we use two different FS models, bothdescribed by the polytropic EoS, P = Kρ Γ , with differentcentral value of the rest-mass density ρ c . We considerthe same scalar-field mass parameter, µ = 1, frequency, ω = 0 .
8, and three different values for the self-interactionparameter Λ = {− , , +30 } . Our model for the bosoniccloud, equation (36), has a couple of free parameters wecan vary, namely the amplitude A and the width of theGaussian profile σ . For all our models we consider σ = 90which corresponds to a bosonic cloud much larger thanthe FS radii. We summarize some of the properties ofour initial models in Table I and Table II.In Fig. 2 we show the evolution of the scalar-field en- ergy contained in spheres of different radii r ∗ calculatedwith Eq. (43), for models MS3, MS4 and MS5 describedin Table I. The growth of the lines E φ , E φ , E φ , and E φ shows that during the evolution the energy of thescalar cloud, which at the initial time is spread over alarge spatial volume, gradually concentrates in a smallervolume, as it is being accreted by the FS. Part of thecloud energy does not fall on to the FS but it is radiatedaway through the gravitational cooling mechanism. Forall three models, from time t (cid:39)
750 the curves start toconverge slowly to each other, indicating that the scalarfield is contained within small radii, radiating the excessenergy to infinity. The remnant energy is confined intoa volume delimited by r (cid:39) E φ (cid:39) .
065 for MS3, E φ (cid:39) .
09 for MS4,and E φ (cid:39) .
11 for MS5.Figure 2 shows some differences between models withand without self-interaction, and also depending on thesign of the self-interaction term. In the case with Λ = −
30 (top panel) the lines E φ , E φ , E φ and E φ slowlyconverge to each other and, at around t (cid:39) E φ ) with a final energyaround E φ (cid:39) . E φ , grows reaching the green one, E φ . This indicates that all the scalar matter aroundthe forming compact object is accreting onto it. Thecase with Λ = 0 (central panel) is an intermediate case,with the lines slowly converging to each other for the en-tire evolution. This result can be understood as follows:a Λ > < Λ =-30 Λ =0 Λ =30 FIG. 2: Scalar-field energy in spheres of different radii, formodel MS3 (top), MS4 (center), and MS5 (bottom). The red,blue, gold and green lines, correspond to r ∗ = 10 , , , finity. Nonetheless the formation process seems to beaccelerated but it is due to the fact that the scalar par-ticles around the formed compact object escape faster toinfinity. We point out that, as | φ | <
1, the self-interactionterm, which is proportional to λ | φ | , gives a lower ordercontribution than the mass term µ | φ | . They are only comparable when the object is compact enough to reachhigh values of φ . This is the reason why the first partof the evolution before the object forms is basically thesame for the three models.In Fig. 3 we depict the late-time radial profiles of thescalar field module | φ | for models MS3 (top panel) andMS5 (bottom panel). For model MS3 we compare threedifferent snapshots during the evolution with an equi-librium configuration of a mixed star with comparablemass and number of bosons and fermions. The compari-son shows that the radial profile of | φ | obtained throughthe dynamical formation process resembles that of thestatic solution.The bottom panel of Fig. 3 shows that for model MS5there are two maxima of the scalar field and there is anode at around r (cid:39) B. Evolutions of the equilibrium configurations
In Section II A we discussed how we identify the re-gion of the parameter space where stable configurationsare found. In this section we intend to verify the resultsobtained by performing numerical evolutions of stableand unstable models. We expect stable mixed stars toshow a combination of the typical behaviour of isolatedstable boson stars and fermion stars. This means thatwe expect the scalar field to oscillate with its character-istic eigenfrequency ω while the fermionic density ρ isexpected to oscillate slightly around its initial state due | φ | r time = 5100time = 5400time = 6075static model | φ | r time = 4800time = 5400time = 6075 FIG. 3: Late-time radial profiles of the scalar field module | φ | for model MS3 with Λ = −
30 (upper panel) and modelMS5 with Λ = 30 (bottom panel). The three snapshots ofmodel MS3 are compared with the radial profile of a staticmixed star model of similar ρ c , φ c , and bosonic and fermionicenergy and number (dashed black line in the plot). ModelMS5 presents a node at r (cid:39) to the numerical truncation errors introduced by the dis-cretization of the differential equations of the continuummodel. All physical quantities of the stable models suchas mass, boson number density or fermion number den-sity are expected to be constant in time. Even under theintroduction of a small perturbation, stable models areexpected to oscillate around their static solutions.For a model in the unstable region, however, we expectthe small-amplitude perturbations induced by the numer-ical errors to grow due to the non-linearity of the sys-tem. The growth of the perturbations can lead to threedifferent outcomes: the migration to the stable region,the gravitational collapse and formation a Schwarzschildblack hole, or the dispersion of the bosonic particles.As our evolution code is based on isotropic coordinates(16) and the mixed-star models are constructed usingSchwarzschild coordinates (9), we must apply a coordi-nate transformation to be able to evolve the initial con-figurations. We follow the procedure proposed in [67] FIG. 4: Evolution of the radial profile of the scalar fieldmodule | φ | for model MS3 with Λ = −
30 (upper panel) andmodel MS5 with Λ = 30 (bottom panel) in the time window t ∈ [3000 , which can be divided in two steps. First, we perform thechange of coordinates noting that from the comparisonbetween the two metrics we have that d ˆ rdr = ˜ a ( r ) ˆ rr . (46)To obtain the coordinate transformation we introduce thefunction β β = ˆ rr . (47)Rewriting equation (46) in terms of ln β we obtain d ln βdr = 1 r (˜ a ( r ) − , (48)which leads to β ( r ) = exp (cid:20) − (cid:90) r max r r (cid:48) (˜ a ( r (cid:48) ) − dr (cid:48) (cid:21) . (49)As initial condition to solve this integral, we imposethat at the outer boundary the spacetime resembles the0 ρ c φ c t ρ c φ c ρ c φ c t ρ c φ c ρ c φ c t ρ c φ c N B N F t N B N F N B N F t N B N F M AH / M AD M t FIG. 5: Time evolution of static models with self-interaction parameter Λ = 30. Left panels depict the central value of the fluiddensity ρ c and of the scalar field φ c (top row) and number of bosons N B and fermions N F (bottom row) for the stable modelMS11. Middle panels show the same physical quantities for the unstable model MS12 without the addition of an artificialperturbation. The right panels show the collapse to a Schwarzschild black hole of model MS12 when a 2% perturbation isinduced in the scalar field. The right bottom plot displays the AH mass in units of the ADM mass (red solid line) and the timeevolution of the ADM mass normalized by its initial value (black dashed line).TABLE III: Static fermion-boson star models. From left to right the columns indicate the model name, its stability, the value ofthe self-interaction parameter Λ, the central value of the fluid density ρ c and of the scalar field φ c , the field frequency obtainedwith the shooting method ω shoot , the normalized frequency ω , the number of bosons to fermions ratio N B /N F , the number ofbosons N B , the radius containing 99% of bosons, fermions and total particles, R B , R F , R T , respectively. All radii are evaluatedusing Schwarzschild coordinates.Model Branch Λ ρ c φ c ω shoot ω N B /N F N B R B R F R T MS9 Stable 0 1 . × − . × − .
199 0 .
732 0 .
191 0 .
208 5 .
51 7 .
22 7 . . × − . × − .
436 0 .
602 0 .
251 0 .
261 3 .
66 5 .
81 5 . . × − . × − .
238 0 .
807 0 .
344 0 .
308 7 .
06 7 .
08 7 . . × − . × − .
629 0 .
815 13 .
03 1 .
134 6 .
93 3 .
06 6 . . × − . × − .
129 0 .
680 0 .
055 0 .
079 5 .
18 7 .
73 7 . . × − . × − .
099 0 .
603 0 .
068 0 .
097 3 .
99 7 .
34 7 . Schwarzschild solution which yieldsˆ r max = (cid:32) (cid:112) ˜ a ( r max )2 (cid:33) r max ˜ a ( r max ) . (50)Once we obtain β , we can finally obtain the conformalfactor which is defined as ψ = (cid:114) r ˆ r = (cid:114) β . (51)We point out that the introduction of the new variable β is necessary to make the integral (49) behave well at theorigin, and to be able to reconstruct the solution in the entire radial domain. The interested reader is addressedto [67] for further details.We perform evolutions of several models for values ofthe self-interaction parameter Λ = {− , , } , both inthe stable and unstable region of the existence surface.These numerical evolutions confirm our analysis aboutthe stability of the models. We summarize their relevantphysical properties in Table III.Figure 5 shows the time evolution of the results ob-tained for the case Λ = 30, in particular models MS11and MS12 of Table III. In the left panels we display theevolution of the central value of the fluid density ρ c andof the scalar field φ c (top row) and the evolution of thenumber of fermions and bosons (bottom row), for the1stable model MS11. As expected all these physical quan-tities remain constant in time confirming that the modelis stable. The middle panels show the time evolution ofthe same physical quantities for model MS12, which is inthe unstable region. We can observe that the central val-ues of the scalar field and the fluid density very rapidlydepart from their initial values, with a large variationwhich is damped in a few cycles. The system settles on anew configuration in the stable branch, oscillating aroundthe new central values ρ c (cid:39) . φ c (cid:39) . φ ( r ) → φ ( r ) (cid:18) A (cid:19) , (52)where A = 2, which corresponds to a 2% level pertur-bation. Despite fairly small, this artificial perturbationis stronger than that introduced by the discretization er-rors alone which triggered the evolution shown in themiddle panels of Fig. 5. We now observe that due tothe stronger perturbation the model does not migrate tothe stable region but rather collapses to a Schwarzschildblack hole, as signalled by the formation of an apparenthorizon (AH). In the top row we show the time evolutionof the central values of the fluid density and of the scalarfield while in the bottom row we show the time evolutionof the mass of the black hole evaluated on the AH in unitsof the ADM mass of the system (which we depict witha dashed black curve). We could not find any modelfor which the bosonic part dispersed, leaving behind apurely FS. The binding energy of the whole configura-tion is never positive and therefore, unstable models canonly either migrate or collapse. VI. CONCLUSIONS
Fermion-boson stars are gravitationally bound struc-tures composed by fermions and scalar particles. Theyare regular and static macroscopic configurations ob-tained by solving the coupled Einstein-Klein-Gordon-Euler system. In this paper we have discussed a possiblescenario through which fermion-boson stars may form as-suming an initial configuration in which an already exist-ing FS (i.e. a neutron star) is surrounded by an accretingdilute cloud (a Gaussian pulse) of a massive, complexscalar field. Our setup has considered positive and nega-tive values of a quartic self-interaction term in the Klein-Gordon potential. We have built constraint-satisfyinginitial data and we have modelled the astrophysical sit-uation by considering different bosonic cloud amplitudesand widths and two different fermion star models. The results of our spherically-symmetric, numerical-relativitysimulations have shown that once part of the initial scalarfield is expelled via gravitational cooling the system oscil-lates around an equilibrium configuration that is asymp-totically consistent with the static solutions of the sys-tem.Existence diagrams of such equilibrium solutions in thecentral-field-amplitude vs central-fermionic-density planehave been constructed to draw such comparisons. Our re-sults are in agreement, in the corresponding limits, withthe work of [33, 36]. The non-linear stability of staticmodels residing in both the stable and unstable regionsof the existence diagrams has been assessed through sim-ulations with a quartic self-interaction potential in thebosonic sector, not attempted in previous works. Thosehave shown that, for stable configurations, all physicalquantities describing the star, such as energy and num-ber of particles, remain constant during the evolution,while unstable models either migrate to the stable regionor collapse to a Schwarzschild black hole.The dynamical formation of fermion-boson stars forlarge positive values of the coupling constant in the quar-tic self-interaction term (namely Λ = 30) has revealed thepresence of a node in the scalar field. This is an intrigu-ing result as purely boson stars with nodes correspondto excited states and are known to be intrinsically un-stable [43, 57]. However, fermion-boson stars with nodesin the bosonic sector can dynamically form and appearlong-term stable. This indicates that an excited state ofthe scalar field in the presence of fermionic matter mayform a stable configuration. This result is akin to thefindings of [66] who found that boson star configurationsin which the ground state and the first excited state ofthe scalar field coexist are stable. In upcoming investi-gations we plan to build equilibrium fermion-boson con-figurations with an excited state of the scalar field andstudy their stability properties to confirm the result re-ported here. Likewise, we will analyze the dynamicalformation of rotating mixed stars as it might as well bepossible that the presence of fermionic matter stabilizedotherwise unstable spinning boson stars [46].
Acknowledgments
We thank Eugen Radu and Carlos Herdeiro foruseful suggestions. This work was supported bythe Spanish Agencia Estatal de Investigaci´on (grantPGC2018-095984-B-I00), by the Generalitat Valenciana(PROMETEO/2019/071 and GRISOLIAP/2019/029),by the European Unions Horizon 2020 RISE programmeH2020-MSCA-RISE-2017 Grant No. FunFiCO-777740,by DGAPA-UNAM through grants No. IN110218,IA103616, IN105920, by the Funda¸c˜ao para aCiˆencia e a Tecnologia (FCT) projects PTDC/FIS-OUT/28407/2017 and UID/FIS/00099/2020 (CEN-TRA), and CERN/FIS-PAR/0027/2019. SF gratefullyacknowledges support by the Erasmus+ International2Credit Mobility Program KA-107 for an academic stayat the University of Valencia. [1] S. Weinberg, Phys. Rev. Lett. , 223 (1978).[2] J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B , 127 (1983).[3] T. Matos and L. A. Urena-Lopez, Phys.Rev. D63 ,063506 (2001), astro-ph/0006024.[4] T. Matos and L. A. Urena-Lopez, Class. Quant. Grav. , L75 (2000), astro-ph/0004332.[5] M. Gasperini and G. Veneziano, Phys. Rev. D , 2519(1994), gr-qc/9403031.[6] P. Svrcek and E. Witten, JHEP , 051 (2006), hep-th/0605206.[7] P. W. Higgs, Phys. Lett. , 132 (1964).[8] G. Aad et al. (ATLAS), Science , 1576 (2012).[9] A. H. Guth, Adv. Ser. Astrophys. Cosmol. , 139 (1987).[10] D. Langlois, in Cargese School of Particle Physics andCosmology: the Interface (2004), pp. 235–278, hep-th/0405053.[11] M. Kawasaki and K. Nakayama, Ann. Rev. Nucl. Part.Sci. , 69 (2013), 1301.1123.[12] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper,and J. March-Russell, Phys.Rev. D81 , 123530 (2010),0905.4720.[13] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys.Rev.
D95 , 043541 (2017), 1610.08297.[14] V. B. Klaer and G. D. Moore, JCAP , 049 (2017),1708.07521.[15] S.-J. Sin, Phys. Rev. D , 3650 (1994), hep-ph/9205208.[16] P.-H. Chavanis and T. Harko, Phys. Rev. D , 064011(2012), 1108.3986.[17] T. Matos, F. S. Guzman, and L. A. Urena-Lopez, Class.Quant. Grav. , 1707 (2000), astro-ph/9908152.[18] W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. , 1158 (2000), astro-ph/0003365.[19] P. Jetzer, Phys. Rept. , 163 (1992).[20] D. J. Kaup, Phys. Rev. , 1331 (1968).[21] R. Ruffini and S. Bonazzola, Phys. Rev. , 1767(1969).[22] M. Colpi, S. L. Shapiro, and I. Wasserman, Phys. Rev.Lett. , 2485 (1986).[23] P. Jetzer and J. van der Bij, Phys. Lett. B , 341(1989).[24] S. Yoshida and Y. Eriguchi, Phys. Rev. D56 , 762 (1997).[25] F. E. Schunck and E. W. Mielke, Phys. Lett.
A249 , 389(1998).[26] E. Seidel and W. Suen, Physical Review Letters , 384(1991).[27] M. Alcubierre, J. Barranco, A. Bernal, J. C. Degol-lado, A. Diez-Tejedor, M. Megevand, D. Nunez, andO. Sarbach, Class. Quant. Grav. , 19LT01 (2018),1805.11488.[28] V. Jaramillo, N. Sanchis-Gual, J. Barranco, A. Bernal,J. C. Degollado, C. Herdeiro, and D. Nez (2020),2004.08459.[29] R. Brito, V. Cardoso, C. A. Herdeiro, and E. Radu,Physics Letters B , 291 (2016).[30] F. E. Schunck and E. W. Mielke, Class. Quant. Grav. , R301 (2003), 0801.0307.[31] S. L. Liebling and C. Palenzuela, Living reviews in rela-tivity , 5 (2017).[32] A. Henriques, A. R. Liddle, and R. Moorhouse,Physics Letters B , 511 (1990), ISSN 0370-2693, URL .[33] S. Valdez-Alvarado, C. Palenzuela, D. Alic, and L. A.Ure˜na L´opez, Physical Review D , 084040 (2013).[34] R. Brito, V. Cardoso, and H. Okawa, Physical reviewletters , 111301 (2015).[35] R. Brito, V. Cardoso, C. F. Macedo, H. Okawa, andC. Palenzuela, Physical Review D , 044045 (2016).[36] S. Valdez-Alvarado, R. Becerril, and L. A. Ure˜na-L´opez,arXiv preprint arXiv:2001.11009 (2020).[37] E. Seidel and W.-M. Suen, Phys. Rev. Lett. , 2516(1994), gr-qc/9309015.[38] F. Di Giovanni, N. Sanchis-Gual, C. A. R. Herdeiro, andJ. A. Font, Phys. Rev. D98 , 064044 (2018), 1803.04802.[39] T. D. Lee and Y. Pang, Nuclear Physics B , 477(1989).[40] S. H. Hawley and M. W. Choptuik, Phys. Rev.
D62 ,104024 (2000), gr-qc/0007039.[41] M. Gleiser, Phys. Rev.
D38 , 2376 (1988), [Erratum:Phys. Rev.D39,no.4,1257(1989)].[42] M. Gleiser and R. Watkins, Nucl. Phys.
B319 , 733(1989), gr-qc/9905067.[43] J. Balakrishna, E. Seidel, and W.-M. Suen, Phys. Rev.
D58 , 104004 (1998), gr-qc/9712064.[44] F. Guzman, Revista Mexicana de Fisica , 321 (2009).[45] N. Sanchis-Gual, C. Herdeiro, E. Radu, J. C. Degollado,and J. A. Font, Phys. Rev. D , 104028 (2017).[46] N. Sanchis-Gual, F. Di Giovanni, M. Zilh˜ao, C. Herdeiro,P. Cerd´a-Dur´an, J. A. Font, and E. Radu, Physical Re-view Letters , 221101 (2019).[47] E. Seidel and W. Suen, Phys. Rev. D42 , 384 (1990).[48] F. S. Guzman and L. A. Urena-Lopez, Phys. Rev. D ,124033 (2004), gr-qc/0404014.[49] F. Guzman and L. Urena-Lopez, Astrophys. J. , 814(2006), astro-ph/0603613.[50] J. D. Brown, Phys. Rev. D , 104029 (2009), URL http://link.aps.org/doi/10.1103/PhysRevD.79.104029 .[51] M. Alcubierre and M. D. Mendez, Gen.Rel.Grav. ,2769 (2011), 1010.4013.[52] T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor.Phys. Suppl. , 1 (1987).[53] M. Shibata and T. Nakamura, Phys. Rev. D , 5428(1995).[54] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D ,024007 (1998).[55] P. J. Montero and I. Cordero-Carrion, Phys.Rev. D85 ,124037 (2012), 1204.5377.[56] N. Sanchis-Gual, J. C. Degollado, P. J. Montero, andJ. A. Font, Phys. Rev. D , 043005 (2015), 1412.8304.[57] T. D. Lee and Y. Pang, Nucl. Phys. B315 , 477 (1989),[,129(1988)]. [58] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astro-phys. J. , 823 (1994).[59] J. L. Friedman, J. R. Ipser, and R. D. Sorkin, Astrophys.J. , 722 (1988).[60] C. Barcel´o and M. Visser, Classical and Quantum Grav-ity , 3843 (2000), URL https://doi.org/10.1088%2F0264-9381%2F17%2F18%2F318 .[61] A. Escorihuela-Tom`as, N. Sanchis-Gual, J. C. Degollado,and J. A. Font, Physical Review D , 024015 (2017).[62] N. Sanchis-Gual, J. C. Degollado, P. J. Montero, J. A.Font, and V. Mewes, Phys. Rev. D , 083001 (2015),1507.08437.[63] N. Sanchis-Gual, J. C. Degollado, P. J. Montero, J. A.Font, and C. Herdeiro, Phys. Rev. Lett. , 141101 (2016), 1512.05358.[64] I. Cordero-Carri´on and P. Cerd´a-Dur´an, ArXiv e-prints(2012), 1211.5930.[65] I. Cordero-Carri´on and P. Cerd´a-Dur´an, Advances inDifferential Equations and Applications , SEMA SIMAISpringer Series Vol. 4 (Springer International PublishingSwitzerland, Switzerland, 2014).[66] A. Bernal, J. Barranco, D. Alic, and C. Palenzuela, Phys.Rev. D , 044031 (2010), URL https://link.aps.org/doi/10.1103/PhysRevD.81.044031 .[67] B. Kleihaus and J. Kunz, Phys. Rev. D ,834 (1998), URL https://link.aps.org/doi/10.1103/PhysRevD.57.834https://link.aps.org/doi/10.1103/PhysRevD.57.834