Neutron Star Mergers Chirp About Vacuum Energy
Csaba Csáki, Cem Eröncel, Jay Hubisz,
Gabriele Rigo, and John Terning
Department of Physics, LEPP, Cornell University, Ithaca, NY 14853, USA
Department of Physics, Syracuse University, Syracuse, NY 13244, USA
Department of Physics, University of California, Davis, CA 95616, USA
, , ,
,
Abstract
Observations of gravitational waves from neutron star mergers open up novel directions for exploring fundamental physics: they offer the first access to the structure of objects with a nonnegligible contribution from vacuum energy to their total mass. The presence of such vacuum energy in the inner cores of neutron stars occurs in new QCD phases at large densities, with the vacuum energy appearing in the equation of state for a new phase. This in turn leads to a change in the internal structure of neutron stars and influences their tidal deformabilities which are measurable in the chirp signals of merging neutron stars. By considering three commonly used neutron star models we show that for large chirp masses the effect of vacuum energy on the tidal deformabilities can be sizable. Measurements of this sort have the potential to provide a first test of the gravitational properties of vacuum energy independent from the acceleration of the Universe, and to determine the size of QCD contributions to the vacuum energy.
1 Introduction
The recent observations of gravitational waves (GW’s) from the merger of neutron stars (NS’s) by LIGO/Virgo [1] along with the corresponding electromagnetic observations of the resulting kilonova have reverberated across most areas of physics and astronomy. From the point of view of particle physics the most important consequence of GW170817 and future merger events is our new ability to directly examine the properties of the QCD matter forming the inner layers of NS’s, allowing us to use NS’s as laboratories for fundamental physics [2, 3, 4]. This might also open up new avenues to testing the gravitational properties of vacuum energy (VE) which may also get at the heart of some of the deepest puzzles in fundamental physics [5].
It has long been speculated that there may be a new phase of nuclear matter at the core of the NS’s [6]. If such a phase indeed exists it is expected to be accompanied by a jump in VE [7] of order (where is the usual QCD scale) making NS’s the only known objects where VE might make up a nonnegligible fraction of the total mass. Therefore studies of the interior structure of NS’s can also probe the gravitational properties of VE, possibly shedding light on some of the most interesting open questions in physics: it could provide verification of the equivalence principle for VE. This would be the first test independent of those obtained from the cosmic acceleration of the Universe. Acceleration of the Universe provides information on VE in the lowtemperature low density phase of the SM, while NS’s could probe a low temperature but very high density phase if it exists in their cores. This could allow us to isolate the QCD contribution to ordinary VE and probe its gravitational properties.
Alongside the exciting advent of the gravitational wave observation era shepherded in by LIGO/Virgo, the Neutron star Interior Composition ExploreR (NICER) mission will soon measure masses and radii of several millisecond pulsars [8]. These measurements as well as the chirps from the inspiral of merging neutron stars can provide information about the equation of state (EoS) of dense nuclear matter. The chirps in particular are sensitive to the tidal deformability of NS’s as they approach each other [9, 10, 11, 12, 13]. There has already been considerable work on constraining the EoS using the new LIGO/Virgo data [14, 15, 16, 17, 18].
There exists an extensive literature focused on trying to put bounds on the nuclear EoS at high densities from neutron star measurements (see for example [19, 20, 21, 22]). Some recent theoretical work has focused on modeling possible new phases at the cores of neutron stars by using quasiparticle quarks rather than neutrons to provide the simplest description of the microscopic physics [23, 24, 25, 26, 27, 28, 29]. Further work has been done using NS’s to constrain “beyond the Standard Model” physics [2, 3, 4].
In this paper we will assume that there is only Standard Model physics involved in the composition of neutron stars, and we will not try to model the microphysics of the putative new phase. Our main goal is to investigate the observable effects of the presence of VE on the GW signal as well as the mass versus radius curve of NS’s, possibly providing new experimental probes of VE. To achieve this we will parameterize the effect of the new phase with a jump in the ground state energy due to a QCD phase transition assumed at the core [5, 29, 30, 31]. This new phase would appear at a critical pressure of order , and is expected to also lead to a change of VE [7] of order . We will follow the conventional models for NS’s where the EoS is divided into 7 layers, but we modify the innermost layer to take the effect of VE at the core into account. We will then evaluate the tidal Love numbers for such models. We stress that the existence of the new phase requires a modification of the traditional assumptions used in NS simulations regarding energy density jumps at the boundary of a new phase [29, 30, 31]. We will explore the effect of a difference in energy densities of the two phases that includes a discontinuous, density independent term reflecting the absence of the low density QCD contribution to VE [5]. We will present several models of NS cores and estimate the effect of VE on tidal Love numbers. We find that VE can have a significant effect on NS merger waveforms with high chirp masses, so that such events serve as a probe of the physics of vacuum energy.
The paper is organized as follows. In Section 2 we present the models we use for nuclear matter in the interior of NS’s, along with a detailed discussion of the treatment of the phase transition at the boundary of the innermost layer. Section 3 contains the description of the tidal deformability of NS’s. The results of our simulations and the effects of VE on the NS observables are given in Section 4: we show the mass versus radius curves and the tidal deformabilities for three different wellstudied NS models and the effects of VE on those observables. Finally we conclude in Section 5.
2 Modeling High Density QCD
The main difference between our work and that of previous studies of tidal deformability of NS’s is that we will fully account for a phase transition to an exotic phase of QCD in the innermost core region of NS’s. Crucially, we take into account the Standard Model expectation that there is a constant shift , independent of baryon number density, in the ground state energy relative to the surrounding layers parametrizing the change in VE due to the phase transition. In the ordinary phase of QCD the nonperturbative condensates of quarks and gluons make contributions [7] of order to the VE. These contributions, along with those from other sectors of the SM, are canceled by the “bare” cosmological constant down to the observed cosmological constant of order :
(2.1) 
The origin of the mechanism leading to this cancelation remains unknown. In an exotic phase of QCD the QCD contributions to the VE will have order one modifications and hence the precise cancelation will no longer apply:
(2.2) 
where is the shift in the QCD vacuum energy due to the phase transition. Hence in the absence of a dynamical adjustment mechanism, Standard Model physics predicts a density independent shift in the energy of the exotic phase compared to the ordinary phase, which will serve as a new effective cosmological constant term for this phase. An estimate of the difference between the VE of the exotic phase and the ordinary vacuum is given by nuclear saturation density: [5]. Such a phase change is strongly suspected to occur at high chemical potential, with theoretical evidence arising from truncated diagrammatic expansions and other approximate methods [24, 25]. The phase change is in fact part of the standard picture of the QCD phase diagram. For many plausible descriptions of the matter in the outer portions of the star, nuclear saturation density is approached near the core of the densest NS’s, making it quite possible that the most massive NS’s contain cores with an exotic phase. In this section, we give a description of how one can model the QCD equation of state at various pressures, with particular attention paid to the phase transition that may occur in the innermost region.
2.1 Modeling the Outer Layers
The physics of neutron stars is an extremely rich field, and there are many details that go into modeling the different regions of NS’s. Such an analysis is well beyond the scope of this work, however there are methods for coarsegraining these complexities to obtain an approximate equation of state for nuclear matter up until the phase transition we are interested in. Such an approximation is sufficient for the purposes of making predictions for gravitational wave signals. The most common methodology for modeling the high density nuclear physics region outside the exotic phase core is to separate the neutron star into multiple layers, with each layer satisfying a nonrelativistic polytropic equation of state. The parameters of the polytrope are fixed either by matching conditions or by fitting results from more detailed studies.
We follow this established methodology and model the nuclear fluid and its corresponding EoS as a piecewise polytrope where the boundaries between each layer are set by a given value of the pressure. Following previous work [21, 22] we will parametrize the EoS with a total of 7 layers. The Israel junction conditions [32] require that the pressure must always be continuous between layers, even if each side of the boundary is separated by a first order phase transition. It is traditional to parameterize the EoS by assuming that the pressure is given by a power of the mass density rather than a power of the energy density (as would be natural for a highdensity, relativistic fluid). Since we want to efficiently compare our results with the existing state of the art simulations (some of which have been used as benchmarks for the LIGO/Virgo analysis) we will bow to this tradition and parametrize the EoS as
(2.3) 
where for and for . The pressures, , dividing the various layers have a one to one correspondence with the boundaries in the mass density: . The Einstein equations contain the energy density, which is related to the mass density via the first law of thermodynamics: . Integrating the first law together with (2.3) yields the corresponding energy density:
(2.4) 
where the are integration constants. Note that the appearance of the parameters is a consequence of using a polytropic ansatz for the mass density. Naively, one would think that using a relativistic polytropic ansatz for the energy density would have led us to a relation with one less free parameter. However another thermodynamical condition, continuity of the chemical potential, would have forced us to reintroduce the baryon number density, and therefore to bring back another parameter. So these simply correspond to different parametrizations of the EoS, and we adopt the one described above in order to follow the traditional approach.
By using 7 layers we have introduced a large number of parameters ( and ). Most of those can be determined by continuity of various quantities at the layer boundaries. For the outer 6 layers we assume the continuity of the energy density at the boundaries, which allows us to determine the ’s:
(2.5) 
If the constant for the outermost layer is known, then the other values (except for the innermost layer) can be determined by the continuity of the pressure:
(2.6) 
For the outermost layer, the “crust”, we have . Requiring that (physically this means that the edge of the star is ordinary nonrelativistic matter) implies that . Thus the parameterization of the EoS of the NS for the outer layers will require us to specify the critical pressures , all the polytropic exponents as well as the outermost polytropic constant , while all other parameters will be determined by the continuity conditions.
2.2 Modeling the Core and the Effect of VE
For the last layer, we use an equation of state that incorporates physics associated with a change in the QCD vacuum state due to high density. There are two effects expected at this phase transition: a vacuum energy term in the fluid that is independent of baryon number density, and a jump in energy density across the boundary. Unlike in the outer layers, in the exotic phase the nature of the baryonic states may be very different from the usual zero temperature baryons. Since QCD conserves baryon number, for the innermost layer it is more natural to use baryon number density in place of the mass density as the variable parametrizing the EoS for the central core (). In this case, the equation of state can be written as:
(2.7)  
(2.8) 
Note that the vacuum energy appears with the opposite sign in the energy density and pressure, just as with the cosmological constant. Our goal is to see how sensitive neutron star observables are to the VE shift .
To keep the form of the EoS unchanged in the various layers we can introduce the density where is the ordinary neutron mass, and use this rescaled number density for the innermost layer. We can easily see that in terms of this rescaled density the EoS will have the same form as for the outer layers:
(2.9)  
(2.10) 
where are just redefinitions of the unknown constants parametrizing the EoS for the inner layer. We adopt this notation in order to stay close to the standard formalism used in the literature.
Let us now examine in detail the continuity (or jump) of the various quantities at the phase boundary between the sixth and the seventh (innermost) layer. The Israel junction conditions [32] still require that the pressure be continuous:
(2.11) 
but due to the appearance of the term this now requires a jump in from to (where ) and consequently also in from to . Since QCD conserves baryon number, another quantity that we need to require to be continuous is the chemical potential (that is we are assuming chemical equilibrium at the phase boundaries with conserved baryon number). The chemical potential at zero temperature is given by
(2.12) 
where is again the baryon number density. This relation holds even if the VE is nonzero. Therefore the jumps from to and from and (in our convention and ) are related to each other by
(2.13) 
The convexity of the free energy can be translated to . This latter form implies that the number density increases with pressure, yielding . This condition together with the continuity of the chemical potential tells us that the jump in energy density should also be positive, i.e. .
A typical phase transition will have both and nonvanishing, and this scenario is the focus of our studies. We choose to parametrize the jump in energy density such that it is proportional to the absolute value of the shift in VE:
(2.14) 
For each value of , , and , this condition, together with continuity of the chemical potential, fixes the values of and . This parametrization of the phase transition has the advantage that the limit reproduces the results obtained in the literature since both the mass density and the energy density become continuous in this case. In principle, could be taken to be zero, isolating the effects of vacuum energy from a jump in the energy density, corresponding to a second order phase transition. Here we will assume that the phase transition is first order with an accompanying jump in most quantities across the phase boundary and take . A final consistency condition is that both the full pressure and the partial pressure of the fluid, , must be positive. This implies that must satisfy .
3 Modeling Neutron Stars
After presenting the relevant physics of the dense QCD matter forming the interior of the NS we are now ready to review the usual method for calculating the structure of the interior of the NS. GW emission observed by LIGO/Virgo originates from the inspiral phase, when the stars are far apart relative to their radii. In this stage of the merger, the NS’s are still well approximated by nearly spherically symmetric static objects, with deviations described by a linear response in an expansion in spherical harmonics. In this paper we will ignore the effects of NS angular momentum but plan to further investigate that in a future publication. First we briefly review the equations relevant for the spherically symmetric solution and then present an overview of the perturbations due to the gravitational field of the other NS.
3.1 Spherically Symmetric Solutions
At lowest order, the stars are spherically symmetric, and their mass distribution is given by the solution to the TolmanOppenheimerVolkoff (TOV) equations [33]. These equations are easily derived by starting with a spherically symmetric metric ansatz
(3.1) 
and using the associated Einstein equations assuming a spherically symmetric fluid distribution with pressure and energy density . The resulting TOV equations are:
(3.2)  
(3.3)  
(3.4) 
where denotes differentiation with respect to the radial coordinate . The TOV metric provides the unperturbed solution around which the gravitational field of the second star will introduce perturbations that can be dealt with using a multipole expansion. From the solution to these equations, one obtains the internal structure of the star: the mass as a function of radius, as well as the thicknesses and masses of the various layers.
3.2 Tidal Distortion and Love Numbers
In a neutron star binary, each neutron star experiences gravitational tidal forces due to the other. This force squeezes the stars along the axis passing through both of their centers, and deforms the stars, inducing a quadrupole moment. The size of this induced quadrupole moment is determined by the structure of each neutron star, which can be characterized by its compactness and the stiffness of the EoS. These in turn depend on the physical properties of the dense QCD matter as described by its EoS discussed in the previous section. The effect of the induced quadrupole on gravitational wave data is to change the power emission as a function of time and frequency. Thus LIGO data on NS inspirals contains information about this tidal deformability, which depends on the equation of state of the matter making up the stars.
A common way to describe the deformability of a star is through the Love number. Love numbers were originally introduced in the study of Newtonian tides [34]. The application of Love numbers to gravitational waves produced in neutron star inspirals was initiated in refs. [9, 10], and further generalized in [35]. Detailed studies of the prospects for gravitational wave detection were provided in [11, 12, 13].
In the local rest frame of one star a small tidal field can be described in terms of a Taylor expansion of the Newtonian gravitational potential, or the timetime component of the metric tensor. There are two contributions, one from the effect of the distant star, and the other from the induced quadrupole moment. At large distances (using Cartesian coordinates, ) takes the form [11]
(3.5) 
Here parametrizes the external tidal gravitational field, and is the induced quadrupole moment. Both matrices are traceless and symmetric. To linear order in the response, the induced quadrupole is determined by the tidal deformability, , defined by
(3.6) 
One can then define a dimensionless quantity by
(3.7) 
where is the radius of the neutron star. This is referred to as the tidal Love number, and is the main physical observable. The advantage of this parametrization is that the Love number does not vary much with the size of the star, with typical Love numbers ranging from to as masses and equations of state are varied.
In order to determine , one performs the perturbative expansion of the solutions to the Einstein equations in the presence of an external gravitational field assuming a multipole expansion. Thus inside and near the star we will write the metric perturbation as an expansion in spherical harmonics . Due to the axial symmetry around the axis connecting the centers of the two stars, is zero, and since the tidal deformation is dominantly quadrupolar, with no dipole, the leading contribution is at [11]. Hence the full perturbed metric (where is the metric from (3.1)) is written as
(3.8) 
where and are the functions in the solution to the unperturbed spherically symmetric metric (3.1):
(3.9)  
(3.10) 
and the Einstein equations relate the functions and :
(3.11) 
Inserting the perturbed metric into Einstein’s equations results in a second order differential equation for :
(3.12) 
To find solutions, one starts with a series expansion of very near the core of the star, at small :
(3.13) 
The linear term drops out since the solution must be regular at . The size of the coefficient is linearly proportional to the size of the external perturbation, , and is not an intrinsic property of the star, as is clear from the fact that it is simply a normalization coefficient for the solution to the linear ODE for . One can thus pick this coefficient arbitrarily in numerically solving for . The tidal Love number, on the other hand is an intrinsic property, and the value for drops out in calculating it. The value for can be calculated once is found, and matched at large onto the metric ansatz in Eq. (3.5). It is given by
(3.14) 
where is the compactness parameter , and is obtained from the solution to evaluated on the surface of the star:
(3.15) 
Another dimensionless quantity, known as the dimensionless tidal deformability, is often found in the literature. It is obtained from the definition of by factoring out the in front:
(3.16) 
4 Results and Fits
We are now ready to present our results for the effects of adding a VE component to the innermost layer. We will use several different benchmark EoS’s for modeling the NS’s and investigate the consequences of the presence of VE in each of those cases. Two of the EoS’s are more “conservative” in the sense that the maximum stable NS mass that can be achieved just barely goes above (the maximum NS mass observed thus far is ). The two more conservative models are the AP4 [36] and SLy [20] EoS’s, which were also used as benchmarks by LIGO/Virgo [1]. We also consider the less restrictive model of Hebeler et al. [22] which permits larger masses, up to nearly . For the AP4 and SLy models we use the piecewise polytropic parametrization for all seven layers provided by Read et al. [21]. We have tabulated the corresponding parameters in Table 1. While the model of Hebeler et al. also uses a piecewise polytropic EoS for the innermost three layers, for the outer four layers corresponding to the crust they use a semianalytic expression. In their parametrization, the outer crust is modeled by the BPS EoS [37] assuming equilibrium^{1}^{1}1equilibrium corresponds to equal rates of neutron decay and proton capture of electrons., followed by a layer for which chiral EFT (valid up to the nuclear saturation density around 0.18 baryons) is used to obtain the EoS. This semianalytic expression is consistent with the piecewise polytropic approach of AP4 and SLy.
SLy  AP4  Hebeler  
See [22]  
See [22]  
See [22]  
Varying the EoS leads to more or less compact NS’s, whose deformability will also change. The compactness of the NS can be characterized by the radius of a NS with a fixed mass. The deformability describes how much the NS reacts to the presence of the gravitational field of the second NS in the binary merger event and is characterized by the tidal deformability. In the first subsection, we present our results for the mass versus radius, , curves of neutron stars with different nuclear EoS’s including the effect of VE, while in the second, we study the tidal deformability and comment on the potential for LIGO/Virgo to discern between models with different assumptions about the change in VE in exotic phases of QCD.
4.1 Results
We first present the results for the mass versus radius curves for the three different benchmark equations of state, whose parameters are displayed in Table 1. These three benchmarks cover a realistic range of possible EoS’s, with a wide variation in the maximum possible stable neutron star mass. We take care not to violate basic constraints on the behavior of high density QCD matter. For example, when pressures near the center of the star become very large, and relativistic effects dominate, one must ensure that the equation of state does not violate causality. Causality requires that the speed of sound in the fluid does not exceed the speed of light:
(4.1) 
However, using simple EoS models this condition is often violated for very large central pressures. Such violation does not imply the instability of the NS, but is rather an indication that the ansatz for the EoS is no longer a good approximation of the underlying nuclear physics in that region. Such causality violation would never arise in the true QCD equation of state at very high densities.
The true stability condition for the central pressures that a neutron star can support is given by
(4.2) 
This constraint arises from considerations of radial pulsation modes of the star, and is directly associated with stability of the fundamental mode of oscillation [38]. The relation in Eq. (4.2) above can be violated when the jump in the energy density is too large [30]:
(4.3) 
Above this bound, the mass of the NS no longer increases with increasing core pressure, and the NS is unstable [27, 28, 29, 30, 31].
We note that the condition in Eq. (4.2) can be satisfied for two stars of the same mass, but different internal pressures [39], corresponding to different phases in the core of the star. In such cases, the critical energy density jump exceeds that in Eq. (4.3) at the transition. However, even with this instability, one sometimes finds for , that there is a disconnected class of solutions that does not exceed the bound in Eq. (4.2). The possibility then arises that some of the exotic, disconnected solutions have the same masses as some of the normal, lower pressure solutions.
Which of the two conditions, causality or monotonicity, will limit the central pressure depends on the EoS. For AP4 and SLy, the limit is set by causality. This bound can be avoided by modifying the EoS via a “causal extension” [22] into the regions where the pressure exceeds the maximal value set by the causality bound. For the models we are working with, we have found that this extension simply flattens out the curves at the point where causality is violated, and hence does not change the value of the maximum mass significantly. For this reason we have chosen not to make this modification and ended the curves at the point where the speed of sound reaches unity.
The curve and the effect of VE for the Hebeler et al. EoS [22] are shown in Fig. 1. Each curve is obtained by varying the pressure at the center of the star but keeping all of the other parameters fixed. We have fixed in this plot, as well as all those that follow.^{2}^{2}2Taking to be small reduces the change in the curves relative the case, however small values of are not representative of most phase transitions, which are typically accompanied by a change in the energy density as well as the vacuum energy. When the central pressure is greater than , the value of becomes relevant and the other curves depart from the behavior of the case. Dotted parts of the curves correspond to unstable regions, i.e. solutions of the TOV equations in which the stability condition (4.2) is violated. The shaded region represents the most massive neutron star measured to date, with a mass of [40]. Notice that for some positive values of , i.e. when the jump in energy density is big enough according to Eq. (4.3), we find a second stable region which is disconnected from the main branch, as discussed above. This means that for a given mass, there are two possible types of stars, one with no exotic phase in the core, and another with a significant portion of the star in the new phase. This gives rise to interesting effects, both for curves and in GW observables. For example, assuming that the curve is the correct one, it would be possible to observe two neutron stars with significantly different radii. That is, there are two stable configurations for stars with the same mass. It is quite interesting that the physics of QCD may allow for a plethora of different compact objects, with population numbers depending on the conditions of their formation.
Our procedure for introducing the VE for this model is the following. In order to make sure that all values of considered here are compatible with a neutron star mass of , we stop the nexttoinnermost polytropic region as soon as the mass reaches . This corresponds to choosing a critical pressure . Once the critical pressure is reached, we transition into the innermost polytropic region with a nonzero VE, and we allow for the central pressure to be as high as possible without violating the causality bound.
Next we present results for the AP4 [36] and SLy [20] EoS models. The curves for the AP4 and SLy models with different values of the VE in the innermost layer are shown in Fig. 2. One can again see that up to a certain critical mass, the curves corresponding to different ’s in the innermost layer do coincide with each other. The reason for this is that below this mass the central pressure is not high enough to enter the exotic high density phase of QCD. The critical pressure for the phase transition to occur is set by which is an input of the model. For the AP4 and SLy models, and respectively which correspond to a critical mass of approximately for both models.
The plots for all three EoS’s show that the maximal mass of the neutron star decreases for both positive and negative values of VE. This is a generic feature of neutron star models with phase transitions with vacuum energy in our study, and is due to the jump in the energy density across the phase transition. This conclusion is similar to that obtained in previous works that study phase transitions without vacuum energy [41].
4.2 Tidal Deformabilities and LIGO/Virgo
Let us now discuss NS observables from GW’s emitted during the merger of NS’s. The frequency versus time behavior of the waveform of the emitted gravitational wave, usually expressed in terms of the “gravitational wave phase” appearing in the Fourier transform of the chirp, can be determined by an expansion in relativistic effects, starting at Newtonian order, and proceeding to postNewtonian corrections in the velocity. At dominant Newtonian order, where the two NS’s are approximated by point masses, the waveform depends only on a particular combination of the masses called the chirp mass:
(4.4) 
where is the reduced mass of the system [42]. For the recently observed merger event, GW170817, the chirp mass was measured to be .
Since this is the dominant contribution to the waveform, the chirp mass can be determined quite accurately. However, the individual masses must be extracted from higher order velocity corrections to the waveform, and are thus more difficult to constrain. At higher order, spincouplings are important as well, and without information about the stars’ rotational speeds and axes, precise extraction of the masses is impossible. This information is in principle retrievable from measurements of the waveform, but is difficult as it relies on data near the end of the inspiral, where current experiments lose sensitivity, and where full numerical simulation of the merger event may be necessary [43].
At present, the individual masses can only be estimated by using the chirp mass and some assumptions for the angular rotation frequency of the stars. For GW170817 in the lowspin case, the estimated mass range is – for the heavy star and – for the light star, while for the highspin case, there is considerably more possible variation in the masses: – for the heavy star and – for its less massive partner.
Similarly, it is not yet possible to measure individual tidal deformabilities. However, it is possible to constrain a weighted combination of the individual masses and deformabilities through their contribution to the gravitational wave phase at order . This socalled “combined dimensionless tidal deformability” is defined as
(4.5) 
For the recent event GW170817, the current constraint placed on is for the lowspin assumption and for the highspin case. In the lowspin case, the neutron star masses are probably too low to contain an exotic QCD phase, and thus event GW170817 would not contain information about VE. Of course, this may not be the case for future merger events, which may involve heavier NS’s. In the highspin case, however, the inner core could be in the exotic phase, and the constraints from GW170817 are relevant for studying VE.
The rest of this section contains our results for the effects of VE on the tidal deformabilities, which will be presented in a series of plots. Each plot will be presented both for the Hebeler et al. EoS, which allows for larger NS masses and hence larger effects from VE, as well as for the AP4 and SLy EoS’s, which cut NS masses off at and thus have smaller VE effects.

Fig. 3 shows the individual tidal deformabilities of both NS’s and the effect of VE using the Hebeler et al. EoS.

Fig. 4 translates the effects of VE into a fractional shift of the combined tidal deformability . This plot shows that the effect of VE can be as large as 70% and is generically sizable for the case of the Hebeler et al. EoS.

Fig. 10 shows the effect of VE on GW170817 (assuming the Hebeler EoS) where it can significantly change the allowed region of NS masses one would infer from the constraint on the tidal deformability.
The EoS parametrization of Hebeler et al. [22] allows for large possible deviations in the tidal deformability when a VE term is added to the central core in the new phase. In Fig. 3, we show the effect of varying the VE term for a selection of 3 different input chirp masses. The curves are obtained by fixing the chirp mass at a few representative values and then scanning over the mass of the heaviest star, . Typically it is found that the heavier the star, the smaller the tidal deformability. This is largely due to the fact that more massive stars typically have smaller radii, and thus respond less to external tidal fields.
We note that the curve in the third plot is composed of two separate branches, corresponding to the two separate stable stars with equal masses but different radii as explained in the previous section. The branch with the highest values of corresponds to only the most massive star laying in the disconnected branch of the curve, while in the other case both stars in the binary would come from the disconnected branch. Part of the reason why the deviations from the curve are significant here can be found directly in Fig. 1. Since there the maximum mass for the curve is close to , the curves that correspond to a nonzero can depart significantly without being excluded by the measurement of the most massive neutron star, .
As we are most interested in the changes brought about by considering nonvanishing VE, it is useful to introduce a variable that quantifies the relative shift in due to VE:
(4.6) 
where is the value of obtained by taking the VE term to zero.
The deviation as a function of the heavy star mass, , for the Hebeler et al. parametrization is shown in Fig. 4. The negative values for mean that introducing a VE term lowers the value of . In order to isolate as much as possible the effects that a nonzero value of has on the internal structure of the stars, we are comparing each point in a given curve with the corresponding event on the curve that has the same neutron star masses. Therefore, any deviation in the value of comes entirely from the change in the tidal deformabilities, .
Even with the more conservative SLy and AP4 models one still finds large deviations in for events with larger chirp masses. The case of the (smaller) chirp mass corresponding to GW170817 is displayed in Fig. 5, and the deviations in deformability are small. This is because the combined deformability is typically dominated by the contribution from the less massive star, which does not contain a core in the new phase where VE plays a role. However for higher chirp masses the effect of vacuum energy can be sizable even for the SLy and AP4 EoS’s, as shown in Figs. 6 and 7. As the chirp mass increases more of the star contains the new phase, and eventually both stars typically contain cores in the new phase, yielding the increased sensitivity to VE. The high chirp mass we have chosen for these figures corresponds, if the stars are of equal mass, to individual masses of , approaching that of the most massive NS observed to date. One can see that for this case the deviation can be as large as , even for these more conservative equations of state.
Since the chirp mass is the most accurately measured property of the NS merger, it is worthwhile to examine the dependence of (characterizing the sensitivity to VE) on the chirp mass. In Figs. 8 (Hebeler) and 9 (AP4 and SLy) we plot and as a function of the chirp mass. The superscript expresses the fact that, when evaluating the quantities in Eqs. (4.5) and (4.6), the mass of the heavy star is kept fixed at the maximal value allowed by the corresponding fixed value of . Fixing one of the stars to have maximal mass will generically (though not always) give the largest VE effect on . The important result of the plots is that the deviation increases substantially above a certain chirp mass denoted by the vertical gray line in the plots. This threshold corresponds to the chirp mass for which the lighter star mass also reaches the critical mass for the phase transition. Therefore, the large deviation can be understood from the fact that both stars are in the new phase.
In our final plot in Fig. 10, we display the limits that GW170817 places on VE assuming the Hebeler et al. parametrization. In particular, we note that including a VE term significantly changes the allowed range of the individual masses for the highspin assumption. The effect is less pronounced for the SLy and AP4 models. As more data on NS mergers are collected with some of those corresponding to mergers of more massive stars, strong limits could be placed on the EoS of dense nuclear matter. This will especially be true once the sensitivities for probing the tidal contributions to the gravitational wave phase further improve.
The future outlook is difficult to extrapolate, but promising. The constraints that are placed in the coming years will depend strongly on currently uncertain characteristics of NS binaries or NSblack hole mergers that will be captured by upcoming datataking runs at LIGO and other GW observatories [44]. Constraints will depend upon masses, spins, and branch populations in cases where there are multiple configurations with the same mass. In addition, the sensitivities of the experiments will evolve, and may be able to better capture higher order contributions to the waveforms. Finally, utilization of neutron stars as laboratories to study very high density physics and VE depends crucially on a precise theoretical calculation of the QCD equation of state at high densities [45]. Taking an optimistic viewpoint on these issues leads us to the conclusion that neutron star mergers can tell us about the interface of gravity and quantum field theory in a regime never before tested.
5 Conclusions
In this paper, we have argued that neutron star mergers can be a valuable tool for testing new phases of QCD at large densities, in particular for finding the contribution of a VE term in exotic high density phases. To study the effects of such a new phase on neutron star observables, we have started with the conventional layer parametrization of the EoS, then assumed a nonzero value for the VE in the innermost layer leading to a jump in the energy density. For the three benchmark models we have chosen, we have calculated the curves and tidal Love numbers for different values of the VE. By using those results, we have obtained individual tidal deformabilities and calculated the combined dimensionless tidal deformability parameter which can be constrained by neutron star mergers observed in gravitational wave observatories. We have found that for larger chirp masses, the nonzero VE at the innermost core can have an effect on the combined dimensionless tidal deformability parameter, hence future observations of neutron star merger chirps can place strong limits on the EoS of dense nuclear matter. We have also shown that for some parameters, introducing a nonzero VE can create a disconnected branch of stable neutron star solutions allowing the possibility of having two neutron stars of the same mass with significantly different radii. This possibility is unique to EoS’s which have a phase transition at the core, hence it is a smoking gun for new phases of QCD.
Acknowledgments
We thank Brando Bellazzini and Javi Serra for many helpful conversations, and for collaboration at early stages of this work. We also thank Mark Alford for a useful discussion. J.H. thanks the Aspen Center for Physics for its hospitality while some of this work was performed. C.C., C.E., and G.R. thank the Galileo Galilei Institute of the INFN in Florence, Italy for its hospitality while this work was in progress. C.E, J.H., and G.R. thank Cornell University for hospitality throughout this work. C.C. is supported in part by the NSF grant PHY1719877. C.E., J.H., and G.R. are supported in part by the DOE under grant DEFG0285ER40237. J.T. was supported in part by the DOE under grant DESC000999.
References
 [1] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. 119 (2017) 161101 grqc/1710.05832.
 [2] M. Baryakhtar, J. Bramante, S. W. Li, T. Linden, and N. Raj, “Dark Kinetic Heating of Neutron Stars and An Infrared Window On WIMPs, SIMPs, and Pure Higgsinos,” Phys. Rev. Lett. 119 (2017) 131801 hepph/1704.01577.
 [3] D. Croon, A. E. Nelson, C. Sun, D. G. E. Walker, and Z. Z. Xianyu, “HiddenSector Spectroscopy with Gravitational Waves from Binary Neutron Stars,” hepph/1711.02096.
 [4] J. Ellis, et al. “Search for Dark Matter Effects on Gravitational Signals from Neutron Star Mergers,” astroph.CO/1710.05540.
 [5] B. Bellazzini, C. Csaki, J. Hubisz, J. Serra, and J. Terning, “Cosmological and Astrophysical Probes of Vacuum Energy,” JHEP 1606 (2016) 104 astroph.CO/1502.04702.
 [6] D. D. Ivanenko and D. F. Kurdgelaidze, “Hypothesis concerning quark stars,” Astrofizika 1 (1965) 479; F. Pacini, “Highenergy Astrophysics and a Possible Subnuclear Energy Source,” Nature 209 (1966) 389; D. Boccaletti, V. de Sabbata, and C. Gualdi, “The red or violetshift of quasars without any source motion,” Nuovo Cim. 45 (1966) 513; N. Itoh, “Hydrostatic Equilibrium of Hypothetical Quark Stars,” Prog. Theor. Phys. 44 (1970) 291.
 [7] J. R. Ellis, J. I. Kapusta, and K. A. Olive, “Phase transition in dense nuclear matter with quark and gluon condensates,” Phys. Lett. B 273 (1991) 123.
 [8] F. Ozel, et al. “Measuring Neutron Star Radii via Pulse Profile Modeling with NICER,” Astrophys. J. 832 (2016) no.1, 92 astroph.HE/1512.03067; S. Bogdanov, et al. 2017, AAS/High Energy Astrophysics Division, 16, 104.04; Gendreau, K and Arzoumanian, Z, “Searching for a Pulse” Nature Astronomy 1 (2017) 895.
 [9] E. E. Flanagan and T. Hinderer, “Constraining neutron star tidal Love numbers with gravitational wave detectors,” Phys. Rev. D 77 (2008) 021502 astroph/0709.1915.
 [10] T. Hinderer, “Tidal Love numbers of neutron stars,” Astrophys. J. 677 (2008) 1216 astroph/0711.2420.
 [11] T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, “Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral,” Phys. Rev. D 81 (2010) 123016 astroph.HE/0911.3535.
 [12] S. Postnikov, M. Prakash, and J. M. Lattimer, “Tidal Love Numbers of Neutron and SelfBound Quark Stars,” Phys. Rev. D 82 (2010) 024016 astroph.SR/1004.5098.
 [13] B. D. Lackey and L. Wade, “Reconstructing the neutronstar equation of state with gravitationalwave detectors from a realistic population of inspiralling binary neutron stars,” grqc/1410.8866.
 [14] B. Margalit and B. Metzger, “Constraining the Maximum Mass of Neutron Stars From MultiMessenger Observations of GW170817,” astroph.HE/1710.05938.
 [15] A. Bauswein, O. Just, N. Stergioulas, and H. T. Janka, “Neutronstar radius constraints from GW170817 and future detections,” astroph.HE/1710.06843.
 [16] E. Zhou, A. Tsokaros, L. Rezzolla, R. Xu, and K. Ury, “Uniformly rotating, axisymmetric and triaxial quark stars in general relativity,” astroph.HE/1711.00198.
 [17] L. Rezzolla, E. R. Most, and L. R. Weih, “Using gravitationalwave observations and quasiuniversal relations to constrain the maximum mass of neutron stars,” astroph.HE/1711.00314.
 [18] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, “Gravitationalwave constraints on the neutronstarmatter Equation of State,” astroph.HE/1711.02644.
 [19] J. M. Lattimer and M. Prakash, “Neutron star structure and the equation of state,” Astrophys. J. 550 (2001) 426 astroph/0002232; K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, “Constraints on neutron star radii based on chiral effective field theory interactions,” Phys. Rev. Lett. 105 (2010) 161102 nuclth/1007.1746; J. M. Lattimer, “The nuclear equation of state and neutron star masses,” Ann. Rev. Nucl. Part. Sci. 62 (2012) 485 nuclth/1305.3510.
 [20] F. Douchin and P. Haensel, “A unified equation of state of dense matter and neutron star structure,” Astron. Astrophys. 380 (2001) 151 astroph/0111092; P. Haensel, A. Y. Potekhin, and D. G. Yakovlev eds., “Neutron Stars 1: Equation of State and Structure” (2007, New York, Springer).
 [21] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, “Constraints on a phenomenologically parameterized neutronstar equation of state,” Phys. Rev. D 79 (2009) 124032 astroph/0812.2163.
 [22] K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, “Equation of state and neutron star properties constrained by nuclear physics and observation,” Astrophys. J. 773 (2013) 11 astroph.SR/1303.4662.
 [23] H. Mueller and B. D. Serot, “Relativistic mean field theory and the high density nuclear equation of state,” Nucl. Phys. A 606 (1996) 508 nuclth/9603037; E. S. Fraga, R. D. Pisarski, and J. SchaffnerBielich, “Small, dense quark stars from perturbative QCD,” Phys. Rev. D 63 (2001) 121702 hepph/0101143; H. Heiselberg, “Neutron star masses, radii and equation of state,” astroph/0201465; S. Banik and D. Bandyopadhyay, “Color superconducting quark matter core in the third family of compact stars,” Phys. Rev. D 67 (2003) 123003 astroph/0212340; J. L. Zdunik and P. Haensel, “Maximum mass of neutron stars and strange neutronstar cores,” Astron. Astrophys. 551 (2013) A61 astroph.SR/1211.1231.
 [24] M. G. Alford, A. Schmitt, K. Rajagopal, and T. Schäfer, “Color superconductivity in dense quark matter,” Rev. Mod. Phys. 80 (2008) 1455 hepph/0709.4635.
 [25] “Equations of state and stability of colorsuperconducting quark matter cores in hybrid stars,” B. K. Agrawal, Phys. Rev. D 81 (2010) 023009 astroph.HE/1001.1584; R. Anglani, et al. , “Crystalline color superconductors,” Rev. Mod. Phys. 86 (2014) 509 hepph/1302.4264; E. S. Fraga, A. Kurkela, and A. Vuorinen, “Interacting quark matter equation of state for compact stars,” Astrophys. J. 781 (2014) 2, L25 nuclth/1311.5154.
 [26] M. Bejger, D. Blaschke, P. Haensel, J. L. Zdunik, and M. Fortin, “Consequences of a strong phase transition in the dense matter equation of state for the rotational evolution of neutron stars,” Astron. Astrophys. 600 (2017) A39 astroph.HE/1608.07049; M. A. R. Kaltenborn, N. U. F. Bastian, and D. B. Blaschke, “Quarknuclear hybrid star equation of state with excluded volume effects,” Phys. Rev. D 96 (2017) 056024 astroph.HE/1701.04400.
 [27] D. E. AlvarezCastillo and D. B. Blaschke, “Highmass twin stars with a multipolytrope EoS,” Phys. Rev. C 96 (2017) 045809 nuclth/1703.02681; J. E. Christian, A. Zacchi, and J. SchaffnerBielich, “Classifications of Twin Star Solutions for a Constant Speed of Sound Parameterized Equation of State,” astroph.HE/1707.07524.
 [28] M. G. Alford and A. Sedrakian, “Compact stars with sequential QCD phase transitions,” Phys. Rev. Lett. 119 (2017) 161104 astroph.HE/1706.01592.
 [29] M. G. Alford, S. Han, and M. Prakash, “Generic conditions for stable hybrid stars,” Phys. Rev. D 88 (2013) 083013 astroph.SR/1302.4732; M. G. Alford and S. Han, “Generic Conditions for Stable Hybrid Stars,” EPJ Web Conf. 80, 00038 (2014); M. G. Alford, G. F. Burgio, S. Han, G. Taranto, and D. Zappalà, “Constraining and applying a generic highdensity equation of state,” Phys. Rev. D 92 (2015) 083002 nuclth/1501.07902.
 [30] Z. F. Seidov, “The Stability of a Star with a Phase Change in General Relativity Theory,” Soviet Astronomy 15 (1971) 347.
 [31] R. Shaeffer, P. Haensel, and L. Zdunik, “Phase transitions in stellar cores,” Astronomy and Astrophysics 126 (1983) 121; L. Lindblom, “Phase transitions and the mass radius curves of relativistic stars,” Phys. Rev. D 58 (1998) 024008 grqc/9802072.
 [32] W. Israel, “Singular hypersurfaces and thin shells in general relativity,” Nuovo Cim. B44, (1966) 1, Erratum: Nuovo Cim. B48 (1967) 463.
 [33] R. C. Tolman, “Static Solutions of Einstein’s Field Equations for Spheres of Fluid,” Phys. Rev. 55 (1939) 364; J. R. Oppenheimer and G. M. Volkoff, “On Massive Neutron Cores,” Phys. Rev. 55 (1939) 374.
 [34] A. E. H. Love, “Some Problems of Geodynamics” (1911, Cambridge University Press, Cambridge).
 [35] T. Damour and A. Nagar, “Relativistic tidal properties of neutron stars,” Phys. Rev. D 80 (2009) 084035 grqc/0906.0096; T. Damour and A. Nagar, “Effective One Body description of tidal effects in inspiralling compact binaries,” Phys. Rev. D 81 (2010) 084016 grqc/0911.5041; T. Binnington and E. Poisson, “Relativistic theory of tidal Love numbers,” Phys. Rev. D 80 (2009) 084018 grqc/0906.1366.
 [36] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, “The Equation of state of nucleon matter and neutron star structure,” Phys. Rev. C 58 (1998) 1804 nuclth/9804027.
 [37] G. Baym, C. Pethick and P. Sutherland, “The Ground state of matter at high densities: Equation of state and stellar models,” Astrophys. J. 170 (1971) 299; J. W. Negele and D. Vautherin, “Neutron star matter at subnuclear densities,” Nucl. Phys. A 207 (1973) 298.
 [38] J. M. Bardeen, K. S. Thorne, and D. W. Meltzer, “A Catalogue of Methods for Studying the Normal Modes of Radial Pulsation of GeneralRelativistic Stellar Models,” Astrophysical Journal, vol. 145, (1966), 505.
 [39] N. Glendenning, C. Kettner, “Possible third family of compact stars more dense than neutron stars,” Astronomy and Astrophysics 353 (2000) L9 astroph/9807155; K. Schertler, C. Greiner, J. SchaffnerBielich, and M. H. Thoma, “Quark phases in neutron stars and a ’third family’ of compact stars as a signature for phase transitions,” Nucl. Phys. A 677 (2000) 463 astroph/0001467; P. Haensel, “Equation of state of dense matter and maximum mass of neutron stars,” EAS Publ. Ser. 7 (2003) 249 astroph/0301073.
 [40] P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, “Shapiro Delay Measurement of a Two Solar Mass Neutron Star,” Nature 467 (2010) 1081 astroph.HE/1010.5788.
 [41] H. Heiselberg and M. HjorthJensen, “Phase transitions in neutron stars and maximum masses,” Astrophys. J. 525, L45 (1999) astroph/9904214.
 [42] C. Cutler and E. E. Flanagan, “Gravitational waves from merging compact binaries: How accurately can one extract the binary’s parameters from the inspiral wave form?,” Phys. Rev. D 49, 2658 (1994) grqc/9402014.
 [43] T. Dietrich, S. Bernuzzi and W. Tichy, “Closedform tidal approximants for binary neutron star gravitational waveforms constructed from highresolution numerical relativity simulations,” Phys. Rev. D 96 (2017) 121501 grqc/1706.02969.
 [44] P. Kumar, M. Pürrer, and H. P. Pfeiffer, “Measuring neutron star tidal deformability with Advanced LIGO: a Bayesian analysis of neutron star  black hole binary observations,” Phys. Rev. D 95 (2017) 044039 grqc/1610.06155.
 [45] L. Rezzolla and K. Takami, “Gravitationalwave signal from binary neutron stars: a systematic analysis of the spectral properties,” Phys. Rev. D 93 (2016) 124051 grqc/1604.00246.