Neutron Star Mergers Chirp About Vacuum Energy

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 non-negligible 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.


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 Λ 4 QCD (where Λ QCD ∼ 200 MeV is the usual QCD scale) making NS's the only known objects where VE might make up a non-negligible 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 low-temperature 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 quasi-particle 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 p c ∝ Λ 4 QCD , and is expected to also lead to a change of VE [7] of order ∆Λ ∝ Λ 4 QCD . 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 well-studied NS models and the effects of VE on those observables. Finally we conclude in Section 5.

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 (100 MeV) 4 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 (meV) 4 : 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: 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: |∆Λ| ∼ Λ 4 QCD [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.

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 coarse-graining 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 non-relativistic 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 ρ(r) = m n n(r) rather than a power of the energy density (as would be natural for a high-density, 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 where the a i are integration constants. Note that the appearance of the a i 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 (γ i , K i and a i ). 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 a i 's: If the K 1 constant for the outermost layer is known, then the other K i values (except for the innermost layer) can be determined by the continuity of the pressure: For the outermost layer, the "crust", we have p 0 = 0. Requiring that lim ρ→0 ρ = 1 (physically this means that the edge of the star is ordinary non-relativistic matter) implies that a 1 = 0. Thus the parameterization of the EoS of the NS for the outer layers will require us to specify the critical pressures p i , all the polytropic exponents γ i as well as the outermost polytropic constant K 1 , while all other parameters will be determined by the continuity conditions.

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 n in place of the mass density as the variable parametrizing the EoS for the central core (p > p 6 ). In this case, the equation of state can be written as: 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 ρ = m n n where m n 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: where K 7 =K 7 /m γ 7 n , (1 + a 7 ) =ã 7 /m n 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: but due to the appearance of the Λ term this now requires a jump in ρ(r) from ρ + to ρ − (where ρ − = ρ 6 ) and consequently also in (r) 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 where n 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 ρ + = m n n + and ρ − = m n n − ) are related to each other by The convexity of the free energy ∂ 2 F ∂V 2 T,N > 0 can be translated to ∂p ∂n T,N > 0. 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 Λ non-vanishing, 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 γ 7 , α, and Λ, this condition, together with continuity of the chemical potential, fixes the values of K 7 and a 7 . This parametrization of the phase transition has the advantage that the Λ = 0 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 α > 0. A final consistency condition is that both the full pressure and the partial pressure of the fluid, K 7 ρ γ 7 , must be positive. This implies that Λ must satisfy −p 6 < Λ.

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.

Spherically Symmetric Solutions
At lowest order, the stars are spherically symmetric, and their mass distribution is given by the solution to the Tolman-Oppenheimer-Volkoff (TOV) equations [33]. These equations are easily derived by starting with a spherically symmetric metric ansatz and using the associated Einstein equations assuming a spherically symmetric fluid distribution with pressure p(r) and energy density (r). The resulting TOV equations are: where denotes differentiation with respect to the radial coordinate r. 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.

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 time-time 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, x i ) g tt takes the form [11] Here E ij parametrizes the external tidal gravitational field, and Q ij 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 One can then define a dimensionless quantity k 2 by where R is the radius of the neutron star. This is referred to as the l = 2 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 k 2 = 0.001 to k 2 = 1 as masses and equations of state are varied.
In order to determine k 2 , 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 Y m l . Due to the axial symmetry around the axis connecting the centers of the two stars, m is zero, and since the tidal deformation is dominantly quadrupolar, with no dipole, the leading contribution is at l = 2 [11]. Hence the full perturbed metric g αβ + h αβ (where g αβ is the metric from (3.1)) is written as where e µ(r) and e ν(r) are the functions in the solution to the unperturbed spherically symmetric metric (3.1):  To find solutions, one starts with a series expansion of H very near the core of the star, at small r: (3.13) The linear term drops out since the solution must be regular at r = 0. The size of the coefficient a is linearly proportional to the size of the external perturbation, E ij , 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 H. One can thus pick this coefficient arbitrarily in numerically solving for H. The l = 2 tidal Love number, on the other hand is an intrinsic property, and the value for a drops out in calculating it. The value for k 2 can be calculated once H is found, and matched at large r onto the metric ansatz in Eq. (3.5). It is given by where C is the compactness parameter GM/R, and y is obtained from the solution to H evaluated on the surface of the star: Another dimensionless quantity, known as the dimensionless tidal deformability, is often found in the literature. It is obtained from the definition of k 2 by factoring out the C 5 in front:λ (3.16)

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 2M (the maximum NS mass observed thus far is M = (2.01 ± 0.04)M ). 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 3M . 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  [37] assuming β equilibrium 1 , followed by a layer for which chiral EFT (valid up to the nuclear saturation density around 0.18 baryons/fm 3 ) is used to obtain the EoS. This semi-analytic expression is consistent with the piecewise polytropic approach of AP4 and SLy.
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, M (R), 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.

M (R) 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: 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 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]: 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 p > p 6 , 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 The M (R) 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 α = 3 in this plot, as well as all those that follow. 2 When the central pressure is greater than p 6 , the value of Λ becomes relevant and the other curves depart from the behavior of the Λ = 0 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 (2.01 ± 0.04)M [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 M (R) curves and in GW observables. For example, assuming that the Λ = (165 MeV) 4 curve is the correct one, it would be possible to observe two 2M 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 (2.01 ± 0.04)M , we stop the next-to-innermost polytropic region as soon as the mass reaches 2.00M . This corresponds to choosing a critical pressure p 6 ≈ (150 MeV) 4 . 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 M (R) 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 p 6 which is an input of the model. For the AP4 and SLy models, p 6 ≈ (179 MeV) 4 and p 6 ≈ (176 MeV) 4 respectively which correspond to a critical mass of approximately 1.6M 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].

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 post-Newtonian 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: where µ is the reduced mass of the system [42]. For the recently observed merger event, GW170817, the chirp mass was measured to be M = 1.188 +0.004 −0.002 M . 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, spin-couplings 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 low-spin case, the estimated mass range is 1.36-1.60M for the heavy star and 1.17-1.36M for the light star, while for the high-spin case, there is considerably more possible variation in the masses: 1.36-2.26M for the heavy star and 0.86-1.36M 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 v 5 . This so-called "combined dimensionless tidal deformability" is defined as For the recent event GW170817, the current constraint placed onΛ is ≤ 800 for the lowspin assumption and ≤ 700 for the high-spin case. In the low-spin 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 high-spin 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 2M 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 de-formabilityΛ. 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.
• Figs. 5, 6, 7 repeat the same analyses for the AP4 and SLy EoS's, where we see that the deviations are generically smaller, but can still reach 25-30% for larger chirp masses.
• Figs. 8, 9 emphasize the role of the chirp mass: they show the maximal achievable effect of VE as the chirp mass is increased. We can see that observing events with chirp masses above 1.6-1.8M will be key to observing the effects of VE.
• 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 M at a few representative values and then scanning over the mass of the heaviest star, M 1 . 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 Λ = (165 MeV) 4 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λ 2 corresponds to only the most massive star laying in the disconnected branch of the M (R) 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 Λ = 0 curve are significant here can be found directly in Fig. 1. Since there the maximum mass for the Λ = 0 curve is close to 3M , the curves that correspond to a nonzero Λ can depart significantly without being excluded by the measurement of the most massive neutron star, (2.01 ± 0.04)M . As we are most interested in the changes brought about by considering non-vanishing VE, it is useful to introduce a variable that quantifies the relative shift inΛ due to VE: whereΛ 0 is the value ofΛ obtained by taking the VE term to zero.
The deviation as a function of the heavy star mass, M 1 , 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 Λ = 0 curve that has the same neutron star masses. Therefore, any deviation in the value ofΛ comes entirely from the change in the tidal deformabilities,λ i .  chirp mass is taken to be M = 1.188M , which is the same value as in GW170817.λ 1 and λ 2 correspond to the dimensionless tidal deformability parameters for the heavy and light stars, respectively. Each curve is obtained by varying the heavy star mass while holding the chirp mass fixed. The α-parameter of (2.14) is chosen to be α = 3.
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 1.9M , approaching that of the most massive NS observed to date. One can see that for this case the deviation can be as large as 37%, 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Λ max 0 and δ max 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 assum-  For the smaller chirp mass the effect is rather small, however for a higher chirp mass the effect can be as large as 38%. The α-parameter of (2.14) is again chosen to be α = 3.
ing 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 high-spin 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 NS-black hole mergers that will be captured by upcoming data-taking 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  For the smaller chirp mass the effect is rather small, however for a higher chirp mass the effect can be as large as 25%.
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.

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 7-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 M (R) 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 O(1) 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.