Cooling of Small and Massive Hyperonic Stars

Our models are based on two new parametrizations of the FSU2 RMF model (Chen & Piekarewicz 2014). The Lagrangian density of the theory reads (Serot & Walecka 1986, 1997; Glendenning 2000; Chen & Piekarewicz 2014)

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \sum _{b}{{ \mathcal L }}_{b}+{{ \mathcal L }}_{m}+\displaystyle \sum _{l=e,\mu }{{ \mathcal L }}_{l},\end{eqnarray} \tag{ 1 }$$

with the baryon (b), meson (m = σ, ω, ρ and ϕ), and lepton (l = e, μ) Lagrangians given by

$$\begin{eqnarray}{{ \mathcal L }}_{b} & = & {\bar{{\rm{\Psi }}}}_{b}(i{\gamma }_{\mu }{\partial }^{\mu }-{m}_{b}+{g}_{\sigma b}\sigma -{g}_{\omega b}{\gamma }_{\mu }{\omega }^{\mu }-{g}_{\phi b}{\gamma }_{\mu }{\phi }^{\mu }\\ & & -\,{g}_{\rho b}{\gamma }_{\mu }{{\boldsymbol{I}}}_{b}{\boldsymbol{\rho }}{}^{\mu }){{\rm{\Psi }}}_{b},\\ {{ \mathcal L }}_{m} & = & \displaystyle \frac{1}{2}{\partial }_{\mu }\sigma {\partial }^{\mu }\sigma -\displaystyle \frac{1}{2}{m}_{\sigma }^{2}{\sigma }^{2}-\displaystyle \frac{\kappa }{3!}{({g}_{\sigma N}\sigma )}^{3}-\displaystyle \frac{\lambda }{4!}{({g}_{\sigma N}\sigma )}^{4}\\ & & -\,\displaystyle \frac{1}{4}{{\rm{\Omega }}}^{\mu \nu }{{\rm{\Omega }}}_{\mu \nu }+\displaystyle \frac{1}{2}{m}_{\omega }^{2}{\omega }_{\mu }{\omega }^{\mu }+\displaystyle \frac{\zeta }{4!}{({g}_{\omega N}{\omega }_{\mu }{\omega }^{\mu })}^{4}\\ & & -\,\displaystyle \frac{1}{4}{{\boldsymbol{R}}}^{\mu \nu }{{\boldsymbol{R}}}_{\mu \nu }+\displaystyle \frac{1}{2}{m}_{\rho }^{2}{{\boldsymbol{\rho }}}_{\mu }{\boldsymbol{\rho }}{}^{\mu }+{{\rm{\Lambda }}}_{\omega }{g}_{\rho N}^{2}{{\boldsymbol{\rho }}}_{\mu }{\boldsymbol{\rho }}{}^{\mu }{g}_{\omega N}^{2}{\omega }_{\mu }{\omega }^{\mu }\\ & & -\,\displaystyle \frac{1}{4}{P}^{\mu \nu }{P}_{\mu \nu }+\displaystyle \frac{1}{2}{m}_{\phi }^{2}{\phi }_{\mu }{\phi }^{\mu },\\ {{ \mathcal L }}_{l} & = & {\bar{\psi }}_{l}(i{\gamma }_{\mu }{\partial }^{\mu }-{m}_{l}){\psi }_{l},\end{eqnarray} \tag{ 2 }$$

where Ψ_b and ψ_l are the baryon and lepton Dirac fields, respectively. The field strength tensors are ${{\rm{\Omega }}}_{\mu \nu }={\partial }_{\mu }{\omega }_{\nu }\,-{\partial }_{\nu }{\omega }_{\mu }$, ${{\boldsymbol{R}}}_{\mu \nu }={\partial }_{\mu }{{\boldsymbol{\rho }}}_{\nu }-{\partial }_{\nu }{{\boldsymbol{\rho }}}_{\mu }$, and ${P}_{\mu \nu }={\partial }_{\mu }{\phi }_{\nu }-{\partial }_{\nu }{\phi }_{\mu }$. The strong interaction couplings of a meson to a certain baryon are denoted by g and the baryon, meson and lepton masses by m. The vector ${{\boldsymbol{I}}}_{b}$ stands for the isospin operator.

The Lagrangian density (2) incorporates scalar and vector meson self-interactions, as well as a mixed quartic vector meson interaction. The nonlinear meson interactions are important to describe nuclear matter and finite nuclei. The scalar self-interactions with coupling constants κ and λ (Boguta & Bodmer 1977) are responsible for softening the EoS of symmetric nuclear matter around saturation density and allow one to obtain a realistic value for the compression modulus of nuclear matter (Boguta & Bodmer 1977; Boguta & Stoecker 1983). The quartic isoscalar–vector self-interaction (with coupling ζ) softens the EoS at high densities (Mueller & Serot 1996), while the mixed quartic isovector–vector interaction (with coupling Λ_ω) is introduced (Horowitz & Piekarewicz 2001a, 2001b) to modify the density dependence of the nuclear symmetry energy. From the Lagrangian density (2), one derives in the standard way the equations of motion for the baryon and meson fields, which are solved in the mean‐field approximation (Serot & Walecka 1986).

To compute the structure of neutron stars, one needs the EoS of matter over a wide range of densities. The core of a neutron star harbours chemically equilibrated (β‐stable), globally neutral, charged matter. We compute the EoS for the core of the star using the Lagrangian of Equations (1)–(2). As in Tolos et al. (2017a, 2017b), for the crust of the star we employ the crustal EoS of Sharma et al. (2015), which has been obtained from microscopic calculations. Once the EoS is known, the solution of the Tolman–Oppenheimer–Volkoff (TOV) equations (Oppenheimer & Volkoff 1939) provides the mass M and radius R of the neutron star.

We first consider the nucleonic RMF parametrization FSU2 of Chen & Piekarewicz (2014). The parameters of the model were fitted by requiring an accurate description of the binding energies, charge radii and monopole response of atomic nuclei across the periodic table, and, at the same time, a limiting mass of ≃2 M_⊙ in neutron stars. The resulting FSU2 parameter set (Chen & Piekarewicz 2014) (the parameters and saturation properties of the models can be found in Table 1) provides a stiff enough EoS in the region of large baryon densities and, as a consequence, reproduces heavy neutron stars. Figure 1 shows the pressure of β‐stable neutron star matter as a function of the baryon density for the different models considered in the present work. It can be seen that the pressure from FSU2—and also the pressure from all the other models that we use—passes through the region constrained by the study of flow data in experiments of energetic HICs (Danielewicz et al. 2002). The predicted mass–radius (M–R) relation of neutron stars is represented in the upper panel of Figure 2. The lower panel of Figure 2 shows the mass of the star as a function of the central density. FSU2 reaches a maximum mass of 2.07 M_⊙ and, thus, accommodates the observed masses of 2 M_⊙ in pulsars PSR J1614–2230 (Demorest et al. 2010) and PSR J0348+0432 (Antoniadis et al. 2013). The radius of the star at maximum mass is of 12.1 km in FSU2 (see Table 2). For a neutron star with a mass of 1.4 M_⊙, FSU2 predicts a stellar radius of 14.1 km.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Parameters of the Models FSU2 (Chen & Piekarewicz 2014), and FSU2R and FSU2H (Tolos et al. 2017a, 2017b) (Here, we use the slightly updated version of the FSU2R and FSU2H parameters given in Tolos et al. 2017b)

Models	FSU2	FSU2R	FSU2H
mσ (MeV)	497.479	497.479	497.479
mω (MeV)	782.500	782.500	782.500
mρ (MeV)	763.000	763.000	763.000
${g}_{\sigma N}^{2}$	108.0943	107.5751	102.7200
${g}_{\omega N}^{2}$	183.7893	182.3949	169.5315
${g}_{\rho N}^{2}$	80.4656	206.4260	197.2692
κ	3.0029	3.0911	4.0014
λ	−0.000533	−0.001680	−0.013298
ζ	0.0256	0.024	0.008
Λω	0.000823	0.045	0.045
n0 (fm−3)	0.1505	0.1505	0.1505
E/A (MeV)	−16.28	−16.28	−16.28
K (MeV)	238.0	238.0	238.0
Esym(n0) (MeV)	37.6	30.7	30.5
L (MeV)	112.8	46.9	44.5

Note. The mass of the nucleon is m = 939 MeV. Also reported are the values in nuclear matter at the saturation density n₀ for the energy per particle (E/A), incompressibility (K), symmetry energy (E_sym), and slope parameter of the symmetry energy (L).

Download table as:  ASCIITypeset image

Recent progress in astrophysical estimates of neutron star radii—see for example Ozel & Freire (2016) for a review—suggests that radii may be smaller than previously thought. Particularly, the extractions of radii from quiescent low‐mass X‐ray binaries (QLMXBs) and X‐ray bursters are pointing toward values no larger than approximately 13 km (Guillot et al. 2013; Guver & Ozel 2013; Guillot & Rutledge 2014; Heinke et al. 2014; Lattimer & Steiner 2014a, 2014b; Lattimer & Prakash 2016; Nättilä et al. 2016; Ozel et al. 2016; Ozel & Freire 2016). The first observation of a binary neutron‐star merger by the LIGO and Virgo collaborations also appears to indicate that neutron stars cannot have excessively large radii (Abbott et al. 2017). Indeed, very recently, by comparing the tidal deformabilities extracted from the gravitational‐wave signal of this merger with the predictions using different types of EoS's, relatively small upper limits for the radius of a 1.4 M_⊙ neutron star have been deduced such as of R < 13.76 km in (Fattoyev et al. 2018), of R < 13.6 km in (Annala et al. 2018; Krastev & Li 2018), and of R < 13.45 km in (Most et al. 2018). It is worth mentioning that the advent of accurate data on neutron star radii should allow one to probe the EoS of neutron‐rich matter in a complementary way to the heavy masses. This is due to the fact that, while the maximum stellar mass depends on the high‐density sector of the EoS, the strongest impact on the radius of the star comes from the EoS in the low‐to‐medium density region of 1–2 times the saturation density n₀ (Lattimer & Prakash 2007; Ozel & Freire 2016). Given that the nuclear symmetry energy, E_sym(n), governs the departure of the energy of neutron matter from symmetric matter, it means that data on neutron star radii pin down the density dependence of the symmetry energy around n₀, and the slope parameter L, defined as $L=3{n}_{0}{\left.\tfrac{\partial {E}_{\mathrm{sym}}(n)}{\partial n}\right|}_{n={n}_{0}}$, which is intimately related with the isospin properties of atomic nuclei, although its value is still uncertain (Li et al. 2014). Larger L values (stiffer symmetry energy) favor larger radii in neutron stars, whereas smaller L values (softer symmetry energy) favor smaller radii. Therefore, astrophysical evidence of small neutron star radii is consistent with a nuclear symmetry energy that is not overly stiff around saturation density.

To account for the possible existence of massive stars with small radii in our theory, yet without compromising the agreement with constraints from the properties of atomic nuclei and from HICs, in Tolos et al. (2017a, 2017b) we developed the FSU2R parametrization, which produces a soft symmetry energy and a soft pressure of neutron matter for densities n ≲ 2n₀. This can be seen in Figure 1 by comparing the EoS's of FSU2R (solid black line) and FSU2 (dashed blue line). For a given neutron star, in the FSU2R EoS up to densities ∼2n₀, there is less pressure to balance gravity, thereby leading to increased compactness of the star and smaller stellar radius. In the high‐density sector of the EoS, the FSU2R and FSU2 EoS's are close to each other (see Figure 1), thus FSU2R also reproduces heavy neutron stars, as can be seen in Figure 2. Smaller radii within the range of 11.5–13 km are obtained in FSU2R for neutron stars between maximum mass and M = 1.4 M_⊙ (see Figure 2 and Table 2), owing to a softer symmetry energy, while reproducing the properties of nuclear matter and nuclei. Our conclusions are in keeping with the results of recent studies with RMF models with a soft symmetry energy (Chen & Piekarewicz 2015a, 2015b). Indeed, FSU2R predicts E_sym(n₀) = 30.7 MeV and L = 46.9 MeV (Tolos et al. 2017b), which are in good accord with the limits of recent determinations (Lattimer & Lim 2013; Li & Han 2013; Hagen et al. 2015; Roca-Maza et al. 2015; Birkhan et al. 2017; Oertel et al. 2017).

As the FSU2 and FSU2R parametrizations assume nucleonic (non‐strange) stellar cores, we shall often use the notation FSU2(nuc) and FSU2R(nuc) to refer to these models. We also analyze in the present work the consequences of the appearance of hyperons inside neutron star cores. As described in Tolos et al. (2017a, 2017b), the couplings of the hyperons to the vector mesons are related to the nucleon couplings by assuming SU(3)‐flavor symmetry, the vector dominance model and ideal mixing for the physical ω and ϕ mesons, as, e.g., employed in recent works (Schaffner & Mishustin 1996; Weissenborn et al. 2012; Colucci & Sedrakian 2013; Miyatsu et al. 2013; Banik et al. 2014). The values of the hyperon couplings to the scalar σ meson field are determined from the available experimental information on hypernuclei, in particular by fitting to the optical potential of hyperons extracted from the data. Finally, the coupling of the ϕ meson to the Λ baryon is reduced by 20% from its SU(3) value in order to reproduce ${\rm{\Lambda }}{\rm{\Lambda }}$ bond energy data (Ahn et al. 2013).

The EoS of neutron star matter and the M–R relation from the FSU2R model with inclusion of hyperons—dubbed as FSU2R(hyp) model—are plotted, respectively, in Figures 1 and 2. As expected, due to the softening of the high‐density EoS with hyperonic degrees of freedom (compare the FSU2R(hyp) and FSU2R(nuc) EoS's in Figure 1), we obtain a reduction of the maximal neutron star mass below 2 M_⊙ in the FSU2R(hyp) calculation. However, we may readjust the parameters of the nuclear model by stiffening the EoS of isospin‐symmetric matter for densities above twice the saturation density, i.e., the region where hyperons set in, simultaneously preserving the properties of the previous EoS for the densities near saturation, which are important for finite nuclei and for stellar radii. The couplings of the hyperons to the different mesons can be determined as before. The parameters of the new interaction (Tolos et al. 2017a, 2017b), denoted as FSU2H, are displayed in Table 1, along with the predicted symmetry energy E_sym(n₀) at saturation density and its slope L, which are safely within current empirical and theoretical bounds (cf. Figure 4 of Tolos et al. 2017a and Figure 1 of Tolos et al. 2017b). The neutron star calculations with the FSU2H model with allowance for hyperons in the stellar core—FSU2H(hyp) model—successfully fulfill the 2 M_⊙ mass limit with moderate radii for the star (see Figure 2 and Table 2), while the base nuclear model FSU2H still reproduces the properties of nuclear matter and nuclei. In isospin‐symmetric nuclear matter, FSU2H leads to a certain overpressure in the EoS at high densities (Tolos et al. 2017a), but the EoS in neutron matter satisfies the constraints from HICs (Danielewicz et al. 2002), as can be seen in Figure 1.

We also draw in Figure 2 the M–R relation predicted by FSU2H if one neglects hyperons—FSU2H(nuc) model. As expected, as the hyperonic FSU2H(hyp) EoS is softer than the FSU2H(nuc) EoS after hyperons appear (see Figure 1), the neutron star calculations with FSU2H(nuc) lead to a higher maximum mass than FSU2H(hyp). In the next section, comparisons between cooling calculations performed with FSU2H(hyp) and FSU2H(nuc) will be used to discuss the influence of the occurrence of hyperons on neutron star cooling.

We first consider the thermal evolution of neutron stars without taking into account any sort of pairing (neither in the core nor in the crust). This is done in Figures 3–7 with the ultimate goal of determining how the different models for the EoS describe the thermal behavior of neutron stars as "benchmark testing" before the inclusion of pairing. However, these results should not be regarded as advocating for the absence of superfluidity/superconductivity in the star, whose effect will be considered in Section 3.2. In Figures 3–7, we depict several cooling curves for low, medium and high neutron star masses for each model, to investigate the cooling behavior for different central density regimes. The figures also display the observed surface temperature versus the age of a set of prominent neutron stars, including that of the remnant in Cassiopeia A (Mereghetti et al. 1996; Possenti et al. 1996; Halpern & Wang 1997; Zavlin et al. 1999; Pavlov et al. 2001, 2002; Gotthelf et al. 2002; Pons et al. 2002; Burwitz et al. 2003; Kaplan et al. 2003; McGowan et al. 2003, 2004, 2006; Zavlin 2007; SafiHarb & Kumar 2008; Zavlin 2009; Ho et al. 2015; Klochkov et al. 2015; Beloin et al. 2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The cooling behavior from each microscopic model is summarized in Table 3 for low‐mass to high‐mass neutron stars. When DU reactions of neutrino production are allowed according to the model, they lead to an enhanced cooling of the star. We summarize here the DU processes that may take place inside the neutron star:

$$\begin{eqnarray*}n & \to & p+{e}^{-}+\bar{\nu },\\ {\rm{\Lambda }} & \to & p+{e}^{-}+\bar{\nu },\\ {{\rm{\Sigma }}}^{-} & \to & n+{e}^{-}+\bar{\nu },\\ {{\rm{\Sigma }}}^{-} & \to & {\rm{\Lambda }}+{e}^{-}+\bar{\nu },\\ {{\rm{\Sigma }}}^{-} & \to & {{\rm{\Sigma }}}^{0}+{e}^{-}+\bar{\nu },\\ {{\rm{\Xi }}}^{-} & \to & {\rm{\Lambda }}+{e}^{-}+\bar{\nu },\\ {{\rm{\Xi }}}^{-} & \to & {{\rm{\Sigma }}}^{0}+{e}^{-}+\bar{\nu },\\ {{\rm{\Xi }}}^{0} & \to & {{\rm{\Sigma }}}^{+}+{e}^{-}+\bar{\nu },\\ {{\rm{\Xi }}}^{-} & \to & {{\rm{\Xi }}}^{0}+{e}^{-}+\bar{\nu }.\end{eqnarray*}$$

We stress that the presence of the particles is not enough for such processes to take place, and one must also account for momentum conservation, as previously discussed. We also recall that all inverse reactions also take place as to maintain (on average) chemical equilibrium. In Table 3, we indicate the density threshold for the nucleonic DU process and, when applicable, the threshold for the hyperonic DU processes, determined by the Λ particle as the first hyperon to appear. We also show the central density n_c for three selected neutron star masses, as well as whether such stars exhibit slow or fast cooling scenarios. The information supplied in Table 3 will help to understand the thermal evolution of neutron stars, displayed in Figures 3–7 for the different models and discussed in the following:

FSU2 (Figure 3). All neutron stars whose microscopic composition is described by this model allow for pervasive DU process, even for low-mass stars with low central densities of around 2n₀, leading to a cooling which is too fast in comparison with the observed data.

FSU2R(nuc) (Figure 4). For this model, most stars (up to M ∼ 1.90 M_⊙) exhibit slow cooling. This model could, in principle, explain most of the observed data. However, it would mean that all colder stars with T ≲ 10⁶ K have masses higher than 1.9 M_⊙, which seems unlikely.

FSU2H(nuc) (Figure 5). For this model, the DU process is absent in all stars studied, from low-mass stars with lower central densities to high-mass stars with higher central densities. Therefore, all stars exhibit slow cooling, which allows for agreement with some of the observed data, but fails to explain most of the observed cold stars.

FSU2R(hyp) (Figure 6). The results of this model exhibit a relatively similar pattern to that of the FSU2R(nuc) model, with the exception that the highest mass possible in this model is 1.765 M_⊙, and that only stars with masses similar to that one exhibit fast cooling.

FSU2H(hyp) (Figure 7). This particular model shows the most promising results, with stars with lower masses exhibiting slow cooling (without DU) and higher masses displaying fast cooling (with DU), while intermediate masses show a behavior in between these extremes. This model explains most of the observed data, including a reasonable agreement with Cas A, without the need of resorting to extensive pairing. Nevertheless, in this case, a large portion of the observational data corresponds to stars with masses between the restricted range of 1.85 and 1.90 M_⊙. This issue becomes less problematic when nucleonic pairing is included, as will be shown in the following subsection.

The different cooling pattern for each model can be understood from the microscopic differences between the models. We start by analyzing the case of low-mass stars of 1.4 M_⊙ with central densities around twice nuclear saturation density. The different cooling pattern of a 1.4 M_⊙ star observed for FSU2 (fast cooling) and FSU2R(nuc) (slow cooling) is mainly due to the different density dependence of the symmetry energy around saturation in these models and, hence, to the different value of the symmetry energy slope parameter (L), which can be read in Table 1 (L ≈ 113 MeV in FSU2 and L ≈ 47 MeV in FSU2R). The larger the value of L is, the more protons are produced, thus the DU process appears at lower densities, making the cooling of the star more efficient. Indeed, even if a M = 1.4 M_⊙ star in the FSU2 model has a smaller central density than in the FSU2R(nuc) model, the cooling is faster in FSU2 as the DU threshold for FSU2 is at a much smaller density, as seen in Table 3. This conclusion is corroborated by the cooling behavior of FSU2R(nuc) and FSU2H(nuc), which show a similar qualitatively slow cooling for neutron stars with M ∼ 1.4 M_⊙, as both models have an alike L value (cf. Table 1). Therefore, one finds that in low-mass stars, where the central density does not go much above 2n₀, large stellar radii (stiff nuclear symmetry energy near n₀) are associated with fast cooling, whereas small stellar radii (soft nuclear symmetry energy near n₀) imply slow cooling, as also seen in Dexheimer et al. (2015) in the framework of the chiral mean-field model.

As for high-mass stars (M ∼ 1.8–2 M_⊙), we find that the different behaviors exhibited by the cooling curves of FSU2 (fast), FSU2R(nuc) (fast), and FSU2H(nuc) (slow) are correlated with the different values of the central densities in these stars, as seen in Table 3. This is understood by considering the fact that the FSU2H(nuc) model produces a stiffer EoS in the region of high densities than the FSU2 and FSU2R(nuc) models (see Figure 1). Therefore, for the same heavy stellar mass, neutron stars obtained within the FSU2H(nuc) model have much lower central densities than in FSU2 and FSU2R(nuc), making the DU process less efficient, thus leading to a slower cooling. For these high stellar masses, the densities reached are much higher than saturation density and, therefore, the slope of the symmetry energy at saturation is not determinant for the cooling behavior.

With regards to the models that include hyperons, as their EoS's are softer than their nucleonic counterparts, they produce stars with higher central density values, which may then overcome the DU threshold, as one can clearly see in Table 3. As a consequence, stars with M ∼ 1.76 M_⊙ that cooled more slowly without hyperons in the FSU2R(nuc) model, change to a faster cooling pattern in the presence of hyperons, as seen for FSU2R(hyp). We also observe in Table 3 that hyperonic models activate the DU process at similar or even lower densities than the models without hyperons. This is mainly due to the fact that the presence of hyperons (mostly Λ particles) reduces the neutron fraction at a given baryon number density and, consequently, the DU constraint ${{\boldsymbol{k}}}_{{Fn}}={{\boldsymbol{k}}}_{{Fp}}+{{\boldsymbol{k}}}_{{Fe}}$, with ${k}_{{Fn}}$, ${k}_{{Fp}}$ and ${k}_{{Fe}}$ being the Fermi momenta of the neutron, proton and electron, respectively, can be fulfilled at a lower density. This is seen in Figure 8, where the particle fractions are shown as functions of the neutron star matter density, for the various interaction models explored in this work. It is clear that the appearance of the Λ hyperons between 0.3 and 0.4 fm⁻³ for the FSU2R(hyp) and FSU2H(hyp) models (middle and lower panels) induces a substantial decrease in the neutron fraction compared to the purely nucleonic models collected in the upper panel. In summary, the presence of hyperons in the cores of medium to heavy mass stars speeds up their cooling pattern. When comparing the cooling curves of the FSU2R(hyp) and FSU2H(hyp) models, we notice that higher stellar masses are needed in the case of FSU2H(hyp) to reach a fast cooling behavior. This is due to the fact that the EoS of the FSU2R(hyp) model is softer and the central densities achieved are larger compared to the ones for FSU2H(hyp), hence the DU threshold is overcome more easily, even if it appears at slightly higher densities. Finally, it may be noticed that the cooling pattern seen in Figure 7 for FSU2H(hyp) is similar to that in Figure 4 for FSU2R(nuc), so that it could argued that it is unnecessary to resort to a hyperonic model such as FSU2H(hyp) to explain most of the observed cooling data. However, purely nucleonic models are forcedly omitting the presence of hyperons, which are physically allowed at intermediate densities when the appropriate chemical equilibrium conditions are fulfilled.

Zoom In Zoom Out Reset image size

We now turn our attention to investigate the neutron superfluidity and proton superconductivity on the cooling of neutron stars. Pairing effects have been considered a key factor for the thermal evolution of neutron stars (Yakovlev et al. 2000; Gnedin et al. 2001; Page et al. 2004; Weber et al. 2007, 2009; Beloin et al. 2018), as well as of extreme importance to find agreement between the theoretical models and the observed data, particularly for Cas A (Ho & Heinke 2009; Heinke & Ho 2010; Page et al. 2011; Shternin et al. 2011; Yang et al. 2011; Blaschke et al. 2012; Negreiros et al. 2013; Ho et al. 2015).

As indicated before, one of the major cooling channels in neutron stars, especially during the first ∼10³ years, is neutrino production reactions. Most of these reactions involve baryons, and chief among those is the DU process, which, if present, is the leading cooling mechanism in neutron stars. The introduction of a superfluidity (conductivity) gap in the energy spectrum of such baryons reduces the reaction rates. One notes that the reduction factor depends on the temperature and its relation to the corresponding superfluidity (conductivity) critical temperature (or the gap), and leads to a sharp drop of neutrino emissivity after the matter temperature drops below the pairing critical temperature. The calculation of the reduction factor for each neutrino emission process is a complicated procedure, which can be obtained by the study of the phase-space of the emission processes (see Yakovlev et al. 2000 for a comprehensive calculation of such factors).

There is, however, a great deal of uncertainties regarding proton and neutron pairing in neutron stars (Yakovlev et al. 2000; Page et al. 2004, 2006), especially at high densities, not to mention the possibility of pairing among hyperons, which is even more uncertain. The inclusion of microscopic calculations for pairing is a challenging task, especially for proton superconductivity. For that reason we have chosen to follow a phenomenological approach like that of Kaminker et al. (2001), which gives us the flexibility to probe different scenarios of proton pairing. Such an approach was used in the work of Shternin et al. (2011) and others (Gusakov et al. 2004, 2005; Yakovlev et al. 2005). Also the work of Beloin et al. (2018) uses a similar phenomenological approach for pairing. Here we explore three different scenarios for proton pairing that reflect three assumptions: (a) shallow proton pairing (limited to densities up to 2–3 n₀), (b) medium proton pairing (extending up to 4 n₀), and (c) deep proton pairing (extending beyond 5 n₀, deep in the inner core). In this manner, we can estimate how extensive the proton pairing must be so as to allow for a good comparison with the observed data. In Figure 9, we show the critical temperature for proton singlet (${}^{1}{S}_{0}$) pairing as a function of density for the three different scenarios. We note that we have used a maximum critical temperature of the order of 10¹⁰ K that is above certain parameterizations used in the literature such as Page et al. (2004) (${T}_{{\rm{c}},{\rm{p}}}^{\max }\lesssim 6\times {10}^{9}$ K) and Kaminker et al. (2001) (${T}_{{\rm{c}},{\rm{p}}}^{\max }\sim 5\,\times {10}^{9}$ K), and closer to Gusakov et al. (2005) (${T}_{{\rm{c}},{\rm{p}}}^{\max }\sim 7\times {10}^{9}$ K) and Beloin et al. (2018) (${T}_{{\rm{c}},{\rm{p}}}^{\max }={7.59}_{-5.81}^{+2.48}\times {10}^{9}$ K). This choice, however, was incidental, as we have used a set of parameters for proton superconductivity that lead to such a high critical temperature. Nevertheless, such a choice makes little difference in the long-term behavior of the cooling of the star, as we will see later.

Zoom In Zoom Out Reset image size

As for neutron pairing, somewhat less uncertain than that of protons (particularly for the neutron singlet pairing in the crust), we chose the standard phenomenological approach: we allow for extensive neutron singlet (${}^{1}{S}_{0}$) pairing at subnuclear densities (crust) and for a limited neutron triplet (${}^{3}{P}_{2}$) pairing in the core, with a maximum critical temperature ∼5 × 10⁸ K, similar to the pairing used to explain the observed temperature of Cas A in Page et al. (2011). The critical temperature of both neutron singlet and triplet pairings as a function of density can be seen in Figure 10.

Zoom In Zoom Out Reset image size

Apart from the suppression of neutrino emission reactions that involve paired particles, a consequence of pairing is the appearance of a new transient neutrino emission process, commonly known as pair breaking-formation (PBF) process. This process can be represented as $b\to b+\nu +\bar{\nu }$, with b denoting paired baryons. It can be also understood as the annihilation of two quasi-baryons with similar anti-parallel momenta into a neutrino pair (Flowers et al. 1976; Voskresensky & Senatorov 1987; Kaminker et al. 1999; Leinson & Perez 2006; Sedrakian et al. 2007; Kolomeitsev & Voskresensky 2008; Steiner & Reddy 2009; Kolomeitsev & Voskresensky 2010). This process, also comprehensively described in Yakovlev et al. (2000), is transient, reaching a maximum near the superfluidity (conductivity) onset and decreasing afterwards.

For the analysis of pairing effects on the cooling of neutron stars, we concentrate our study on the two most relevant cases, the FSU2R(nuc) and the FSU2H(hyp) models (Tolos et al. 2017a, 2017b), which are the ones that best reproduce the observed data on cooling in the previous section. We also show the predictions from the FSU2(nuc) model (Chen & Piekarewicz 2014) for comparison. The main features of the results are as follows:

FSU2 (Figure 11). The original FSU2 model, as discussed before, has a small DU threshold, leading to fast cooling for all stars studied. The inclusion of shallow and medium proton pairing (in addition to neutron pairing, common to all simulations) is ineffective in slowing down the thermal evolution. Although part of the DU (among other processes) is suppressed, there is still a relatively large region at high densities in which the DU takes place, thus leading to a cooling behavior very similar to that of Figure 3. The deep proton pairing, on the other hand, is effective in slowing the cooling, and leads to a slow cooling behavior, as shown in Figure 11. While the first "knee" in the cooling curves, which usually happens between 50 and 150 years, is associated with the core-crust thermal coupling, the second "knee" around ∼5 × 10³ years is linked to the onset of neutron pairing and the subsequent production of neutrinos coming from the PBF process.

Zoom In Zoom Out Reset image size

FSU2R(nuc) (Figure 12). As previously seen, this model exhibits slow cooling only for stars below 1.9 M_⊙ as the DU appears at high densities of n₀ ∼ 0.61 fm⁻³. Shallow to medium proton pairing, plus neutron pairing, allows for a satisfactory agreement with data, with the caveat that all observations with T ≲ 10⁶ K would need to be stars of relatively high mass (above 1.9 M_⊙) within this model, as seen in Figure 12. Deep proton pairing leads to the complete absence of DU processes, thus leading to a slow cooling scenario for neutron stars of all masses.

Zoom In Zoom Out Reset image size

FSU2H(hyp) (Figure 13). This model, which was in good agreement with the observed thermal data without the inclusion of pairing, also benefits from pairing. The inclusion of shallow to medium proton pairing, together with neutron superfluidity, leads to an overall improved agreement with the observed data. This can be seen in Figure 13, where we plot the results of the FSU2H(hyp) model in the case of the medium proton pairing scenario (the red curve in Figure 9). We note that within this model the best fit to Cas A is obtained for a 1.88 M_⊙ star. This is due to the change of slope of the cooling curve by the implementation of the onset of neutron triplet pairing and the PBF process within the FSU2H(hyp) model, as seen in the inset of Figure 13. The results for a 1.4 M_⊙ star are also close to the Cas A observations, but the onset of neutron superfluidity, and thus of the PBF processes, happens a little later (∼10³ years), so that the agreement with the slope of the cooling curve of the Cas A data is not as good. In any case, we recall that the Cas A thermal states are still under debate in the literature (Posselt et al. 2013), as new stars, contrary to old stars, need carbon atmospheres for their description. Therefore, we have to note that, so far, the agreement for Cas A of a certain EoS model and the employed pairing gaps is to be considered with some caution.

Zoom In Zoom Out Reset image size

One should also notice that at the moment, there exist no theoretical calculations of proton pairing in neutron-star matter that incorporate consistently all media effects, which subjects the predicted pairing gaps to uncertainties in size and in extension. We observe in Figure 13 for a 1.88 M_⊙ star that the use of a maximum critical temperature for the proton pairing of 1.41 × 10⁹ K, considerably below the maximum T_c ∼ 10¹⁰ K of our assumed parametrizations in Figure 9, makes little difference for the cooling behavior of the star. Most microscopic models of proton singlet pairing, such as those detailed in Figure 9 of Page et al. (2004), lead to proton superconductivity extending up to proton Fermi momenta of about 0.8–1.25 fm⁻¹. Our phenomenological pairing model is in agreement with several other publications (Gusakov et al. 2004, 2005; Yakovlev et al. 2005; Shternin et al. 2011) that have used a similar approach. We obtain that the medium proton pairing of the cooling calculations of Figure 13 extends up to k_Fp ∼ 1.3 fm⁻¹, as in the pairing model CCDK—but larger than in the other models—of the aforementioned reference (Page et al. 2004). In the shallow pairing scenario (black curve in Figure 9), we obtain a lower value k_Fp ∼ 1.19 fm⁻¹ for the maximum extension of proton pairing. The major effect in the cooling simulations is that neutrino emission processes are less suppressed, thus making the temperature drop sharper and stronger, as can be seen for a 1.88 M_⊙ star in Figure 13. The modification of the cooling curve, however, is not found to be dramatic and does not largely alter our main conclusions.

In summary, while including neutron pairing together with shallow or medium proton pairing is inefficient in slowing the thermal evolution of the stars predicted by the FSU2 model that has a stiff nuclear symmetry energy, deep proton pairing does improve the agreement with data, but still does not explain the cooling of intermediate and low temperature stars. The FSU2R(nuc) and FSU2H(hyp) models, which are characterized by a softer symmetry energy than FSU2 and postpone the onset of DU to higher baryon density, already perform well without consideration of pairing. Nevertheless, we know that a complete neglect of pairing in stellar cooling is not realistic, especially for the well-established neutron ${}^{1}{S}_{0}$ and ${}^{3}{P}_{2}$ pairing (Page et al. 2004). We find that including medium proton pairing, in addition to the neutron pairing, improves the agreement of the cooling curves of the FSU2R(nuc) and FSU2H(hyp) models with data. This stems from the fact that the proposed pairing scheme suppresses DU processes until about 4 times saturation density. This pairing scheme does not block the nucleonic DU processes in massive stars of the FSU2R(nuc) and FSU2H(hyp) models (as it takes place at higher densities; see Table 3), but suppresses other less efficient cooling mechanisms, such as the modified Urca processes that take place in intermediate-mass stars. In other words, we have found that a shallow/medium proton superconductivity does not lead to an over-suppression of the DU processes, a fact that, combined with the underlying properties of the microscopic model, leads to an optimum agreement with observed data, especially concerning the slope for the cooling of Cas A in the case of the FSU2H(hyp) model, as can be seen in the inset of Figure 13. This is our preferred model as the colder stars are described with masses a little below those for the FSU2R(nuc) model, and also because it allows for the presence of hyperons when the chemical equilibrium conditions are fulfilled, being at the same time able of supporting neutron stars with 2 M_⊙.

We note that in our work, hyperon pairing has not been taken into account. This was addressed recently by Raduta et al. (2018), where the cooling of hyperonic stars was studied, although following a somewhat different perspective than ours. Raduta et al. (2018), employing different relativistic density functional models, focused their investigation on the role of hyperonic DU processes on cooling, as well as on the effect of pairing on such processes. They have found that the ${\rm{\Lambda }}\to p+{e}^{-}+\bar{\nu }$ is the dominant hyperonic DU process. Furthermore, they have found that observed data for colder objects can be explained by stars with M ≤ 1.85 M_⊙ whereas the hotter ones are well explained by stars with M ≤ 1.6 M_⊙. Our results are in agreement with their conclusions regarding the dominant hyperonic DU process. However, in our preferred FSU2H(hyp) model, we have found that hotter stars can be described by objects with masses of up to M ≈ 1.85 M_⊙, as opposed to M ≤ 1.6 M_⊙. One must note that for the aforementioned results Raduta et al. (2018) have considered pervasive pairing such that the nucleonic DU process is completely suppressed. The possibility of nucleonic DU was also considered, but did not produce agreement with observed data. In this regard, our studies differ considerably from those of Raduta et al. (2018), as we have not completely excluded the nucleonic DU process and we have not considered hyperonic pairing. For a future investigation, we wish to expand our study to also include hyperonic pairing.