On the sound velocity bound in neutron stars

It has been suggested in the literature that the sound velocity of the nuclear matter $v_s$ violates the so-called sound velocity bound $v_s \le c/\sqrt{3}$ at high density, where $c$ is the speed of light. In this paper, we revisit this issue and confront the current measurements of mass, radius, and tidal deformability of neutron stars with $10^5$ different equations of state which are parametrized at low density and saturates the sound velocity bound beyond twice the saturation density where the equation of state has not been constrained yet, by which we can conservatively obtain the maximum mass of the neutron stars compatible both with the observed properties of neutron stars and the sound velocity bound. We find that majority of the models are eliminated by the incompatibility with the observations and, especially, the recently detected massive pulsar ($2.35\pm 0.17 M_\odot$) is hardly realized by our simulations. Our study strongly supports the violation of the sound velocity bound.


Introduction
Even though the quantum chromodynamics (QCD) has long been established as a fundamental theory of strong interaction, properties of nuclear matter beyond the saturation density are poorly understood both theoretically and experimentally.Since neutron stars (NSs) are very compact stellar objects consisting of nuclear matter and the matter density inside can exceed the saturation density, observed properties of NSs such as mass, radius, and tidal deformability enable us to constrain the equation of state (EoS) of nuclear matter at such large density.At zero temperature, which is a good approximation for NSs, EoS refers to the pressure as a function of the energy density: P = P ( ).
Observations of massive NSs exceeding 2 M [1] suggest that pressure rapidly increases as the energy density increases to support the star against gravity.The slope of the EoS is related to the sound velocity of the nuclear matter as v 2 s = c 2 ∂P ∂ .In this respect, there has been a growing interest in whether the sound velocity v 2 s exceeds the so-called conformal limit c 2 3 at some density or not (i.e., [2]).Obviously, the bound v 2 s ≤ 1 3 c 2 is satisfied for systems consisting of non-interacting particles or interacting non-relativistic particles.At ultra-density where interactions between ultra-relativistic quarks become negligible due to asymptotic freedom, v s → c √ 3 and the bound is also satisfied there.It is known that the bound is satisfied for a wide class of field theories (e.g., [3]).Meanwhile, recent studies show that the bound can be largely violated in theories close to QCD [4][5][6][7] although it is still not theoretically established that the bound is violated in QCD.Thus, there is a strong motivation to investigate if the bound is satisfied or violated by studying NSs.
By making comparison between various possible EoSs and observations of massive pulsars and measurements of tidal deformability of NSs by gravitational-wave experiments, recent papers strongly indicate that the bound is actually violated and the sound velocity takes a local maximum (larger than 1 √ 3 c) at intermediate nuclear density [8][9][10][11][12][13].This has been done by adopting the known EoS at low density and connecting it to higher density regime where EoS is given by piecewise form (linear in the chemical potential or in the energy density) or some specific functional form characterized with variable parameters.In particular, the recent study [13] incorporating the heaviest NS (M = 2.35 ± 0.17 M ) recently discovered as a pulsar shows that the bound is violated even if the lowest value in the allowed mass range (i.e.M = 2.18 M ) is adopted.
Although the existing studies already provide a strong support in favor of the violation of the bound, we would like to address this issue by adopting the type of EoS used in [8].This EoS matches the EoS of nuclear matter below twice of saturation density where the EoS is thought to be well understood.Above the critical density, the EoS sharply turns into stiff one given by v s = 1 √ 3 c at any density.Obviously, this EoS is artificial and unrealistic.Yet, the benefit from using it is that we can obtain the robust maximum mass of the NS compatible with the bound.In other words, observation of a NS heavier than this maximum mass provides an absolute evidence that the sound velocity exceeds 1 √ 3 c irrespective of the uncertainties of EoS at large density.Notice that our approach is the same philosophy as the one used to obtain the robust maximum mass of NS consistent with causality v s ≤ c [14] by adopting the EoS which becomes the stiffest one given by v s = c above a certain density .
In the next section, we first define the EoS used in this paper based on some relevant papers.We then briefly review how to compute the tidal deformability and the effect of rotation on NS masses.In Sec. 3, we present our results, followed by discussions about the implications and comparison with previous studies.Final section is devoted to the conclusion.

β-equilibrated matter
We use the Skyrme based parametrization [15] refined by Hebeler et.al [16] under which the energy density including the major contributions from rest masses(electron mass neglected) can be written as where x = n p /n represents the proton fraction present, n 0 = 0.16 ± 0.01 fm −3 , T 0 = 2 /2m.The set of parameters {α, η, α L , η L , γ}, control how the EoS will act and vary as a function of n = n/n 0 .The pressure corresponding to Eq. (2.1) is given by p = n 2 ∂( /n) ∂n .Due to the symmetric nature of matter at saturation, we have Also one may define the quantity B such that, which is also computed at the saturation density as 16 ± 0.1 MeV [17].Additionally the EoS can be described in terms of some properties : (a) The nuclear incompressibility coefficient K which measures the stiffness of symmetric nuclear matter can be described as K(n) = 9 ∂p ∂n .At saturation density in the symmetric phase, it becomes From study of giant resonances [18,19], K 0 has been determined to be 235 The symmetry energy S(n) is defined by Additionally we may also define the slope of this symmetric energy curve at saturation density by Using experimental data from heavy ion collisions, giant dipole resonances, and dipole polarizabilities [20,21], the values of S and L were taken to be 32 ± 2 MeV and 50 ± 15 MeV.
The set of these previous 5 equations can be used to solve for the values of the parameters {α, η, α L , η L , γ} characterizing the EoS.Furthermore, the proton fraction x can be computed by the condition for the β-equillibrium of matter [22]: Now, the EoS obtained from Eq. (2.1) may be correct at lower densities, but fails at higher densities.As we have discussed in the Introduction, our aim is to determine the maximum mass consistent not only with the EoS compatible with experimental data at low density but also with the bound v 2 s ≤ c 2 /3 in a conservative manner.Thus we adopt the following form of EoS similar to the one used in [8]: where (n, x) is obtained after using x determined from the minimization condition Eq. (2.7).The latter term in the region above twice the saturation density is added to preserve the continuity of the EoS [23].
In order to incorporate the experimental uncertainties of n 0 , B, K, S, L in the determination of the maximum mass of the neutron star, we randomly generate the values of these five parameters by treating them as probabilistic variables obeying the normal distribution with each mean and variance given by the central value and the error of experiments.We exclude any sample that are outside the 2σ range of the normal distribution in order not to include the outliers which inevitably appear when the sample size becomes large.In our analysis in the next section, we generate a dataset of 10 5 different EoS which are further processed according to the observational constraints as we will discuss in 3.2.

Mass, radius, and tidal deformability of static neutron stars
The metric sourced by static and spherically symmetric neutron star can be written in the form ds 2 = −e ν(r) (cdt) 2 + e Φ(r) dr 2 + r 2 (dθ 2 + sin 2 θdφ 2 ). (2.9) Outside the star, ν(r) and Φ(r) are given by where M is the mass of the neutron star.Solving Einstein's field equations and using the energy-momentum tensor for this metric, we obtain the Tolman-Oppenheimer-Volkoff (TOV) equation [24,25] where the mass M (r) contained in the sphere of radius r is related to the energy density as (2.12) By solving the TOV equation starting at r = 0 under a given central pressure up to the radius R where p drops to zero, we can obtain the mass and the radius of the neutron star.
The tidal deformability is defined as the ratio of mass quadrupole moment of a star Q ij to the external tidal field ε ij as [26,27] where k 2 is the Tidal Love number whose analytical formula is given as follows [27]: where C = GM (R)/c 2 R and y = Rβ(R)/H(R).The latter two unknown quantities can be calculated by solving the set of equations [28]: where X(r) and Y (r) are represented by ) (2.17) The dimensionless tidal deformability Λ is written as (2/3)k 2 [(c 2 /G)(R/M )] 5 .Using the above one can systematically find the mass, radius and Λ of the neutron star.

Formalism to compute the mass of the rotating neutron star
The rotation deforms the neutron star: it flattens the star to some extent, depending upon its angular velocity Ω. Formalism to compute the structure of the slowly rotating neutron The metric of a slowly rotating, stationary, and axially symmetric system may be defined as [31]: where H, Q, K, L are functions of r and θ alone.L can be recognised as the angular velocity (dφ/dt) acquired by an observer who falls from infinity to the point (r, θ).The Einstein equations are then solved under this metric to derive various properties.Components of the Einstein equations at zeroth-order in Ω gives the TOV equation (Eq.(2.11)).We assume that the star's rotation is uniform, i.e., the angular velocity does not depend on the position.At first order in Ω, the Einstein equations boil down to the following differential equation for ω(r, θ) = Ω − ω(r, θ) which is the relative angular velocity of fluid with respect to the local inertial frame as: where j = e −(ν+Φ)/2 .By solving this equation, we can determine the metric up to first order in Ω.
In order to compute the mass increase due to rotation, we need to consider components of the Einstein equations at second order in Ω.The mass increase δM can be decomposed as [29] where J represents the angular momentum of star given by c 2 R 4 (dω/dr)/6G and m 0 (R) is obtained by solving the following coupled differential equations: ) (2.22) Figure 2: The logarithm of the ratio of mass of rotating neutron stars to that of non-rotating one as a function of the rotation velocity Ω. Uniform rotation is assumed.

Mass-Radius plot and computing the Maximum Mass
For our analysis here, we would need the maximum mass arising out of a single parameter set {n 0 , B, K, S, L}.For this, we vary the initial starting pressure for a range of values, and see where the maximum appears.The effect of variation of initial central pressure on the mass of neutron stars has been shown in Fig. 1a, from which the maximum mass is extracted upto the required accuracy.The Mass-radius plot is also depicted for the same set of parameters in Fig. 1b.
3 Results and Discussions

Effect of rotation
Firstly let us investigate the quantitative effect of rotation.To this end, we numerically solved the set of coupled equations (2.19)-(2.22) to obtain the ratio of mass of rotating neutron star to the one of non-rotating neutron star for various values of Ω.The result is plotted in Fig. 2, in which at a single rotation value, we have composed 50 different random neutron stars simulations all markedly scattered very close to one point.best fit line that passes through the mean of all these values has a slope of 2.007, which verifies the fact that δM is proportional to Ω 2 [29].In our calculations, we restricted ourselves to a maximum rotation frequency of 10 3 Hz given the observations of 3 ms and 40 ms rotating massive neutron stars [32,33].From Fig. 2, we find that the relative mass increase is still O(0.1)% even at the maximum rotation velocity we consider.Since this level of tiny mass variation is irrelevant to our study, we neglect rotation in the following analysis.1.892 M .We also notice that there are non-negligible fraction of EoS that achieve the maximum mass exceeding 2 M .This may be ascribed to our choice of EoS (2.1) whose pressure sharply increases beyond 2n 0 to realize the highest sound speed v s = 1/ √ 3c.We mention that our histogram does not exactly coincide with the one given in [8] which also derived the histogram based on the same procedure as ours.The origin of the discrepancy remains unsettled.We then proceed to place the following observational constraints on the mass-radius curve for each EoS: i) NICER experiment which recently measured mass and radius of the two pulsars PSR J0030+0451 and PSR J0740+6620 as 1.34 +0.15 −0.16 M and 2.072 +0.067 −0.066 M for the mass respectively and 12.71 +1.14 −1.19 km and 12.39 +1.30 −0.98 km for the radius respectively [4,34] #1 , ii) tidal deformability Λ constrained from the observation of the neutron star merger event GW170817 by the LIGO-Virgo Collaboration [35].We impose Λ < 1400 for the neutron star of mass 1.4 M (picked from the M-R diagram).As an illustration, Fig. 4 shows the mass-radius relation of an EoS that passes the first condition i) and subsequently also the second ii).
#1 Although there are some other pulsars, such as PSR J1614-2230 [33], for which the mass ranges are known, we limit ourselves to the NICER data because it also contains the radius ranges, allowing for tighter constraints.Fig. 5 is the histogram of the maximum mass of the neutron star which has been obtained by applying the two observational constraints to Fig. 3. First of all, we find that many models have been eliminated after the selection is applied, only 7.05 % models remained.Intriguingly, in the histogram, the maximum mass is at most 2.178 M which is very close to the lower end of the mass range of the recently detected pulsar (M = 2.35 ± 0.17 M ) [13].Even if it is just one such set amongst 10 5 parameter sets, we emphasize that the sound velocity bound is not ruled out robustly because our sample is limited.Even after we begin with 10 5 different EoS, we end up with a few EoS that achieve M max 2.1 M after imposing the observational conditions.In such a situation, different runs will give different tail in the histogram at large M max .Unfortunately, the numerical computations with larger number of sample are not possible within reasonable time with our computers.Thus, we leave the question as to whether the sound velocity bound is violated in nature an open issue.Nevertheless, since the adopted EoS shows the sudden jump of its slope across 2n 0 , which is highly artificial, more realistic EoS satisfying the sound velocity bound should be smoother.This implies that the histogram based on such EoS shifts toward smaller M max , making the sound velocity bound more incompatible with observations.Thus, our analysis strong supports the idea that the sound velocity bound is violated inside the heavy neutron stars observed in nature.

Conclusion
Observations of the massive neutron stars indicate that the pressure of nuclear matter rises rapidly as the energy density is increased at extremely large density.In particular, it has been suggested in the literature that the sound velocity bound v s ≤ 1/ √ 3c is violated at such high density.In this paper, we have studied to what extent the current measurements of mass, radius, and tidal deformability of neutron stars reinforces the violation of the sound velocity bound.In order to minimize the impact of the uncertainties of the equation of state of nuclear matter and derive a conservative consequence, we employed the equation of state that is parameterized by a set of parameters below twice the saturation density and saturates the sound velocity bound above that density and obtained the mass-radius relation for 10 5 different sets of parameters.We found majority of the models are eliminated by the incompatibility with the observations.Especially, the recently detected massive pulsar (2.35 ± 0.17M ) is hardly realized by our simulations.Our study strongly supports the violation of the sound velocity bound.

Figure 1 :
Figure 1: (a) The Mass-Initial Central Pressure plots for calculating maximum mass and (b) The Mass-Radius plot of a Neutron Star

Figure 3 :
Figure 3: The maximum mass histogram without any constraints

Figure 4 :
Figure 4: The successful case of passing through the regions of uncertainty of both the neutron stars PSR J0030+0451 and PSR J0740+6620

Figure 5 :
Figure 5: The final histogram after applying the constraints from observed neutron masses and tidal deformability