Neutron star cooling in modified gravity theories

We study thermal evolution of isolated neutron stars in scalar-tensor theories for the first time. Whether the rapid cooling due to the direct Urca process occurs or not is an interesting question in the viewpoint of the temperature observation of isolated neutron stars. Moreover, investigation of the cooling effect of nucleon superfluidity also has the large uncertainties though it is important in modern cooling theory. The cooling effect is typically influenced by the proton fraction and the central density. If a fifth force is mediated due to modification of gravity, the relation between the central density and mass of neutron stars differs from one in general relativity, and the cooling curve is also naively expected to be varied. We find that an unscreened fifth force near the surface of neutron stars changes mass-central density relation, and the direct Urca process can be triggered even for neutron stars with smaller mass. We also present cooling curves including nucleon superfluidity under the scalar-tensor theory. These results show that it might be useful to test gravitational theories with cooling observations of neutron stars.

We study thermal evolution of isolated neutron stars in scalar-tensor theories for the first time. Whether the rapid cooling due to the direct Urca process occurs or not is an interesting question in the viewpoint of the temperature observation of isolated neutron stars. Moreover, investigation of the cooling effect of nucleon superfluidity also has the large uncertainties though it is important in modern cooling theory. The cooling effect is typically influenced by the proton fraction and the central density. If a fifth force is mediated due to modification of gravity, the relation between the central density and mass of neutron stars differs from one in general relativity, and the cooling curve is also naively expected to be varied. We find that an unscreened fifth force near the surface of neutron stars changes mass-central density relation, and the direct Urca process can be triggered even for neutron stars with smaller mass. We also present cooling curves including nucleon superfluidity under the scalar-tensor theory. These results show that it might be useful to test gravitational theories with cooling observations of neutron stars.

Introduction
The recent observations of gravitational waves from binary system [1] play an important role in astrophysics and cosmology for testing the relativistic nature of gravity. One of the most significant results of LIGO and VIRGO is the detection of gravitational wave event (GW170817) from a neutron star (NS) binary and its electromagnetic counterpart (GRB170817A) [2], which provides a stringent constraint on the propagation speed of the gravitational waves [3]. This observation enables us to exclude a large theoretical space of gravitational theories beyond general relativity, which could be responsible for explaining the accelerated expansion of the universe [4][5][6][7][8][9]. Among a vast class of models proposed so far, the remaining scalar-tensor theories in the Horndeski theory [10] after imposing this constraint consist of a non-minimal coupling to gravity with the kinetic gravity braiding (K-essence plus cubic galileon) [11], including f (R) theories [12][13][14]. In such theories, the modification of gravity can be significant at cosmological scales in general, and hence observations such as type Ia supernovae [15,16], cosmic microwave background [17], large scale structures [18][19][20] and clusters of galaxies [21,22] are able to pin down parameter space in modified gravity theories. In addition, the solar-system experiments allow us to test gravity at short distances even for extremely small deviations from general relativity, and it is therefore possible to put a stringent constraint on a certain class of modified gravity theories [23]. During the last few decades, these tests are well-formulated in a various context of observations for modified gravity theories, and the future observations are expected to improve these observational constraints (See for review e.g., [25]). These cosmological observations and solar-system experiments can, however, test gravity only at low energy scales; therefore, testing strong gravity regime could provide extra information of modified gravity theories, and non-trivial effects from relativistic objects, especially NSs, should be extremely important in future observations, and such tests in strong gravity environments could be a powerful tool for measuring deviation from general relativity.
NSs, which are born after supernova explosion, cool down by emitting a large number of neutrinos and photons. In the absence of companion stars, isolated neutron stars (INS) cool by the neutrino emission for the time t 10 5 yr after their birth, since then the photonemission cooling follows. In the era of the neutrino cooling, thermal evolution of INS depends on the interior physics of the equation of state (EOS) including nucleon superfluidity. The most important problem for INS cooling is whether the direct Urca (DU) process, which is the strong neutrino emission through β and inverse β decay, occurs or not. The DU process is unnecessary for most of the observed neutron stars, which can be explained by the slow cooling processes with an additional cooling effect due to nucleon superfluidity, minimal cooling scenario [26]. The minimal cooling scenario is therefore widely recognized to be the most successful standard cooling theory for most of temperature observations. However, some of neutron stars are thought to have enhanced cooling such as the DU process in order to explain their observed low temperature. For example, the pulsar J0205+6449 in supernova remnant 3C58 and RX J0007.0+7302 are cold enough to need the fast cooling processes [26]. For accreting neutron stars in steady state, SAX J1808.4−3658 and 1H 1905+000 have low quiescent luminosity for high mass accretion rate [27] and therefore they are thought to be the strong evidence of occurrence of fast cooling processes. Moreover, some cold compact objects such as, G127.1+0.5, G084.2+0.8, G074.0−8.5, G065.3+5.7 [28], and G043.3−0.2 [29], have been detected though it's unknown whether they are NSs or black holes, respectively. This implies that, if any one of them is NS, a fast cooling process should occur in this star (e.g., see Figure. 11 in Ref. [30]). Focusing on the occurrence of the DU process, the thermal evolution of neutron stars is determined by the Y p -density relation [31] and the mass as well as the EOS. Therefore, the EOS could be constrained by observations of thermal evolution of INSs (e.g., [32]). The recent observations indeed have rejected some of the EOSs [33][34][35][36][37][38][39], though there are uncertainties in the mass-radius relation as well as the possible effects of exotic particles. Moreover, the change of the total neutrino emissivities by nucleon superfluidity is unclear for especially triplet gap of neutrons because of uncertain nuclear potentials [40][41][42][43]. Due to these uncertainties, there can be many scenarios to account for the observations of the age and surface temperature [26,44].
The observations of NS structure are rapidly advancing. For example, a few heavy NSs with M ∼ 2 M have been recently discovered by relativistic Shapiro delay combined with the 2/34 mass function of Kepler motion [33,35,37], where M is the solar mass. Especially, the millisecond pulsar PSR J0740+6620 discovered using the Green Bank Telescope has the most heaviest mass so far, which is estimated to be M = 2.14 +0.10 −0.09 M within 1σ region [37]. Such constraints about massive NSs enable us to reject many soft EOSs. Furthermore, the NICER has reported the observations of another millisecond pulsar PSR J0030+0451 , where the observational mass M and the equatorial radius R eq are estimated to be M = 1.34 +0. 15 −0.16 M , R eq = 12.71 +1.14 −1.19 km [38] and M = 1.44 +0.15 −0.14 M , R eq = 13.02 +1.24 −1.06 km [39] within 1σ regions, respectively. 2 Going back to modified gravity theories, an extra degree of freedom in scalar-tensor theories typically modifies the law of gravity even for small scales, and therefore screening mechanisms such as the Chameleon mechanism [45,46] or the Vainshtein mechanism [47] are crucial to pass stringent tests in the solar system. On the other hand, a screening mechanism might not completely hide the fifth force in NSs as reported in [48]. Once the mass-radius relation (or the mass-central density relation) are modified, one would expect that thermal evolution of NSs can be also different from general relativity. If this is true, there should be deviation in cooling curves compared with that of general relativity, and one might be able to test modified gravity theories by using cooling curves of NSs. With this in mind, for the first time, we investigate how the cooling curve of NSs changes in the scalar-tensor theory proposed in [48]. Although this model can be tightly constrained by the local gravity experiments [49], it could potentially change the structure of neutron stars from the one in general relativity as shown in [48] and is therefore suitable to see the effect on the cooling curve. The structure of NSs in other modified gravity theories such as f (R) and galileon theories has been intensively investigated in many literatures [50][51][52][53][54][55][56][57][58][59][60][61][62][63]. In the present paper, we show how the attractive nature of the fifth force affects to the cooling curve, which can be applicable for other Chameleon theories.
This paper is organized as follows. In Sec. II, we review the model and structure of NSs studied in [48] and the EOSs adopted in the present paper. In Sec. III, we give an overview of NS cooling, including observations of isolated NSs. In Sec. IV. we show cooling curves of NSs in the scalar-tensor theory. Sec. V is devoted to discussions and conclusion. Details on the relation between Jordan and Einstein frame are summarized in Appendix A. In Appendix B, we briefly review the Chameleon mechanism. In Appendix C, we present the cooling curves of NSs in the massless Brans-Dicke theory as a complimentary material. Throughout this paper, we use the reduced Planck mass M pl = 1/ √ 8πG, the densityρ 0 = m n n 0 = 1.6749 × 10 14 g · cm −3 with n 0 = 0.1 (fm) −3 , where m n = 1.6749 × 10 −24 g is the neutron mass, and the distance r 0 = c/ √ Gρ 0 = 89.664 km for normalization.

Theory
In this section, we give a brief summary of theoretical set-up, basic equations, and numerical solutions of NSs. We consider the following model in scalar-tensor theories [64], 2 Although the same observational data are used in the Refs. [38,39], the discrepancy between these constraints stems from independent analysis of these ground.

3/34
where g is the determinant of metric tensor g µν , R is the Ricci scalar, Q is a constant, Ψ m describes the matter field minimally coupled with g µν , and the non-minimal coupling F (φ) is chosen to be We focus on the model with the following potential Here, the value of Q in (3) is chosen such that φ becomes a canonically normalized field in the Einstein frame without a filed redefinition [64], and the infinite limit λ → ∞ corresponds to general relativity. In this model, the canonical kinetic term for the scalar field is apparently absent due to the choice of Q, however, the scalar field regains a kinetic term from the nonminimal coupling with gravity, as we will see below. The structural features of the action (1) with Q = −1/ √ 6 is the well-known relation with the metric f (R) gravity [65] given by through an appropriate choice of the function F and the potential V , The φ 4 potential given in (3) arises from the f (R) theories described by f (R) = R + aR 4/3 where a is a constant [48]. Although the theories (1) and (4) are equivalent, for simplicity, we use the scalar-tensor action (1) throughout this paper. We keep the parameter Q in this section for the sake of completeness, and the case in the massless Brans-Dicke theory is discussed in Appendix C. The model with the φ 4 potential potentially possesses a screening mechanism for the fifth force in dense environments [45,46]. Thanks to this screening effect called the Chameleon mechanism, the model with the φ 4 potential can be consistent with the local gravity experiments (see e.g., [49]). In Appendix B, we summarized a quick introduction to the Chameleon mechanism. In this paper, we choose relatively larger λ where the corresponding Compton wavelength is roughly the size of NSs. In such case, the scalar field φ can be ignored in cosmological dynamics although one needs the cosmological constant or dark energy to explain the accelerated expansion of the universe. For this reason, a constraint on λ is irrelevant to cosmological observations in our case. The variation with respect to the metric gives the Einstein equation, where the energy-momentum tensor is defined as in usual, T µν = −(2/ √ −g)(δS m /δg µν ), F ,φ = ∂F/∂φ, and F ,φφ = ∂ 2 F/∂φ 2 . The Euler-Lagrange equation for the scalar field φ is given by The Ricci scalar in (7), can be removed by taking the trace of the Einstein equation. As one can see, with the choice of Q specified by (3), the kinetic term in (7) vanishes. However, substituting the trace of Einstein equation back into the scalar field equation, it can be easily seen that the scalar field sourced by the trace of the energy-momentum tensor of the matter field, and the kinetic term of the scalar field appears after this substitution. As in general relativity, the Bianchi identity provides the well-known form of the energymomentum conservation,

Modified Tolman-Oppenheimer-Volkoff (TOV) equation
In this subsection, we summarize the modified Tolman-Oppenheimer-Volkoff (TOV) equations and boundary conditions. We adopt the following static and spherically symmetric background metric, and consider a perfect fluid, which is given by T µ ν = diag(−ρ(r), P (r), P (r), P (r)). Here, ρ is the energy density and P is the pressure for the matter. The energy-momentum conservation (8) gives the continuity equation where a prime represents the derivative with respect to r. After solving the Einstein equations and the scalar field equation in terms of f , M, and φ , we get [48] where we introduced the mass function M(r),

Equation of state
We adopt the three EOSs constructed on the Brussels-Montreal-Skyrme (BSk) functionals [66], which are the Hatree-Fock-Bogoliubov (HFB) mass models with 'usual' Skyrme effective interaction including extra two terms of density-dependent generalizations, BSk19, BSk20 and BSk21 [67]. These functionals are fitted to the experimental nuclear masses from Atomic Mass Evaluation (AME) in 2010, and the experimental constraints of the symmetryenergy parameters, including the slope parameter L, from measurements of heavy-ion collisions and neutron-skin thickness. Moreover, these EOSs are accurately fitted to realistic 5/34 pure-neutron matter (PNM) EOS; BSk19, BSk20, and BSk21 EOSs are fitted to FP [68], APR [69], and LS2 [70] PNM EOSs, respectively, which are based on a variational method with the realistic 2-body interaction with 3-body force. Thus, the EOSs constructed with BSk functionals are successful to describe properties of both measured nuclei in the ground state and the infinite nuclear matter. To obtain the basic properties of cold NS matter, we use the public code 3 which fully reproduces the three EOSs [71].
In Fig. 1, we show the EOS softness and the distribution of proton fraction Y p as functions of the density. As one can see, the softest EOS among three EOSs is BSk19, while the stiffest EOS is BSk21. The value of Y p of the BSk21 is much higher than that of the other EOSs for ρ 2ρ 0 . We note that the values of Y p of the BSk19 and BSk20 are clearly lower than DU threshold for all region, which predicts non-rapid cooling of NS. These profiles of pressure and Y p distribution can be explained by the difference of the slope parameter of the symmetric energy defined by L = ρ nuc , where S(ρ) is the symmetry energy and ρ nuc is nuclear saturation density. The value of L calculated in Ref. [67] is L = 31.9 MeV for the BSk19, L = 37.4 MeV for the BSk20, and L = 46.6 MeV for the BSk21, whose values are consistent with the behaviors of the pressure and Y p distributions in Fig. 1. Thus L is an indicator of the relevancy for the EOS softness and the threshold mass of the nucleon DU process (hereafter M DU ). For more precise predictions, it is necessary to know the higherderivative terms of the symmetry energy with respect to the finite density and Y p , where the interactions between nucleons are still unclear. Therefore, the EOS itself should have some uncertainties of the order of 10% in ρ-P relation. As we will see in the next section, the 3 http://www.ioffe.ru/astro/NSG/BSk/index.html 6/34 fifth force is sensitive to ρ − 3P , and this will therefore change the mass-radius relations. However, our method in the present paper should be still valid for other EOSs. , and E φ (lower panel) as a function of radius r in the unit of km. The purple, blue, and green line correspond toλ = 10 3 , 10 4 , 10 5 , respectively. The central density is given by ρ 0 = 8.1 × 10 14 g/cm 3 and the EOS is adopted BSk21.

Numerical solutions
Here, we give an overview of the structure of NSs by solving the modified TOV equations with the EOSs, BSk19, BSk20, and BSk21. In order to numerically solve these equations, we impose the following boundary conditions such that the regularities at the center of NSs are satisfied, For simplicity, we assume the asymptotic flatness at spatial infinity, 4 Hereafter, we use the rescaled dimensionless parameter, with r 0 M pl = 1.107 × 10 39 . By imposing the boundary conditions (15), we numerically solve Eqs. (10)-(13) with the EOS outward from the center of star (r = 0) in order to clarify the profiles of f, M, φ, ρ, and P inside the star. The star radius r s is determined by the condition P (r s ) = 0, and we identify the mass of star M s = M(r s ). Outside the star, we simply set ρ = P = 0 and solve Eqs. (11)-(13) for f, M, φ up to the sufficient large distance, e.g., r = 10 5 r s , at which both φ and φ sufficiently approach 0. Here we note that the asymptotic values of f and φ at r → ∞ do not satisfy the other boundary conditions (16) in general. However, we can shift the asymptotic value of f to 1 by virtue of time reparametrization invariance. Regarding the scalar field, we resort to a shooting method in order to find the initial value φ(r = 0) which eventually satisfies the boundary condition φ(r → ∞) = 0. See [48] for more details of numerical procedure. In this work, following the Ref. [48], we fix the dimensionless parameterλ to be 0, 10 3 , 10 4 , 10 5 . By choosing these values, the corresponding Compton wavelength of the scalar field can be comparable with the size of NSs, implying that the fifth force effect can significantly affect the structure of NSs. However, these parameters are unfortunately inconsistent with the constraint by the Eöt-Wash torsion balance experiment [101], λ O(1) (e.g., [102]), while it is still consistent with other astrophysical and laboratory tests [24]. Nonetheless, this choice of parameters enables us to capture the effect of the fifth force in the cooling curve even for other modified gravity theories possessing the Chameleon mechanism as we will discuss in Sec. IV and V.
Let us first discuss the behavior of the scalar field. Since the scalar field equation, (7) or (13), is complicated due to the non-minimal coupling to gravity, we work in the Einstein frame for convenience. By using the conformal transformation of the metric (g µν ) E = F (φ)g µν , one can remove the non-minimal coupling of the scalar field to gravity, and the resultant theory consists of the Einstein-Hilbert term, the canonical scalar field, and the matter field coupled to the scalar field. In Appendix A, we have summarized the equation of motions for gravity, the scalar field, and the matter component, and the relation between the Jordan frame and the Einstein frame. We note that all the boundary conditions (15) and (16) and the central density ρ 0 are given in the Jordan frame in numerical computation, although we discuss the behavior of the scalar field in the Einstein frame in the below.
As summarized in Appendix A, the scalar field obeys the Klein-Gordon equation with the effective potential given by (A7). Outside NSs, the effective potential is nothing but V (φ). As ρ E − 3P E increase, the minimum of the effective potential is no longer located at φ = 0, and the effective mass m 2 eff = d 2 V eff /dφ 2 , being dependent on φ, becomes large inside NSs. When the scalar field acquires sufficiently large mass in such situation, the fifth force due to the scalar field is in general screened, and deviation from general relativity becomes relatively small. Therefore, we naively expect that a screening of the fifth force could be effective inside NSs. To see the fifth force effect we introduce the following parameters in the Einstein frame as defined in (A13) and (A14) where r E , (T µ ν (φ) ) E , (G µ ν ) E , ρ E , and P E are the radius, the energy-momentum tensor of the scalar field, the Einstein tensor, the energy density, and the pressure in the Einstein frame. Ω scalar field, normalized by pressure-gradient, i.e., the first term in the continuity equation (A11). Namely, E φ characterizes the ratio of the fifth force to the total force. Roughly, we may regard the effective gravitational constant as G eff (1 + E φ )G. In Fig. 2, we plot these dimensionless parameters as a function of the radius in the Jordan frame, r. For anyλ, the parameter Ω (h) φ is approximately zero well inside NSs and rapidly increases a little inside the surface of NS. On the other hand, the plot of Ω (f ) φ shows that it can be large, for instance Ω (f ) φ 0.2 at r = 10 km forλ = 10 3 . In addition, the effect of the fifth force is still active even well inside NSs forλ = 10 3 , as one can see the plot of E φ . On the other hand, for λ = 10 5 , at sufficiently small r, the fifth force is negligible, thus this implies that Chameleon 10/34 mechanism effectively works for sufficiently largeλ well inside NSs 5 . However, since the energy density and pressure for the matter is zero for r ≥ r s , the effective mass m φ tends to increase for r < r s . This implies that the effect of the fifth force is still significant near the surface even for largeλ, which can be also seen in the Fig. 2, and this fact affects the structure of NSs.
These arguments in the Einstein frame agree with mass-radius relation in our numerical calculation. Fig. 3 shows the mass M s (normalized by the solar mass M = 1.9884 × 10 33 g) as functions of the radius r s for BSk19 (upper left panel), BSk20 (upper right panel), and BSk21 (lower panel) EOSs. Even forλ = 10 5 , the fifth force effect near the surface of NS significantly affects MR relation. As one can see, the maximum mass of NSs decreases as the parameterλ decreases. Since the fifth force due to the scalar field is attractive and the net gravitational force from the metric g µν and the scalar field φ also becomes larger. Based on this fact, NSs tend to be compact for smallλ, which can be also seen from the Fig. 4, and this is why the maximum mass of general relativity is the largest value. This feature is similar to the gas density profile of galaxy clusters in an f (R) gravity model, as found in Refs. [21,22].

SETUP OF NEUTRON STAR COOLING
INS with the low magnetic field has no heating source, and after being born by a supernova explosion, INS cools down by the emission process of neutrino at the early stage (t 10 4−5 yr) and by that of photons at the late stage (t 10 5 yr). The emissivities of neutrinos and photons are therefore important to describe the thermal evolution of INS. In this section, we give a brief review of physics about the neutrino and photon energy loss.
The most important factor for comparing with the observations of INS is the neutrino emission process, which has been studied by many works (for review, see Ref. [72]). We solve the thermal evolution of INS in the Jordan frame, then we assume that the radiative processes through photons and neutrinos are not affected by the modification of the gravity. In general, the neutrino emission processes are classified into fast and slow cooling processes. Within the standard matter, n, p, e, µ, the fast cooling process is the Direct Urca (DU) process, which is the forward and inverse β decay. The neutrino emissivity of the DU process is given by [31,73] where m * n /m n and m * p /m p are effective mass ratio of neutrons and protons, respectively, which depend on the EOS. Here ρ B is the baryon density, T 9 is the temperature in units of 10 9 K, and Y l is the lepton fraction, Θ npl denotes the step function: Θ npl = 1 when the momentum conversation is satisfied k Fn = k Fp + k Fl , while Θ npl = 0 when k Fn = k Fp + k Fl , where k Fn , k Fp , and k Fl are the wavenumber of neutrons, protons, and leptons, respectively. From the energy-momentum conservation of k Fn , k Fp , and k Fe , the proton fraction Y p for the DU process being effective satisfies the following condition: where Y DU p is the net threshold Y p of DU and x e = Y e / (Y e + Y µ ). In the very high-density regions, (19) becomes Y DU p 0.1477 because of the condition Y e Y µ . Whether the DU process occurs or not is determined by the proton fraction in the central region of INS, which may be changed by gravity and EOS models. Thus, the gravity and the EOS models are significant for the cooling of INS.
Generally, the DU process is prohibited for INS with light masses, in which the thermal evolution of INS is determined by the slow cooling processes, the Modified Urca (MU) process and the neutrino bremsstrahlung radiation in nucleon-nucleon collisions [74,75]. The MU process is the reaction of n + N → p + N + l +ν l , p + N + l → n + N + ν l , where N is neutron or proton, although the DU process is n → p + l +ν l , p + l → n + ν l . The bremsstrahlung is the neutrino radiation emitted when electrons or muons are braked by the Coulomb field of protons. These emissivities can be approximately expressed as Here, we ignore the dependences of effective mass ratio and fraction of particles in (20) though we consider them in numerical calculation. Compared to (18), we see that the slow cooling processes are much weaker than the DU process, but most of the reactions always occur above the threshold density of neutron drip 6 . Therefore, the slow cooling processes, especially the modified Urca process which is stronger than the bremsstrahlung [75], are important for describing thermal evolution of all INSs. In modern cooling theory, however, the above cooling processes are not sufficient because neutrons or protons are highly expected to be in a superfluid state due to low-temperature environment. Specifically, the 1 S 0 neutrons superfluidity is contained in the inner crust, and the 3 P 2 neutrons superfluidity and 1 S 0 protons superconductivity exist in the core. According to the BCS theory [76], such nucleons make their pairs which can be expressed as following reaction: Here N is a neutron or proton. [N N ] denotes the N − N Cooper pair and behaves as a bosonic particle. Each nucleon superfluidity has their superfluid transition temperature T cr , but T cr dependence of the density still has large uncertainties above all for the 3 P 2 neutrons superfluidity. In this work, we adopt the superfluid models as shown in Fig. 5; T84 model for neutrons 1 S 0 [40], 'a" model for neutrons 3 P 2 [26], and CCDK model for protons 1 S 0 [78].
The important superfluid effect on NS cooling curve is the pair breaking formation (PBF) process. By making or breaking Cooper pairs, additional neutrino emissions occur with above 6 The MU process with the proton branch does not occur for Y p ≤ 1/65. all T ≈ 0.5T cr , and the emissivities are roughly expressed as follows [75]: whereF i is the control function with T /T cr , which depends on whether the state of nucleons is singlet (i = s) or triplet (i = t), andF s andF t are given in Ref. [79]. Compared with the emissivity of slow cooling processes formulated as Eq. (20), the PBF process is stronger by 1∼2 order of magnitude, depending on the temperature. Without any fast cooling processes, therefore, the PBF process is dominant for INS cooling curves. Thus, nucleon superfluidity effect is though to be an essential factor to consider cooling curves in modern cooling theory. The photon cooling is dominant for the thermal evolution of INS at the late stage. Assuming the blackbody emission, the photon luminosity is given by where a is the Stefan-Boltzmann constant, and T eff,7 is the effective temperature in units of 10 7 K. L γ depends on the surface compositions of the envelope of INS. The surface compositions reflect the relation between surface temperature and interior temperature at ρ = 10 10 g cm −3 . The reason why the surface compositions affect the cooling curves is ascribed to the opacities, which include the radiative opacity and conductive opacity mainly due to electrons. We simply adopt the GPE relation [80], which assumes the envelope only with 56 Fe. We have assumed that the emissivities by photons and neutrinos as well as the transport equation are not affected by the modification of the gravity in the Jordan frame because the equation of motion of the matter is not directly coupled to the scalar field, which is clearly seen in Eq. (1) or Eq. (4). Then, we use the following cooling transport equations for thermal 13/34 Fig. 6 Age and effective temperature for a distant observer of INSs adopted from Ref. [91]. The difference of the color indicates that of models of the atmosphere; carbon atmosphere for the star with red color, hydrogen atmosphere with blue, magnetized hydrogen atmosphere with green, and black-body models with purple. The insets in this upper-left corners indicate magnification of the data with Cassiopeia A and the continuous data are plotted in linear axis about the age. evolution of NS without heating source [81] where C V is the specific heat, and ν is the total neutrino emission including Eq. (18) and Eq. (20). The boundary conditions for T and L γ are given by L γ (r = 0) = 0 from Eq. (22) and T (t = 0) = 10 10 K, but the choice of the initial temperature with T (t = 0) 10 8.5 K does not affect cooling curves with t 10 1−2 yr at all. After giving the initial temperature, we calculate the cooling curves by the low surface temperature of T ∞ eff = 10 5 K. In our numerical computation to find the solution to Eq. (23), we use the Henyey scheme [82] after solving the NS structure from Eqs. (10)- (13). Then, we obtain theoretical cooling curves which describes the thermal evolution of INS. For numerical calculation of INS cooling, we use the public code NSCool 7 [83] which has already included the neutrino and photon emissions.
To examine the validity of the models for the thermal evolution of INS, it is useful to compare the theoretical cooling curves with observational data of the age and the surface temperature. The age is mainly estimated by two method: One is based on a measurement of spin-down rate of a pulsar due to the braking by a rotating magnetic filed. The other is based on estimating the "kinetic age" by connecting the transverse velocity to a distance from the geometric central position of a supernova remnant. Generally, though the "kinetic 7 http://www.astroscu.unam.mx/neutrones/NSCool/ 14/34 ages" are not always known, it is thought to be more reliable than the method of measuring the spin-down rate. The surface temperature is estimated by measurements of X-ray flux. For the estimation, the uncertainty comes from modeling of the atmosphere of NS, which significantly affects the net X-ray flux. We adopt the data of Ref. [84] for Cassiopeia A and other 18 INS in Ref. [91]. The data of 19 INS are plotted in Fig. 6. Cassiopeia A, a young supernova remnant with the age around t ≈ 340 yr, is especially important because the uncertainties of the age and the surface temperature are much smaller than those of the other observations. The surface temperature of Cassiopeia A has been measured for about 20 years by Chandra X-ray detectors. However, we should note that the results may contain systematic errors because the observed decay rate depends on the Chandra X-ray detectors and modes [85,86]. Constraint from the observational data for Cassiopeia A has been investigated in Refs. [87][88][89][90].

Cooling Curves without Nucleon Superfluidity
Let us now discuss thermal evolution of INS in the scalar-tensor theory. Figs. 7∼9 present the cooling curves with various mass for the EOSs, BSk19, BSk20, and BSk21, respectively. As expected from Y p distribution in Fig. 1, we find that the DU process occurs with the BSk21 EOS, while it does not occur with the other EOSs. This is because the slope parameter L is large enough only in the BSk21 EOS to derive the DU process. To see the cooling effect with the BSk21 EOS in more detail, Fig. 10 shows the cooling curves with fixing the mass in steps of 0.01M . This figure shows that the threshold mass that the DU process begins to occur, M DU , tends to be small asλ decrease. We find the threshold mass is in the range 1.25 < M DU /M ≤ 1.26 forλ = 10 3 , 1.36 < M DU /M ≤ 1.37 forλ = 10 4 , and 1.46 < M DU /M ≤ 1.47 forλ = 10 5 . Hence, the DU process occurs even for smaller mass with smallerλ, where the fifth force effect is more active. Furthermore, this figure also indicates that the occurrence of the DU process is determined by the threshold mass in a very narrow range smaller than 0.01M , which agrees with Ref. [91]. The threshold value of M DU obtained from the numerical results is slightly higher than M DU expected from Y p distribution in Fig. 1 and the mass-central density relation in Fig. 4. This is explained by the fact that the occurrence of the DU process is not determined by the value of Y p at the center of a NS but the value of Y p at some finite radius.
As one can see from Figs. 7 and 8 with the BSk19 and BSk20 EOSs, respectively, the effect of the fifth force looks small because the DU process does not occur for any mass. Comparing the theoretical curves and the observational data of Cassiopeia A, we see that the theoretical cooling curves depend on the value ofλ. This is because the central density is lower asλ is larger. Therefore the effective region of the neutrino emissions becomes narrower. On the other hand, the cooling curves at the photon cooling stage can be basically regarded as independent of neutrino cooling. This explains that modification of gravity does not affect the cooling curves at the late stage t 10 5 yr.
In order to compare the dependence ofλ on the cooling curves, Fig. 11 shows the cooling curves with fixing the mass M s = 1.4 M for the three EOSs. One can check that the curves approach to the curve of GR asλ increases. Furthermore, one can see that the theoretical curves with the small valueλ fail to account for the observational data of Cassiopeia A. This feature is clear in the model with the EOS BSk21 because the DU process occurs. This figure exemplifies how the modification of gravity changes the cooling curves. By comparing 15/34 Fig. 7 Cooling curves using the BSk19 EOS.  the curve withλ = 0 and the curve of the general relativity in Fig. 9, the change of M DU is at most 0.4 M with the BSk21 EOS. This panel shows that the cooling of NSs can be a tool for testing the scalar-tensor theories in the strong regime. However, uncertainties of the modeling of NSs, e.g., the parameters of mass and the EOS, must be well fixed.

Cooling Curves with Nucleon Superfluidity
The nucleon superfluidity, which is well known to be one of the important physics in the NS cooling theory was not taken into account in the results of the previous section. For completeness, we here calculate cooling curves with nucleon superfluidity based on the models described in the Sec. 3. Throughout this section, we assume the mass of NSs to be 1.4 M . We present the results with different superfluid models in Figs. 12 and 14 which focus on the dependence of nucleon superfluid models. Moreover, we also show theλ dependence of realistic cooling curves in which all kinds of nucleon superfluidity are considered in Figs. 13 and 15. Since the cooling behaviors with the BSk19 and BSk20 EOSs are similar as we see in the previous section, we do not show the case of BSk19 EOS.
First of all, let us see the case of BSk20 EOS (See Fig. 12), where the DU process is prohibited. Without nucleon superfluidity, the surface temperature gets lower around t ∼ 10 2 yr because the compactness of NSs due to the fifth force effect makes the neutrino luminosity higher. Such a difference from the GR case tends to appear as a knee, when the cooling rate are highly changed in a moment. Meanwhile, the surface temperature at photon cooling stage, whose timescale is approximately proportional to M   The curve of GR in each panel corresponds to the infinity ofλ according to [48]. neutrons 3 P 2 superfluidity promotes, and finally the former makes the difference from the GR case smaller while the latter makes larger with the age when the knee appears. Thus, surface temperature near the knee is slightly different between GR and scalar-tensor theory. This signature due to the fifth force effect might enable us to discriminate between GR and the scalar-tensor theory if a rapid cooling at t ∼ 10 2 yr is observed. At photon cooling stage, on the other hand, surface temperature is increased with scalar-tensor theory because the specific heat is reduced due to the pairing effect.
For the 1.4 M stars with the BSk21 EOS, the DU process occurs withλ 10 4 as we see in Fig. 11. In such models with scalar-tensor theory, the surface temperature is lower by 60 ∼ 80% for t 10 2 yr compared with the GR case (see left-and middle-bottom panels of Fig. 12). If the nucleon superfluidity is considered, the cooling curves take values in highertemperature regions because of the cooling suppression effect. Especially, the suppression effect of 3 P 2 neutrons superfluidity is larger than that of other superfluidity in Fig. 14. However, such a suppression effect is not large enough to cancel the effect of the DU process at least within the models adopted here. Actually, the cooling curves without superfluidity (right panel in Fig. 11) and with superfluidity (left panel in Fig. 15) are not significantly different. Furthermore, the change of the surface temperature between the GR and the scalar-tensor theory within minimal cooling scenario is at most 20% (see right panels of Fig. 13 and Fig. 15 withλ 10 5 ). This implies that the superfluid effect is not so important when the fast cooling processes occur. Although the effect of the cooling suppression due to the superfluid depends on the models, the actual threshold mass of the DU process is useful to distinguish between the theories of gravity even with nucleon superfluidity as suggested in previous section.
Summarizing the above, the effect of fifth force basically tends to make the surface temperature around the knee highly decreased with neutrons superfluidity though the difference is not as large as the case of the DU process. Meanwhile, proton conductivity makes the difference of cooling curves from the GR case smaller. At the photon cooling stage, the surface temperature under scalar-tensor theory is increased. This tendency is similar with the slow cooling scenario. Therefore, we might distinguish between the GR and the scalar-tensor theory from the behavior of the knee, although the dominant cooling effect comes from the DU process depending on the models. Moreover, if there are cold or rapidly cooled observations in the moment, such as the Cassiopeia A [93], it would be possible to put a constrain on the scalar-tensor theory.

Discussion & Conclusion
We have considered the cooling effect of NSs in the scalar-tensor theory with φ 4 potential as well as massless Brans-Dicke theory discussed in Appendix C, adopting BSk19, BSk20, and BSk21 as the EOS for NSs. In the scalar-tensor theory, an additional scalar degree of freedom introduces so-called fifth force and modifies the gravitational law predicted by general relativity. If such fifth force is present, the structure of NSs can also differ from the one predicted in general relativity. We first took a closer look at a fifth force effect in NSs and discussed how it modifies, e.g., the mass-radius relation. We have confirmed that a screening of the fifth force is sufficiently active at high-density region, i.e., near the center of NSs. On the other hand, the fifth force effect cannot be still neglected near the surface of  . Dotted curves do not include superfluidity while the yellow curves all kinds of neutron( 1 S 0 , 3 P 2 ) and proton( 1 S 0 ) superfluidity. Except these cases, some of nucleon superfluidity are considered; curves with "n( 1 S 0 , 3 P 2 )" do not consider proton conductivity, for example.
NS even for largeλ. Such an incomplete screening near the surface significantly affects the mass-radius relation. Since the scalar field always invokes an attractive force in our case, the net gravitational force is effectively stronger and NSs tend to be compact compared with general relativity, as found for the cluster's gas density profile in Refs. [21,22]. This implies that, if we fix the central density, the smallerλ case has smaller mass as confirmed in Fig. 4, and one can draw the constant mass curve inλ-ρ 0 plane 8 as in Fig. 16.
For the nucleon superfluid effect on cooling curves, the behavior of the knee, where the cooling rate is drastically changed, is affected by the fifth force. Hence, the properties of Cooling curves (Left panel) and the derivation of surface temperature from that with GR (Right panel) with the BSk20 EOS, assuming 1.4M stars. All kinds of neutron( 1 S 0 , 3 P 2 ) and proton( 1 S 0 ) superfluidity are considered. The difference of color denotes that ofλ values.
the knee can distinguish between the GR and the scalar-tensor theory, although the clear distinction between them would not be easy due to many model parameters such as EOSs, superfluidity, surface compositions. Observations of cold NSs or rapid-cooling NSs like Cassiopeia A might test the theories of gravity. The observations of Cassiopeia provide us with a precise cooling rate, which would be useful to test the scalar-tensor theory. This is left for a future work. Now let us discuss the cooling effect. Since the rapid cooling of NSs by the DU process does not occur for the EOSs, BSk19 and BSk20, as discussed in Sec. 2.2, we hereafter focus on the case for BSk21. In this case, whether the DU process occurs or not is solely determined by the central density of NSs as shown in Fig. 1. Thus, if the central density of NSs is above a certain value denoted as a threshold density ρ DU (dashed line in Fig. 16), the DU process can be active in NS with the masses above the corresponding M DU for each parameterλ. As clearly seen from Fig. 16, compared with general relativity, even smaller-mass NSs can rapidly cool by the DU process in the scalar-tensor theory, and M DU tends to be small as λ decreases, that is, as the fifth force effect is enhanced. The inverse of the parameterλ typically represents the strength of the fifth force. Therefore, we stress that these arguments on the rapid cooling of NSs can be also applied to other modified gravity theories as long as a fifth force has an attractive nature 9 . On the other hand, scalar-tensor theories that possess  Fig. 14 The same as Fig. 12, but with the BSk21 EOS.
the Vainshtein mechanism might be difficult to see differences in a cooling curve since the interior structure of NSs is almost the same as one in general relativity as reported in [63].
Thus the information about masses and cooling curves of NSs could provide a signature of the Chameleon force. It would be also interesting to investigate cooling effects in other modified gravity theories having the symmetron mechanism [99] and the spontaneous scalarization [100], and we leave these for future work.
In this work, we assumed the specific scalar-tensor gravity models with the non-minimal coupling which mediates the fifth force as F (φ) = e 2φ/( √ 6Mpl) and the self-scalar coupling potential as V (φ) = 1 4 λφ 4 . If we consider other kinds of them, the structure and cooling of NSs naturally differ from those with the current model. However, the property to make NSs compact by the breaking down of the screening mechanism should be universal for the scalar-tensor models because the mediated fifth force is attractive. This implies some universal NS properties for gravity models; If the screening mechanism effectively works, the maximum mass and M DU are higher, and the effect of nucleon superfluidity on cooling curves is lower. Hence, this work will provide a basis for these forthcoming discussions because the 22 Fig. 15 The same as Fig. 13, but with the BSk21 EOS.
qualitative behaviors of NS structure and cooling curves are not affected by changing the law of gravity. We fixed the dimensionless self-coupling constant asλ = 0, 10 3 , 10 4 , 10 5 in this study. Then, the current models are inconsistent with the Eöt-Wash torsion balance experiment demanding λ O(1) orλ O(10 78 ) (e.g., [102]). Since the screening mechanism effectively works withλ = 10 5 . the NS structure and cooling curves should be almost the same as those with GR in the experimentally good λ value. For calculation with such an extremely large λ, however, there is numerical difficulty as presented in [48] (such problems are mentioned in other literatures as well). By tackling this numerical problem, we hope that the above point will be investigated in detail as a future work.
Here the second term and third term in the continuity equation (A11) respectively represent the force from the Einstein frame metric and the scalar field. In analogy with the density parameters in the Friedmann equation, we introduce the following dimensionless parameters which measures the fraction of the energy-momentum tensor for the scalar field in the Einstein equation. In addition to these parameters, we also introduce the force fraction parameter For ρ E − 3P E > 0, the scalar field is an attractive force, and this parameter then satisfies 0 ≤ E φ ≤ 1. Since the derivative of the pressure vanishes at the origin and outside of NSs, we define this parameter within 0 < r < r s . 25/34 Finally, by using the fact that ds 2 E = F ds 2 , we can relate quantities in the Jordan and Einstein frame as With these relation, one can compute the Einstein frame quantities from one in the Jordan frame obtained in numerical calculation.

B. Chameleon mechanism
In this appendix, we review the Chameleon mechanism in a non-relativistic setup (see e.g., [49] for the detail.). Hereafter, we work in the Einstein frame summarized in the appendix A. Let us first consider the geodesic equation in the Einstein frame, where τ represents a proper time in the Einstein frame, and the Christoffel symbol Γ µ αβ is evaluated by g µν = F −1 (φ)(g µν ) E . Then taking the non-relativistic limit, we arrive at the Newtonian equation of motion [103], where Φ N is the Newtonian gravitational potential, and Φ φ is the scalar potential defined as Φ φ = Qφ/M pl . Therefore, the fifth force due to the scalar field given by Now we would like to find a profile of the scalar field. To this end, we, for simplicity, consider a point source of mass M in a medium of homogeneous background density ρ 0 and split the scalar field as the background value φ 0 and perturbation δφ. Assuming the scalar field is located at the minimum of the effective potential (A7), the background value φ 0 is then simply determined by Now we consider the scalar field equation (A5) around this background, which is given by where m eff is the effective mass defined as 26/34 and the energy density of a point source is given by ρ = M δ 3 (x). Then, the solution of this equation can be written as When m eff r 1, the fifth force is comparable to the Newtonian one, −∇Φ N , if Q ∼ O(1). On the other hand, when m eff r 1, the scalar force is short ranged, and the fifth force can be negligible due to the exponential factor, that is, the fifth force is screened. Thus, whether the Chameleon screening works or not is determined by the Compton wavelength of the scalar field λ C = m −1 eff . As one can see the schematic plot of the effective potential in Fig. B1, the shape of the effective potential depends on the density of environments. As the density ρ 0 increase, the effective mass m eff becomes larger, and the Compton wavelength λ C becomes shorter. This tells us that the Chameleon screening becomes more effective in a high density environment. Note that the effective potential also depends on pressure in the case of neutron stars, and also the profile of the scalar field including relativistic effects becomes more complicated.

C. Massless Brans-Dicke theory
In this appendix, we consider the massless Brans-Dicke theory whose solutions of NSs are discussed in [48]. The massless Brans-Dicke theory is given by the action (1) with V (φ) = 0, and the relation between the parameter Q and the original Brans-Dicke parameter is given by 2Q 2 = 1/(3 + 2ω BD ) [104]. This model is indistinguishable from general relativity with the canonical scalar field in the limit of Q → 0. Although this model for Q ∼ O(1) is not completely consistent with the solar-system experiments [23] due to a significant modification of gravitational law at all scales, we, for demonstration, consider the effect of NS cooling in 27/34 this model. Although we adopt the different EOSs (BSk19, BSk20, and BSk21) from ones used in [48], the details of the structure of NSs are qualitatively similar here. Therefore, we only plot mass-radius relations for the EOSs, BSk19, BSk20, and BSk21 in Fig. C1. Furthermore, we plot the cooling curves with massless Brans-Dicke models in Figs. C2∼7. The change of the cooling curves under scalar-tensor theory can be explained in case of Figs. C2∼C5; With small |Q| the central density with a fixed mass becomes higher. In BSk21 EOS, M DU becomes drastically lower with higher-|Q| and in case of Q = −1/ √ 2 (massless dilaton), the cooling curves cannot explain all observations though some uncertain parameters such as superfluid effect are not considered. The effect of cooling curves by changing scalar-tensor theories through |Q| can be qualitatively understood by replacingλ with 1/|Q| in Fig. 16. These figures enable us to clearly see the positive correlation between |Q| and the central density, related to M DU , from a fixed mass.
We have also calculated the cooling curves including the effect of nucleon superfluidity in the massless Brans-Dicke model as in Figs. C6 and C7. As in the case of the DU process, the |Q| dependence on the cooling curves is qualitatively the same as the tendency of 1/λ in Figs. 13 and 15. In Q = −0.10, the cooling curves are not almost changed with the GR case, but in Q = −1/ √ 6, −1/ √ 2, the cooling curves basically tend to locate in lowertemperature regions compared with the GR case. This is because the larger-|Q| values have higher-central density, that is, stronger neutrino emissions, as with the case of cooling curves without superfluidity.