Gravitational-wave Signature of a First-order Quantum Chromodynamics Phase Transition in Core-Collapse Supernovae

A first-order quantum chromodynamics (QCD) phase transition (PT) may take place in the protocompact star (PCS) produced by a core-collapse supernova (CCSN). In this work, we study the consequences of such a PT in a non-rotating CCSN with axisymmetric hydrodynamic simulations. We find that the PT leads to the collapse of the PCS and results in a loud burst of gravitational waves (GWs). The amplitude of this GW burst is $\sim30$ times larger than the post-bounce GW signal normally found for non-rotating CCSN. It shows a broad peak at high frequencies ($\sim2500-4000$ Hz) in the spectrum, has a duration of $\lesssim5 {\rm ms}$, and carries $\sim3$ orders of magnitude more energy than the other episodes. Also, the peak frequency of the PCS oscillation increases dramatically after the PT-induced collapse. In addition to a second neutrino burst, the GW signal, if detected by the ground-based GW detectors, is decisive evidence of the first-order QCD PT inside CCSNe and provides key information about the structure and dynamics of the PCS.


INTRODUCTION
Quarks are confined in hadrons such as protons and neutrons at low temperatures and densities. Nonetheless, free quarks should exist in the early universe when the temperature is extremely high (k B T > ∼ 150 MeV) [1]. They may also exist in the cold and superdense interior of compact stars with a density above the nuclear saturation density (ρ sat 2.6 × 10 14 g cm −3 ) [see e.g., [2][3][4]. Moreover, a first-order quantum chromodynamics (QCD) phase transition (PT), i.e., hadron-quark PT, may take place in the protocompact stars (PCS) produced by a core-collapse supernova (CCSN) [5][6][7] or binary neutronstar merger [8,9]. Such a PT can result in a more compact PCS and even collapse of the PCS to a black hole (BH). For CCSNe, this can provide an additional energy source for the explosion [7], and leads to interesting observational consequences, such as a second neutrino burst [5] and the production of rare r−process elements [10].
A galactic CCSN is a yet-undiscovered candidate gravitational-wave (GW) source for ground-based GW detectors [11]. The GWs from a CCSN, in combination with the neutrino and electromagnetic signals, will boost our understanding of the CCSN explosion mechanism [12]. Sophisticated multi-dimensional simulations have predicted the GW signals emitted by rotational CCSNe, the oscillations of the proto-neutron stars (PNS), and the standing accretion shock instability [see, e.g., [13][14][15][16][17]. In the meantime, the collapse of a neutron star (NS) to a quark star has been studied in [18][19][20] and GW emission is also found in this scenario. However, in these studies the collapse is triggered artificially by using different equations of state (EOS) for the construction of NS and the hydrodynamic simulation. It is unclear whether the PT-induced collapse of PCS in CCSNe can leave an imprint in the GW signal. In this Letter, we demonstrate the effects of a first-order QCD PT on the GW signal from a non-rotating CCSN, with two-dimensional simulations and a simplified hybrid EOS including hadrons and quarks.

Equation of state
To study a first-order QCD PT in CCSNe, we use a hybrid EOS from Ref. [5,21,22]. This EOS employs the STOS EOS [23] for the hadronic phase and the MIT bag model EOS [2] for the quark phase, with the Gibbs con-arXiv:2007.04716v2 [astro-ph.HE] 29 Oct 2021 struction for the mixed phase [24]. The Bag constant has a value of B = (165 MeV) 4 . The Gibbs construction enforces global charge neutrality; therefore, it allows different charge fractions for the hadronic and quark portions in the mixed phase. This leads to a smooth transition from the hadronic phase to the mixed phase (at around ρ sat for the composition of a CCSN core), and the pressure is continuously increasing. A pure quark phase is realized at densities higher than ∼ 3.5ρ sat , which is similar to the values of those more sophisticated hybrid EOSs [7,9]. More information about the hybrid EOS can be found in Ref. [5,21,23].
For this hybrid EOS, there exists a stable branch of the hot third family of compact stars [4] with a pure quark core (see Fig. 2 in [25]) for matter properties similar to the CCSN core (entropy 3 k B /baryon and lepton number fraction ∼ 0.4). The maximum mass for the third family (∼ 1.50 M ) is larger than that of the second family whose core is in the mixed phase. As we will see, this unique property is important for the dynamics of the PCS.

FLASH simulation
We carry out CCSN simulations in two dimensions with the assumption of axisymmetry, using the FLASH code [26] with an "M1" scheme for the neutrino transport [27]. We take the 12-M , solar-metallicity, presupernova progenitor s12 from [28] as the initial conditions. To apply general-relativistic approximations, gravity is calculated with the Case A formulation of [29]. Unlike previous simulations with the FLASH code, we include the lapse function in the Euler equations to mimic the time-dilation effect in general relativity (see the modified Euler equations and some numerical tests in [30], also see [31]). This is found to affect the GW frequency significantly after the PT-induced collapse. A cylindrical grid with adaptive mesh refinement is used. It extends out to 2 × 10 9 cm in radius and ±2 × 10 9 cm along the cylindrical axis, with a finest resolution of 150 m. We extract the plus GW strain h + from our Newtonian simulation using the standard quadrupole formula [32].

Dynamics
To show the consequences of the PT, we run two simulations with the same settings, one using the hybrid EOS and the other using the STOS EOS. The resulting dynamics are shown in the upper panel of Fig. 1. The iron core of the s12 model collapses to above ρ sat and bounces at t b 151 ms for both EOSs. At t b , the core of the hybrid EOS has already entered the mixed phase with a central quark mass fraction X q = 18%. However, because the hybrid EOS transitions smoothly from the pure hadronic phase to the mixed phase, the PCS remains in the mixed phase with a low X q shortly after t b . The bounce shock turns into an accretion shock and stalls at ∼ 150 km at t b + 50 ms and begins receding inward.
During the accretion phase, the central density ρ c of the PCS with the hybrid EOS is always larger than that of the PNS with the STOS EOS and X q continuously increases. The mass of the PCS grows and reaches the maximum of the second family for the hybrid EOS at ∼ t b + 286 ms. The PCS becomes unstable against gravity and experiences a second dynamical collapse. The central density ρ c grows to 1.5 × 10 15 g cm −3 (∼ 6ρ sat ) and the PCS core enters the pure quark phase (X q = 1).
The pure quark core bounces in less than 1 ms at t 2b as the PCS enters the new stable branch of the third family and ρ c drops to ∼ 5ρ sat . This bounce shock expands quickly to explode the outer envelope. At the end of the simulation, the mean shock radius extends out to ∼1500 km with an explosion energy of ∼ 2.0 × 10 50 erg. The PT-induced collapse is associated with a second neutrino burst with more electron anti-neutrinos than electron neutrinos (lower panel of Fig. 1), which is consistent with the results of the spherically-symmetric simulation in [5].

Gravitational waves
In Fig. 2 we show the GW waveforms h + (t) up to 400 ms after the first bounce (∼ t 2b + 113 ms) extracted from both simulations. We assume that the distance from the source is 10 kpc. The signal from t b to ∼ t b + 50 ms comes from the prompt convection behind the stalling accretion shock. It is followed by an episode of continuous emission from the oscillations of the PCSs [33]. There is no qualitative difference between the two waveforms until t 2b and the cumulative emitted GW energies are quantitatively similar (bottom of Fig. 2). In accord with the more compact PCS, the peak GW frequency for the hybrid EOS is always higher than that for the STOS EOS.
Around t 2b , the PT-induced collapse results in a burst of GW emission with a much larger amplitude than those of earlier episodes. In Fig. 2, the 10 ms window around t 2b is stretched in time to show clearly this burst, which is associated with the PT-induced collapse and bounce. The maximum amplitude of h + reaches 10 −20 and ∼ 30 times larger than those of the other episodes. The energy carried by this burst is ∼ 4.6 × 10 −7 M c 2 , which is ∼ 3 orders of magnitudes more than the GW energy of the other episodes (and also that of the signal for the STOS EOS). Our numerical test shows that this GW burst results from asphericities developed between t b and t 2b [30]. After this burst, the amplitude damps quickly to the same level as before t 2b . This part of signal should come from the oscillations of the PCS with a pure quark core.
A time-dependent spectrum (or spectrogram) is useful for understanding the emission mechanisms of GWs, as well as designing efficient detection strategies. shows the spectrogram of the GW signal extracted from the simulation using the hybrid EOS. We use a Kaiser window with a width of 25 ms for the short-time Fourier transform except for around t 2b where a width of 10 ms is used. Before t 2b , the spectral evolution is similar to that of the STOS EOS (see [30]). The GW peak frequency is continuously increasing, in accord with the evolution of the Brunt-Väisälä frequencies f BV (Eq. (3) in [30]) at densities between 10 11 and 10 12 g cm −3 (blue band in Fig. 3), which is approximately the PCS surface [34]. Around t 2b , the GW burst has a much higher frequency (∼ 2500 − 4000 Hz). This is related to the change of the dominant GW emission region from ∼ 10 − 20 km to ∼ 5−10 km (see [30]), which is inside the quasi-static core during the second collapse and bounce. During this time, f BV peaks at 2900 Hz near a radius of 10 km (ρ 8×10 13 g cm −3 ) and is closer to the observed GW frequency.
Shortly after t 2b , the GW peak frequency drops back to ∼ 1000 Hz and continues to increase afterwards, albeit at a much faster rate. We find that f BV has a much larger spread inside the PCS and the GW spectral evolution does not match the track of f BV . After t b + 300 ms, the peak frequency of the dominant GW emission is closer to f BV at densities ∼ 5 × 10 12 g cm −3 . Nevertheless, due to the much larger ρ c ( > ∼ 4 times) and compactness of the PCS for the hybrid EOS (Fig. 1), the peak GW frequency is 2 − 3 times higher than that for the STOS EOS. The GW spectral evolution after t 2b contains information about the structure and evolution of the PCS with a pure quark core, from which one may infer the properties of the quark EOS (e.g. Bag constant). of the GW signals from t b to t b +400 ms with the STOS (blue) and hybrid (orange) EOSs. Also shown are the h char (f ) of the GW signal for the hybrid EOS in the time interval of: from t b to t b +280 ms (green) and between t 2b −3ms and t 2b +7ms ms (red). The black line is the sensitivity spectrum of Advanced LIGO [39].

Detection prospect
To estimate the detectability of the GW signals, we calculate the dimensionless characteristic GW strain (h char ) [35] assuming a distance of 10 kpc, and compare it with the sensitivity of Advanced LIGO in Fig. 4. Below ∼ 1000 Hz, h char are quantitatively similar for the hybrid and STOS EOSs. At higher frequencies, h char for the hybrid EOS shows a broad peak between ∼ 2500−4000 Hz, which is also above the detector's sensitivity curve. This part is mainly contributed by the burst associated with the PT-induced collapse, seen from the comparison between the entire h char and that between t 2b − 3ms and t 2b + 7ms.
We calculate the single-detector signal-to-noise ratio (SNR) of the GW waveforms assuming the optimal orientation using Eq. (1.1) in [35]. If a confident detection requires an SNR of 8, then for the hybrid EOS, inclusion (exclusion) of the burst yields a detection radius of 22 (12) kpc. The detectability of the burst is not significantly better because the current detectors are optimized for GW signals at ∼ 10 − 1000 Hz. The amplitude of h char and detector noise increase together by a factor of 10 from 100 Hz to 3000 Hz (Fig. 4). Future-generation detectors, such as the Einstein Telescope [36] and the Cosmic Explorer [37], may consider improving the sensitivity at several kHz if such signals are targeted (also for BH forming CCSNe [38]).
In this work, the GW waveforms h + (t) are extracted from 2D simulations with the assumption of axisymmetry. In various studies [e.g., 16,17,40], the amplitude of h + in 3D simulations can be 5 − 10 times smaller than those in their 2D counterparts, which lowers the expectation for the GW detection. Future study is needed to explore 3D effects for the PT-induced collapse and observables. We expect that the burst associated with the PT-induced collapse would still be present in 3D simulations but with smaller amplitudes. DISCUSSION We present here a specific case in which the PTinduced collapse results in a bounce shock that successfully explodes the mantle. However, for other progenitors [41] or other hybrid EOSs [5], the star may fail to explode and collapse into a black hole (BH) in two scenarios. First, the PCS at the onset of the PT-induced collapse may exceed the maximum mass which the hybrid EOS permits and it directly collapses into a BH. In this case, the GW (and neutrino) burst reported here will be absent. Nevertheless, the existence of free quarks in the PCS might be inferred from the shortening of the BH formation time [6,42], though it is subjected to the uncertainties of the pure hadronic EOS.
In the second scenario, the second bounce shock is launched but the PCS still collapses into a BH at a later time. In this case, the burst of GWs and neutrinos associated with the PT-induced collapse will be present, followed by the shut-off of both signals at BH formation. This is an interesting possibility to be explored. Moreover, in both cases of BH formation, if the iron core is rapidly rotating, the inclusion of PT can produce different BH ring-down signals compared to those for a hadronic EOS due to the different free-fall time of the PCS, which is found in binary NS-merger simulations [8].

CONCLUSIONS
In this Letter, we demonstrate the effects of a firstorder QCD PT on the GW signals from a non-rotating CCSN. We find that the PT results in the collapse of the PCS at ρ c ∼ 3.5ρ sat , and the core radiates a loud GW burst in < ∼ 5 ms. The amplitude of this burst reaches h + = 10 −20 assuming a source distance of 10 kpc and is larger by a factor of ∼ 30 than other episodes of GW emission (and generally those using a hadronic EOS). The spectrum of this burst shows a broad peak at ∼ 2500 − 4000 Hz, which is higher than that generally found for CCSNe without the PT-induced collapse. The peak GW frequency following this burst is also much higher (>1 kHz) than that for the hadronic EOS due to the large compactness of the PCS with a pure quark core. Therefore, the PT inside a CCSN can be inferred from the GW detection. However, the louder burst is not necessarily easier to detect because of the increasing noise level of high frequencies for current ground based GW detectors. Nevertheless, the loud, high-frequency burst of GW radiation over a short period of time may be a prime target for future searches of coherent wave burst signals [43].
The hybrid EOS transitions from the hadronic phase to the mixed phase at a low density (∼ ρ sat ). Ref. [7] simulated CCSNe in spherical symmetry using a more physical and complex hybrid EOS (DD2F-SF, transition density ∼ 2.4ρ sat ), and the dynamics of the second collapse are similar to our results. Therefore, we expect that the properties of the GW burst (i.e. the frequency and amplitude) associated with the PT-induced collapse should still be present with a more physical EOS, such as DD2F-SF [7,9] or a Chiral Mean Field model [8], which are consistent with the maximum NS mass measurement [44]. A natural extension is to employ such EOSs in multi-dimensional CCSN simulations. Moreover, progenitor dependence, such as the initial mass and rotation, should be studied to acquire a more comprehensive picture for the effects of a first-order QCD PT on the GWs from CCSNe. Particularly, we expect that if the iron core is rapidly rotating before collapse, the GW burst associated with the PT-induced collapse will be much louder and it may allow the detection for sources farther away. The FLASH code solves the Newtonian Euler equations. To mimic the deeper gravitational well in general relativity (GR), we use the effective gravitational potential with the Case A formalism of [29], which has been tested in corecollapse supernova (CCSN) simulations routinely [27,42,45].
In this work, we extend this by including the GR time-dilation effect directly in the Euler equations. The modified Euler equations read: where ρ, v, P and τ are the rest mass density, velocity, pressure and kinetic plus internal energy density of the fluid; Φ is the effective GR gravitational potential; α = exp(Φ) is the lapse function. The lapse function is included in the fluxes and source terms. The additional source term in the momentum equation, αP/c 2 ∇Φ, is verified by the derivation of the GR hydrodynamic equations with the metric g µν = [−α 2 , 1, r 2 , r 2 sin θ] in spherical symmetry [46]. This source term is important for maintaining the mechanical equilibrium in the hydrostatic regions. The energy equation is derived through the same procedure and has no additional source term. The consistency of the equations can be checked with a polytropic equation of state in which the pressure and specific internal energy are analytic functions of the rest mass density. To test the performance of the GR approximations, we simulated the oscillations of a compact star constructed using the Case A potential and hybrid EOS. We choose an initial central density of 1.5 × 10 15 g cm −3 to mimic the protocompact star (PCS) after the second bounce in our CCSN simulations. For reference, we also performed a simulation with the GR1D code using a fully relativistic TOV star with the same central density as the initial conditions. The initial conditions are different for the FLASH and GR1D simulations because the fully relativistic 0.5 TOV solution is not a stable configuration for the FLASH simulations (also see the Appendix A of [27]). The results of the central density evolution for 5 ms are shown in Fig. S1. The small-amplitude oscillations originate from the imperfect mapping of the initial conditions to the computational grids of the hydrodynamic simulations. This mapping leads to slightly different equilibrium states (central densities) in different codes. The frequencies of the PCS oscillations are ∼ 2500 Hz, 4500 Hz and 3300 Hz for the simulations using GR1D, FLASH and FLASH with the inclusion of the lapse function. This test shows that our implementation of the lapse function can approximate the GR time-dilation effect to some extent.
For the CCSN simulations in the main text, we find that the lapse function reduces the gravitational-wave (GW) frequency significantly, especially after the second collapse when the PCS is extremely compact. For example, the frequency of the burst around the second collapse is ∼ 2500 − 4000 Hz (∼ 4000 − 5000 Hz) for the simulation with (without) the lapse function. We expect full GR simulations will further reduce the GW frequency.

Test of potential numerical artifacts
We perform a test simulation to evaluate the contribution of potential numerical artifacts from the computation domain to the GW signal during the PT-induced collapse and bounce. The test simulation starts from a sphericallysymmetric compact star constructed using the STOS EOS, while the hybrid EOS is used for the 2D hydrodynamic simulation. The sudden reduction of pressure results in the collapse of the compact star. We compare the results of this test simulation to those of the PT-induced collapse in the CCSN simulation in Fig. S2. In the test simulation, the amplitude of the GW strain remains less than 0.2 × 10 −21 until 2 ms after bounce, which is at least an order of magnitude smaller than that of the burst in the CCSN simulation. This indicates that the GW burst in the CCSN simulation results from the asphericity already developed during the episode between the first bounce and the second collapse, but not artifacts from the dynamical collapse itself.

Resolution dependence
We perform a set of simulations with different resolutions, for the episode of the second collapse and bounce in the main text. The simulations start from ∼ 10 ms before the second collapse and the finest resolutions are 300 m, 150 m and 75 m, respectively. We find the PCS structure after the second bounce agrees well for the different resolutions. In Fig. S3, we plot the GW waveform and spectrum for the loud burst around the second collapse. Although the GW signals do not match exactly in phase, the amplitude and frequency of the burst agrees quantitatively well for the different resolutions. Because the high-resolution simulation has less numerical dissipation, the amplitude of h + damps more slowly in the high-resolution runs (75 m and 150 m) than in the low-resolution run (300 m) after t 2b + 2 ms.
Spatial distribution of GW strain Fig. S4 shows the spatial distribution of the tangential velocity (v θ , left half) and GW strain (right half): at 2 ms before (left) and 0.5 ms after (right) the second bounce. Here we assume the distance from the source is D = 10 kpc. The results are from the simulation with a finest resolution of 150 m. The distribution of v θ roughly shows the convective regions inside the PCS, which changes from ∼ 10 − 20 km to ∼ 5 − 10 km. The distribution of h + is correlated to that of v θ , which suggests the connection between the GW emission and convective motions inside the PCS. The contribution to h + from the regions outside 15 km is generally less than ∼ 10% during the burst around the second collapse.
Spectrogram for the STOS EOS Fig. S5 shows the spectrogram of the GW signal extracted from the simulation using the STOS EOS. This spectrogram is similar to those of CCSN simulations using a hadronic EOS in the literature, e.g. [38]. The green bands shows the time-dependent Brunt-Väisälä frequency (f BV ) at densities between 10 11 and 10 12 g cm −3 , where f BV is independent of the specific position. We estimate f BV by the Newtonian Brunt-Väisälä frequency multiplied by the lapse function: where c s is the speed of sound and others are the same as in Eq. 1. The evolution of the peak GW frequency roughly follows that of f BV for the 400 ms after bounce.