A New Approach to the Mass and Radius of Neutron Stars with Supernova Neutrinos

The mass and radius of neutron stars provide clues to revealing the properties of high-density nuclear matter that are characterized by the equation of state (EOS; Lattimer & Prakash 2016; Özel & Freire 2016; Oertel et al. 2017; Li et al. 2019). Recently, constraints on the mass–radius relation have been obtained from observations of X-ray bursts from neutron star binaries with low-mass companions (Özel et al. 2010; Steiner et al. 2010) and the gravitational wave from a binary neutron star merger, GW 170817 (Abbott et al. 2018). As a result, EOSs that predict large (≳13.5 km) radii of neutron stars are disfavored. Furthermore, very recently, Riley et al. (2019) and Miller et al. (2019) have obtained new insights for the mass and radius of PSR J0030+0451 using NICER, which is a soft X-ray telescope launched in 2017 June.

As is well known, neutron stars are formed in core-collapse supernova explosions. In this process, a large number of neutrinos are emitted, which were actually detected for SN 1987A (Bionta et al. 1987; Hirata et al. 1987). At that time, mass estimations for the newly born neutron star were attempted by using the neutrino event number, which should correspond to the binding energy of the neutron star emitted as neutrinos (Sato & Suzuki 1987; Lattimer & Yahil 1989). Although there were large uncertainties owing to low statistics for SN 1987A, a statistically significant constraint will be provided by the next Galactic core-collapse supernova (Suwa et al. 2019).

In this paper, we investigate how the supernova neutrino signal depends on not only the mass but also the radius of a neutron star. For this purpose, we perform comprehensive simulations of the neutrino emission from a protoneutron star (PNS), which is a nascent compact remnant of a core-collapse supernova. About half of the neutrinos are emitted within the dynamical phase of a supernova core, and the others originate from the cooling stage of a PNS (Nakazato et al. 2013; Mirizzi et al. 2016; Müller 2019). In the latter phase, a PNS is almost hydrostatic but it is hot and lepton-rich. Since the entropy and lepton number are reduced by neutrino diffusion, a PNS evolves into a cold neutron star, which takes about a few minutes.

The cooling process and neutrino emission of a PNS have been studied for a few decades (Burrows & Lattimer 1986; Suzuki 1994; Pons et al. 1999; Roberts 2012). In particular, the effect of the EOS is one of the central issues in PNS cooling (Sumiyoshi et al. 1995; Pons et al. 1999; Roberts et al. 2012; Camelio et al. 2017; Nakazato & Suzuki 2019), and even the phase transitions that lead to black hole formation have been examined (Keil & Janka 1995; Pons et al. 2001). The impacts of inhomogeneous matter at subnuclear densities have also been discussed and it was reported that the neutrino luminosity is insensitive to inhomogeneous matter and is mainly determined by the high-density EOS (Nakazato et al. 2018). Nevertheless, there have been no systematic studies on the EOS dependence except for our previous paper (Nakazato & Suzuki 2019). In this study, we extend our investigations to further variations of the EOS. Moreover, we present a formulation of the cooling timescale for PNSs, which depends on their mass and radius, in order to compare it with our numerical results.

The numerical methods employed in our simulations of PNS cooling are basically described in our previous paper (Nakazato & Suzuki 2019). Quasi-static evolutions of the PNS are computed by solving the Tolman–Oppenheimer–Volkoff and neutrino transport equations. The neutrino transfer, for which we adopt the multigroup flux-limited diffusion scheme (Suzuki 1994), is responsible for the time evolutions of the entropy and lepton-number profiles. We treat the three species of neutrinos, ν_e, ${\bar{\nu }}_{e}$, and ν_x ($={\nu }_{\mu }={\bar{\nu }}_{\mu }={\nu }_{\tau }={\bar{\nu }}_{\tau }$), that interact with matter. For the neutrino interactions, while most of the rates are taken from Bruenn (1985), we include neutrino pair processes via nucleon bremsstrahlung (Suzuki 1993) and plasmon decay (Kohyama et al. 1986). Further sophistication of the medium modifications of neutrino opacities (Fischer 2016; Camelio et al. 2017), such as the nucleon dispersion-relation dependence of the interaction rates, is deferred to a future work. In this study, we do not consider effects of additional mass accretion (Burrows 1988) and convection (Roberts et al. 2012) for simplicity.

The initial conditions of our PNS cooling are adopted from the numerical results of hydrodynamical simulations as in Nakazato et al. (2013). We compute the core collapse, bounce, and shock propagation for the progenitor models of 15M_⊙ and 40M_⊙ in Woosley & Weaver (1995), utilizing the numerical code of general relativistic neutrino-radiation hydrodynamics in spherical symmetry (Sumiyoshi et al. 2005). In these computations, we employ the Togashi EOS (Togashi et al. 2017). Then, we obtain profiles of the entropy and electron fraction inside the shock wave as functions of the baryon mass coordinate, which are adopted as the initial conditions of our PNS cooling. Note that, although spherically symmetric models do not yield explosions in general, the region behind the stalled shock wave at some point in time can be regarded as a PNS model after shock revival. For the progenitor model of 15M_⊙, snapshots when the shock wave is stalled at the baryon mass coordinates of 1.47M_⊙ and 1.62M_⊙ are adopted and referred to as models 147a and 162a, respectively. For the 40M_⊙ progenitor, snapshots when the shock wave is at 1.62M_⊙ and 1.78M_⊙ are referred to as models 162b and 178b, respectively. While the PNS models 162a and 162b have the same baryon mass, the initial conditions are different.

For our PNS cooling, we employ a series of phenomenological EOSs. For uniform nuclear matter at zero temperature, we express the energy per baryon as

$$\begin{eqnarray}&&w({n}_{b},{Y}_{p})={w}_{0}+\displaystyle \frac{{K}_{0}}{18{n}_{0}^{2}}{\left({n}_{b}-{n}_{0}\right)}^{2}+S({n}_{b}){\left(1-2{Y}_{p}\right)}^{2},\end{eqnarray} \tag{ 1 }$$

with the baryon number density n_b and proton fraction Y_p. Here, the saturation density and saturation energy are set to n₀ = 0.16 fm⁻³ and w₀ = −16 MeV, respectively. The stiffness of symmetric nuclear matter, for which Y_p = 0.5, is characterized by the incompressibility K₀. The density-dependent symmetry energy is written as

$$\begin{eqnarray}&&S({n}_{b})={S}_{0}+\displaystyle \frac{L}{3{n}_{0}}({n}_{b}-{n}_{0})+\displaystyle \frac{1}{{n}_{0}^{2}}\left({S}_{00}-{S}_{0}-\displaystyle \frac{L}{3}\right){\left({n}_{b}-{n}_{0}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

with the coefficients of the symmetry energy S₀ and its density derivative L at the saturation density. Here, the symmetry energy at the density of 2n₀ becomes S₀₀: $S(2{n}_{0})={S}_{00}$. To construct the finite-temperature EOS, we include the contribution of thermal effects evaluated using the ideal Fermi gas model, where the nucleon effective mass is a free parameter. Furthermore, our EOS shares the subsaturation-density region containing inhomogeneous nuclear matter with the Shen EOS (Shen et al. 2011). Note that our EOS is described in detail in our previous paper (Nakazato & Suzuki 2019).

So far, the properties of nuclear matter in the vicinity of the saturation density have been probed in several terrestrial experiments. In particular, for symmetric nuclear matter, experimental data on isoscalar giant monopole resonance suggest that K₀ lies around 230 or 250 MeV (Shlomo et al. 2006; Garg & Colò 2018). On the other hand, for asymmetric nuclear matter, the values of S₀ and L inferred from nuclear binding energies are correlated (Kortelainen et al. 2010; Lattimer & Lim 2013). They have been indicated to be 30 $\mathrm{MeV}\leqslant {S}_{0}\leqslant 32\,\mathrm{MeV}$ and 40 $\mathrm{MeV}\leqslant L\leqslant 60\,\mathrm{MeV}$ if other experimental constraints are also considered (Tews et al. 2017). To approach the symmetry energy at supranuclear densities, observations of neutron stars are advantageous. According to Zhang & Li (2019), the symmetry energy at 2n₀ is constrained to 46.9 ± 10.1 MeV using the data of GW 170817. On the basis of these constraints, we calculate the structure of cold neutron stars using our EOS models and select the models that have a reasonable mass–radius relation. Incidentally, some models are inappropriate owing to the maximum mass of neutron stars being considerably smaller than the observations (Arzoumanian et al. 2018; Cromartie et al. 2020).

In Figure 1, we show the mass–radius relations of neutron stars for the EOS models selected in this study. Here, we examine the values of K₀ = 220, 245, and 270 MeV. For the symmetry energy S(n_b), we employ the models with $({S}_{0},L,{S}_{00})=(30,35,35)$, (30, 35, 40), (30, 35, 45), (30, 35, 55), (31, 50, 40), (31, 50, 45), (31, 50, 55), (32, 65, 45), (32, 65, 55), and (33, 80, 55), in units of MeV. Thus, for all possible combinations of K₀ and (S₀, L, S₀₀), we consider 30 models in total for the zero-temperature EOS. The neutron stars constructed in our EOS models have different radii ranging from 11 to 13 km for typical masses. Nevertheless, the gravitational mass is insensitive to the EOS and is about 1.33M_⊙, 1.44M_⊙, and 1.57M_⊙ for neutron stars with baryon masses of 1.47M_⊙, 1.62M_⊙, and 1.78M_⊙, respectively (Figure 1).

Zoom In Zoom Out Reset image size

The effective mass of nucleons is a key ingredient in characterizing the finite-temperature EOS. It was reported that the likelihood of supernova explosion is increased for EOSs with a large effective mass (Schneider et al. 2019). In our EOS, the effective mass in units of nucleon rest mass is denoted by u and we assume that neutrons and protons have the same value of u. Furthermore, we do not deal with the density and temperature dependences of the effective mass. Nevertheless, our EOS is advantageous for investigating the thermal contribution due to the effective mass because the dependencies on the effective mass are separated from the variation of the zero-temperature EOS. In this study, we examine the cases with effective masses of u = 1, 0.75, and 0.5 for each zero-temperature EOS. Thus, we employ 90 models of the finite-temperature EOS to perform the simulations of PNS cooling.

The light curves of ${\bar{\nu }}_{e}$, which are the time evolutions of neutrino luminosity, are shown in Figure 2 for the models with effective masses of u = 1 and u = 0.5. The PNS models with a larger mass have a longer timescale of neutrino emission. This trend is qualitatively consistent with a previous study (Camelio et al. 2017). The resultant neutrino light curves become similar for the models of 162a and 162b when the same EOS is adopted. This fact means that the neutrino signal in the late phase is insensitive to the initial profiles of the entropy and electron fraction (Suwa et al. 2019). As reported in our previous paper (Nakazato & Suzuki 2019), the decay timescale of the neutrino luminosity is longer for the models with larger effective masses because the thermal energy stored in the PNS is larger.

Zoom In Zoom Out Reset image size

As shown in the neutrino light curves obtained from our computations (Figure 2), the PNS cooling is divided into three phases. First, the neutrino luminosity decreases steeply due to the rapid contraction of the PNS because the neutrino luminosity is determined by the surface area. Second, since the structure of the PNS becomes almost stationary, the decrease of the neutrino luminosity becomes slower, which is identified as a shallow decay phase in the neutrino light curves. In this phase, the neutrinos trapped inside the PNS leak out from the surface gradually. The neutrino light curve of this phase is sensitive to the EOS (Pons et al. 1999; Nakazato & Suzuki 2019) and we focus on this phase hereafter. Finally, the matter in the PNS becomes neutrinoless β-equilibrium and the neutrino luminosity reduces steeply again.

In this study, the e-folding time of the ${\bar{\nu }}_{e}$ luminosity, ${\tau }_{{\bar{\nu }}_{e}}(t)$, is introduced as

$$\begin{eqnarray}&&{L}_{{\bar{\nu }}_{e}}(t+{\tau }_{{\bar{\nu }}_{e}}(t))=\displaystyle \frac{{L}_{{\bar{\nu }}_{e}}(t)}{e},\end{eqnarray} \tag{ 3 }$$

where ${L}_{{\bar{\nu }}_{e}}$ is the ${\bar{\nu }}_{e}$ luminosity and e is the base of the natural logarithm. Furthermore, we define the cooling timescale of a PNS as the maximum value of ${\tau }_{{\bar{\nu }}_{e}}$:

$$\begin{eqnarray}&&{\tau }_{\mathrm{cool}}\equiv \mathop{\max }\limits_{t}{\tau }_{{\bar{\nu }}_{e}}(t).\end{eqnarray} \tag{ 4 }$$

In Figure 3, τ_cool is plotted as a function of the radius of the cold neutron star for each EOS model. Here, we can confirm that, whereas the PNS cooling is computed for various EOS models with different incompressibilities and symmetry energies, the cooling timescale is significantly correlated with the radius of the cold neutron star, and it is longer for an EOS model with a smaller neutron star radius.

Zoom In Zoom Out Reset image size

As described above, the cooling timescale is longer for a model with a larger neutron star mass and a smaller neutron star radius. This feature can be formulated as below. The Kelvin–Helmholtz timescale is given by ${\tau }_{\mathrm{KH}}\approx | {E}_{g}| /{L}_{* }$, where E_g is the stellar gravitational binding energy and L_* is the stellar luminosity (Kippenhahn & Weigert 1990). Since the luminosity is proportional to the stellar surface area, we can suppose

$$\begin{eqnarray}&&{\tau }_{\mathrm{KH}}\propto \displaystyle \frac{| {E}_{g}| }{{R}^{2}},\end{eqnarray} \tag{ 5 }$$

where R is the stellar radius. In the Newtonian case, it is rewritten as τ_KH ∝ M²/R³ because of $| {E}_{g}| \propto {M}^{2}/R$, where M is the stellar mass. Thus, the timescale is longer for stars with larger mass and smaller radius. However, for extending it to the neutron star with the mass of m and radius of r, we consider the following two relativistic effects. First, a factor of $1/\sqrt{1-2\beta }$ is multiplied due to the time dilation being β = Gm/rc², where c and G are the velocity of light and the gravitational constant, respectively. Second, $| {E}_{g}| \propto {M}^{2}/R$ is replaced by the relativistic binding energy of the neutron star E_b. Therefore, assuming that the cooling timescale of PNSs is evaluated by the Kelvin–Helmholtz timescale, we obtain

$$\begin{eqnarray}&&{\tau }_{\mathrm{cool}}\propto \displaystyle \frac{{E}_{b}}{{r}^{2}\sqrt{1-2\beta }},\end{eqnarray} \tag{ 6 }$$

instead of Equation (5). According to Lattimer & Prakash (2001), the binding energy of neutron stars is approximated for a large class of EOSs as

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{b}}{m}=\displaystyle \frac{0.6\beta }{1-0.5\beta }.\end{eqnarray} \tag{ 7 }$$

Substituting Equation (7) into Equation (6), we can express the cooling timescale as

$$\begin{eqnarray}&&{\tau }_{\mathrm{cool}}={\tau }^{* }{\left(\displaystyle \frac{m}{1.4{M}_{\odot }}\right)}^{2}{\left(\displaystyle \frac{r}{10\mathrm{km}}\right)}^{-3}\displaystyle \frac{1}{(1-0.5\beta )\sqrt{1-2\beta }},\end{eqnarray} \tag{ 8 }$$

with a coefficient τ*. Actually, as shown in Figure 4, we find that the cooling timescale of PNSs obtained from our simulations can be faithfully described by Equation (8) with τ* = 37.0 s, 35.2 s, and 33.7 s in the cases with effective masses of u = 1, 0.75, and 0.5, respectively.

Zoom In Zoom Out Reset image size

If the cooling timescale is measured by the future detection of supernova neutrinos, it will be useful to probe the mass and radius of a newly formed neutron star. For example, in our previous paper (Nakazato & Suzuki 2019), we have computed the PNS cooling of the model with a baryon mass of 1.47M_⊙ for the EOS model proposed by Lattimer & Swesty (1991) with an incompressibility of 220 MeV (LS220 EOS), and we have obtained τ_cool = 21.8 s. The relation between m and r given by Equation (8) for this case is shown in Figure 5 with the mass–radius relation of neutron stars for the LS220 EOS. The intersection of them is consistent with the mass and radius of this model (shown by the red point). Furthermore, this model has E_b = 0.136M_⊙ and the relation between m and r given by Equation (7) is also shown in Figure 5. Since E_b corresponds to the total emission energy of supernova neutrinos, both of τ_cool and E_b are measurements of supernova neutrinos. Thus, by combining Equations (7) and (8), we may be able to obtain constraints on the mass and radius of the neutron star from the neutrino observation.

Zoom In Zoom Out Reset image size

In practice, some uncertainties should be considered to estimate the mass and radius of a neutron star applying our method. In this study, the effective mass of nucleons, which also affects neutrino opacities in the medium, is treated as an unknown parameter. As seen in the difference between the models 162a and 162b (Figure 3), the cooling timescale has a slight dependence on the initial condition for the case with the large effective mass. Note that the initial condition of the model 162b has a higher entropy than that of the model 162a. While the entropy is originally generated by the shock propagation in the supernova core, the convection smooths the entropy gradient especially in the early phase (Roberts et al. 2012). Therefore, the convection would affect the cooling timescale as well as the initial entropy. While we define τ_cool with the luminosity of ${\bar{\nu }}_{e}$, E_b is the total emission energy of all neutrino flavors. Since energy equipartition among different flavors may not be achieved (Müller 2019), different types of detectors are needed to evaluate the emission energy for each flavor (Laha & Beacom 2014; Li et al. 2018; Nikrant et al. 2018). In any case, statistical uncertainties would be included when the supernova neutrinos are actually detected. Then, for reducing uncertainties, other constraints on the mass and radius of neutron stars, such as a mass–radius relation predicted by nuclear physics, are desired to be available. Therefore, various approaches are certainly worth investigation.

In this paper, we have carried out a comprehensive simulation study of PNS cooling using 90 EOS models with different incompressibilities, symmetry energies, and nucleon effective masses. We have found that if the PNS mass is fixed, the cooling timescale depends on the radius of the cold neutron star in the final state and the nucleon effective mass, which is introduced to characterize the thermal properties of the EOS. Furthermore, we have presented a theoretical expression of the cooling timescale that describes our numerical results faithfully.

The findings in this study suggest that, as well as the total energy of emitted neutrinos, the decay timescale will be useful to probe the mass and radius of a newly formed neutron star. Actually, in our results for the PNS cooling timescale, the dependence on the zero-temperature EOS is encapsulated in a single parameter, the radius of the cold neutron star r, while we have examined various EOS models. In general, since the integration and slope of a neutrino light curve should have different dependences on m and r, they can provide complementary information, as demonstrated in Figure 5. For predicting m and r more accurately than this study, some sources of uncertainties, such as nucleon effective masses, neutrino opacities, and convection, need to be considered. However, this paper will provide an important and reliable basis for future work because Equation (8) does not depend on details of the cooling model.

The authors are grateful to Ken'ichi Sugiura, Kohsuke Sumiyoshi, Yudai Suwa, and Shoichi Yamada for valuable comments. In this work, numerical computations were partially performed on the supercomputers at Research Center for Nuclear Physics (RCNP) in Osaka University. This work was partially supported by JSPS KAKENHI grant Nos. JP26104006, JP17H05203, JP19H05802, and JP19H05811.