Thermal luminosity degeneracy of magnetized neutron stars with and without hyperon cores

The internal temperature evolves via the heat diffusion equation (Aguilera et al., 2008; Pons et al., 2009)

In Eq. (1), the heat capacity per unit volume of nucleons, leptons and hyperons is denoted by $c_{\textrm{V}}$. The internal, local temperature and the dimensionless gravitational potential are $T$ and $\Phi$ respectively, and the differential operator $\nabla$ includes the metric factors. The heat flux $\bm{F}$ reads $\bm{F}=-e^{-\Phi}\hat{k}\cdot\bm{\nabla}(e^{\Phi}T)$, where $\hat{k}$ denotes the thermal conductivity tensor (Potekhin, A. Y. & Yakovlev, D. G., 2001; Potekhin et al., 2003), while $Q_{\textrm{J}}$ and $Q_{\nu}$ denote respectively the Joule heating rate per unit volume and neutrino emissivity per unit volume.

Given the high thermal conductivity of the core, the latter becomes isothermal a few decades after the neutron star birth. The crust relaxes slower thermally, during a period that typically lasts $t\sim 10^{2}$ yr, depending on the thermal conductivity, heat capacity of the crust layers and whether neutrons are superfluid (Lattimer et al., 1994; Gnedin et al., 2001). During the thermal relaxation stage the crust temperature is higher than the core temperature, and the thermal luminosity of the star does not reflect the thermal evolution of the core. The relaxation stage ends when the “cooling wave” (Gnedin et al., 2001) propagating from the core reaches the stellar surface, and the thermal luminosity drops by orders of magnitude (depending, among other factors, on the presence of superfluid phases).

We solve the heat-diffusion equation everywhere in the stellar interior, except in the outer envelope. The typical time-scales in the envelope are shorter than in the deeper layers, requiring a smaller time-step and increasing the computational cost. Instead, we rely on an effective relation between the internal temperature at the bottom of the outer envelope $T_{\textrm{b}}$ and the surface temperature $T_{\textrm{s}}$. The latter is obtained from the $T_{\textrm{s}}$ – $T_{\textrm{b}}$ relation employed in Potekhin et al. (2015), Viganò et al. (2021) and Anzuini et al. (2021). The $T_{\textrm{s}}$ – $T_{\textrm{b}}$ relation depends on the magnetic field (Potekhin et al., 2015); see also the discussion in Anzuini et al. (2021). In the following, we assume that the outer envelope is composed of iron.

The magnetic field $\mathbf{B}$ evolves according to the magnetic induction equation, which in the crust reads

where $c$ is the speed of light, $\sigma_{e}$ is the temperature- and density-dependent electrical conductivity, $e$ is the elementary electric charge and $n_{e}$ is the electron number density. The first term is the Ohmic (dissipative) term and the second is the nonlinear Hall term.

As the temperature drops due to neutrino emission, the thermal and electric conductivities increase and become temperature-independent for sufficiently low temperatures (Aguilera et al., 2008), gradually decreasing the Ohmic dissipation rate. At the same time, the decay of the magnetic field is enhanced by the Hall term in the magnetic induction equation, because although the Hall term does not directly dissipate magnetic energy, it produces small-scale magnetic structures, where Ohmic dissipation is enhanced. Furthermore, the Hall drift tends to push the electric currents toward the crust-core boundary, where they may be dissipated more efficiently by the presence of impurities and pasta phases (Pons et al., 2013; Viganò et al., 2013), producing higher Joule heating rates. As the magnetic field evolves, the thermal conductivity along and across the magnetic field lines changes, affecting the local temperature. Hence, the magnetic evolution influences the thermal evolution and vice versa.

The evolution of the magnetic field in the core is more uncertain due to its multifluid nature and the occurrence of proton superconductivity. There may be regions with protons in the normal phase or in the superconducting phase, the latter being of type-II (Baym et al., 1969; Sedrakian & Clark, 2019) or type-I (leading to magnetic field expulsion due to the Meissner effect). Recent calculations (Wood et al., 2020) predict phase coexistence in mesoscopic regions (larger than the flux tubes, but smaller than the macroscopic length-scales), introducing several length and time-scales into the problem. When superconductivity is neglected, the typical time-scales for Ohmic dissipation and Hall advection exceed the cooling time-scales, so that the magnetic field undergoes little change in the stellar core (Elfritz et al., 2016; Dehman et al., 2020; Viganò et al., 2021). The inclusion of ambipolar diffusion could partially speed up the dynamics under certain conditions (Castillo et al., 2020). A more consistent approach including hydrodynamic effects in the superfluid/superconducting core could substantially accelerate the evolution, as recently discussed in terms of estimated time-scales by Gusakov et al. (2020) (see also references therein). Moreover, as noted above, an initial complex magnetic topology can also reduce the typical length- and time-scales, compared to the usually assumed purely dipolar fields.

In the simulations reported in this work, the induction equation in the crust includes both the Ohmic and Hall terms. In the core, only the Ohmic term is included, so that the core magnetic field is frozen, since the typical diffusion timescale in the core is larger than the typical age of isolated neutron stars studied here.

The evolution is independent of the initial internal temperature (if the latter is sufficiently high), and we typically adopt a temperature of $10^{10}$ K (Yakovlev et al., 1999; Page et al., 2004; Page et al., 2006; Potekhin et al., 2015).

The magnetic fields of neutron stars may be sustained by electric currents both in the crust and in the core. In particular, in the crust the currents can produce small-scale magnetic fields that enhance Joule heating via the Hall cascade (Gourgouliatos & Pons, 2019; Brandenburg, 2020). In this work we consider various possible initial magnetic field configurations (listed in Tables 1 and 2). The two main categories are the following. (i) Crust-confined fields. The radial magnetic field component vanishes at the crust-core interface, while the latitudinal ($B_{\theta}$) and toroidal ($B_{\phi}$) components are different from zero. (ii) Core-threading fields. At the crust-core interface the radial component of the magnetic field is $B_{r}\neq 0$, and the magnetic field lines penetrate into the core. In both cases, at the surface the magnetic field is matched continuously with the potential solution of a force-free field (i.e. the electric currents do not leak into the magnetosphere).

In Table 1 we list the crust-confined magnetic field configurations considered in this work. From a computational point of view, there is limited capability to follow numerically the rich dynamics of small-scale magnetic fields in the crust; however, the configurations studied here may reproduce typical Joule heating rates of more realistic, small-scale crustal fields. We vary the ratio of the magnetic energy stored in the toroidal component ($E^{\textrm{T}}_{\textrm{mag}}$) and the total magnetic energy ($E_{\textrm{mag}}$), as well as the number of poloidal multipoles $l_{\textrm{pol}}$. Crust-confined magnetic fields generate typically higher Joule heating rates than core-extended fields, for similar values of the total magnetic energy. As a matter of fact, in the crust-confined case all currents are forced to circulate in the crust, where the resistivity is orders of magnitude larger than in the core.

We consider two families of core-extended configurations (listed in Table 2). The first family includes the B1 and B2 configurations (studied for example by Dehman et al. (2020) and Viganò et al. (2021)), where the electric currents that sustain the magnetic field reside exclusively in the core. The corresponding Joule heating is typically lower than crust-confined configurations for two reasons. First, the resistivity in the core is low, so that Joule heating is lower than in the crust. Second, any additional heat produced via Joule heating in the core is carried away by neutrinos. The second family includes configurations with both crustal and core electric currents. To mimic the total currents in both the crust and core in neutron stars, in the C1, C2 and C1m_2,0 configurations we assume the existence of large-scale poloidal-dipolar fields threading the core, plus crustal fields. For example, in the C1 configuration, a large-scale dipolar-poloidal field threads the star, sustained by currents in the core. Additionally, there is a crustal dipolar-poloidal field, sustained by crustal currents. The azimuthal field is almost entirely confined to a torus in the core. In the C1m_2,0 configuration, there is a core-threading large-scale poloidal dipole sustained by electric currents in the core, plus an additional poloidal field with two multipoles in the crust. The core-extended configurations with crustal and core electric currents reproduce qualitatively similar Joule heating rates to the ones expected from small-scale, crustal fields, attained by the crust-confined configurations in Table 1.

Neutron star models calculated with the GM1A EoS include concentrations of nucleons, leptons and hyperons ($npe\mu Y$ matter, where $Y$ denotes hyperonic species). In particular, the EoS is fitted to modern hypernuclear data (e.g. Millener et al. (1988); Schaffner et al. (1994); Takahashi et al. (2001); Weissenborn et al. (2012), see Gusakov et al. (2014)), and predicts that only the $\Lambda$ and $\Xi^{-}$ hyperons appear in dense matter, while $\Sigma^{-}$ hyperons are absent because their potential in dense nuclear matter is repulsive. Below we list concisely the microphysics input in our simulations, such as heat capacity, neutrino emissivity, and thermal and electric conductivities, emphasizing the novelties compared to the last version of the code (Viganò et al., 2021; Anzuini et al., 2021).

• •

  Heat capacity. We include the contribution to the heat capacity of $npe\mu Y$ matter (Yakovlev et al., 1999) as well as the contribution of ions in the crustal rigid lattice.

• •

  Neutrino emission. In the core, neutrinos are produced via nucleonic and hyperon direct Urca reactions, Cooper pair breaking and formation processes, neutrino bremsstrahlung and modified Urca. We implement the in-medium corrections to the modified Urca process emission rates in Shternin et al. (2018), where the enhancement factors are calculated only for the neutron branch (Eqs. (8) and (9) in Shternin et al. (2018)). We apply the same formulae to the proton branch as well, which we interpret as upper limits to the in-medium corrections¹¹1We check that the correction factors implemented in our code reproduce the results in Shternin et al. (2018) for the BCPM EoS (Sharma et al., 2015). . Such corrections make a negligible impact on the cooling curves in our case, given the superfluid model adopted (see below) and the activation of direct Urca cooling processes. For the crust, we match the GM1A EoS with the SLy4 EoS (Douchin & Haensel, 2001), including the crustal neutrino emission processes considered in Anzuini et al. (2021). We note that neutron star models obtained with the GM1A EoS cool down fast via nucleonic direct Urca (which is active for $M\geq 1.1\ M_{\odot}$). For $M\gtrsim 1.49\ M_{\odot}$ the hyperon direct Urca involving protons and $\Lambda$ hyperons is triggered, and for $M\gtrsim 1.67\ M_{\odot}$ the hyperon direct Urca involving $\Xi^{-}$ and $\Lambda$ hyperons activates.

• •

  Electric and thermal conductivities. The conductivities depend on density and temperature and vary by orders of magnitude in the crust and core regions. In our simulations, the thermal conductivity is a tensor because of anisotropic heat transport caused by the magnetic field (Potekhin, A. Y. & Yakovlev, D. G., 2001; Potekhin et al., 2003), with components parallel and perpendicular to the magnetic field lines. We do not include contributions of nucleons, muons or hyperons to the electrical conductivities due to their lower mobility with respect to electrons.

• •

  Superfluid phases. Nucleons and hyperons can be superfluid, suppressing both the heat capacity and most channels of neutrino production (only partially compensated by the Cooper pair breaking and formation neutrino channel). As in Anzuini et al. (2021), we assume that neutrons pair in the singlet channel in the crust (“SFB” model in Ho et al. (2015)) and in the triplet channel in the core (“c” model in Page et al. (2004), see also Yanagi et al. (2020)). Protons pair in the singlet channel throughout the stellar core (“CCDK” model (Ho et al., 2015)). We use the parameters reported in Appendix A in Anzuini et al. (2021) for singlet pairing of hyperon species, reproducing similar gaps to the ones calculated by Raduta et al. (2018). We also neglect the occurrence of nucleon-hyperon superfluid phases arising from the interaction of nucleons and hyperons (Zhou et al., 2005; Nemura et al., 2009; Haidenbauer et al., 2020; Sasaki et al., 2020; Kamiya et al., 2022). The study of the magneto-thermal evolution with different hyperon superfluid gaps, obtained for example from lattice quantum chromodynamics simulations (Aoki et al., 2008, 2012; Hatsuda, 2018; Sasaki et al., 2020; Kamiya et al., 2022) is left for future work.

We study the thermal luminosity of hyperon and non-hyperon stars by comparing the magneto-thermal evolution of light-mass models ($M=1.3\ M_{\odot}$) and massive models with hyperon concentrations in the core ($M=1.8\ M_{\odot}$). Given the large nucleon and hyperon gaps, the thermal luminosity of low-mass stars with $M=1.3\ M_{\odot}$ is similar to models with masses in the range $1.1\ M_{\odot}\lesssim M\lesssim 1.4\ M_{\odot}$ (see Anzuini et al. (2021)). The same applies to the thermal luminosities of high-mass stars with $M=1.8\ M_{\odot}$, which are similar to models with masses in the range with $1.5\ M_{\odot}\lesssim M\lesssim 1.8\ M_{\odot}$ (Anzuini et al., 2021).

We emphasize that neutron stars are commonly found with masses of the order of $M\approx 1.3\ M_{\odot}$, while heavy stars are less common. Massive stars may form in merger events (if the remnant does not collapse into a black hole) (Fryer et al., 2015; Mandel & Müller, 2020; Ruiz et al., 2021), or due to matter accreted over long time-scales (up to $\approx 0.1\ M_{\odot}$ in roughly $10$ Gyr, in the optimal scenario) (Chevalier, 1989; Kiziltan et al., 2013) for example. Some mass measurements obtained via X-ray and optical observations of neutron stars in binary systems with white dwarfs fall in the range $M\gtrsim 1.6\ M_{\odot}$ (Kiziltan et al., 2013; Alsing et al., 2018). Heavy stars formed via merging or accretion may have inhomogeneous internal temperatures and complex magnetic field configurations, far from the initial conditions commonly employed in the literature of neutron star cooling. For the purpose of this work (i.e. the study of the $L_{\gamma}$ degeneracy between low-mass and high-mass models), we adopt the standard initial conditions employed by several authors (Yakovlev et al., 1999; Page et al., 2004; Potekhin & Chabrier, 2018; Raduta et al., 2018; Raduta et al., 2019), bearing in mind that the magneto-thermal evolution of heavy stars likely requires more realistic temperature and magnetic field configurations initially.

In principle, it should be possible to constrain the internal composition of a neutron star and hence the EoS of dense matter by comparing the output of cooling simulations, specifically $L_{\gamma}$ as a function of the stellar age, with optical and X-ray measurements of $L_{\gamma}$ (Viganò et al., 2013; Potekhin et al., 2020). In practice, there are several scenarios where the task is complicated by internal heating.

Consider an internal heating mechanism in the crust. When the thermal power supplied by the source is high enough, the local temperature increases and becomes almost independent of the influence of the neutrino emission processes deep in the core, so that the crust and core are thermally decoupled (Kaminker et al., 2006; Kaminker et al., 2007; Kaminker et al., 2009, 2014). We emphasize that decoupling in this context means that the crust and core temperatures evolve approximately independently. It does not mean that the crust and core are thermally insulated; there is still a heat flux between the crust and core. One key role is played by the location where the additional thermal power is supplied (Kaminker et al., 2006; Anzuini et al., 2021). If the additional heat is supplied at the bottom of the crust, it is transported via thermal conduction to the core, where it is easily lost via neutrino emission processes. As a result, the crust and core are thermally coupled, and the star cools down faster. On the other hand, if the heater deposits heat close to the outer envelope, it increases the local temperature, a fraction of the heat is transported via thermal conduction to the surface (increasing the surface temperature and hence $L_{\gamma}$), and a fraction makes its way into the core, where it is lost by neutrino emission processes. If the heating rate is sufficiently high, the local temperature in the crust is dominated by the heater, and the thermal evolution of the crust decouples from the core. In this scenario, it is challenging to ascertain whether the observed $L_{\gamma}$ is the result of fast cooling counteracted by internal heating, or if fast cooling is not active at all.

In this paper we focus on Joule heating. Joule heating rates may be sufficiently high to hide the cooling effect of nucleonic direct Urca, as discussed in Aguilera et al. (2008). Here we extend those results to include hyperon species and to consider stars with $B_{\textrm{dip}}\lesssim 10^{14}$ G and strong internal fields. Although $B_{\textrm{dip}}$ can be inferred from timing properties, there is no direct method to infer the strength of the internal field, which may decay and keep the star hot via Joule heating. Other internal heating mechanisms are also plausible but are not modelled in this paper, such as vortex creep (Shibazaki & Lamb, 1989; Page et al., 2006) or rotochemical heating (Reisenegger, 1995; Hamaguchi et al., 2019) (see Gonzalez & Reisenegger (2010) for a concise review).

We first focus on magnetars, which typically have an inferred dipolar poloidal field strengths at the polar surface in the range $10^{14}\ \textrm{G}\lesssim B_{\textrm{dip}}\lesssim 10^{15}$ G and ages $t\lesssim 10^{5}$ yr.

In Figure 1 we display the cooling curves corresponding to the A4, A5 and C2 configurations (red, blue and orange curves respectively). Magnetars are represented by red dots; stars with lower $B_{\textrm{dip}}$ are represented by black dots. Figure 1(a) studies the model with $M=1.3\ M_{\odot}$, which cools mainly via nucleonic direct Urca, and Figure 1(b) reports the cooling curves corresponding to the model with $M=1.8\ M_{\odot}$, cooling via both nucleonic and hyperonic direct Urca.

In Figure 1(a), the A5 configuration (blue, dotted curve) maintains higher $L_{\gamma}$ with respect to the A4 configuration (red, solid curve) and the C2 configuration (orange, dotted-dashed curve). The blue curve matches some of the most luminous sources with $L_{\gamma}\gtrsim 10^{35}$ erg s^-1. The model maintains $L_{\gamma}\gtrsim 10^{34}$ erg s^-1 up to $t\approx 2\times 10^{5}$ yr. The red cooling curve attains lower values of $L_{\gamma}$, and at later times ($t\gtrsim 10^{4}$ yr) it maintains a similar thermal luminosity to the blue curve ($L_{\gamma}\gtrsim 10^{34}$ erg s^-1). The C2 configuration matches lower thermal luminosities of both magnetars and stars with lower fields.

In Figure 1(b) we study the same magnetic configurations as in Figure 1(a), but for a star with $M=1.8\ M_{\odot}$ and superfluid hyperons. The blue dotted and red solid lines (A5 and A4 configurations respectively) attain similar thermal luminosities to the corresponding low-mass models in Figure 1(a) for $10^{3}\ \textrm{yr}\lesssim t\lesssim 10^{5}$ yr. For $10^{5}\ \textrm{yr}\lesssim t\lesssim 2\times 10^{5}$ yr, both curves fall below $L_{\gamma}\approx 10^{34}$ erg s^-1, contrarily to the corresponding curves in Figure 1(a). The orange curve (C2 configuration) attains lower $L_{\gamma}$ with respect to the corresponding curve in Figure 1(a) for $t\lesssim 10^{4}$ yr. At later times, the orange curves in Figure 1(a) and Figure 1(b) reach similar $L_{\gamma}$.

The comparison between the cooling curves displayed in Figure 1(a) and Figure 1(b) shows that the thermal luminosity observations and age estimates of magnetars can be explained equally by stars cooling mainly via nucleonic direct Urca emission ($M=1.3\ M_{\odot}$) and stars cooling mainly via both nucleonic and hyperonic direct Urca emission ($M=1.8\ M_{\odot}$). The decay of the strong magnetic field produces sufficient thermal power to decouple the thermal evolution of the crust and the core, and the measured $L_{\gamma}$ of magnetars is dominated by Joule heating in the crust, regardless of the neutrino emission mechanisms active in the core involving hyperons.

Similar conclusions hold if hyperons are not superfluid, and hyperon direct Urca operates without being suppressed by superfluid effects (Figure 2). The cooling curves in Figure 2(a) correspond to models with the A4 initial magnetic configuration. The red curve shows the case $M=1.3\ M_{\odot}$, and the orange dotted-dashed curve the case $M=1.8\ M_{\odot}$ with hyperon superfluidity. The blue, dotted curve is obtained for $M=1.8\ M_{\odot}$ without hyperon superfluidity. Up to $t\approx 10^{4}$ yr, the red, orange and blue cooling curves are similar, and are compatible with the $L_{\gamma}$ data of the same magnetars. The orange and blue curves are degenerate up to $t\approx 10^{4}$ yr. Only at later times Joule heating becomes weaker (e.g. $t\gtrsim 10^{4}$ yr), and one can distinguish between low-mass and high-mass stars, with or without hyperon superfluidity. Figure 2(b) reports a scenario similar to Figure 2(a), but for the A5 initial configuration. The red and orange curves ($M=1.3\ M_{\odot}$ and $M=1.8\ M_{\odot}$ models respectively, the latter including hyperon superfluidity) and the blue curve ($M=1.8\ M_{\odot}$, without hyperon superfluidity) are similar up to $t\approx 2\times 10^{4}$ yr. As in Figure 2(a), the three cooling curves are distinguishable only for $t\gtrsim 2\times 10^{4}$ yr, which exceeds most of the age estimates of the magnetar population (Viganò et al., 2013).

Some magnetar sources lie above the blue curves in Figure 1. We remind the reader that in this work we consider magnetized, iron-only outer envelopes. The data points with $10^{35}\lesssim L_{\gamma}/\textrm{erg s}^{-1}\lesssim 10^{36}$ may be explained by invoking accreted envelopes and/or higher initial magnetic fields (Potekhin, A. Y. & Yakovlev, D. G., 2001; Potekhin et al., 2003; Viganò et al., 2013). Furthermore, the thermal luminosity of these magnetars may be higher because of the presence of small hot spots, produced by the inflow of magnetospheric currents on the stellar surface. The latter process is not modeled in our simulations. We also note that our simulations do not include Joule heating in the highly resistive layer of the outer envelope, which may help to increase the thermal luminosity and match the data of the brightest magnetars.

Above we consider low-mass and high-mass models obtained with the same EoS. However, there are several more scenarios in which the cooling curves become nearly indistinguishable. Consider for example two models with the same (high) mass, the first obtained with the GM1A EoS (hosting $npe\mu Y$ matter and cooling via nucleonic and hyperonic direct Urca) and the second with a different EoS, for example hosting only $npe\mu$ matter and cooling only via nucleonic direct Urca. If both stars are born with strong magnetic fields, their magneto-thermal evolution can lead to similar observed thermal luminosities, making it hard to infer the presence of hyperons in the stellar core. We also emphasize that our results depend unavoidably on the superfluid model adopted, in particular on the nucleon energy gaps in the core. Smaller nucleon gaps lead to higher nucleon direct Urca emissivity, widening the $L_{\gamma}$ difference between light and heavy models. In this case, stronger Joule heating may be required to obtain the $L_{\gamma}$ degeneracy discussed above.

In summary, thermal luminosity data of magnetars with $L_{\gamma}\gtrsim 10^{34}$ erg s^-1 are not suitable to infer whether the core contains hyperons or not and hence infer the internal composition. We note also that it is not possible to constrain the properties of hyperon superfluid phases, since $L_{\gamma}$ is degenerate for stars with cores hosting normal or superfluid hyperons with large energy gaps. Our results show that low-mass models composed of $npe\mu$ matter and high-mass models composed of $npe\mu Y$ matter have a similar magneto-thermal evolution due to crust-core decoupling.

Magnetars are only a subset of the observable neutron star population. Timing measurements reveal that most neutron stars have inferred fields satisfying $B_{\textrm{dip}}\lesssim 10^{14}$ G. We study the evolution of such stars in Figure 3, where we display the cooling curves corresponding to the A2, A1m_2,2 and A1m_3,3 configurations (blue dotted, red solid and orange dotted-dashed curves respectively). The superfluid energy gaps are the same as in Figure 1.

Figure 3(a) reports the case $M=1.3\ M_{\odot}$. For $t\lesssim 10^{4}$ yr, the blue, red and orange lines attain similar thermal luminosities. The A1m_2,2 and A1m_3,3 configurations produce almost degenerate cooling curves, which are compatible with X-ray emitting isolated neutron stars (XINS), such as RX J$1856.5$–$3754$ and RX J$1605.3$+$3249$, and ordinary pulsars such as PSR J$0357$+$3205$ (Potekhin et al., 2020).

The case $M=1.8\ M_{\odot}$ with superfluid hyperons is presented in Figure 3(b). The blue, red and orange curves correspond to the A2, A1m_2,2 and A1m_3,3 initial magnetic configurations respectively. They differ from the curves in panel (a) for $t\lesssim 10^{5}$ yr, attaining lower values of $L_{\gamma}$. However, at later times they match the same sources as in Figure 3(a), e.g. RX J$1605.3$+$3249$ and PSR J$0357$+$3205$. It is not possible to distinguish at $t\gtrsim 10^{5}$ yr between the cooling curves of low-mass stars cooling via nucleonic direct Urca (Figure 3(a)) and high-mass stars cooling via both nucleonic and hyperonic direct Urca (Figure 3(b)). The cooling due to the appearance of hyperons is masked by the high Joule heating rate caused by the decay of the unobserved internal field. One may ask why the cooling curves in Figure 3(a) differ from the corresponding ones in Figure 3(b) for $t\lesssim 10^{5}$ yr, and attain similar values of $L_{\gamma}$ for $t\gtrsim 10^{5}$ yr. The reason is that the thermal power supplied by Joule heating for the A2, A1m_2,2 and A1m_3,3 initial configurations is insufficient to counterbalance the power lost due to nucleonic and hyperonic direct Urca emissivity when the star is relatively hot. The crust-core decoupling is incomplete. However, at later times the direct Urca emissivity is weaker due to the lower internal temperature, and Joule heating dominates the thermal evolution of the star. Consequently, the cooling curves of low-mass and high-mass stars are similar for $t\gtrsim 10^{5}$ yr.

Below we investigate further the consequences of the incomplete crust-core thermal decoupling in stars with $B_{\textrm{dip}}\lesssim 10^{14}$ G. We show that if hyperons are in the normal phase, the cooling curves of high-mass hyperon stars are clearly distinguishable from the curves of stars without hyperons in their core. Figure 4 displays the cooling curves corresponding to the A1m_3,3 and A2 magnetic configurations. In Figure 4(a) (A1m_3,3 initial configuration) the red (solid) and orange (dotted-dashed) lines correspond to models with $M=1.3\ M_{\odot}$ and $M=1.8\ M_{\odot}$ (the latter assuming that hyperons are superfluid). The blue, dotted line corresponds to a model with $M=1.8\ M_{\odot}$, but with hyperons in the normal phase. The blue curve falls below $L_{\gamma}=10^{33}$ erg s^-1 already for $t\lesssim 10^{3}$ yr, and matches the data corresponding to PSR J$0357$+$3205$ for example. There is no degeneracy between the orange and the blue curve.

Similar results are found in Figure 4(b) (A2 initial configuration), where the red and orange curves become almost degenerate for $t\gtrsim 10^{5}$ yr. However, the blue curve ($M=1.8\ M_{\odot}$ model without hyperon superfluidity) is clearly distinguishable, attaining $L_{\gamma}\lesssim 10^{32}$ erg s^-1 for $t\gtrsim 10^{5}$ yr.

In summary, we find two trends. If the star is born with a magnetar-like magnetic field with $B_{\textrm{dip}}\gtrsim 10^{14}$ G and/or a strong internal field that stores a large fraction of the total magnetic energy, the crustal temperature is regulated by Joule heating and is almost independent of neutrino cooling in the core, causing crust-core thermal decoupling. Magnetar data can be explained by both low-mass and high-mass models, regardless of the presence of hyperons. If the star has a lower field at birth ($B_{\textrm{dip}}\lesssim 10^{14}$ G) but the internal field stores most of the magnetic energy, the crust-core thermal decoupling is incomplete. In the latter scenario, if hyperons are superfluid, the magneto-thermal evolution is similar for low-mass and high-mass stars for $t\gtrsim 10^{5}$ yr. This raises the question of how to distinguish between stars that cool via nucleonic and hyperonic direct Urca heated by magnetic field decay, and stars that cool down only via nucleonic direct Urca (or even stars where direct Urca is not active at all), given that the internal field configuration and strength are unknown. On the contrary, if hyperons are not superfluid, the cooling curves are clearly distinguishable also for $t\gtrsim 10^{5}$ yr. We emphasize that we do not claim that neutron stars with low inferred values of $B_{\textrm{dip}}$ always have strong internal fields, nor that the available data of thermal emitters with low $B_{\textrm{dip}}$ must be interpreted in terms of strong internal heating. Such magnetic configurations may be characteristic of a subset of the neutron star population, rather than a common feature of thermally emitting stars.