On the formation of"supermassive"neutron stars and dynamical transition to spontaneous scalarization

It is well known that neutron stars can undergo a phase transition under a certain class of Scalar Tensor Theories of gravity (STT's) where a new order parameter, the {\it scalar charge}, appears within the star. This is the well known phenomenon of spontaneous scalarization (SC) discovered by Damour and Esposito-Far\`ese in 1993. Under such mechanism neutron stars can afford in principle a maximum mass larger than in general relativity (GR) for a given equation of state without taking into account additional observational constraints (e.g. binary systems). This opens the possibility that neutron stars might be formed with masses as large as $\sim 2 M_\odot$ without the need of stiff, or more exotic, equations of state for the nuclear matter. Thus, STT's through SC may account for compact objects with large masses observed recently in the sky in the form of pulsars (PSR J0348+0432 with $M= 2.01 M_\pm 0.04\odot$ observed in 2013, PSR J1614-2230 with $M= 1.97\pm 0.04 M_\odot$ observed in 2010 or J0740+6620 $M= 2.14^{+0.10}_{-0.09} M_\odot$ observed in 2019). However, we argue that even if that was possible such maximum mass models within STT cannot be formed solely from the dynamic transition of an initial"isolated"unscalarized neutron star whose mass cannot exceed the maximum mass in GR. This is because SC, being an energetic-preferred configuration, produces a final static star with a mass lower that the initial one with a fixed baryon mass. The mass difference between the initial and final configurations is radiated away in the form of a scalar-field. Thus, maximum mass models of scalarized neutron stars, if present in nature, must have formed by a different process, perhaps of cosmological origin or by the subsequent accretion of additional scalar charge and mass.

Introduction. Relativistic theories of gravity different from general relativity (GR) [1][2][3][4][5] have received a major interest in recent years for various reasons. For instance, several phenomena in our universe can be explained by the presence of some unknown form of matter and energy, dubbed dark [6]. The anisotropies of temperature variations in the cosmic microwave background, the flat rotation curves of galaxies, different clustering properties of large scale structures in the cosmos, and the accelerated expansion in the universe, are perhaps the most conspicuous phenomena that require dark components. However, the failure of a direct detection of some of these dark substances by several dark-matter laboratories, as well as the unsuccessful efforts to accommodate the cosmological constant (the simplest form of dark energy) in a consistent manner within the framework of a robust quantum field theory (see [7] for a thorough review), have led to an alternate direction consisting in modifying the gravitational sector, i.e., modifying GR [8][9][10][11][12][13]. Nonetheless, the amazing success of GR at different scales makes also difficult for alternative theories of gravity (ATG) to reproduce in a fully consistent fashion the standard tests of GR while providing the required features without dark sectors. Scalar tensor theories of gravity (STT) in its simplest form (see below) [9,[14][15][16] are one of the few ATG that remains as a viable candidate to generalize Einstein's GR while making clear cut predictions that can be falsified by several kind of experiments. At cosmological level, they can produce an EOS for the dark energy (in the form of a scalar-field) which can vary in cosmic time, unlike the cosmological constant Λ. This variation will be tested in the forthcoming years by numerous projects [17][18][19][20] which will help to validate or rule out a simple ΛCDM model. Since STT are a generalization of the famous Brans-Dicke theory, these theories predict small variations of the effective gravitational coupling G eff , which can be further tested observationally. Furthermore, STT predicts the existence of an additional polarization mode of gravitational waves, the breathing mode, due to the presence of a scalar-field that is non-minimally coupled (NMC) to gravity which can be detected or constrained in the future [21]. This type of polarization mode can be revealed itself in the strong gravity regime even in spherical symmetry, unlike the usual tensor polarizations (h + , h × ) that require the presence of sources with relatively large quadrupolar time variations. The natural place to look for that kind of gravitational radiation is in neutron stars, notably in pulsars. It is well known since the early 1990's, from the pioneer analysis by Damour and Esposito-Farése (DEF) [22,23] that neutron stars can undergo a phase transition in the form of spontaneous scalarization (SC) where a new order parameter, dubbed scalar charge, can appear within the star if the latter is sufficiently compact. Therefore, at the threshold of this transition corresponding to a critical central energy-density (or alternatively a critical total baryon mass) the star transits to a configuration of lower energy (i.e. lower gravitational mass) while keeping the total number of baryons fixed. The difference of gravitational mass between the initial (unstable) configuration and the final static scalarized config-uration is radiated away in the form of scalar gravitational waves of the sort described above. Following the DEF discovery, subsequent analyses showed that SC is robust in that it is independent of the equation of state (EOS) for the nuclear matter adopted to model a neutron star [24][25][26][27]. Finally, a notable feature of SC is that maximum mass models of neutron stars constructed within some specific classes of SST can be much larger than the corresponding models in GR while keeping fixed the EOS [22,24]. The recent observations of massive neutrons stars in the form of pulsars (PSR J0348+0432 with M = 2.01 ± 0.04M ⊙ [28], PSR J1614-2230 with M = 1.97 ± 0.04M ⊙ [29] or J0740+6620 M = 2.14 +0. 10 −0.09 M ⊙ [30]; if confirmed PSR J2215+5135 might host the largest neutronstar mass observed to date M = 2.27 +0.17 −0.15 M ⊙ [31]) put stringent constraints on the current EOS of nuclear matter used to model neutron stars in GR [32][33][34]. This opens the door for the large observed masses of neutron stars to be explained by the phenomenon of SC in STT without the requirement of very stiff or more exotic EOS as it would be the case if one works under the framework of pure GR .
Nevertheless, and this is the most important contribution of this letter, we argue that if such massive neutron stars are explained by the SC transition, the latter must be accompanied by an additional process, maybe of cosmological origin or by accretion, because SC, at least in its standard conception, requires that the initial configuration (one which coincides with a static configuration in GR) has a mass that cannot be larger that the maximum mass allowed by the unscalarized neutron star which maybe lower than 2M ⊙ if the EOS considered is not very stiff [32,33]. In order to support this conclusion, we perform a thorough numerical analysis by evolving numerically the full non-linear system of equations in STT in the Jordan frame (as opposed to the Einstein frame which is the most frequently used [25,35]), including the relativistic hydrodynamical sector that represents the star (modeled by a perfect fluid), under the assumption of spherical symmetry. As initial data we select a large sample of possible initial configurations of unscalarized neutron stars in hydrostatic equilibrium. By definition, the initial data corresponds to an initial data in GR where the scalar-field is null. Near the threshold of instability towards SC the initial configuration is perturbed by a Gaussian scalar-field profile, and then a numerical evolution is performed until reaching the final state of the system which corresponds to a static scalarized neutron star with gravitational mass lower than the initial one but with the same initial baryon mass. This analysis has several consequences, both theoretical and observational. From the theoretical point of view it seems impossible that a dynamical process like the one described here may lead by its own to a scalarized static star with a gravitational mass larger than the initial one since by definition, spontaneous scalarization corresponds to a static configuration with energy (i.e. gravitational mass) lower than the energy in GR while keeping the total baryon mass fixed (see Fig.1 in [36], and Fig. 6 in [24]). That is, a scalarized neutron star is the energetically preferred configuration above a critical energy-density. As we stressed above, scalarized configurations with gravitational masses larger than the maximum mass in GR were predicted in the past by constructing spherically symmetric configurations by solving directly the field equations of STT under the assumption of strict staticity [24,36]. Thus, the scalarized neutron star models with masses larger than the maximum mass in GR, with fixed EOS, cannot be the result from the purely dynamical transition of initial lower mass configurations. A fortiori such massive scalarized stars must be formed dynamically by a more complicated process that is worth exploring in the future. From the observational point of view, it is still unknown if (large or low mass) scalarized neutron stars really exist in nature. On one hand, detailed observations in binary pulsars put constraints on the possible values of the NMC between the scalarfield and the Ricci curvature. For certain class of STT, the weaker the NMC the lower the scalar charge, in which case, dynamical SC might not be a viable mechanism to explain the recently detected massive neutron stars without appealing to stiff or to more exotic EOS. On the other hand, the current observations of gravitational waves by the LIGO-VIRGO collaboration, notably, by the event GW170817 involving the collision of two neutron stars [37], together with the simultaneous detection of gamma-ray bursts [38] have allowed to rule out several ATG [39][40][41]. Notwithstanding, such event did not have enough precision to constraint by its own the existence of an additional polarization mode, like the breathing scalar mode alluded above. One can avoid, in principle, the constraints by the binary pulsars in STT by including a mass term for the scalar field [42] as in this case, the field becomes of shorter range depending on the magnitude of the mass. In the following sections we describe the formalism and the numerical results on the dynamic transition to SC using ab initio the Jordan frame where the physics is better understood. We conclude the paper by discussing several directions of study like considering a mass term and speculate also about possible scenarios leading to the formation of massive scalarized neutron stars other that the dynamical transition explored here.
Scalar Tensor Theories and Spontaneous Scalarization. We consider the action for STT in the Jordan frame given by where S matt represents the action for the matter part that in the present case corresponds to a perfect fluid, f (φ) is a NMC function f (φ) = 1 κ (1 + κξφ 2 ) , with κ = 8πG 0 , G 0 is Newton's gravitational constant, ξ a positive dimensionless constant, and V (φ) is a scalar potential. Variation of the action (1) with respect to the metric and with respect to the scalar field provide, respectively, the following two field equations: where a prime indicates derivatives with respect to the scalar field. The effective energy-momentum tensor (EMT) is given by the contribution of three parts: µ stands for the total mass density of matter in the rest frame of the fluid, and p is the pressure as measured in the same frame. The diffeomorphism invariance of STT leads to the conservation of the EMT of matter alone: ∇ a T ab fluid = 0, which in turn provides the hydrodynamic equations for the fluid. For simplicity we take V (φ) ≡ 0 (i.e., a massless and non-self interacting scalar field) like in the original SC scenario [22,24]. Notwithstanding, the massive case has been analyzed recently in static situations as a mechanism to avoid the constraints imposed by pulsars [42]. In the future we plan to study the dynamical transition to SC by taking into account a mass term as well. This paper focuses on solving the field equations of STT under the assumption of spherical symmetry. The details concerning the system of equations under this symmetry will be reported elsewhere. The interested reader is urged to consult [43] in order to have a glimpse of this system, but suffice is to say that in this analysis we use radial coordinates of area-type as opposed to isotropic. In [43] we considered a complex-valued boson field (i.e. a boson star) as matter instead of a perfect fluid. As we mentioned before, the dynamical transition to SC have been analyzed by several authors in the past, but in the Einstein frame [25,35] where the field is coupled minimally to the curvature but non-minimally to the matter sector.
Stellar models. We consider a neutron-star model that initially is in hydrostatic equilibrium. The conservation of the EMT leads to an equation which is formally identical to the Tolman-Oppenheimer-Volkoff (TOV) equation, except that the metric components take a different form due to the contribution of the scalar-field. For simplicity and in order to avoid interpolation of data of realistic EOS, we model the internal structure by a simple polytropic EOS parametrized by the baryon density n b : where m b = 1.66 × 10 −24 g and n 0 = 0.1 fm −3 . The values of the parameters γ andκ are adjusted to fit a stiff and a soft EOS. In particular, we use γ = (2.34, 2.46), and κ = (0.0195, 0.0093) like in Ref. [36]. Given this explicit EOS, we solved the TOV equation and constructed an initial static configuration where the star is not yet scalarized, and thus, the scalar field is absent. Thus, initially all the equations reduce to the field equations in GR. However, as the scalar-field evolves, automatically the field equations of STT are taken fully into account. In order to trigger easily the dynamical transition to SC the initial static configuration must be near the threshold of instability towards SC and then add a small Gaussian perturbation to the scalar field. If we departed from a static configuration not very near that threshold, while keeping the perturbation small, the latter is simply radiated away and the star is not scalarized. Under those initial conditions the star evolves towards the scalarized state, and the scalar field grows at the center of the star until the star stabilizes again into a different "hydrostatic" equilibrium configuration with a lower central energy density. Part of the total energy is radiated away in the form of scalar gravitational waves, and during this process the star develops a new global quantity, dubbed the scalar charge ω, defined asymptotically in terms of a surface integral as follows: The factor 1/4π √ G 0 as well as the minus sign are a matter of convention, and ω has mass units. Since at the end of the transition the scalar field behaves asymptotically as φ ∼ √ G 0 ω/r, we can simply extract ω from this expression in the asymptotic region.
We choose ξ ≥ 15 in order to enhance the appearing of SC [44]. However, according to [23,45] the curvature parameter β 0 = −2ξ associated with STT is such that −5 β 0 in order to satisfy the constraints imposed by several pulsars. This constraint translates into ξ 2.5. Intriguingly, for sufficiently high ξ the star properties become universal (cf. Figure 1). For values 0 < ξ 15 with the polytropic EOS described before, SC is basically unobserved. That is, the resulting star configurations are practically the same as in GR.
Maximal masses. The mass of a neutron star is computed using the definition of the Misner-Sharp mass function, where r out is the radial coordinate at the outer boundary of the numerical domain. This mass is basically the ADM mass M ADM minus the energy radiated away in the form of scalar gravitational waves. At the end of the evolution M MS < M ADM , but initially, M t=0 MS = M ADM . We checked that during the evolution the total baryon mass M bar remains the same, corroborating that the preferred static configuration with fixed M bar is the one with lower gravitational mass M MS , here t f stands for the end of the transition to SC.
Representative sequences of equilibrium models of neutron stars are plotted in Figure 1 which shows different mass profiles as a function of µ 0 , the density at the center of the star. The location of the maximum mass indicates the critical point separating the stable and the unstable branches towards gravitational collapse. In GR, static configurations located at the left of that point are stable, while those on the right are unstable. However, in the current analysis there also exists a point of instability (i.e. a critical density) towards SC. This critical density is marked with a green box in Fig. 1. The configurations between that density and the density leading to the maximum mass model, scalarize spontaneously and end up in one of the configurations associated with the colored lines below the black one. Figure 2 shows some examples of the metric component g rr (r) (upper panel) and the lapse function α(r) (lower panel) prior scalarization (dashed lines) and after scalarization (solid lines) when ξ = 25. Figure 3 (lower panel) depicts the corresponding densities prior (dashed line) and after scalarization (solid line). The solid lines of Figures  2 and 3 correspond to a configuration associated to the point (marked with a * ) in Figure 1. The upper panel of Figure  3 shows the scalar field profile after the star becomes scalarized. Before the spontaneous scalarization initiates the scalar field is null. Figure 4 depicts the product P := −g rr × g tt , where g tt = −α 2 , as a function of the area coordinate r. As we stressed, prior scalarization the scalar field is absent, thus, outside the star there is strict vacuum and the solution there is given by the Schwarzschild solution for which the product P ≡ 1 (solid line). However, at the end of scalarization the scalar field is present outside the star, and the solution is not longer the Schwarzschild solution. This is reflected by the fact that the product P is not unity (dashed line), but only approaches that value asymptotically as the scalar field vanishes.

Dynamical evolution.
We performed non-linear numerical evolutions in spherical symmetry of the STT-perfect-fluid system using as initial data a static TOV configuration with the OllinSphere code, a numerical relativity finite-difference code for spherical symmetry.
The code implements a variant of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) evolution equations for the geometry [46,47], which includes the contribution of the NMC scalar field, and for the fluid the Valencia formulation for relativistic hydrodynamics is implemented [48]. Explicit details about our numerical implementation have been reported for instance in [43,49]. We also note that the convergence properties, gauge choices, and boundary conditions of our numerical code have been extensively tested before in various physical systems [50,51].
Under the assumption of spherical symmetry we adopt a conformal decomposition for the 3-metric as follows, where dΩ 2 = sin 2 θdϕ 2 + dθ 2 is the solid angle element. We adopt a BSSN formulation of the STT field equations as written in the form of Einstein's equations (under spherical symmetry) with an effective EMT given by T ab as described below Eq.(1), which is composed by the three contributions Eqs.
(2)-(4). Thus, the following quantities are evolved in time from the initial data: the metric functions a(t, r) and b(t, r), the conformal factor ψ, the trace of the extrinsic curvature K, the traceless part of the conformal extrinsic curvature and the radial component of the conformal connection functions. However, during the evolution ψ and b remain constant to their initial values ψ = 1 and b = 1.
The matter fields include the contribution of both the fluid and the scalar field. The dynamical equations for the fluid are the relativistic Euler equation for the velocity field and for the total energy density measured by Eulerian observers, both resulting from ∇ a T ab fluid = 0. These equations can be written in conservative form following the the Valencia formulation [48], which are complemented by the equation for the conservation of the baryon density (i.e. rest mass density). For the scalar field φ we transform the second order "Klein Gordon" equation included below Eq.(1) into a system of two first-order evolution equations [52]. At each time step we compute the matter sources for the effective Einstein equations. In addition to the BSSN spacetime variables and matter content, there are two more variables left undetermined, the lapse function α, and the shift vector β r . For our simulations we choose for simplicity a vanishing shift throughout the evolution, while the lapse function is evolved from the initial data using the standard 1+log slicing condition: ∂ t α = −2αK. We perform a free evolution, and so we monitor at each time step the Hamiltonian and the momentum constraints to check the accuracy of the numerical solution, without solving the constraints at every time step, but only initially. A detailed description of the specific equations in spherical symmetry and their numerical implementation will be provided in a forthcoming investigation where a thorough study about this subject matter will be reported, which include the collapse of a neutron star into a black hole and the radiation in the form of scalar gravitational waves associated with the NMC field φ [53]. Figure 5 depicts the evolution of the central value of the density of the star (lower panel) and the central value of the scalar field (upper panel) during the scalarzation process. Notice that as the scalarization ends, both values provide the central values of the scalarized profiles of Figure 3. Several physical quantities evolve during the scalarization process but for brevity we report only these two. In a more detailed report we plan to provide a more complete set of plots showing the evolution of other interesting variables. Remarkably, the scalarized star is less compact than the initial configuration as one can appreciate from the value of the lapse at r = 0, which is larger after the scalarization finishes. This is also reflected in the height of the peak in grr. The black dot indicates the values at the surface of the star (i.e. where the pressure and density vanish). The radial coordinate is given in units of rs = G0M⊙/c 2 .

Discussion and Outlook
The phenomenon of spontaneous scalarization in STT allows to describe neutron stars with maximum masses that are larger than those in GR [22,24]. Beyond some critical baryon mass configurations with a non vanishing scalar field are energetically more favorable than the corresponding configurations at the same baryon mass with zero scalar field. However, scalarized neutron stars with masses larger than the maximum mass models in absence of the scalar field (i.e. those in GR) cannot be produced dynamically while keeping the baryonic mass of the star fixed since a dynamical scalarization process produces configurations with energy (i.e. gravitational mass) lower that the non scalarized counter part. Thus, the stationary scalarized neutron stars with masses larger that the maximum mass models in GR described in the past (for a given equation of state) [22,24], if existing in nature, must have formed following a dynamical process different from the pure spontaneous scalarization as described in this letter. Such process may include a cosmological origin or accretion of a surrounding scalar field. No doubt that the formation of scalarized neutron stars with masses larger than the maximum mass models of GR deserves further investigation as such supermassive scalarized neutron stars may explain the existence of the recently observed supermassive  Unlike the initial (GR) configuration where this product is unity outside the star due to the absence of the scalar field (dashed line), for the scalarized solution (solid line) this product is not unity outside the star, illustrating that the spacetime is not given by the Schwarzschild solution there due to the contribution of the scalar field (the Birkhoff theorem does not apply in this case). Notwithstanding, in the asymptotic region the product does tend to unity.
pulsars with masses larger that 2M ⊙ without necessarily appealing to exotic or very stiff equations of state for the nuclear matter at high densities. While understanding this issue is a matter of principle, from the observational point of view it may well happen that such supermassive scalarized stars are already ruled out by observations in binary systems for requiring a large value for the NMC constant ξ. Such constraints emerge, however, when the scalar field is massless, and thus of large range, so the inclusion of a mass term for φ may help to avoid those constraints [42] while possibly leading to super- "tov_rho0.0028k1527/A.rl" "tov_rho0.0028k1527/fluid_rho0.rl" u 1:($2*1000) "superficie.txt" "tov_rho0.0028k1527/psi.rl" "tov_rho0.0028k1527/B.rl" "tov_rho0.0028k1527/alpha.rl" u 1:($2*$2)