Constraints on the maximum mass of neutron stars with a quark core from GW170817 and NICER PSR J0030+0451 data

The QMF (Zhu et al., 2018) and DD2 (Fortin et al., 2016) models are employed to represent fiducial soft and stiff hadronic matter, respectively. To avoid uncertainties due to the crust, we utilize unified EOSs with a consistent treatment of NS crusts and the crust-core transition properties along with their cores. For NSs without a quark core, the soft QMF EOS leads to a maximum mass $M_{\rm TOV}=2.07\,{\rm M}_{\odot}$ and a typical radius $R_{\rm 1.4}=11.77\,\rm km$, whereas for the stiff DD2 EOS $M_{\rm TOV}=2.42\,{\rm M}_{\odot}$ and $R_{\rm 1.4}=13.17\,\rm km$. Both EOSs are constructed within the widely used relativistic mean-field (RMF) approach, based on an effective Lagrangian with meson fields mediating strong interactions between quarks or hadrons. Compatibility with available experimental as well as observational constraints at sub-nuclear and higher densities is ensured (e.g., Drischler et al., 2020).

For high-density quark matter the “constant-speed-of-sound” (CSS) parameterization (Alford et al., 2013) is applied, and the EOS from the onset of the phase transition up to the maximum central density of a star is determined by three dimensionless parameters: the transition density $n_{\rm trans}/n_{0}$, the transition strength $\Delta\varepsilon/\varepsilon_{\rm trans}$, and the sound speed squared in quark matter $c^{2}_{\rm QM}$ (we work in units where $\hbar=c=1$). The full EOS is

where $n_{\rm trans}\equiv n_{\rm HM}(p_{\rm trans})$ and $\varepsilon_{\rm trans}\equiv\varepsilon_{\rm HM}(p_{\rm trans})$, and $\Delta\varepsilon/\varepsilon_{\rm trans}$ is essentially the finite discontinuity in the energy density at the phase boundary. In practice, our study is limited to first-order phase transitions with a sharp interface (Maxwell construction).

Since there exists a one-to-one correspondence between the underlying EOS and global properties of NSs such as ($M,R,\Lambda$), predicted values for these quantities can therefore be characterized by the three CSS parameters ($n_{\rm trans}/n_{0},\Delta\varepsilon/\varepsilon_{\rm trans},c_{\rm QM}^{2}$) in the hybrid EOSs. In Fig. 1 we show four parameter sets as examples for the high-density quark matter together with the unified QMF and DD2 EOSs for the inner crust and the outer core. Before reporting our results, we first introduce the scheme of the Bayesian inference approach that allows direct and consistent utilization of the LIGO/Virgo and NICER observational data.

In the present analysis, the likelihood is defined as

where $\mathcal{L}_{M_{s}}$ is the likelihood to encapsulate the mass measurement of a pulsar source. $\mathcal{L}_{\rm GW}$ and $\mathcal{L}_{\rm PSR}$ are the likelihoods of the GW170817 data and the mass-radius measurement of PSR J0030+0451 by NICER, respectively. If we only focus on the GW (PSR) data, we will neglect the PSR (GW) data.

Lower bound on $M_{\rm TOV}$ from MSP J0740+6620 (catalog ). The mass measurements of massive pulsars establish a firm lower bound on the NS maximum mass. Only the EOSs that support a $M_{\rm TOV}$ larger than this lower bound can pass this constraint, while others will be rejected. Assuming the mass measurement is a Gaussian distribution N$(\mu,\sigma^{2})$, the likelihood can be written as

where $\Phi(x)\equiv\int_{-\infty}^{x}(2\pi)^{-\frac{1}{2}}e^{-\frac{x^{2}}{2}}dx$ is the standard Gaussian cumulative distribution function. In the present analysis, we choose the highest mass measured through Shapiro delay, $M=2.14_{-0.09}^{+0.10}\,{\rm M}_{\odot}$ (68% confidence level) of MSP J0740+6620 (catalog ) (Cromartie et al., 2020), to place the $M_{\rm TOV}$ constraint. We fit the mass distribution of this source by a Gaussian distribution with $\mu=2.14\,\,{\rm M}_{\odot}$ and $\sigma=0.1\,\,{\rm M}_{\odot}$.

GW170817. The first binary NS merger GW170817 detected by the LIGO/Virgo collaboration (Abbott et al., 2017a) provided significant insight into the matter EOS and the NS internal structure because of their effects on the gravitational waveform. Assuming that the noise in the LIGO/Virgo detectors is Gaussian and stationary, the functional form of likelihood is often expressed as

where $S_{\rm n}(f)$, $\tilde{d}(f)$, and $\tilde{h}(f,\vec{\theta}_{\rm GW})$ denote the power spectral density (PSD), the Fourier transform of the measured strain data, and the frequency domain waveform generated using the parameter set $\vec{\theta}_{\rm GW}$, respectively. In the present analysis we adopt a 128s strain data¹¹1https://www.gw-openscience.org/eventapi from GPS time 1187008682s to 1187008890s and the PSDs²²2https://doi.org/10.7935/KSX7-QQ51 of GW170817 (Abbott et al., 2019b) released by the LIGO/Virgo collaboration. We also take the waveform model IMRPhenomD_NRTidal (Dietrich et al., 2017) for illustrative purposes.

PSR J0030+0451 (catalog ). The recent simultaneous mass-radius measurements of PSR J0030+0451 from NICER, $M=1.44^{+0.15}_{-0.14}\,{\rm M}_{\odot}$, $R=13.02^{+1.24}_{-1.06}~{}\rm km$ (Miller et al., 2019) and $M=1.34^{+0.15}_{-0.16}\,{\rm M}_{\odot}$, $R=12.71^{+1.14}_{-1.19}~{}\rm km$ (Riley et al., 2019), to the $68.3\%$ credibility interval, provide an additional useful constraint on the NS EOS. Since results in Riley et al. (2019) and Miller et al. (2019) are consistent with each other, we select the best fit scenario within the ST+PST model of Riley et al. (2019). In this case, the likelihood of generating the mass-radius distribution of PSR J0030+0451 (catalog ) is expressed as

where the right-hand side is a Gaussian Kernel Density Estimation (KDE) of the posterior samples $\vec{S}$ of the mass and radius given by Riley et al. (2019).

The parameters related to the hybrid star EOSs are $\vec{\theta}_{\rm EOS}=\{n_{\rm trans}/n_{0},\Delta\varepsilon/\varepsilon_{\rm trans},c_{\rm QM}^{2}\}$ within the CSS parameterization. Given that a normal NS reaches the maximum-mass configuration with its central density below $\approx 7.0\,n_{0}$/$6.0\,n_{0}$ for the QMF/DD2 hadronic EOS, we choose a uniform distribution of $n_{\rm trans}/n_{0}$ for hybrid stars in the range of $[1,7]$/$[1,6]$, respectively. For the sound speed squared in quark matter $c_{\rm QM}^{2}$, we vary its value from $c_{\rm QM}^{2}=1/3$ (the conformal limit in perturbative QCD matter) to $c_{\rm QM}^{2}=1$ (the causal limit), with a uniform distribution. We also assign a reasonably wide range for the energy density discontinuities as $\Delta\varepsilon/\varepsilon_{\rm trans}\in[0,2]$ with a uniform distribution.

For GW170817, the GW parameters are chirp mass in the detector frame, mass ratio, tidal deformabilities, aligned spins, inclination angle, geocentric time, polarization of GW, orbital phase at coalescence, red-shift of the source, and its position in the sky, i.e., $\vec{\theta}_{\rm GW}=\{\mathcal{M}^{\rm det},q,\Lambda_{1},\Lambda_{2},\chi_{1z},\chi_{2z},\theta_{\rm jn},t_{\rm c},\Psi,\varphi,z,\alpha,\delta\}$. We fix the location of the source to the position determined by EM observations (Abbott et al., 2017b; Levan et al., 2017) with $\alpha\,({\rm J}2000.0)=13^{\rm h}09^{\rm m}48^{\rm s}.085$, $\delta\,({\rm J}2000.0)=-23^{\circ}22^{\prime}53^{\prime\prime}.343$ and $z=0.0099$. $\mathcal{M}^{\rm det}$, $q$, $\chi_{1z}(\chi_{2z})$, $t_{c}$ and $\Psi$ are chosen to distribute uniformly in the range $[1.18,1.21]M_{\odot}$, $[0.5,1.0]$, $[-0.05,0.05]$, $[1187008882,1187008883]$s, and $[0,2\pi]$. For the inclination angle $\theta_{\rm jn}$, we assume a uniform prior distribution.³³3Both the GW data (Abbott et al., 2017a) and the EM data (Mooley et al., 2018) as well as structured jet modeling (Troja et al., 2020) can give some constraints on the viewing angle for GW170817. However, since the constrained viewing angle in these studies gives a very wide distribution, the uniform prior distribution can serve for the purpose of our study. The tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$ are related to the component-star masses through the hybrid star EOS, i.e.,

where the component-star masses $m_{1}$, $m_{2}$ can be obtained with

To improve the convergence rate of the nest sampler, we marginalize the likelihood over the orbital phase at coalescence.

The mass and radius entering the likelihood of PSR data are related to the hybrid star EOS via the central pressure of the star, $p_{c}$, i.e.,

Therefore we only need one additional parameter $p_{c}$ to perform the analysis if we take account of the PSR data.

We carry out four main tests to investigate how each data set (see details in Sec. 3.1) would affect the result, namely:
(i) $2.14\,{\rm M}_{\odot}$ pulsar: where we apply the minimum $M_{\rm TOV}$ constraint from MSP J0740+6620 (catalog );
(ii) +GW170817: where we consider both the constraint in (i) and the GW data of GW170817;
(iii) +NICER: where we consider both the constraint in (i) and the NICER measurement of PSR J0030+0451;
(iv) +GW170817+NICER: where we combine the data of PSR J0030+0451 and GW170817 all together with the maximum-mass constraint in (i).

Figure 3 illustrates the posterior distributions of the hybrid EOS, and the most probable values for the phase transition parameters and their $90\%$ confidence boundaries are summarized in Table 1, $2^{\rm nd}$ to $4^{\rm th}$ columns. The $5^{\rm th}-7^{\rm th}$ columns of Table 1 list the corresponding results for three properties of maximum-mass hybrid stars ($M_{\rm TOV},R_{\rm TOV},n_{\rm TOV}^{c}$). In Fig. 4, we report in details the posterior PDFs of five stellar configuration parameters ($M_{\rm TOV},R_{\rm TOV},n_{\rm TOV}^{c},M^{\rm core}_{\rm TOV}/M_{\rm TOV},R^{\rm core}_{\rm TOV}/R_{\rm TOV}$). See also Appendix B for similar predictions for $1.4\,{\rm M}_{\odot}$ and $2.0\,{\rm M}_{\odot}$ hybrid stars.

It is commonly accepted that the inclusion of extra degrees of freedom (such as deconfined quarks discussed in the present work) tends to soften the NS matter EOS. However, as long as the high-density phase is sufficiently stiff i.e. with a large enough sound speed, through the introduction of a hadron-quark phase transition the parameter space for NS EOSs could be extended significantly (Miao et al., 2020) (see e.g. colored solid curves and pink shaded regions in Fig. 3), especially for soft hadronic EOSs such as QMF. In fact, the QMF hadronic EOS is quite soft due to the inclusion of the quark interaction for describing the inner structure of a nucleon in vacuum; see details in a recent review (Li et al., 2020). As a result, $M_{\rm TOV}$ of hybrid NSs with a quark core is larger/smaller (thus the corresponding central density is relatively smaller/larger) than that of normal NSs based on the soft QMF/stiff DD2 hadronic EOS, falling in a relatively robust range of $\sim 2.10-2.85\,{\rm M}_{\odot}$ (90% confidence level). The probability distribution of quark matter for the maximum-mass star peaks within a sphere of $R^{\rm core}_{\rm TOV}/R_{\rm TOV}\approx 0.90$ with its mass fraction larger than $90\%$, corresponding to the scenario of a big quark core.

From Table 1 one can see that the joint analysis of the GW170817 and NICER data prefers the hadron-quark phase transition taking place at not-too-high densities ($n_{\rm trans}/n_{0}\approx 2.04$ for QMF and $1.61$ for DD2, respectively). Such a low threshold is consistent with previous studies (Al-Mamun et al., 2021; Xie & Li, 2021), while the sound speed squared $c_{\rm QM}$ above the transition should be larger than $0.53$ ($0.51$) for QMF (DD2), indicating a rather stiff quark EOS generically. The current statistical analysis is also consistent with previous model calculations which confirm that to ensure a large mass for hybrid stars above two solar masses, the core sound speed needs to be necessarily large (Xia et al., 2019), even when the assumption of a first-order phase transition is relaxed (Drischler et al., 2020).

We now turn to examine the nature and properties of several merger events based on our results of the hybrid star maximum mass. Even though there is a wide distribution of $M_{\rm TOV}$ from our results, the peak value $\sim 2.4\,{\rm M}_{\odot}$ suggests that the possibility that the secondary component of GW190814 with a mass $M\sim 2.6\,{\rm M}_{\odot}$ (Abbott et al., 2020b) could be a spinning supermassive NS is not ruled out, even if one takes into account the limitation of the intrinsic r-mode instability (Zhou et al., 2020). Within this scenario, however, the supramassive NS needs to carry a low dipolar magnetic field so that it has not been spun down yet at the time of merger. A more likely situation is that this secondary component is a black hole, since the typical spin-down timescale is quote short (Sarin et al., 2020). Considering the maximal spin observed in Galactic NSs, Abbott et al. (2019a) inferred that the 90% credible intervals for the component masses of the GW170817 event lie between $1.16$ and $1.60\,{\rm M}_{\odot}$ (with the most probable total mass around $M_{\rm tot}=2.73\,{\rm M}_{\odot}$), and Abbott et al. (2020a) reported the corresponding component masses ranging from $1.46$ to $1.87\,{\rm M}_{\odot}$ (with the most probable total mass around $M_{\rm tot}=3.4\,{\rm M}_{\odot}$) for the GW190425 event. Our results on the maximum mass of NSs (up to $\sim 2.85\,{\rm M}_{\odot}$ to the 90% confidence level) suggest that the remnant of GW170817 could be a massive rotating NS (Ai et al., 2020). In comparison, the GW190425 event may be more compatible with a BH remnant, especially given that it has no detection of a luminous counterpart, and the inferred total mass is $0.45\,{\rm M}_{\odot}$ above the most massive observed binary NS system (PSR J1913+1102). Further studies on the stellar evolution are warranted to clarify these important issues. See Appendix C for our results using the GW190425 data.

The scenario that GW170817 left behind a massive NS has implications in interpreting the EM counterparts of binary NS mergers. First, since GW170817 was associated with the short GRB 170817A, it suggests that the launch of a GRB jet does not require the formation of a black hole to begin with. The argument that the existence of a short GRB necessitates the formation of a black hole, and hence, a small $M_{\rm TOV}$ (e.g. Margalit & Metzger, 2017; Rezzolla et al., 2018; Ruiz et al., 2018), is flawed. Rather, it suggests that there might exist an underlying long-lived engine that powers the late EM counterpart, including the kilonova (Li et al., 2018; Yu et al., 2018) and late X-ray emission (Piro et al., 2019; Troja et al., 2020). The 1.7 s delay between GRB 170817A and GW170817 was not due to the delay time to form a black hole, but rather due to the time delay for a promptly launched jet to reach the energy dissipation radius (Zhang et al., 2018; Zhang, 2019). The lack of significant energy injection in the GW170817/GRB 170817A remnant required that the dipolar magnetic field strength is below $10^{13}$ G (Ai et al., 2018, 2020). Alternatively, most of the spin energy of the post-merger massive NS could have been carried away via secular gravitational waves (Fan et al., 2013; Gao et al., 2016) or deconfined quarks (Drago et al., 2016; Li et al., 2016; Sarin et al., 2020). For a review on the post-merger remnant of binary NSs and GW170817 in particular, see Bernuzzi (2020); Sarin & Lasky (2020). The existence of a long-lived NS engine for short GRBs is consistent with the interpretation of the X-ray plateau commonly existing in short GRB afterglows (Rowlinson et al., 2010, 2013; Lü et al., 2015) and the X-ray transients CDF-S XT2 and XT1 (Xue et al., 2019; Sun et al., 2019), which pointed towards a $M_{\rm TOV}$ in the range of $\sim 2.4\,{\rm M}_{\odot}$ (Lasky et al., 2014; Gao et al., 2016; Li et al., 2016; Radice et al., 2018). According to this picture, a large fraction of binary NS mergers will be accompanied by a bright X-ray emission component (Zhang, 2013). Future regular monitoring of X-ray counterparts from binary NS merger GW sources with e.g. Einstein Probe (Yuan et al., 2018), will test the scenario of a relatively large $M_{\rm TOV}$ suggested here.