Metastable hypermassive hybrid stars as neutron-star merger remnants

Hypermassive hybrid stars (HMHS) are extreme astrophysical objects that could be produced in the merger of a binary system of compact stars. In contrast to their purely hadronic counterparts, hypermassive neutron stars (HMNS), these highly differentially rotating objects contain deconfined strange quark matter in their slowly rotating inner region. HMHS and HMNS are both mestastable configurations and can survive only shortly after the merger before collapsing to rotating black holes. The appearance of the phase transition from hadronic to quark matter in the interior region of the HMHS and its conjunction with the emitted GW will be addressed in this article by focussing on a specific case study of the delayed phase-transition scenario that takes place during the post-merger evolution of the remnant. The complicated dynamics of the collapse from the HMNS to the more compact HMHS will be analysed in detail. In particular, we will show that the interplay between the spatial density/temperature distributions and the rotational profiles in the interior of the wobbling HMHS after the collapse generates a high-temperature shell within the hadron-quark mixed phase region of the remnant.


Introduction
On September 14, 2015, almost exactly a hundred years after Albert Einstein developed the field equations of General Relativity and predicted the existence of GWs, these curious spacetime-ripples have been observed from a pair of merging black holes by the LIGO detectors (GW150914, [61]). Using the GW detectors LIGO and Virgo, two signals so far have been associated with the collision of two neutron stars: GW170817 [2] and GW190425 [1]. For GW170817 electromagnetic radiation in all frequency ranges was also detected during this event [62] and an emitted gamma-ray burst (GRB 170817A, [37]) hit the gamma-ray satellite telescopes with a delay of 1.7 s. Space-based gamma-ray telescopes (e.g., the Fermi's gamma-ray burst monitor or the Swift gamma-ray burst mission) detect on average approximately one gamma-ray burst per dayhowever, the gamma-ray burst [37] that had been associated with GW170817 is an outstanding event and, in addition with the observations of the electromagnetic counterparts of the associated kilonova, provides a conclusive picture of the whole merger event. This coincidence of the direct detection of a GW from a neutrona e-mail: hanauske@th.physik.uni-frankfurt.de (corresponding author) star collision with the emitted short gamma-ray burst was the first observational proof that binary neutron star (BNS) mergers generate short gamma-ray bursts. Finally, with the use of the observed tidal deformabilities of the two neutron stars from the late inspiral phase and other properties of GW170817, the equation of state (EOS) of dense matter could be severely constrained [3,4,13,19,22,40,44,47,48,53].
Taking the results of various BNS merger simulations into account, it is apparent that the overall picture of the GW events GW170817 and GW190425 matches well with the results obtained in numerical simulations. These simulations show that right after the neutron stars merge, a metastable hypermassive neutron star (HMNS) is formed. This remnant either collapses promptly and leaves no electromagnetic counterpart, which is the most likely scenario for GW190425 given the high binary mass and the absence of an electromagnetic counterpart or it collapses delayed within less than a second to a rotating Kerr black hole (see [27] for the case of GW170817). Alternatively, after redistributing its angular momentum it could turn into a supramassive neutron star with smaller gravitational mass and which is stable for longer times. In the case of a collapse, a short gamma-ray burst is emitted, releasing in less than 1 s the energy emitted by our Galaxy over one year [52]. To understand the post-merger evo- lution of a BNS merger, numerous numerical-relativity simulations have been performed in the last decades by various groups (see [11] for a review). The emitted GW spectrum [20,59], the impact of initial spin [15,24,35], eccentricity [23,46], and mass ratio [25] have been analyzed in detail, however, the possibility of a hadron-to-quark phase transition (HQPT) and its effect on the GW signal have been addressed only recently [12,43,65]. The temperatures and densities reached during the merger and post-merger evolution of the remnant are extreme and a phase transition to deconfined quark matter is most likely to happen in the HMNS, which then becomes a hypermassive hybrid star (HMHS [5]).
A natural question that emerges from these considerations is: Can we detect the HQPT? In future GW detections of BNS merger systems it might be possible to detect the QCD phase transition by analysing the spectrum of the emitted GW during the merger and post-merger phase. The post-merger GW signatures of phase transitions in binary mergers of compact stars have been recently systematically analysed [65]. While in the prompt phase transition (PPT) scenario [12] a strong first-order phase transition leads, immediately after merger, to the formation of a gravitationally stable extended pure-quark matter core in the merger remnant, in the case of a phase-transition triggered collapse (PTTC), the phase transition takes place in the postmerger evolution and leads to the rather rapid collapse of the HMHS to a Kerr black hole. We here discuss in more detail the third possibility of a delayed phase transition (DPT) scenario that was recently presented in [65]. Like in the PTTC case, the phase transition within the DPT scenario takes place at late post-merger times, however, the collapse of the HMNS is halted by the formation of a metastable HMHS and the corresponding GW signatures of a HQPT are more distinct than in the PTTC case.
In particular, the properties of the HMHS produced in a DPT scenario will be discussed here. We focus on the specific case study of an irrotational, equal mass BNS simulation using the FSU2H-PT EOS. In contrast to [65], where we mostly described the gravitational wave signal, we here provide more details on the density, temperature and rotation profiles inside the inner area of the HMHS. We also provide a detailed explanation on how these thermodynamic properties create a large m = 1 oscillation mode connected to the postmerger GW signal and for the first time show the strong correlation between the central density and the GW frequency. The paper is structured as follows: Sect. 2 summarizes briefly the essential aspects of the microphysical matter treatment and shows the static properties of hybrid stars obtained within the FSU2H-PT EOS. In Sect. 3, the numerical setup of the BNS merger simulation will be briefly summarized. In Sect. 4 the results of the DPT simulation will be discussed in detail. The collapse from the HMNS to the HMHS, the properties of the produced HMHS and the GWs emitted during its post-merger evolution will be compared to a purely hadronic evolution. Finally, a summary and outlook will be presented in Sect. 5.

Microphysical description
For the cold, purely hadronic part of the EOS, a relativistic mean-field model (FSU2H, for details see [63,64]) has been used. To account for a first-order transition from hadronic to quark matter, a mixed phase and pure-quark matter region has been included. The overall cold EOS (FSU2H-PT) used within our simulation consists of a piecewise polytrope representation [54] of the FSU2H model (purely hadronic part), a soft mixed phase region with polytropic index Γ = 1.04 in the rest-mass density range ρ/ρ 0 ∈ [2.085, 4.072] and a stiff pure-quark matter part (Γ = 5.1) for densities above 4.072 ρ 0 . With these special choices of the parameters determining the phase transition, a strong HQPT has been modeled and, additionally, the recent constraints from electromagnetic and GW detections, namely the constraints on the maximum mass [9,40,53,56,58], radius [8,18,22,39,44,51,60] (which are also compatible with recent constraints from the NICER mission [41,55]) and tidal deformability of neutron stars [3,21,22,36], are fullfilled. To guarantee that the sound speed of pure-quark matter is always below the speed of light, three additional piecewise polytropes have been added for ultra-high densities (ρ/ρ 0 = 4.823, 4.969, 5.289; Γ = 4.7, 4.1, 3.1). To describe the low-density regime (ρ < 0.005 ρ 0 ), namely the crust consisting of leptons and nuclei, the standard prescription of Baym, Pethick and Sutherland has been used [14]. For the intermediatelow density regime (0.005 < ρ/ρ 0 < 0.65), the results of Negele and Vautherin have instead been employed [45]. The special choice of parameters in the FSU2H-PT EOS produces a mass-radius relation with a small twin-star branch (see the supplemental material of [65]), belonging formally to the twin-star Category III (see [42] for details).
To account for additional shock heating during the merger and post-merger phase, thermal effects are included by adding a "thermal" ideal-fluid component (p th = ρ th (Γ th − 1)) to the cold EOS where ρ is the rest-mass density, and Γ th = 1.75. The effective temperature obtained within this ideal-gas approach can be roughly approximated as T = (m n p th )/(k B ρ), where m n is the nucleonic mass and k B the Boltzmann constant. It should be stressed that the estimated temperature within this simple approximation only accounts for contributions of the ideal gas of neutrons and protons (mass differences have been neglected). Especially at the low-density regions of the outer crust of the hybrid star (which is composed of gas of leptons and nuclei at high values of electron fraction Y e ), the underlying temperature estimates deduced from the thermal pressure should be decreased by a factor (1 + Y e ) and account for the contribution due to the radiation pressure [see Eq.
(2) in [26]]. Finally, the thermal component of the FSU2H-PT EOS does not account for any contri-bution resulting from a high-temperature hadron-quark crossover transition and, as a result, the structure of the phase diagram for T > 80 MeV is different from that resulting from fully temperature-dependent EOSs (see e.g., [43]).

Mathematical and numerical setup
The simulations were performed using the generalrelativistic hydrodynamics code WhiskyTHC [49,50] and the dynamical spacetime solver McLachlan [17], which is part of the Einstein Toolkit [38]. These codes solve the combined set of the Einstein within the CCZ4formulation [6,7] and of the general-relativistic hydrodynamic equations in a conservative form [54].
The initial configuration of the irrotational equalmass binary is computed with the LORENE code [28], where we have used an initial separation of 45 km and a total binary mass of 2.64 M . The grid, provided by the Carpet mesh-refinement driver [57], consists of six refinement levels with the finest level having a resolution of dx = 0.16 M ( 237 m) and uses a reflection symmetry across the z = 0 plane to reduce the computational cost, but does not employ a π-symmetry to allow for the appearance of m = 1 features triggered by the phase transition (see discussion below). We have performed the simulations also at a higher (dx 198 m) and lower (dx 368 m) resolution and verified the qualitative consistency of our results.

Dynamics of the hypermassive hybrid star
The signals detected from GW170817 and GW190425 can only give information about the late-inspiral phase of the merging systems, when the stellar matter is still rather cold and the rest-mass densities are slightly below the central rest-mass densities of the initial configurations. In the late-inspiral phase, the maximum rest-mass density value of the tidally deformed purely hadronic neutron stars is slightly decreasing (see green area in Fig. 1). During the last few orbits, the tidal deformations increase and the temperature in the outer stellar layers increases rapidly as the two stars collide [29,32]. The time of the merger (t = 0) has been defined at the moment when the amplitude of the emitted GW reaches its first maximum.
Shortly after the merger, the maximum rest-mass density of the remnant increases, reaching a maximum at t ≈ 0.5 ms (see the pronounced first local maximum marked by a star symbol in Fig. 1), that is high enough so that a mixed phase including quark matter is created in the core. For this reason we refer to this object as a HMHS. The temperature of the mixed-phase matter within the central region of the HMHS is at moderate values (T < 20 MeV), whereas in the inner and outer shocked regions, spots of high-temperature values can increase up to 100 MeV at moderate densities ρ < ρ 0 . Next, the density quickly drops again and for t ∈ [0.5, 2] ms the density everywhere in the remnant is low enough so that there is no mixed phase present, i.e. the HMHS turns into an HMNS. Then for t > 2 ms the evolution of the FSU2H-PT simulation (red line in Fig. 1) starts to differ significantly from the simulation without phase transition (NPT; black line in Fig. 1) and for t > 3.5 ms, the evolution of the DPT scenario differs qualitatively from the NPT case. This can be seen in the left column of Fig. 2, which reports the distributions of the density/temperature at t = 3.46 ms (top and middle panels) and visualizes these quantities in a (T, ρ) phase diagram (bottom panel). The DPT occurs in the transient post-merger phase and as a result, the density profile shows a pronounced doublecore structure (for details see [30]). In the time window t ∈ [3.5, 4.0] ms the HMHS, which at this point still consists mostly of hadronic matter, collapses rapidly to a smaller and more compact configuration with a significant pure-quark core, whose overall angular velocity increases. The pure-quark core now comprises ≈ 20-30% of the total binary mass.
Before we analyse the properties and dynamics of the HMHS generated from the DPT collapse, the results of the NPT simulation will be briefly illustrated. After the violent, early post-merger phase of the HMNS, two high-temperature hot-spots emerge. The high-temperature values (T ≈ 40 MeV) are reached in regions where the rest-mass density is in a range ∼ (1 − 2)ρ 0 , while the maximum-density values are always located in the center at moderate temperatures T < 5 MeV. The dynamics of the underlying fluid can be visualized using tracer particle (see [16] for details). Some of the tracer particles circle around the high-temperature hot-spots, others populate the low temperature-high density inner region, and some are in the outer surface of the HMNS within the lowdensity region. At later post-merger times, the HMNS reaches a quasi-stationary configuration and the temperature hot-spots are smeared into a ring-like structure in the equatorial plane (this structure is effectively a thick torus in the three dimensions). The "peanut" shape of the density distribution will also dissolve and the populated area in the (T, ρ) plane will consist of only a small region that is quasi-stationary. The left column in Fig. 2 shows the configuration of the HMHS at a time right before the collapse to the more compact star. The small asymmetry in the density profile and especially the double-core structure is amplified by the collapse, resulting in a large one-sided asymmetry (i.e., an m = 1 asymmetry in a sphericalharmonics decomposition [10]), which triggers a sizeable h 21 + GW strain (see supplemental material of [65] for details). It is reasonable that this asymmetry is present also in the vertical direction but we cannot measure it because of the imposed symmetry across the z = 0 plane.
The configuration following this collapse is shown in the middle column of Fig.2, which corresponds to a time near the second-density maximum at t = 6.25 ms (see  in Fig. 1). The large m = 1 contribution can be seen by looking at the asymmetry of the spatial location of the pure-quark core, which is marked with the blue contour line in the top row of Fig. 2. As a result of this asymmetry, the location of the two temperature hot-spots (see second row, middle column) are at different radial distances from the grid center. Additionally, the temperature maximum of the hot-spot which is further away from the center is higher (T ≈ 95 MeV) than the more central hot-spot (T ≈ 65 MeV), as seen from the double-peaked structure near ρ/ρ 0 ∼ 2 in the bottom panel, middle column, of Fig. 2.
The following evolution of the HMHS is characterized by large density oscillations (see Fig. 1 for t > 4 ms), during which the density maxima reached in the center of the HMHS (ρ max ) and the emitted instantaneous GW frequency f GW are highly (but not exactly) correlated. In particular, the time interval between consecutive peaks of the density oscillations, i.e., Δt ρ , and of the instantaneous graviational-wave frequency, i.e., Δt GW , differs by less than 5% as can be seen from Fig. 3. These oscillations are accompanied by a change of the differential rotation profile of the HMHS. More specifically, when the rest-mass density of the HMHS reaches its maximum during each oscillation, the HMHS is more compact and spins faster and, as a result, the emitted GW frequency also reaches a maximum. Following these large macroscopic oscillations strong shocks are created within the soft mixed-phase region of the HMHS, which are accompanied by local increases of the temperature within the mixed-phase region. The interesting aspect of this correlation between the density and the GW oscillations is that it is in principle possible to determine with good precision from the evolution of f GW when the hadron-to-quark PT actually takes place.
The right-most panels of Fig. 2 report the HMHS properties at t = 13.15 ms and shows that in addition to the two temperature hot-spots, a new high-temperature shell surrounding a cold-core appears within the mixed phase region of the remnant. For subsequent postmerger times, the two temperature hot-spots (and the maxima of the angular-velocity profile; cf. Fig. 4) will be smeared out to become a ring-like structure on the equatorial plane. Furthermore, the high-temperature sphere in the mixed phase will increase in value until the m = 1 density asymmetry has dissolved and the HMHS has reached a metastable configuration.
Finally, as discussed in the case of BNS simulations with purely hadronic EOSs (see [31,[33][34][35]), the spatial location of the temperature hot-spots and the subsequent torus-like structure are closely connected with the angular-velocity profile Ω of the differentially rotating HMHS. The distribution of the latter on the equatorial plane is depicted in Fig. 4 at two representative times, and shows that the two maxima of the angularvelocity profile occur at the same spatial location as the two temperature maxima. Interestingly, however, the hot-spot with the highest temperature coincides with the location of the smaller rotation maxima (cf. mid-  dle and right columns in the central row of Fig. 2); this behaviour is likely due to the conservation of the Bernoulli constant, which anticorrelates the evolution of the angular velocity to that of the density [33]. Last but not least, it should be mentioned that the hightemperature shell, which is present in the mixed phase region of the HMHS has no counterpart in the angular velocity profile, as it is caused by the softness of the EOS.

Summary and outlook
Gravitational waves emitted from merging neutronstar binaries open a new observable window to explore the high-density/temperature regime of nuclear matter. The post-merger GW signatures of a HQPT in a binary merger of compact stars have recently been analyzed systematically in [65]. This paper presents an extension to this study and focuses on the DPT scenario, which has been possible when considering the specific case of the FSU2H-PT EOS. In particular, we have shown that the collapse of the merger remnant to a black hole can be stopped by the formation of a differentially rotating HMHS. Due to the large stiffness of the pure-quark phase, the collapse is hindered and the late post-merger dynamics (i.e., for t > 4 ms) is characterized by large macroscopic oscillations of the maximum density. The amplitude and frequency of these oscillations are mainly determined by the properties of the mixed-and pure-quark matter phases (e.g., they depend on the total extent of the mixed-phase region and on the sound speed in the pure-quark phase) and will be studied systematically in a subsequent work. Overall, the scenario presented here highlights that with future GW detections of binary compact-star mergers and by analysing the post-merger GW emission of the differentially rotating HMHS, it might be possible to detect the QCD phase transition and hence the formation of a quarkgluon plasma in the present universe. Theoretische Physik at Goethe Universitaet Frankfurt. The simulations were performed on the SuperMUC and SuperMUC-NG clusters at the LRZ in Garching, on the LOEWE/GOETHE cluster in CSC in Frankfurt, and on the HazelHen cluster at the HLRS in Stuttgart.
Funding Information Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.