Spontaneous scalarization with an extremely massive field
and heavy neutron stars
Abstract
We investigate the internal structure and the massradius relation of neutron stars in a recently proposed scalartensor theory dubbed asymmetron in which a massive scalar field undergoes spontaneous scalarization inside neutron stars. We focus on the case where the Compton wavelength is shorter than , which has not been investigated in the literature. By solving the modified Einstein equations, either purely numerically or by partially using a semianalytic method, we find that not only the weakening of gravity by spontaneous scalarization but also the scalar force affect the internal structure significantly in the massive case. We also find that the maximum mass of neutron stars is larger for certain parameter sets than that in general relativity and reaches even if the effect of strange hadrons is taken into account. There is even a range of parameters where the maximum mass of neutron stars largely exceeds the threshold that violates the causality bound in general relativity.
I Introduction
General relativity (GR) has been tested by various experiments Will (2014). Especially, the precise tests in the Solar system put stringent constraints on the deviation from GR. In fact, GR is a nonrenormalizable theory in quantum field theory and cannot be applied beyond the Planck scale. In this sense, it should be regarded as an effective field theory of an asyetunknown quantum gravity theory. However, even at the classical level of the effective field theory, the correct theory may not be GR. Actually, the existence of dark matter and dark energy may suggest the modification of GR in the extreme regime not probed by the terrestrial or solarsystem experiments, and various theories which pass the solarsystem experiments and accommodate the dark sector have been proposed Joyce et al. (2015). In order to test these possibilities, it is important to test gravity beyond the currently accessible scale.
Recently, gravitational waves emitted from a binary black hole coalescence were detected for the first time Abbott et al. (2016a) and two signals as well as a less significant candidate followed it Abbott et al. (2016b, 2017). Needless to say, gravitational waves will become a new tool to test gravity in the near future. In fact, the detected signals put constraints on the gravity theory in the highly dynamical and strongfield regime for the first time Abbott et al. (2016c, b); Yunes et al. (2016). Although these new constraints are weaker than the existing ones for most of the theories, more precise tests will be possible as the detector sensitivity will be improved.
One of the main targets of the gravitationalwave observations, which is also relevant to our present study, is neutron stars. A neutron star is a highly dense star, whose density reaches the nuclear density. There is not enough evidence to show that GR is correct in such a highly dense region. On the other hand, the structure of neutron stars in various theories of gravity has been investigated by many authors (See Berti et al. (2015) and references therein. Babichev et al. (2016); Sakstein et al. (2017); Minamitsuji and Silva (2016); Maselli et al. (2016); Aoki et al. (2016); Cisterna et al. (2015, 2016) are recent work on this subject not included in Ref. Berti et al. (2015).). Since the structure varies from that in GR, it is possible to distinguish these theories by observations of neutron stars in principle. However, things are not so trivial because of the degeneracy between the uncertainties in the equation of state for nuclear matter and gravitational physics. One way to break this degeneracy is to use equationofstateindependent relations Postnikov et al. (2010); Hinderer et al. (2010); Yagi and Yunes (2013). The ILoveQ relation Yagi and Yunes (2013) is one such relation, which holds between the moment of inertia, I, the Love number and the quadrupole moment, Q. The Love number can be measured by gravitationalwave observations in the near future Flanagan and Hinderer (2008). Therefore, tests of gravity with neutron stars will become feasible in the coming era.
Another topic related to neutron stars is the existence of massive neutron stars, whose masses are about 2 Demorest et al. (2010); Antoniadis et al. (2013). Especially, PSR J16142230 Demorest et al. (2010) shows a strong Shapiro delay signature, and the estimated mass of the pulsar turned out to be . This measurement does not assume any models of emission mechanism of the stars and the result is robust. On the other hand, although the state of nuclear matter inside the highly dense core of the neutron stars is unknown, it is natural that strange hadrons appear there since the chemical potential of the neutrons is large enough. For example, it is pointed out that hyperons appear when the baryon number density, , surpasses the threshold density of Nishizaki et al. (2002), where is the saturation density. Such strange hadrons soften the equation of state and significantly reduces the maximum mass of the neutron stars that can be supported against gravity. A 2 neutron star is difficult to realize with these strange hadrons (see Fig. 3 of Demorest et al. (2010)). In order to solve this potential inconsistency, many ways to improve the equations of state for nuclear matter have been suggested, such as taking into account threebaryon interactions Nishizaki et al. (2002); Gandolfi et al. (2015); RikovskaStone et al. (2007) or a quark matter core Bonanno and Sedrakian (2012); Masuda et al. (2013). On the other hand, it is also possible that modification of GR solves this inconsistency Astashenok et al. (2014).
One simple extension of GR is scalartensor theory Fujii and Maeda (2007). In this theory, a new scalar degree of freedom is added to GR. This additional degree mediates a new force called the scalar force and it affects the motion of particles. For example, it modifies the trajectory of light around the Sun from that in GR. Currently, the measurement of the time delay of light with the Cassini spacecraft Bertotti et al. (2003) puts stringent constraints on the coupling constant between the scalar field and matter, and the scalar force is negligible in the Solar System.
Damour and EspositoFarese proposed an interesting scalartensor theory which passes the constraints and modifies the structure of neutron stars significantly Damour and EspositoFarese (1993). In this theory the value of the scalar field in the Solar System is so small that the coupling constant, which is zero when the scalar field is vanishing, is negligibly small. Therefore, this theory safely passes Solar System experiments. On the other hand, inside neutron stars, the scalar field gets a much larger value than in the Solar System and the deviation from GR is significant. This phenomena is called spontaneous scalarization Damour and EspositoFarese (1996). Especially, the effective gravitational constant decreases and gravity is weakened in the spontaneousscalarization phase. Therefore, this theory may allow more massive neutron stars than in GR and solve the aforementioned problem of the observed existence of massive neutron stars.
The model proposed by Damour and EspositoFarese has two problems. The first is that this model is already severely constrained by binary pulsar observations Antoniadis et al. (2013); Damour and EspositoFarese (1996). Especially, from the observation of PSR J03480432 Antoniadis et al. (2013), it was found that the effect of spontaneous scalarization must be negligibly small even inside the massive neutron star of the binary, whose mass is . The second is that spontaneous scalarization can occur during the inflation and matter dominated eras. Therefore, the scalar field gets a large value in the present universe for a broad range of its initial value, which means finetuning for the initial value is necessary for this theory to pass solarsystem experiments Sampson et al. (2014); Damour and Nordtvedt (1993a, b).
These problems do not exist if the scalar field is massive Chen et al. (2015). First, if the Compton wavelength of the scalar field, , is much smaller than the periapse of the orbit of PSR J03480432, which is the order of Antoniadis et al. (2013), the scalar field does not affect the orbital motion of the binary and the test with this binary is safely passed Ramazanoğlu and Pretorius (2016). In addition, even if spontaneous scalarization occurs in the early universe, the scalar field begins damped oscillation when the Hubble parameter becomes equal to the mass and converges to zero without fine tuning Chen et al. (2015). In other words, GR is a cosmological attractor in the late time universe. Therefore, models of spontaneous scalarization with a massive scalar field have recently attracted attention and been investigated. In Ramazanoğlu and Pretorius (2016), the structure of neutron stars in the case where was investigated. On the other hand, the structure of neutron stars in the case of the shorter Compton wavelength has not been investigated. This region is interesting since the oscillating component of the scalar field may account for the overall amount of dark matter. This scenario is dubbed asymmetron scenario Chen et al. (2015).
In this paper, we investigate the structure of neutron stars in the model of spontaneous scalarization with a massive scalar field introduced in Chen et al. (2015), focusing on the case where . We also investigate the maximum mass of neutron stars and identify the parameter space where massive neutron stars exceeding are allowed with strange hadrons. The rest of the paper is organized as follows. In Sec. II, we explain our model and the basic equations. In Sec. III, we explain the numerical method we employ to obtain the structure of neutron stars. In Sec. IV, we show the results of our analysis, that is, the internal structure, the massradius relation, and the maximum mass of neutron stars in our model. The last section is devoted to the conclusion.
Ii Basic equations
We briefly review the model of the massive scalar field proposed in Chen et al. (2015) and the basic equations we use in our analysis.
ii.1 Model
As in Chen et al. (2015), we introduce a real massive scalar field in the gravitational sector, whose Compton wavelength , where is the mass of , is smaller than . We assume that this scalar field couples with ordinary matter fields universally through the physical metric , which is called the Jordan metric (compared to this metric, is called the Einstein metric.). This coupling guarantees that the weak equivalence principle holds true. We consider the action given by
(1)  
(2) 
where is the gravitational constant measured in laboratory experiments and is the Lagrangian of all the other matter fields. The equations of motion derived from this action are the following:
(3)  
(4) 
where is defined by . and are the stressenergy tensor and its trace defined by the Jordan metric respectively, that is,
(5) 
According to Eq. (3), the effective gravitational constant, , changes when the value of deviates from . The equation of motion for the scalar field can be rewritten as
(6) 
This shows that dynamics of is described by an effective potential , which depends on the density of surrounding matter. As is clear from the expression of , if the second derivative of at the origin is negative, the spontaneous scalarization happens when exceeds a critical value Chen et al. (2015).
In this study, the scalar force plays an important role. The strength of the scalar force is characterized by . As an example, we consider two point particles interacting with each other only through in vacuum. We assume that the background value of is , which is the stable point of the potential . In the Newtonian limit, the equation of motion for the particles is
(7) 
where
(8) 
and is the mass of the other particle. The first term in the righthand side of Eq. (7) corresponds to gravitational force and the second corresponds to the scalar force. It can be easily seen that the absolute value of controls the strength of the scalar force.
Following Chen et al. (2015), we use
(9) 
throughout our study. In our study we assume that the matter is perfect fluid and , where is energy density and is pressure in the Jordan frame. Therefore the effective potential is
(10)  
(11) 
It can be easily seen that the point becomes unstable and spontaneous symmetry breaking occurs if , where the critical density is defined by
(12) 
The shape of the effective potential is given in Fig. 1. The value of is assumed to be somewhat smaller than mass density of a nucleon, . In this case, spontaneous symmetry breaking occurs only inside neutron stars, as we explicitly demonstrate later. In the Solar System, where , vanishes everywhere. Accordingly, also vanishes and the general relativity is classically recovered. Therefore, our model safely passes the solarsystem tests. On the other hand, inside a neutron star, where , has stable points at
(13) 
In this region, stays near the nontrivial stable point, which is , and . Therefore, the effective gravitational constant, , becomes smaller than that in the low density region and gravitational force is weakened inside neutron stars.
ii.2 Modified Einstein equations for a static and spherical star
Throughout our study we neglect spins of neutron stars and consider a static and spherically symmetric configuration, that is,
(14) 
The field equations Eq. (3) and Eq. (4) lead to
(15)  
(16)  
(17)  
(18)  
(19) 
where denotes differentiation with respect to the radial coordinate . We have introduced a new variable in order to reduce the equations to the first order differential equations. Equation (17) is the hydrostatic equilibrium condition for nuclear matter. In addition to the standard gravitational force, the scalar force described by the second term also contributes to the equilibrium condition. As we will see later, this scalar force significantly changes the internal structure of neutron stars. The surface of the star, , is defined as the surface on which becomes , and is set to be outside the star.
Next, we discuss boundary conditions. Nonsingular solutions satisfy the following conditions at the origin:
(20) 
Outside the star, has an exponentially growing solution and an exponentially decaying solution. Since the former one is physically unacceptable, we impose another boundary condition at infinity as
(21) 
This condition is equivalent to choosing a suitable value of at the origin. For , becomes almost identical to the physical metric . Since we are interested in is much shorter than the typical orbital distance of binary pulsars, the mass of the neutron star estimated by the binary pulsars measurement, , is equal to .
Finally, we explain roughly how changes in . As explained in II. C of Chen et al. (2015), ignoring the curvature of the spacetime and performing the change of variables as
(22) 
we can reduce Eq. (18) and Eq. (19) to
(23) 
It is the same as the equation for a motion of a particle under the potential with time dependent friction, and it helps us to understand the profile of Khoury and Weltman (2004). In the core of the star, where , the particle, , rolls down from (the left plot of Fig. 2). The potential around can be approximated by
(24) 
with . For , is
(25) 
and of the order of . ^{*1}^{*1}*1For , Eq. (25) can not be applied. However can be written as and is of the order of even in this case. Therefore, the deviation from increases exponentially with the scale of . After the density decreases and becomes lower than , the potential U is convex around . According to the boundary condition equation (21), climbs this convex potential and reaches at infinity (the right plot of Fig 2).


ii.3 Equations of state
In addition to the equations from Eq. (15) to Eq. (19), we need an equation of state. We adopt an equation of state defined in the Jordan frame,
(26) 
since it is the one that is directly provided by the nuclear physics experiments and theory. In our study, we consider two types of equations of state, that for hadronic matter and for strange quark matter. The latter is the matter composed of deconfined up, down and strange quarks. Although whether stars consisting of strange quark matter exist or not has not been clarified, their formation path and properties have been discussed in the literature Itoh (1970); Witten (1984); Haensel et al. (1986); Alcock et al. (1986); Olinto (1987).
ii.3.1 Hadronic matter
To model equations of state for hadronic matter, we use piecewise polytropic equations, which are constructed by a few polytropic equations Read et al. (2009). A polytropic equation is in the form of
(27) 
where is the restmass density of the matter and is the adiabatic index. Integrating the first law of thermodynamics
(28) 
we obtain the relation in the form of Eq. (26) as
(29) 
where is an integration constant. Piecewise polytropic equations are constructed by connecting polytropic equations at each density contained in a set of dividing densities . For each interval the equation of state is
(30) 
and
(31) 
with . In the limit the matter becomes nonrelativistic, that is, . It leads to
(32) 
The continuity condition for and at leads to
(33)  
(34) 
Therefore, the equation of state is characterized by , and . In our study we use two piecewise polytropes, which approximate AP4 (which contains only as nuclear matter) Akmal et al. (1998) and GS1 (which contains kaons, which are strange hadrons, in addition to ) Glendenning and SchaffnerBielich (1999). The GS1 is one of the equations of state taking into account the effect of strange hadrons and the maximum mass for it in GR is , which does not reach the measured mass of the neutron star PSR J16142230. It is one of the goals of this work to clarify whether the maximum mass for the GS1 reaches this measured value in our model. Following Read et al. (2009), we use the piecewise polytrope approximating the SLy equation of state at low density and three polytropes at high density for both cases. The parameters for both of the equations of state at low density and at high density are separately listed in Tables 1 and 2.
AP4  2.830  3.445  3.348  

GS1  2.350  1.267  2.421 
ii.3.2 Strange quark matter
Iii Method to solve the equations
In this section, we explain the numerical methods we use in this study. We adopt different methods for (mildly massive case) and (very massive case).
iii.1 case (mildly massive case)
We numerically integrate the equations (15)  (19) outwards from the center using the RungeKuttaFehlberg method, imposing the initial conditions Eq. (20) and
(36) 
Because some terms in the equations are apparently singular at the origin and can not be handled with by numerical computations, we use the Taylor expansion of the solution around to start the computations at .
In order to satisfy the boundary condition Eq. (21), the value of the scalar field at the origin must be tuned to a certain value. We use shooting method to find such physical solutions. We look for the appropriate value for in the interval between and a sufficiently large value by means of bisection search. The detail of the bisection search is as follows.

If we observe grows exponentially to positive direction, we take smaller value of for the next trial.

If we observe becomes negative, we take larger value of for the next trial^{*2}^{*2}*2We assume that tracks the stable point of and is positive everywhere for correct .

After repeating this procedure many times, we obtain which shows convergent behavior at sufficiently large r. We stop the numerical integration if the integration reaches the point in the vacuum at which the following conditions are satisfied.
(37) The third condition guarantees that the contribution from the energy of the scalar field at to is small enough to be ignored. Around that point, the following WKB solution can be applied (The derivation of the WKB solution is summarized in Appendix A).
(38) We connect the numerical solution with the WKB solution and calculate and . If , the deviation from the correct configuration of arising from the incorrectness of the value of is small all over the region, . Therefore, we stop the bisection search and adopt the resultant configuration as a good approximate solution. We approximate the mass of the neutron star as . On the other hand, if , we take larger value of for the next trial if and vice versa.
iii.2 case (very massive case)
If , it is very difficult to solve the equations numerically all over the space because the required accuracy for the numerical integration is extremely high. For example, if the numerical error in arises at , it increases by a factor of
(39) 
up to , where is a coefficient. represents the difference between and , and for . Since , the numerical error arising near the origin increases by a factor of
(40) 
Therefore, the necessary precision is too high to achieve, and we have to develop an alternative method to obtain an approximate solution in the limit, , which we will explain below.
We can model the profile of as shown in Fig. 3. As pointed out in II. C of Chen et al. (2015), stays extremely close to up to a certain radius in this regime. Then, leaves and transitions toward somewhere inside the star with the scale of . In the transition region, the scalar force becomes significant and surpasses gravitational force, that is,
(41) 
This inequality can be easily confirmed as follows. According to Eq. (16), noticing is larger than the other terms in the right side because it is the integrated value of the energy density from the center, we can estimate by
(42) 
Since inside neutron stars, we have
(43) 
On the other hand, the scalar force can be estimated by
(44) 
since increases by with the scale of . Since we now consider extremely short Compton wavelength, we can assume in this region. Therefore, Eq. (41) holds. Notice that is typically in the present case and the effect of the scalar force becomes significant when . Since the scalar force compresses the star significantly, the pressure decreases significantly in this region. For later convenience, we define as the radial coordinate of the point where the scalar force becomes comparable to the scalar force, that is,
(45) 
where the subscript i means the value at . After the rolldown phase, climbs up the mountain of and reaches at infinity.
Next, we explain the detail of each region. In the region, , remains extremely close to . Therefore, we can use the approximation,
(46) 
In addition, we have
(47) 
since changes with much longer scale than . The justification of Eq. (46) and Eq. (47) is explained in Appendix B. Therefore, we can approximate Eq. (15), Eq. (16) and Eq. (17) as
(48)  
(49)  
(50) 
By integrating these equations numerically, we can obtain an approximate solution for this region. The important point is that we do not have to solve Eq. (18) and Eq. (19). Therefore, we have no numerical difficulties in this integration.
In the region where rolls down on the length scale , we can approximate Eq. (17) as
(51) 
since the scalar force becomes dominant. In addition, since
(52) 
the region where the scalar force is dominant is an extremely thin shell and we can treat as constants, that is,
(53) 
With Eq. (52) and Eq. (53), we can approximate Eq. (18) and Eq. (19) as
(54) 
Equation (51) and Eq. (54) can be analytically integrated once. We first consider the hadronic equations of state and explain the case of the strange quark matter after that. In this case, the integration from to leads to
(55)  
(56) 
The derivations of Eq. (54), Eq. (55) and Eq. (56) are explained in Appendix C in detail.
Since decreases significantly in this region, there are two cases to consider. If
(57) 
is satisfied, becomes while is rolling down. In this case, Eq. (56) leads to the following condition satisfied by the physical quantities at and on the surface:
(58) 
where the subscript f means the value on the surface in this case. On the other hand, if
(59) 
is satisfied, even after has rolled down. In this case, Eq. (55) and Eq. (56) lead to
(60)  
(61) 
where the subscript f means the value at the position where is satisfied and the scalar force becomes negligible.
For both cases, climbs the mountain afterwards, and in this phase can be approximated by the following WKB solution:
(62) 
for the case of Eq. (57) and for the case of Eq. (59). Substituting Eq. (62) into Eq. (58) and Eq. (61), we obtain the equations to determine . As a result, we find that can be obtained as a zero point of , which is defined as follows:
(63) 
where is determined by Eq. (60). For the MIT bag model, is defined as follows:
(64) 
Obtaining , we can obtain an approximate solution as follows. First, we integrate Eq. (48), Eq. (49) and Eq. (50) from to , where . Then if satisfies Eq. (57), we stop the integration since becomes suddenly at that point. and can be obtained as and respectively. On the other hand, if satisfies Eq. (59), we have to continue the numerical integration from . Since we can neglect the contribution of , we just have to integrate the Einstein equations of GR up to the point where . and can be obtained as and at that point respectively.
Iv Result
In this section, we explain the results of our analysis. We show the internal structure, massradius relationship and the maximum mass of neutron stars.
iv.1 The internal structure of neutron stars
The profiles of for four different wavelengths and two different central densities are shown in Fig. 4. For each case, we observe that as decreases the profile approaches what we obtain by means of the semianalytical method which we discuss in Sec. III.2. This demonstrates the validity of our semianalytical method. We also observe that grows in the core for the higher central density while is monotonic decreasing function of for the lower central density. It is because is not a monotonically increasing function of . For the hadronic equations of state and , can be written as
(65) 
Since typical hadronic equations have , the second term becomes dominant and turns to decrease as , or , increases. Such a behavior is specific to the hadronic equations of state since and is always positive for the MIT Bag model. For the central density of the right plot of Fig. 4, at . For such high central densities and short Compton wavelengths, sits extremely close to in the core and starts to grow near the spherical surface, , on which (See the green line of the right plot of Fig. 4). Strictly speaking, the approximation, , is not valid at , and this seems to imply that the semianalytic method explained in the previous section is no longer valid. However, this is not true, as we explain below. Since the behavior of the effective mass, , near can be written as
(66) 
becomes large close to for large . On the other hand, must return back to before becomes much larger than since otherwise rolls down quickly and the boundary condition is not satisfied. Therefore, the region where the approximation, , is invalid is so thin that its effect on the internal structure is negligibly small, which means our semianalytical method is still applicable in this case.


Next, we show the profiles of , and for different central densities in Fig. 5. For the lower central density, decreases more gradually than in GR because and the gravitational force is weakened. Near the surface, decreases drastically due to the compression by the scalar force. For the higher central density, the deviation of GR is extremely small in the core because . In the region where grows as increases the decrease of becomes more gradual because the scalar force pushes the matter toward the outside. This effect also alters the internal structure significantly and such an effect of the scalar force is incorporated in the second term of the denominator of the right side in Eq. (50) in the semianalytical approach. The gravitational force is weakened in the outer part of the star and the scalar force compresses the star near the