Possibility of rapid neutron star cooling with the realistic equation of state

Whether fast cooling processes occur or not is crucial for the thermal evolution of neutron stars. In particular, the threshold of the direct Urca process, which is one of the fast cooling processes, is determined by the interior proton fraction $Y_p$, or the nuclear symmetry energy. Since recent observations indicate the small radius of neutron stars, a low value is preferred for the symmetry energy. In this study, simulations of neutron star cooling are performed adopting three models for equation of state (EoS): Togashi, Shen, and LS220 EoSs. The Togashi EoS has been recently constructed with realistic nuclear potentials under finite temperature, and found to account for the small radius of neutron stars. As a result, we find that, since the direct Urca process is forbidden, the neutron star cooling is slow with use of the Togashi EoS. This is because symmetry energy of Togashi EoS is lower than those of other EoSs. Hence, in order to account for observed age and surface temperature of isolated neutron stars (INS) with use of the Togashi EoS, other fast cooling processes are needed regardless of the surface composition.

Whether fast cooling processes occur or not is crucial for the thermal evolution of neutron stars. In particular, the threshold of the direct Urca process, which is one of the fast cooling processes, is determined by the interior proton fraction Y p , or the nuclear symmetry energy. Since recent observations indicate the small radius of neutron stars, a low value is preferred for the symmetry energy. In this study, simulations of neutron star cooling are performed adopting three models for equation of state (EoS): Togashi, Shen, and LS220 EoSs. The Togashi EoS has been recently constructed with realistic nuclear potentials under finite temperature, and found to account for the small radius of neutron stars. As a result, we find that, since the direct Urca process is forbidden, the neutron star cooling is slow with use of the Togashi EoS. This is because symmetry energy of Togashi EoS is lower than those of other EoSs. Hence, in order to account for observed age and surface temperature of isolated neutron stars (INS) with use of the Togashi EoS, other fast cooling processes are needed regardless of the surface composition.

Introduction
Neutron stars are formed as a remnant of core-collapse supernovae. Nascent neutron stars are hot and their temperature decreases with their age. Since the emissivity of neutrinos is much larger than that of photons, the neutrino emission determines the thermal evolution of neutron stars until 10 5−6 yr after the birth. Since the ages of observed neutron stars are mostly younger than 10 6−7 yr [1], the observed temperature of neutron stars depends mainly on neutrino emission processes.
The thermal evolutions of neutron stars are described by cooling curves; the relation between the age t and the effective temperature T ∞ eff of neutron stars. So far, many theoretical models of cooling curves have been compared to the observed temperature and age of isolated neutron stars (INS) [2][3][4][5][6][7][8][9][10]. The cooling process of neutron stars depends sensitively on the state of matter, which is described by the equation of state (EoS) beyond the nuclear saturation density. One of the most important ingredients on the thermal evolution of neutron stars is whether fast cooling processes occur or not. Especially, the direct Urca (DU) process is the strongest process in the nucleon cooling processes [11]. The threshold of DU process is determined by the proton fraction Y p , which depends on EoS. Therefore, the EoS dependence of the cooling curves has been studied actively [12,13]. Another significant effect on the neutron star cooling comes from the superfluidity [8,[14][15][16]. While the superfluid critical temperature is a crucial parameter to govern the efficiency of the cooling, it is uncertain especially for the triplet gap (e.g. [17][18][19][20]). Due to the superfluidity, thermodynamical quantities are affected severely through the neutrino loss rates. Therefore, to calculate cooling curves we adopt models of EoS, neutrino loss rates, and superfluid critical temperature, as indispensable physical quantities.
In this paper, we focus on the EoS dependence of the cooling curves. Recently, several constraints on the EoS are considered. For example, the discovery of neutron stars with ≃ 2 M ⊙ [21,22] provides a strong constraint because the EoSs with the maximum mass of neutron stars with <2 M ⊙ are inconsistent with them (see also Ref. [23]). Another examples are from the analysis of low mass X-ray binaries and gravitational wave from neutron star merger (GW170817) [24,25]. They provide constraints on the mass-radius relation of neutron stars and the neutron star radius is indicated to be less than 13 km. In our computations of neutron star cooling, we adopt three models of EoS with different radius of neutron stars. Then, for the EoS model preferred by the observations of neutron star radius, we investigate whether the resultant cooling curves are consistent with the temperature observations of INS. Furthermore, we examine the possibility of rapid cooling process with different surface compositions.
Contents in this paper are as follows: In section 2, we explain basic equations of our cooling simulations, and neutrino processes related to superfluid models. In section 3, we introduce the models of EoS adopted in our computations. In section 4, we show the results of cooling simulation comparing with the observations of INS, and examine the EoS dependence and the effect of superfluidity. In section 5, we discuss the relation between the threshold of DU process and symmetry energy and examine the influence of surface composition. In section 6, we give a concluding remark and present future investigations.

Setup of cooling simulation 2.1. Basic equations
When nuclear burning does not occur, fundamental equations of stellar structure are described as follows [26], with the gravitational constant G and velocity of light c. Here, ρ is the mass energy density, ρ r is the rest mass density, P is the pressure, M tr is the gravitational mass, M r is the rest mass inside the radius r, T is the local temperature, ε ν is the energy loss rate by neutrino emission, C V is the specific heat per baryon, ∇ rad is the radiative gradient, and φ is the gravitational potential in unit mass. From the latest neutron star mass measurements [27], observed minimum and maximum masses are 1.174 ± 0.004 M ⊙ [28] and 2.01 ± 0.04 M ⊙ [22], respectively. This range of the masses observed is consistent with simulations of supernova explosions [29]. In this study, we perform the cooling simulations of neutron stars with the masses from 1.1 M ⊙ to 2.1 M ⊙ . For the age and the surface temperature of INS, we adopt observational data in Ref. [12] for 18 INS and Ref. [30] for Cassiopeia A. Note that, while the temperature variation is observed by Chandra for Cassiopeia A [31] (but see also Ref. [32] ), we do not intend to account for it in this paper. Although the surface composition affects the photon emission rate due to the differences in the opacity κ [33,34], the surface compositions are fixed to be light elements, which are 73% of 1 H, 25% of 4 He, and 2% of 56 Ni. The mass in the envelope is supposed to be 5.0 × 10 −13 M ⊙ , ignoring the accretion. In Section 5.2, however, we vary the surface composition and consider the pure Ni case. Furthermore, ∇ rad depends on mainly κ, and is calculated by using public conductivity codes [35,36]. We adopt the evolution code of a spherically symmetric neutron star [37,38]. As a result, we construct cooling curves by solving the structure equations Eqs. (1) -(6) from the core to the surface of neutron stars.

Neutrino emission processes
So as to evaluate the neutrino energy loss rate ε ν in Eq. (3), we adopt the neutrino emission processes listed in Table 1. The modified Urca and the nucleon pair bremsstrahlung processes always turn on and their loss rate is approximately 10 19−21 erg cm −3 s −1 [39]. While they are slow cooling processes, the DU process is fast cooling process. The emissivity is around 10 27 erg cm −3 s −1 [11]. Since the DU process occurs satisfying the momentum conservation, the threshold of proton fraction Y DU p for the DU process to turn on is written as [11]: where x e = Y e /(Y e + Y µ ), defining Y e and Y µ as electron and muon fractions, respectively. In the very high-density regions, Y DU p ≃ 0.1477 because of Y e ≃ Y µ . When muons are absent, Y DU p is 1/9. Since Y p increases with density in high-density regions, the DU process occurs generally inside massive neutron stars.

Nucleon superfluidity
Nucleon superfluid effect is also important for neutron star cooling. When the local temperature T is lower than the critical temperature T cr , proton-proton pairs and neutron-neutron pairs are formed. Then, the neutrinos and antineutrinos are produced, which is called the Table 1 Nucleon neutrino emission processes adopted in this study. l means electron or muon.

Name
Process Efficiency Modified Urca (neutron branch) nucleon pair breaking formation (PBF) process [8,40,41]. As a consequence both protons and neutrons become 1 S 0 states for low-density regions and neutrons become 3 P 2 state for high-density regions. As a result, neutrino emissivity is suppressed in proportional to exp (−aT cr /T ), where a is 1.76 for 1 S 0 state and 8.40 for 3 P 2 state [42]. However, the relation between the superfluid energy gap (∆ = aT cr ) and the critical temperature is rather uncertain because of the uncertainty of the nuclear force. While many studies about nucleon superfluidity have been done [5,20,43], T cr dependence on the Fermi wave number is still uncertain. We adopt CLS model as the superfluidity of 1 S 0 neutrons, and four models as the other superfluidities: Combination of two 1 S 0 proton superfluid models (AO & CCDK) and two 3 P 2 neutron superfluid models (BEEHS & EEHO) [16] as shown in Fig.1.

Equation of state
As already noted, the EoS is important to calculate the thermal evolution of neutron stars. In this study, so as to include the temperature-dependent structure of neutron stars, we adopt three models of EoS constructed under finite temperatures: Shen, LS220, and Togashi EoSs. While the Shen EoS is based on a relativistic mean field model [44][45][46], the LS220 EoS is based on a Skyrme energy-density functional [47]. These EoSs have been widely used to simulate various astrophysical phenomena [48][49][50][51][52][53][54][55]. The Togashi EoS [56] is constructed with use of realistic two-body potential and phenomenological three-body potential [57][58][59] under the finite temperature. In this section, we describe the properties of these EoSs. Note that, in our calculations, while these EoSs are adopted for the pressure P > 10 30 dyne cm −2 , we adopt BPS EoS [60] in the low-density regions.
In the left panel of Fig. 2, we show total pressures against the baryon density for the EoS models adopted in this study. We find that the Togashi EoS is softer than the other two EoSs for ρ 10 15 g cm Superfluid transition temperature vs. the Fermi wave number of a particle x = n or p [16]. Abbreviations of superfluid models are inserted inside the panel. of heavy ion collisions [61] is also shown in the left panel of Fig. 2 and we can recognize that the EoSs are consistent with the experimental results of heavy ion collisions. Defining u = ρ B /ρ 0 with the baryon mass density ρ B , the energy per nucleon w(u, Y p ) is approximately expanded around the saturation density [62], where w 0 , K, S 0 , and L are, respectively, the energy per nucleon, incompressibility, symmetry energy, and symmetry energy slope at the nuclear saturation density ρ 0 . The term inside the square bracket of Eq. (8) is a density-dependent symmetry energy S(u) = S 0 + L 3 (u − 1) + · · · . In the vicinity of the saturation density, the properties of nuclear matter are characterized by the saturation parameters, ρ 0 , w 0 , K 0 , S 0 , and L, and they are tabulated in Table 2 for the EoSs adopted in this study. Recent terrestrial experiments suggest 220 MeV K 260 MeV, S 0 36 MeV, and L 80 MeV [63,64]. Therefore, the Togashi and LS220 EoSs are consistent with the experimental values while the Shen EoS has large values of S 0 and L.
Properties for three nuclear EoSs: Togashi (green), Shen (pink), and LS220 (orange). Left panel: Relation between baryon density and total pressure. The red vertical line indicates nuclear saturation density ρ 0 of the LS220 EoS. Blue region is from the flow in experiment of heavy ion collisions [61]. Right panel: Mass-radius relation of neutron stars. Light blue bands show two measurements of 1.908 ± 0.016 M ⊙ of pulsar J0348+0432 [21] and 2.01 ± 0.04 M ⊙ of pulsar J1618-2230 [22]. Many dark orange dots indicate the results from the observation of GW170817 [25], and blue regions indicate the results based on the low mass X-ray binary observations [24] (Solid curve: 1σ, Dashed curve: 2σ).
The mass-radius relations of neutron stars are plotted in the right panel of Fig. 2 for the EoSs adopted in this study. In the present paper, we ignore the effects of rotation and magnetic field of neutron stars. Therefore, provided the relation between the total mass energy density ρ and the pressure P , the gravitational mass M and the radius R of neutron stars are obtained from TOV equation as seen from Eqs. (1) and (2) [65]. In the right panel of Fig. 2, we also show the suggestions from the observations: masses of heavy neutron stars observed so far [21,22], gravitational wave GW170817 emitted from neutron star merger [25], and flux from low mass X-ray binaries [24]. The maximum mass of neutron stars is 2.21 M ⊙ for the Togashi EoS, 2.17 M ⊙ for the Shen EoS, and 2.04 M ⊙ for the LS220 EoS. Hence, the three EoSs are consistent with the mass measurements of the heavy neutron stars. On the other hand, these EoSs are different with respect to the radius of neutron stars. The observations of the neutron star merger and the low mass X-ray binaries suggest the smaller radius of neutron stars less than around 13 km. Thus, the Togashi EoS is the most preferred model among the EoSs adopted in this study. Note that, the Togashi EoS has the smallest value for L and it is consistent with the fact that R and L have a positive correlation (e.g. R ∝ L 1/4 [12,66]) in general.

Results
We describe our numerical results of the neutron star cooling. The temperature profiles as a function of the baryon density are shown in Fig. 3. We adopt isothermal models for the initial conditions of the calculation except for the outer layer where ρ B 10 5 g cm −3 .
6/16 Without superfluidity, the thermal structure of the models with the Shen EoS (1.8M ⊙ ) and LS220 EoS (1.7M ⊙ ) show that the DU process works in the core of the star, and the temperatures of the cores drop rapidly. The others show that the DU process does not work and keep the thermal structure almost isothermal. It represents the slow cooling. We present the cooling curves of neutron stars obtained from our calculations with the Togashi EoS, Shen EoS, and LS220 EoS in Figs. 4, 5, and 6, respectively. We focus on the dotted curves corresponding the cooling curves without superfluidity. In Fig. 4, the cooling curves of the Togashi EoS locate at high-temperature regions due to the slow cooling as shown in Fig. 3. This model does not account for the observational data below the curves. In Fig. 5, the cooling curves show that the models with the Shen EoS cools rapidly at any masses, and all of them do not cross observational data of INS. In Fig. 6, the cooling curves are shown for the models with the LS220 EoS. The heavy neutron stars (M ≥ 1.4M ⊙ ) cool rapidly, while the light neutron stars (M < 1.4M ⊙ ) keep warm. By varying the mass for the models with the LS220 EoS, the cooling curves cover the range of the temperature 7/16 For the cooling data of Cassiopeia A (red symbols), we adopt the data in Ref. [30]. The others are taken from Ref. [12] (black dots with errors).
observations. Therefore, the models with the LS220 EoS are suitable to account for the observational data. We examine the models with superfluid effect. First, we notice the models with the Shen EoS, whose cooling curves locate lower-temperature region than all observational data without considering superfluid effect. We show the cooling curves of the models with the Shen EoS with (without) superfluid effect in the solid (dotted) curves in Fig. 5. The effective temperature for the models with the Shen EoS with superfluidity is higher than that without superfluidity, and cooling curves with superfluidity locate in higher-temperature region but are insufficient for the observations with any superfluid models. For instance, any cooling curves with the Shen EoS are inconsistent with the two observations with t 10 4 yr and T ∞ eff 10 6.2 K. Next, we focus on the models with the LS220 EoS, whose cooling curves match to the observational data without considering the superfluid effect as shown in Fig. 6. The cooling curves are sensitive to the difference of the superfluid models. If the superfluid model for 3 P 2 neutrons is the EEHO, the cooling curves are consistent with the observations. In contrast, the models with BEEHS become inconsistent due to the quantitative difference 8  of the superfluid effect. At high-density region, the neutron star with the LS220 EoS causes the strong DU process, which is greatly suppressed by nucleon superfluidity. Therefore, with the DU process, the superfluid effect on the cooling curves is large and the cooling becomes slow. However, without the DU process, the effect from the suppression of neutrino emission is clearly weak because the PBF process works. The PBF process is about 2-3 order of magnitude weaker than the DU process, and 1-2 order stronger than the slow cooling processes. The EoS models which do not induce the DU process such as the Togashi EoS are therefore insusceptible from the superfluidity as seen in Fig. 4. The cooling curves of the Togashi EoS locate in higher-temperature regions with considering superfluid effect, but compared to those of the Shen and LS220 EoSs, the effect of superfluidity is clearly ineffective. With and without the superfluidity, the Togashi EoS does not suit the observations for t > 10 2 yr. We also examine the superfluid-model dependence of the cooling curves. Without the DU process, the superfluid contribution to cooling is small as we noted. If the DU process occurs, the suppression of the neutrino emission by superfluidity is important for the neutron star cooling. In Fig. 7, the central temperatures of each model are plotted with T cr as a function of the Fermi wave number. From this figure, we can see which model has the superfluid effect as a result of the cooling (below 10 8 K). For instance, considering the superfluid effect of 1 S 0 protons for the Shen and LS220 EoSs, the threshold mass against the CCDK is about 1.1M ⊙ and that against the AO is less than 1.1M ⊙ . This implies that the suppression of neutrino emission by the CCDK is stronger than that by AO. Figures 5 and 6 indicate such difference, but the difference is smaller than that of 3 P 2 neutrons. As seen in Fig. 7, the threshold mass for 3 P 2 neutrons against the EEHO for the Shen EoS is about 2.1M ⊙ and that for the LS220 EoS is around 1.6M ⊙ . Therefore the difference of the cooling curves between 1.4M ⊙ to 1.7M ⊙ (1.8M ⊙ to 2.1M ⊙ ) for the LS220 (Shen) EoS are large as shown in Fig. 6 (Fig. 5). Similarly, the cooling is slow for BEEHS as in Figs. 5 and 6 since the superfluid range of the BEEHS is wide enough to easily exceed the central Fermi wave number with 2.1M ⊙ for any EoSs. We recognize that the best models in comparison with the INS observations are with the LS220 EoS and the EEHO 3 P 2 neutron superfluidity. 11

S(u)
[MeV] Shen LS220 Fig. 9 Symmetry energies of three EoSs: Togashi (green), Shen (pink), and LS220 (orange). We adopt the saturation number density n 0 = 0.155 fm −3 for LS220 EoS. The DU process thresholds of npe and npeµ matter correspond to the dotted and solid black curves, respectively. As in Fig. 8 5. Discussion 5.1. The DU process and the symmetry energy at high densities As described in the previous section, the EoS determines whether the DU process turns on or not. This is because Y p depends on the EoS. In our calculations, EoSs are assumed to be in the β equilibrium and charge neutrality: where µ n , µ p , µ e , and µ µ are the chemical potentials of neutrons, protons, electrons, and muons, respectively. In Fig. 8, we show Y p as a function of the baryon mass density for the EoSs adopted in this study. In this figure, adopting the LS220 EoS, we also plot the DU threshold, Y DU p , obtained by substituting Y e and Y µ in Eq. (7). The DU process turns on for Y p ≥ Y DU p . We can recognize that, with use of the LS220 EoS, the DU process turns on in neutron stars with 1.4M ⊙ but this is not the case for neutron stars with 1.1M ⊙ . In case of the Togashi EoS, Y p is lower than the DU threshold and, therefore, the DU process is disallowed even for the maximum-mass neutron stars. In contrast, since Y p is higher than the DU threshold, the DU process turns on for all the models of the Shen EoS considered in our calculations. Note that, whereas Y DU p depends on the EoS, the difference of Y DU p among the EoS models is insignificant compared to that of Y p . 12/16 For the EoS of neutron star matter, Y p is related to the symmetry energy, S(u). When the symmetry energy is defined as the energy difference between the pure neutron matter and the symmetric nuclear matter as in Eq. (8), the chemical potentials of neutrons and protons satisfy the relation: Then, Y p is determined by Eqs. (9), (10) and (11). In general, Y p is smaller for the EoS with lower symmetry energy, provided that the baryon number density is the same. This trend is confirmed for the EoSs adopted in this study from Fig. 9, where the symmetry energy is shown as a function of the baryon number density. At high densities, the Togashi EoS has the lowest symmetry energy while the Shen EoS has the highest. This is consistent with the fact that the symmetry energy slope parameter, L, is small for the Togashi EoS but is large for the Shen EoS. Whereas the value of L is defined at the saturation density, it can be used as an indicator of the symmetry energy at high densities. As a result of the discussion presented so far, the DU process is disallowed even for heavy neutron stars if the EoS adopted for the cooling calculations has low symmetry energy at high densities, or, small L value. Here, we consider the critical symmetry energy as a function of density, S DU (u), for the DU process. When muons are absent, the charge neutrality requires Y e = Y p and, therefore, S(u) and Y p satisfy the relation [11,67]: with the Planck constant h and the nuclear saturation number density n 0 . Then, the condition for the DU process is Y p ≥ Y DU p = 1/9, or, Similarly, taking into account the muon contributions, we can obtain Y p , Y e , and Y µ from Eqs. (9), (10) and (11) if the symmetry energy S(u) is provided as a function of the density u. Then, we can evaluate the DU threshold Y DU p by Eq. (7) and judge whether the DU process turns on. The critical symmetry energy S DU (u) derived in this way is plotted in Fig. 9. As seen in this figure, the difference of S DU (u) between npe and npeµ matter is a few MeV.
Since the DU process is disallowed for any masses, the models with the Togashi EoS are inconsistent with the observations of low-temperature INS as described in the previous section. It stems from the fact that the Togashi EoS has low symmetry energy at high densities, or, small L value. In contrast, a small L value is preferred from the observations of neutron star radius, which indicates the small value for the radius. Therefore, we emphasize that we should consider not only the radius but also the cooling curve so as to examine the neutron star EoS. Note that the fast cooling processes due to exotic particles are omitted in the above discussion.

Influence of the surface composition on cooling curves
In our cooling calculations as described in the previous section, the neutron star surface is assumed to have the light compositions, which are 73% of 1 H, 25% of 4 He, and 2% of 56 Ni. Many previous studies show that the surface compositions affect the effective surface 13 Fig. 10 Same as Fig. 4 but with pure-Ni surface temperature of INS [33,34]. Therefore, using the Togashi EoS, we calculate cooling curves with the pure-Ni surface, which is the heaviest element in conceivable compositions. The results are shown in Fig. 10. Compared to Fig. 4, the overall effective temperature of INS with pure-Ni surface is lower than that of INS with light compositions. In addition, when the age of INS becomes 10 5−6 yr at the photon cooling stage, INS with light surface compositions cool rapidly. We see however that these modifications of surface compositions on cooling curves are insufficient to account for the observed data of INS with low temperature. This is similar to the result of Ref. [68] (see Figure 4 in it). Therefore, the cooling curves with use of the Togashi EoS are inconsistent with the INS with low surface temperature, even if the INS surface is composed of heavy elements.

Concluding remark
We have performed cooling simulation of neutron stars with use of the Togashi, Shen and LS220 EoSs. As a result, the LS220 EoS can account for the temperature observations of INS through the DU process by considering appropriate superfluid model. On the other hand, the Togashi EoS is inappropriate to account for the observations of INS regardless of surface composition because the DU process is prohibited. The fundamental reason comes from rather low symmetry energy, which is appropriate for the neutron star radius. Therefore, other fast cooling processes than the DU process may work in INS with low temperature. For example, the neutrino emissions concerned with exotic particles such as hyperons, pions, kaons, and quarks lead to rapid cooling (e.g. [69,70]). Since the value of the symmetry energy is predicted to be small and the DU process is hard to occur, we can imply the possibility 14/16 that some exotic particles exist in neutron stars. We will investigate their impacts on the cooling curve, as well as on the EoS, elsewhere. Furthermore, the critical temperature of the superfluid transition also depends on the EoS. Calculations taking into account the consistency between the EoS and neutrino emissivity are worth performing.