Thermal evolution of neo-neutron stars. I: envelopes, Eddington luminosity phase and implications for GW170817

TL;DR

This paper investigates the early thermal evolution of neo-neutron stars, focusing on their luminosity, contraction, and observational signatures following supernovae, mergers, or white dwarf collapses.

Contribution

It provides the first spherically symmetric model of neo-neutron stars, analyzing their luminosity limits, energy sources, and potential observability post-merger or collapse.

Findings

Neo-neutron stars can sustain near-Eddington luminosities for a significant time.

Positrons and contraction energy play crucial roles in the star's early evolution.

Potential observability of neo-neutron stars in merger events and white dwarf collapses.

Abstract

A neo-neutron star is a hot neutron star that has just become transparent to neutrinos. In a core collapse supernova or accretion induced collapse of a white dwarf the neo-neutron star phase directly follows the proto-neutron star phase, about 30 to 60 seconds after the initial collapse. It will also be present in a binary neutron star merger in the case the "born-again" hot massive compact star does not immediately collapse into a black hole. Eddington or even super-Eddington luminosities are present for some time. A neo-neutron star produced in a core collapse supernova is not directly observable but the one produced by a binary merger, likely associated with an off-axis short gamma-ray burst, may be observable for some time as well as when produced in the accretion induced collapse of a white dwarf. We present a first step in the study of this neo-neutron star phase in a spherically symmetric configuration, thus neglecting fast rotation, and also neglecting the effect of strong magnetic fields. We put particular emphasis on determining how long the star can sustain a near-Eddington luminosity and also show the importance of positrons and contraction energy during neo-neutron star phase. We finally discuss the observational prospects for neutron star mergers triggered by LIGO and for accretion-induced collapse transients.

Equations

$$\tilde{L}=-K{\left(4\pi r^2\right)}^{2}n{\mathrm{e}}^{\varphi }\frac{\partial \tilde{T}}{\partial a},$$

$$\begin{array}{rl}& {\mathrm{e}}^{\varphi }\frac{\partial (\tilde{T}{\mathrm{e}}^{-\varphi })}{\partial t}=-\frac{1}{C_{\mathrm{V}}}({\tilde{Q}}_{\mathrm{L}}+{\tilde{Q}}_{\nu }+{\tilde{Q}}_{\mathrm{V}})\end{array}$$

$${\tilde{Q}}_{\mathrm{L}}\equiv n\frac{\partial \tilde{L}}{\partial a}$$

$${\tilde{Q}}_{\mathrm{V}}\equiv -\tilde{T}{\left(\frac{\partial P}{\partial T}\right)|}_n\frac{\partial \ln n}{\partial t}$$

$$\tilde{L}(0)=0,r(0)=0,m(0)=0,$$

$$L_{\mathrm{s}}=4\pi {\sigma }_{\mathrm{SB}}{R}^2{T}_{\mathrm{s}}^4$$

$$P_{\mathrm{s}}=\frac{2}{3}\frac{g_{\mathrm{s}}}{{\kappa }_{\mathrm{s}}}\left(1+\frac{L_{\mathrm{s}}}{L_{\mathrm{Edd}}}\right)$$

$$g_{\mathrm{s}}={\mathrm{e}}^{\lambda }GM/R^2$$

$$L_{\mathrm{Edd}}(R)=\frac{4\pi cGM{\mathrm{e}}^{\lambda }}{{\kappa }_{\mathrm{s}}}=\frac{4\pi R^{2}c{g}_{\mathrm{s}}}{{\kappa }_{\mathrm{s}}}$$

$$M=m_{\mathrm{s}}\text{and}R=r_{\mathrm{s}},$$

$${\mathrm{e}}^{\varphi (R)}={\mathrm{e}}^{-\lambda (R)}=\sqrt{1-2GM/R{c}^2}.$$

$$\frac{\mathrm{d}P}{\mathrm{d}l}=-g_{\mathrm{s}}\rho$$

$$P_{\mathrm{a}}=y_{\mathrm{a}}{g}_{\mathrm{a}}\mathrm{with}{y}_{\mathrm{a}}\equiv {\int }_{\mathrm{a}}^{\infty }\rho (r)\mathrm{d}l$$

$$T_{\mathrm{s}}(T_{\mathrm{b}},g_{\mathrm{s},14})=g_{\mathrm{s},14}^{1/4}{T}_{\mathrm{s}}(T_{\mathrm{b}},g_{\mathrm{s},14}=1)$$

$$T_{\mathrm{eff},\mathrm{Edd}}=g_{\mathrm{s}}^{1/4}{\left(\frac{c}{{\sigma }_{\mathrm{B}}{\kappa }_{\mathrm{s}}}\right)}^{1/4}$$

$$t_{\tau =1}\approx 0.5{\left(\frac{E_{\mathrm{ej}}}{10^{51}\mathrm{erg}}\right)}^{-1/2}\left(\frac{M_{\mathrm{ej}}}{5{\mathrm{M}}_{\odot }}\right)\mathrm{yrs},$$

$$t_{\tau =1}\approx 1.1{\left(\frac{E_{\mathrm{ej}}}{10^{50}\mathrm{erg}}\right)}^{-1/2}\left(\frac{M_{\mathrm{ej}}}{10^{-2}{\mathrm{M}}_{\odot }}\right)\mathrm{days},$$

$$L(t)\simeq 3\times 10^{37}{\left(t/10^4\text{s}\right)}^{-3/4}\mathrm{erg}\text{leavevmode}{\mathrm{s}}^{-1}.$$

$$\kappa =\frac{4ac{T}^3}{3K\rho }.$$

$$C_V=k_{\mathrm{B}}T\left(2\frac{Z^{\prime }}{Z}-T{\left(\frac{Z^{\prime }}{Z}\right)}^2+T\frac{Z^{\prime \prime }}{Z}\right),$$

$${\Omega }_{\mathrm{DS}}(E)=\frac{\sqrt{\pi }\exp \left(2\sqrt{a(E-E_1)}\right)}{12a^{1/4}(E-E_1{)}^{5/4}}=\sqrt{2\pi }\sigma {\rho }_{\mathrm{LD}}(E)$$

$$Z(T)=\sum _{i=0}^{i=i_{\mathrm{cf}}}{g}_i\exp \left(-\frac{E_i}{k_{\mathrm{B}}T}\right)+{\int }_{E_{\mathrm{cf}}}^{\infty }{\Omega }_{\mathrm{DS}}(E)\exp \left(-\frac{E}{k_{\mathrm{B}}T}\right)\mathrm{d}E,$$

$$T_0(\rho )=T_{\mathrm{c},0}-\Delta T{\left(\frac{\log [{\rho }_{\mathrm{c}}/\rho ]}{\log [{\rho }_{\mathrm{c}}/{\rho }_{\mathrm{b}}]}\right)}^{\gamma }$$

Full text

Neutron stars — Type II supernovae — Gamma-ray bursts — Accretion induced collapse

††journal: ApJ††software: NSCool (Page, 1989, 2016)

1 Introduction

Neutron stars are by far the most intriguing objects in the Universe. They are superdense, can be superfast rotators, may have superstrong magnetic fields, and are surrounded by the strongest gravitational fields (see, e.g. Haensel et al., 2007). They are born in core collapse supernova events (Baade & Zwicky, 1934) or in accretion induced collapse of white dwarves (Canal & Schatzman, 1976) and start their life as proto-neutron stars (Burrows & Lattimer, 1986). Moreover, a hot born-again massive neutron star may also be produced in the merging of a binary neutron star system and survive as such (Kluźniak & Ruderman, 1998), or collapse into a black hole. During the first hot phase, neutrinos are copiously produced but are trapped in the stellar interior and only escape by slowly diffusing outward. This early evolution, lasting less than a minute, has been extensively studied theoretically, in large part because a Galactic core-collapse supernova would allow us to follow it observationally through the detection of the emitted neutrinos. The subsequent phase, which we will call the neo-neutron star phase, from an age of a minute after the birth/re-birth to a few hours/days, has, however, never been carefully considered. Later phases have been the object of numerous studies (see, e.g., Yakovlev & Pethick 2004; Page, Geppert, & Weber 2006).

After the supernova, it may take decades till the ejecta become transparent to electromagnetic radiation from the central object (Bahcall et al., 1970). In the case of the supernova SN 1987A it is only recently, after more than thirty years, that credible evidence of the presence of a compact object has been found (Cigan et al., 2019). The youngest observed neutron star is the compact object in the center of the Cassiopeia A supernova remnant (Tananbaum, 1999), with an age of about 340 years (Fesen et al., 2006). It is thus doubtful we will have, in the near future, valuable observational data on the very early cooling history of a neutron star, and even less of a neo-neutron star. The neo-neutron phase is, however, the phase during which the neutron star crust is formed and it is, thus, establishing the basic structure for a large amount of neutron star phenomenology.

A complementary set-up is provided by binary neutron star mergers (Lee & Ramirez-Ruiz, 2007; Faber & Rasio, 2012). Although it is often considered that the outcome of such event would be the formation of a low-mass black hole (Eichler et al., 1989; Rezzolla et al., 2011; Murguia-Berthier et al., 2014), there is a possibility that the merged object survives as a massive neutron star (Usov, 1992; Kluźniak & Ruderman, 1998; Metzger et al., 2008). In such a scenario we would have a “born-again” neutron star, with trapped neutrinos because of its high temperature, followed by a massive neo-neutron star. This possibility is real only if the high density equation of state (EOS) is stiff enough to have a high maximum mass, $M_{\mathrm{max}}$. The maximum mass of a cold slowly- or non-rotating neutron star is at least $2\,\mathrm{M}_{\odot}$, from the masses of the pulsars PSR J1614-2230 (Demorest et al., 2010) and PSR J0348+0432 (Antoniadis et al., 2013) and possibly higher than $2.3\,\mathrm{M}_{\odot}$ from the upper value of the mass of PSR J0740+6620 (Cromartie et al., 2019). Analyses of the GW170817/GRB170817A gravitational wave/gamma-ray burst event have also provided new constraints on $M_{\mathrm{max}}$ based on the delayed collapse of the merged object into a black hole. Margalit & Metzger (2017) obtain $M_{\mathrm{max}}\leq 2.17\,\mathrm{M}_{\odot}$ (90%) and Rezzolla et al. (2018) find $2\,\mathrm{M}_{\odot}\leq M_{\mathrm{max}}\leq 2.3\,\mathrm{M}_{\odot}$ while the more detailed study of Shibata et al. (2019) conclude that $M_{\mathrm{max}}\leq 2.3\,\mathrm{M}_{\odot}$. If such is the case only mergers of binaries containing low mass neutron stars could produce a stable merged object so that our neo-neutron star description would be of interest.

As a first step, in the present paper, we consider the evolution of the outer layer of the neo-neutron star, its envelope, just after the formation of nuclei when the surface temperature is high enough to be of the order of the Eddington luminosity, i.e., of the order of $10^{38}$ erg s*-1*. An important question we tackle is the duration of a possible Eddington or super-Eddington phase, and then consider the subsequent evolution. Based on previous studies of proto-neutron star we explore the impact of different possible initial temperature/luminosity profiles in the envelope on the Eddington phase.

This paper is structured as follows. In Sect. 2 and 3 we setup the problem and present our results on the low-density regime of an Eddington envelope. In Sect. 4 we define what a neo-netron star phase is and contrast the method for study of long term cooling of isolated neutron stars versus the neo-neutron star case. Our resulted are described in Sect. 5 and their observational relevance is discussed in Sect. 6. A summary and conclusions are presented in Sect. 7. Finally in Appendix A and B we describe the physical properties of hot neutron star envelopes and our numerical scheme is detailed in Appendix C.

2 Thermal Evolution Equations

We consider a spherically symmetric problem, neglecting the effects of rotation and magnetic fields. Since the structure of the outer layers of our stars will expand or contract we employ the enclosed baryon number $a$ as a (Lagrange) radial variable instead of the circumferential radius $r$. The full set of general relativistic structure and mechanical evolution equations can be found, e.g., in Potekhin & Chabrier (2018). The thermal evolution equations we will solve are:

[TABLE]

where $\widetilde{L}=L\mathrm{e}^{2\phi}$ and $\widetilde{T}=T\mathrm{e}^{\phi}$ are red-shifted luminosity and temperature, respectively, $\mathrm{e}^{2\phi}$ being the time component of the metric. $C_{\mathrm{V}}$ is the heat capacity and $K$ the thermal conductivity. The energy sources/sinks are

[TABLE]

which gives the heat loss/injection from the luminosity gradient, $\widetilde{Q}_{\nu}=\mathrm{e}^{2\phi}Q_{\nu}$ which gives the neutrino energy loss and

[TABLE]

that will be called the contraction energy. This last term comes from the “$P\,\mathrm{d}V$” work and the volume dependent part of the internal energy and gives the gravitational and internal energy release owing to the contraction of the star during its cooling.

We follow the evolution of the star with a sequence of models in hydrostatic equilibrium, i.e., we neglect acceleration compared to gravity. We have explicitly checked that our results are consistent with this approximation.

The microphysics we employ is standard and we describe it in Appendix A and B.

3 Outer Boundary and Envelope

Boundary conditions at the center, where $a=0$, are obvious:

[TABLE]

while $P(a=0)\equiv P_{\mathrm{c}}$ is an arbitrary parameter that will determine the mass of the star. Boundary conditions at the surface are more delicate and “surface” must be properly defined. The simplest and naive condition is the “zero condition”: $P_{\mathrm{s}}=\rho_{\mathrm{s}}=T_{\mathrm{s}}=0$, and one takes $R\equiv r_{\mathrm{s}}$ and $M\equiv m_{\mathrm{s}}$. Subscript “s” refers to the quantities at the surface. However, this is too naive since $P$, $\rho$, and $T$, likely never really reach zero and, more likely, there is a smooth transition from the stellar interior to the surrounding magnetosphere (or the interstellar medium in the case of a non magnetized star). It is more appropriate, when studying the thermal evolution of the star, to define its surface as located at the photosphere, i.e., the layer where the outflowing thermal radiation is produced. We adopt the commonly used Eddington, or photospheric, condition (see, e.g., Hansen, Kawaler, & Trimble 2004) in which detailed radiative transfer (where the energy dependence of the opacity is wholly taken into account) is replaced by a diffusion approximation (where the energy dependent opacity is replaced by its Rosseland mean), i.e., the same equation (1), and the photosphere is defined as the layer where the optical depth is $2/3$. This lead to the conditions

[TABLE]

and

[TABLE]

where

[TABLE]

is the free-fall acceleration at the surface, $\mathrm{e}^{2\lambda}$ being the radial component of the metric, $\sigma_{\mathrm{SB}}$ the Stephan-Boltzmann constant, $\kappa_{\mathrm{s}}$ the Rosseland mean opacity, at the surface, and

[TABLE]

the Eddington luminosity at the stellar surface.

These are then complemented by the obvious relations that define the mass $M$ and radius $R$ of the star

[TABLE]

and the continuity of the metric coefficient with the external Schwarzschild solution

[TABLE]

It is numerically inconvenient to directly apply the outer boundary conditions of Eq. (6) and (7) and we rather apply the standard scheme of separating out an envelope (Gudmundsson, Pethick, & Epstein, 1982) as described below.

In these outer layers the equation of hydrostatic equilibrium is simply

[TABLE]

when written in terms of the proper radial length $l$ which is defined through $\mathrm{d}l\equiv\mathrm{e}^{-\lambda(r)}\mathrm{d}r$. From a given layer at $\rho_{\mathrm{a}}$ and $P_{\mathrm{a}}$, hydrostatic equilibrium can be integrated outward giving the well known classical result

[TABLE]

being the proper column density of matter above point “a” and where it has been assumed that the upper layers are sufficiently concentrated that $g$ can be considered constant. In practice we can replace $g_{\mathrm{a}}$ by $g_{\mathrm{s}}$.

3.1 High Luminosity Envelopes

The envelope is defined as the outer layers, from the surface down to a bottom layer at some pressure $P_{\mathrm{b}}$ or, equivalently, density $\rho_{\mathrm{b}}$, in which the EOS is temperature dependent and thus require a special treatment compared to the highly degenerate interior. For our present purpose we extended the previous models of Beznogov, Potekhin, & Yakovlev (2016) to higher temperatures by adding radiation pressure to the EOS of fully ionized plasma of Potekhin & Chabrier (2010) 111The corresponding Fortran code is available at http://www.ioffe.ru/astro/EIP/ and by adding the $L_{\mathrm{s}}/L_{\mathrm{Edd}}$ term to the surface condition as in Eq. (7). As we justify below in Sect. 4, we restrict ourselves to envelopes made of pure iron.

Fig. 1 shows six envelope temperature profiles labelled by $\log$ $T_{\mathrm{s}}$ and the locations of several critical loci: the melting curve (i.e., ions form a Coulomb crystal at densities above this curve), the appearance of electron-positron pairs, the onset of electron degeneracy (labelled as “$T=T_{\mathrm{F}}$”), the transition from matter to radiation pressure dominated regimes (labelled as “$P_{\mathrm{m}}=P_{\mathrm{r}}$”), and the temperature below which nuclei are formed (Lattimer & Swesty, 1991). One can also see that the profiles with the highest temperatures ($\log$ $T_{\mathrm{s}}$ [K] $>7$) cross the electron-positron pairs curve at the density about $10^{5}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$ and also the nuclear formation/dissociation line when $T>10^{10}$ K. Neither pairs nor nuclei dissociation are included in our envelope models and these parts of the profiles should not be trusted. However, as discussed in Sect. 4, we will locate our outer boundary density $\rho_{\mathrm{b}}$ at $10^{5}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$ and this regime of high density inaccurately modeled envelopes will not be actually used and is plotted here only for illustrative purpose.

Notice that in the matter dominated regime, and with opacity dominated by free-free-absorption, $\rho_{\mathrm{s}}$ increases with $T_{\mathrm{s}}$ while in the radiation dominated regime, and opacity dominated by electron scattering, the relationship is inverted. One sees from Fig. 1 that the transition between these two regimes occurs just above $T_{\mathrm{s}}\sim 10^{7}$ K.

These envelope models provide us with a relationship between the temperature at the bottom of the envelope, $T_{\mathrm{b}}$, and at its surface, $T_{\mathrm{s}}$, the so-called “$T_{\mathrm{b}}-T_{\mathrm{s}}$ relationship”: $T_{\mathrm{b}}=T_{\mathrm{b}}(T_{\mathrm{s}})$. Since energy sources and sinks are neglected within the envelope, the luminosity at its bottom, $L_{\mathrm{b}}$, is equal to the surface luminosity and thus we obtain a relationship between the two searched for solutions of Eq. (1) and (2): $L_{\mathrm{b}}=L_{\mathrm{b}}(T_{\mathrm{b}})$. This allows to replace the outer boundary condition $L_{\mathrm{s}}=L_{\mathrm{s}}(T_{\mathrm{s}})$ of Eq. (6) at $P_{\mathrm{s}}$ by a new one applied deeper at $P_{\mathrm{b}}$.

It was shown by Gudmundsson et al. (1982) that in the resulting relationship $T_{\mathrm{s}}=T_{\mathrm{s}}(T_{\mathrm{b}})$ the dependence on $M$ and $R$ is only through $g_{\mathrm{s}}$ in the form

[TABLE]

where $g_{\mathrm{s,14}}\equiv g_{\mathrm{s}}/(10^{14}\,\mathrm{cm\leavevmode\nobreak\ s^{-2}})$. We have explicitly checked that this result is still valid for our hot envelopes with a lower density inner boundary at $\rho_{\mathrm{b}}=10^{5}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$.

We notice that the approximations that lead to the Eddington boundary condition of Eq. (7) are actually self-inconsistent (Hansen et al., 2004) and the identification of a “surface” layer at temperature $T_{\mathrm{s}}$ has to be rather seen as a convenient ansatz for a more realistic atmospheric boundary condition. It is however well-known (see, e.g., Kippenhahn, Weigert, & Weiss 2012) that envelope models will converge toward the “zero condition” and the exact definition of the “surface” is not important when studying the deeper layers. We will henceforth adopt the common notation of writing the outflowing luminosity in terms of an effective temperature $T_{\mathrm{eff}}$ as $L=4\pi\sigma_{\mathrm{SB}}R^{2}T_{\mathrm{eff}}^{4}$ and use red-shifted quantities as $L^{\infty}=4\pi\sigma_{\mathrm{SB}}R_{\infty}^{2}T_{\mathrm{eff}}^{\infty\,4}$ with $L^{\infty}\equiv\mathrm{e}^{2\phi}L=\widetilde{L}_{\mathrm{s}}$, $T_{\mathrm{eff}}^{\infty}\equiv\mathrm{e}^{\phi}T_{\mathrm{eff}}$ and $R_{\infty}\equiv\mathrm{e}^{-\phi}R$, and in our case $T_{\mathrm{eff}}=T_{\mathrm{s}}$.

4 Neo-Neutron Stars

The early evolution of a newly-born, or a born-again, neutron star can be divided in two separate phases:

• •

proto-neutron star phase, $0\leq t\lesssim 30-60$ s. The star is opaque to neutrinos, $T\gg 10^{10}$ K. The chemical composition of the core slowly evolves toward the zero-temperature one as neutrinos leak out and the star’s lepton number decreases.

• •

neo-neutron star phase, $30-60\leavevmode\nobreak\ \text{s}\lesssim t\lesssim 1$ day. The star becomes transparent to neutrinos, $T\ll 10^{11}$ K. The crust is being formed.

The standard approach used in long term cooling studies (Yakovlev & Pethick 2004; Page et al. 2006) needs adjustments to study neo-neutron stars since we are now interested in very short term evolutions and very high temperatures. To be able to resolve short timescales it becomes necessary to push the outer boundary to much lower densities and we will typically use $\rho_{\mathrm{b}}=10^{5}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$ resulting in an envelope with a thermal time of the order of a second. A direct consequence of this is that the outer layers of the interior have an EOS that becomes temperature dependent. We thus distinguish three regions

• •

Outer (heat blanketing) envelope at densities $\rho_{\mathrm{s}}\leq\rho\leq\rho_{\mathrm{b}}$, treated separately in a time independent way (see Sect. 3.1). It has Eq. (6) and (7) as a surface boundary condition that defines $\rho_{\mathrm{s}}$, $P_{\mathrm{s}}$, and $T_{\mathrm{s}}$, for every given $L_{\mathrm{s}}$.

• •

Inner envelope in the regime $\rho_{\mathrm{b}}\leq\rho\leq\rho_{\mathrm{c}}$ in which the EOS is still temperature dependent and where both structure and thermal equations have to be solved simultaneously. The outer envelope provides the outer boundary condition $L_{\mathrm{b}}=L_{\mathrm{b}}(T_{\mathrm{b}})$ for $T$ and $L$ while for $P$ and $\rho$ we use Eq. (13) to write $P_{\mathrm{b}}(t)=g_{\mathrm{b}}(t)y_{\mathrm{b}}$, the time dependence coming from the contraction of this inner envelope. With the EOS $T_{\mathrm{b}}(t)$ and $P_{\mathrm{b}}(t)$ give us $\rho_{\mathrm{b}}(t)$ [and even if $P_{\mathrm{b}}$ is constant $\rho_{\mathrm{b}}$ will still change as long as $T_{\mathrm{b}}$ does]. In the absence of mass loss $y_{\mathrm{b}}$ is constant, which is what we will assume in the present work. The fact that both $P$ and $\rho$, and consequently the radius $r$, change with time in the inner envelope is the reason we prefer to use the baryon number $a$, a conserved quantity, as radial variable.

• •

Stellar interior at $\rho\geq\rho_{\mathrm{c}}$ where the EOS is temperature independent and only the thermal equations have to be solved at each time step.

Models of both proto-neutron stars (see, e.g., Burrows & Lattimer 1986) and neutron star mergers (see, e.g., Rosswog & Liebendörfer 2003) show that the star relaxes to a temperature of a few times $10^{10}$ K in less than a minute and so we will take as initial temperature $(2-3)\times 10^{10}$ K above $\rho_{\mathrm{c}}$. Taking $\rho_{\mathrm{c}}=10^{11}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$ as the inner boundary of the inner envelope is sufficient to guarantee that the stellar interior EOS can be considered as temperature independent. The microphysics we apply in the interior is the same as in long term cooling models and was described in Page et al. (2004, 2011) while the microphysics of the inner envelope is described in Appendix A.

Numerically, solving the equations of mechanical structure and thermal evolution at very high temperatures where the radiation and pair pressure are significant in the inner envelope is much more challenging that in the later evolution. We describe in Appendix C the details of our solver.

5 Results

5.1 Initial Configurations

We begin our modeling once the star is transparent to neutrino, i.e. after the $\sim$ 30 s long proto-neutron star phase in the case of a core-collapse supernova or after a similar duration after the fusion of the two stars in the case of a neutron star-neutron star merger (in the case the merged object survives instead of having collapsed into a black hole). In both cases the interior temperatures are of the order of $2-3\times 10^{10}$ K (see, e.g., Pons et al. 1999 and Rosswog & Liebendörfer 2003). We will thus take as an initial temperature $T_{0}\simeq 2.5\times 10^{10}$ K at all densities above $\rho_{\mathrm{c}}=10^{11}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$. At such densities this $T_{0}$ is just below the transition temperature where nuclei are formed in the crust (Lattimer & Swesty 1991; see, e.g., Nakazato et al. 2018 for a proto-neutron star evolution study with formation of nuclei). We then introduce a temperature gradient in the inner envelope, from $\rho_{\mathrm{c}}$ down to $\rho_{\mathrm{b}}=10^{5}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$, our initial outer boundary point. Numerical simulations of neither proto-neutron stars nor mergers resolve the temperature profile at low densities (outside the neutrinosphere) and we have thus no information about this outer layer temperature gradient. (Simulation of core-collapse supernovae do model the lower density layers but they typically follow the evolution of the system for less than a second, see, e.g., Janka 2012.) We want to start with a star emitting at the Eddington limit at its surface and this uniquely fixes the initial temperature $T_{\mathrm{b},0}$ at $\rho_{\mathrm{b}}$: the “Eddington effective temperature” $T_{\mathrm{eff,Edd}}$ is obtained, see Eq. (9), from $L_{\mathrm{Edd}}=4\pi c\,GM\mathrm{e}^{\lambda}/\kappa_{\mathrm{s}}=4\pi R^{2}cg_{s}/\kappa_{\mathrm{s}}\equiv 4\pi R^{2}\sigma_{\mathrm{B}}T_{\mathrm{eff,Edd}}^{4}$ or

[TABLE]

while envelope models relate $T_{\mathrm{eff}}$ to $T_{\mathrm{b}}$ with the same $g_{\mathrm{s}}^{1/4}$ scaling, see Eq. (14), implying that the $T_{\mathrm{b}}$ resulting in an Eddington luminosity is a unique temperature, $T_{\mathrm{b,Edd}}$, determined by the boundary density $\rho_{\mathrm{b}}$ and the chemical composition of the envelope, but independent of $M$ and $R$. For a pure iron envelope we find that $T_{\mathrm{b,Edd}}=1.07\times 10^{9}$ K at $\rho_{\mathrm{b}}=10^{5}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$.

We will consider three series of stellar models, with three different surface gravities, and implement in them different initial inner envelope luminosity or temperature profiles:

• •

models A, A′, B1, B2, B3, B4, and E, F: $M=1.4\,\mathrm{M}_{\odot}$ and $R\simeq 11.6-11.8$ km with $g_{\mathrm{s,14}}\simeq 1.6-1.7$.

• •

model C: $M=2\,\mathrm{M}_{\odot}$ and $R\simeq$ 11 km with $g_{\mathrm{s,14}}\simeq 3.2$.

• •

models D and D′: $M=0.25\,\mathrm{M}_{\odot}$ and $R\simeq 17-19$ km with $g_{\mathrm{s,14}}\simeq 0.1$.

The quoted values of the radii come from our specific choice of the core EOS from Akmal, Pandharipande, & Ravenhall (1998). The last two models D and D′ are aimed at mimicking the effect of fast rotation where centrifugal acceleration can be seen as resulting in a small effective surface gravity: a complete treatment of rotations would need a 2D code and our results are only intended to give a first approximation to the possible effects of fast rotation. In model C we do not include the fast neutrino emission by the direct Urca process (Boguta, 1981; Lattimer et al., 1991) acting deep in the inner core since it has no effect on the evolution of the outer parts of the star at early stages.

As explained in Appendix C we find it more convenient numerically to define the initial luminosity profile, $L_{0}(\rho)$, in the envelope rather that defining directly $T$. We show in Fig. 2a our choices: models A, B1, B2, B3, C, and D, have $L=L_{\mathrm{Edd}}$ at $\rho_{\mathrm{b}}$, with the value of $L_{\mathrm{Edd}}$ for their corresponding $M$, and the variation of $L$ with increasing density is constrained so that $T$ reaches $T_{0}\simeq 2.5\times 10^{10}$ K at $\rho_{\mathrm{c}}=10^{11}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$. The models E is, in contradistinction, defined by the temperature profile, following Eq. (C1), and results in a super-Eddington $L$ at the surface. Model B4 is obtained from the $L$ profile of model E scaled down so that its resulting surface luminosity is again the Eddington one (but then it cannot reach $T_{0}$ at $\rho_{\mathrm{c}}$). In Fig. 2b we plot the corresponding temperature profiles. For the reason discussed above, all models with $L(\rho_{\mathrm{b}})=L_{\mathrm{Edd}}$ start at the same $T_{\mathrm{b,0}}=T_{\mathrm{b,Edd}}=1.07\times 10^{9}$ K, while the super-Eddington model E has a higher $T_{\mathrm{b,0}}$.

Notice that since the total opacity $\kappa$ is much smaller in the inner envelope than at the surface, the local Eddington luminosity $L_{\mathrm{Edd}}(r)=4\pi r^{2}\,c\,g_{\mathrm{s}}/\kappa$ is much larger than $L_{\mathrm{Edd}}(R)$ and, hence, in all models the luminosity in the inner envelope is always below $L_{\mathrm{Edd}}(r)$.

Finally, model F represents a cold start with an initial $T_{\mathrm{eff}}$ about twice lower and hence an initial surface luminosity about 15 times below $L_{\mathrm{Edd}}(R)$. To avoid saturating the figure the initial $L$ and $T$ profiles of this model are not displayed in Fig. 2.

5.2 Evolution of a $1.4\,\mathrm{M}_{\odot}$ Star

The cooling curves resulting from our initial $L$ and $T$ profiles are presented in Fig. 3 for our $1.4\,\mathrm{M}_{\odot}$ case. One sees that all models converge, i.e., forget their initial conditions, in about $10^{4}$ s (except model A′, see below) and this initial relaxation phase is denoted as phase “1”. After this, during phase “2”, the cooling is driven by neutrino emission from the pair annihilation process and, after the knee222This knee can already be seen in the results of Nomoto & Tsuruta (1987), but with no interpretation provided, and in Page (1989)., at age $\sim 3\times 10^{5}$ s, by neutrino from the plasmon decay process, the phase “3”. The model A′ has the same initial temperature and luminosity as model A but the neutrino emission by the pair annihilation process has been arbitrarily turned off: this model confirms that pair annihilation is responsible for the evolution during phase “2”, while during phase “3” (driven by plasmon decay) model A′ converges toward model A. At an age of about one year the luminosity, and the surface temperature, reach a stagnation phase, “4”: this is the “early plateau” already well-known in neutron star cooling studies, that will last for a few decades and corresponds to thermal relaxation of the whole neutron star crust which will eventually reach thermal equilibrium with the core (see Nomoto & Tsuruta 1987, Page 1989, Lattimer et al. 1994, and Gnedin, Yakovlev, & Potekhin 2001). The shift from phase “3” to “4” is due to the inner envelope temperature dropping below the plasma temperature and the consequent exponential suppression of plasmon formation and decay: the main neutrino process available is then the very inefficient electron-ion bremsstrahlung resulting in a significant slow-down of the cooling.

In the right panel of Fig. 3 we show a close up of the early evolution of $T_{\mathrm{eff}}^{\infty}$. It is interesting to notice here that these cooling curves map their initial temperature profiles that were displayed in Fig. 2b: it results from a mapping of $T_{0}(\rho)$ into $T_{\mathrm{eff}}(t)$. This mapping is good up to time $\sim 10^{3}$ s with $\rho$ up to $10^{8}\,\mathrm{g\leavevmode\nobreak\ cm^{-3}}$: it was shown by Brown & Cumming (2009) that as long as $T_{\mathrm{eff}}^{\infty}$ is controlled by heat transport from deeper layers its value is determined mostly by the initial $T_{0}(\rho)$ at a depth whose thermal diffusion time scale to the surface is equal to the time elapsed from when this initial $T_{0}(\rho)$ was set. In our case, the mapping ends when $t-\rho$ reaches a density where the evolution is driven by neutrinos more than by heat diffusion toward the surface and this happens when approaching the phase “2” dominated by pair-annihilation neutrinos. In Fig. 4 we show the time evolution of the dominant energy term in the energy balance Eq. (2) as a function of density for model A: the details of such a plot are dependent on the assumed initial $T$ profile but that $Q_{\nu}$ eventually dominates at high densities (which, as one can see, turn out to be above $\sim 10^{8}\,\mathrm{g\leavevmode\nobreak\ cm^{-3}}$) is a simple result of the high $T$ dependence of neutrino processes and the strongly raising $T$ profile as $\rho$ increases.

In the Fig. 5 we show the evolution of the boundary radius $R_{\mathrm{b}}$ of our $1.4\,\mathrm{M}_{\odot}$ models. The different radii at early times are a direct reflection of the inner envelope temperature profiles: hotter envelopes are naturally more expanded. Excluding the model E we find contractions of $R_{\mathrm{b}}$ of the order of 50 to 100 meters. Similarly, in Fig. 6, we show the evolution of the outer boundary density $\rho_{\mathrm{b}}$ of the same models. Since $\rho_{\mathrm{b}}$ evolves with $T_{\mathrm{b}}$ in such a way that $P_{\mathrm{b}}$ remains almost constant 333 There is a small time evolution of $P$ in the outer layers because of contraction and the resulting small change in $g$, as seen, e.g., from Eq. (13)., and since $T_{\mathrm{b}}$ is instantaneously correlated with $T_{\mathrm{eff}}$ through the outer envelope, one sees that $\rho_{\mathrm{b}}$ is directly anti-correlated with $T_{\mathrm{eff}}$ shown in Fig. 3b. On the contrary, $R_{\mathrm{b}}$ results form the integral of the thickness of underlying layers and its evolution is not directly correlated with the detailed evolution of $\rho_{\mathrm{b}}$ or $T_{\mathrm{eff}}$ during phase “1”.

Considering, again, our model A in more detail as an illustrative case, we present in Fig. 7 a series of envelope temperature profiles. Since we have no mass-loss in our models, the column density $y_{\mathrm{a}}$ of any layer is constant during the evolutions and hence each layer evolves at (almost) constant pressure 3. We display in the background of Fig. 7 the pressure of the medium and a series of isobars: matter evolves along these isobars during the cooling. The profile at 600 s corresponds to the end of the early plateau during which $T_{\mathrm{eff}}$ is locked to $T_{\mathrm{eff,Edd}}$: we can divide the inner envelope in two regions, layer “a” at densities above $\sim 10^{8}\,\mathrm{g\leavevmode\nobreak\ cm^{-3}}$ where neutrino losses have had a significant effect (compare with Fig. 4) and layer “b” below $\sim 10^{8}\,\mathrm{g\leavevmode\nobreak\ cm^{-3}}$ where the temperature profile has almost not evolved. During this phase $R_{\mathrm{b}}$ has decreased by some 40 meters (Fig. 5) but this is due to the contraction of layer “a” while layer “b” has not contracted but has rather been slowly sinking, keeping its initial density and temperature profile. After this first phase $T_{\mathrm{b}}$ and the whole layer “b” are thermally connected to the temperature in layer “a”: $T_{\mathrm{b}}$, and $T_{\mathrm{eff}}$, begin to drop following the cooling of “a”. As a result the layer “b” begins to contract and $\rho_{\mathrm{b}}$ to increase as exhibited in Fig. 6. The evolution of the temperature profile from 600 up to $10^{6}$ s shows a clear difference between the region dominated by pair neutrinos, layer a2, versus plasmon neutrinos, layer a3, (separated in the figure by the (white) dotted line). As time runs the layer a2 encompasses less and less, while layer a3 encompass more and more, mass. (At these phases layer “b”, whose energetics is dominated by either $Q_{\mathrm{V}}$ or $Q_{\mathrm{L}}$, always start at densities around $10^{8}$ as seen in Fig. 4.) This different evolution of layers a2 versus a3 is easily understood by considering the difference in temperature dependence of these two neutrino processes (see Fig. 17) that result in the strongly different cooling time scales displayed in Fig. 18a. As long as part of the envelope is in the pair neutrino regime this layer a2 will drive the evolution of the outer layers and we are in phase “2” while after $\sim 10^{6}$ the layer a2 has disappeared , the cooling of the outer layers is driven by layer a3 and we entered phase “3”. It is interesting to see in Fig. 4 that at age $\sim 10^{5.5}$ s, when the cooling curve passes through the “knee”, the decrease in the pair neutrino emission is so strong that the energetics of layers that were previously dominated by $Q_{\nu}$ are now dominated by $Q_{\mathrm{L}}$ up to densities of $10^{9.5}\,\mathrm{g\leavevmode\nobreak\ cm^{-3}}$.

5.3 Robustness of our $1.4\,\mathrm{M}_{\odot}$ Star Results.

After this thorough study of our model A let us have a look at models B1 to B4. They are all based on the same two starting points, an interior initially at a temperature $T_{0}\simeq 2.5\times 10^{10}$ K as implied by studies of proto-neutron stars and binary mergers, and a surface luminosity initially at $L_{\mathrm{Edd}}$, but they have different $L$ and $T$ profiles in-between. These four scenarios have different evolutions only during the early relaxation phase “1” as clearly seen in Figs. 3, 5, and 6, a phase where $T_{\mathrm{eff}}$ is driven by the heat diffusion in the low density part of the inner envelope. However, once the cooling is controlled by neutrino emission, phases “2” and later ones, their evolutions are identical to scenario A: neutrinos are so efficient that they rapidly erase any remembrance of the initial conditions. Nevertheless, during the first $10^{3}$ s these different scenarios only span a range of $T_{\mathrm{eff}}^{\infty}$ between $1.6$ to $2.1\times 10^{7}$ K and a surface luminosity $L^{\infty}$ between $1$ to $3\times 10^{38}$ erg s*-1*. Hence, we have very similar luminosity evolutions during the first half an hour and then a basically universal evolution for the first year. Notice that the neutrino processes from either pair annihilation and plasmon decays depend only on the temperature and the electron density and do not depend on the type of nuclei present in the medium. It is only later, during phase “4” controlled by neutrino emission from electron-ion bremsstrahlung that the actual chemical composition of the medium becomes important.

On the other side, it is well known that the chemical composition of the outer envelope has a strong effect on the $T_{\mathrm{b}}-T_{\mathrm{s}}$ relationship, lighter elements having a larger thermal conductivity and resulting in higher $T_{\mathrm{s}}$ for a given $T_{\mathrm{b}}$. How large is this effect and how likely is the presence of light elements in the high temperature envelope we employ is an open question. Notice that at densities $\sim 10^{5}\,\mathrm{g\leavevmode\nobreak\ cm^{-3}}$ and temperatures $\sim 10^{9}$ K thermonuclear rates are enormous and the survival of light elements is doubtful. We intent to tackle these issues in a forthcoming work.

As a distinct initial configuration let us consider our model F which started with the same $T_{\mathrm{c},0}=2.5\times 10^{10}$ K at density $\rho_{\mathrm{c}}=10^{11}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$ but a much lower outer boundary temperature $T_{\mathrm{b},0}\simeq 0.4\times 10^{9}$ K at $\rho_{\mathrm{b}}=10^{5}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$. This model started with a clearly sub-Eddington $L^{\infty}$ as seen in Fig. 3 but after a few hundreds seconds its surface layers heat up because of the high flux coming from the inner envelope. In Fig. 8 we show the temperature profile evolution in the inner envelope: compared to model A the initial profile has no choice but to have a stronger gradient in the inner part in order to reach a lower $T_{\mathrm{b},0}$ and this is the cause of the rise in $T_{\mathrm{eff}}$ at later times when this enhanced flux reaches the surface. As in the other models, after $\sim 10^{4}$ seconds the temperature profiles have converged to the universal profiles and are indistinguishable from the ones of model A.

As described in Appendix B there is some uncertainty regarding the nuclei contribution to the specific heat, but it is only relevant at densities above $10^{10}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$ and temperatures well above $10^{9}$ K. This implies that this uncertainty has no effect on the duration of the initial Eddington luminosity phase: this phase terminates when neutrino cooling in region a2 (see Sect. 5.2) drives the evolution of $T_{\mathrm{eff}}$ and in this region the nuclei specific heat is negligible.

5.4 Evolution of High and Low Gravity Stars

In the Fig. 9 we show the cooling curves for different gravities. Model C with $g_{\mathrm{s,14}}\simeq 3.2$ turns out to be very similar to model A with $g_{\mathrm{s,14}}\simeq 1.6$: we are plotting the red-shifted luminosity and its intrinsically higher luminosity is in large part compensated by a higher red-shift. For the low gravity models, D and D′, the lower Eddington luminosity clearly shows and moreover the relaxation time is much longer: the initial relaxation phase “1” last much longer and an Eddington luminosity can be sustained for more than $10^{4}$ s versus less than $10^{3}$ s for models A and C. As a curiosity, in model D′ we have arbitrarily switched-off the contraction energy of Eq. (4): as a result during the initial relaxation phase the luminosity slightly decreases instead of staying almost constant as in model D. Nevertheless, at ages between $10^{4}$ up to $10^{7}$ s (i.e., between three hours up to three months), luminosities of these three models with very different surface gravities are still very similar and only actually differ in details (as, e.g., the time at which the “knee” occurs).

5.5 At the crossroads of different physical regimes

We finally describe our model E which has a surface luminosity twice higher that $L_{\mathrm{Edd}}$ implying a significant mass-loss. However, in the inner envelope the luminosity in this model is still sub-Eddington due to the fact that the opacity $K$ is much lower in this region than at the photosphere. We can thus still model the inner envelope within our quasi-static formalism. As seen in Fig. 3 this super-Eddington phase can last longer than the Eddington phase of our other models: about 2,600 s after which time $T_{\mathrm{eff}}$ suddenly drops.

In Fig. 10 we illustrate the evolution of this model through it $T$-profiles. Notice that at early times the low-density part of the inner envelope is clearly in the radiation/pair dominated regime. This regime corresponds in this figure to the region where the isobars are horizontal, i.e., $\rho$ independent, with $P\propto T^{4}$. In contradistinction, in all our other Eddington luminosity models the inner envelopes were always in a regime where matter made a strong contribution fo the pressure as can be seen, e.g., in Fig. 7. This super-hot model E results in a strongly puffed-up envelope, because of radiation pressure, as seen from the larger radius $R_{\mathrm{b}}$ in Fig. 5. The first four $T$-profiles displayed in Fig. 10 show a rapid contraction at the lowest densities. This contraction occurs at (nearly) constant pressure 3, hence at (nearly) constant temperature maintaining a (nearly) constant $T_{\mathrm{eff}}$, and the gravitational energy released by this contraction is used to power the super-Eddington surface luminosity. This contraction wave propagates inward until it reaches the cooling wave from the neutrino cooling propagating outward from the denser regions. After $\sim 2,600$ s further evolution along the isobars implies a significant temperature drop, the inner envelope entering a different pressure regime, and the end of the super-Eddington phase. After $\sim 10^{4}$ s this model has forgotten its initial configuration and follows the same evolution as all our other $1.4\,\mathrm{M}_{\odot}$ Eddington luminosity models.

5.6 Effects of a Strong Magnetic Field

Our neo-neutron star study assumes spherical symmetry and neglect the effect of both fast rotation and the presence of a strong magnetic field and is, thus, only a very first step in this direction. Since most realistic scenarios where a neo-neutron star may be observable are likely to produce a fast rotating and strongly magnetized star and it is imperative to estimate expectable deviations of realistic models from our idealized ones. A simple model intended to mimick fast rotation with a low gravity model was presented in subsection 5.4. The case of a strong magnetic field is more involved as it not only breaks spherical symmetry, but introduces strong anisotropies in almost every physical ingredient of our models at both the micro- and macroscopic levels, and modifies the equation of state, opacities, thermal conductivities and also possibly neutrino emission (for reviews about magnetic field effects we refer the reader to Yakovlev & Kaminker 1994 and Potekhin et al. 2015).

The heat transport anisotropy in the presence of a strong magnetic field and its effects on the thermal emission have been amply studied for a long time in neutron star envelopes (starting from, e.g., Greenstein & Hartke 1983 and Page 1995) and also in deeper layers of the crust (see, e.g., Geppert et al. 2004, 2006 and Pérez-Azorín et al. 2006). The overall effect is that regions of the surface where the magnetic field is tangential to it will be much colder that regions where the field is radial. Due to quantum effects on the thermal conductivity the hot regions with a radial field are moreover slightly hotter than they would be in the absence of a magnetic field. As a result, when integrating the outcoming flux over the whole stellar surface, one generally obtains luminosities similar to the non magnetic case with uniform surface temperature (see, e.g., Page & Sarmiento 1996 for examples with dipolar+quadrupolar magnetic field geometries). We, thus, do no expect that the magnetically induced anisotropy in heat transport results in significant deviations from our results.

A second point of interest is the duration of the early Eddington luminosity phase. It is controlled by neutrino emission from the pair annihilation process which occurs at densities above $\sim 10^{8}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$, the layer “a” in Fig. 7 and clearly illustrated by the difference between the models A and A′ in panel (a) of Fig. 3. This process is not affected by magnetic field at least up to a strength $\sim 10^{14}$ G (Kaminker & Yakovlev, 1994). It is only for much stronger fields, which become strongly quantizing even at temperature above $10^{10}$ K in this density range, that the pair annihilation process may increase the neutrino losses and reduce the duration of an Eddington phase, but unfortunately there is no complete calculation available in this density-temperature range to confirm this statement. Synchrotron neutrino emission cannot compete at very high temperatures with the pair annihilation process, at least for fields below $10^{15}$ G, and thus is not expected to affect the Eddington luminosity phase for these magnetic field intensities. However, in later phases dominated in the non magnetized cases by the plasmon decay process (phase “3” shown in panel (a) of Fig. 3) synchrotron neutrinos will increase losses when the field is above $\sim 10^{14}$ G and accelerate the cooling.

The third interesting effect of the presence of a strong magnetic field is the strong suppression of the opacity for extraordinary mode (XO) photons, $\kappa_{\mathrm{XO}}(\omega)\approx(\omega/\omega_{c})^{2}\kappa_{\mathrm{Th}}$ when $\omega\ll\omega_{c}$, where $\omega_{c}$ is the electron cyclotron frequency and $\kappa_{\mathrm{Th}}$ the electron-scattering (Thomson) opacity. As a results, when mode switching between the XO and the ordinary (O) mode is taken into account, the critical luminosity $L_{c}$ (i.e., the Eddington luminosity in the presence of the magnetic field) is increased compared to the zero field $L_{\mathrm{Edd}}$ of Eq. (9) by $L_{c}\approx(\omega_{c}/\omega)L_{\mathrm{Edd}}$ and can easily reach $10^{41-42}$ erg s*-1* (Miller, 1995). Even higher luminosities, $\sim 10^{44}$ erg s*-1*, have been observed during the magnetar Giant Flares (see, e.g., Kaspi & Beloborodov 2017) where a hot plasma is likely confined in the magnetosphere (Thompson & Duncan, 1995). How long could a strongly magnetized neo-neutron star sustain a high thermal luminosity close to $L_{c}\sim 10^{41-42}$ erg s*-1* is the crucial question. In terms of energetics, our model A sustained $L\simeq 3\times 10^{38}$ erg s*-1* for about 1000 second and emitted a total of $\sim 3\times 10^{41}$ ergs while our extreme model E kept $L\simeq 6\times 10^{38}$ erg s*-1* for about 2500 second emitting a total of $\sim 2\times 10^{42}$ ergs: this second model had a much hotter and bloated envelope that contained much more energy, both gravitational and thermal, allowing it to sustain a higher luminosity for a long time. Whether a strongly magnetized neo-neutron star envelope will be hot enough and contains sufficient energy to sustain a luminosity close to its possible $L_{c}$ for a long time (a significant fraction of an hour) is difficult to assess at present time and definitely requires a detailed calculation with all the appropriate physics included. We hope to address this issue in a future paper.

6 Observational Prospects

6.1 Core Collapse Supernova and Supernova Remnants

Neutron stars are great laboratories for studying the equation of state of nuclear-density matter. The study of supernova remnants, on the other hand, help us elucidate the composition and structure of their stellar progenitors (e.g. Lopez et al., 2011). By associating neutron stars with supernova remnants, we can obtain unique information about these systems that is unavailable when we study them separately. What is more, supernova associations provide a way to independently constrain the age of the neutron star as well as searching for former binary surviving companions.

There are, however, clear limitations that prevent us from uncovering young systems; the most noticeable being that the ejecta gas needs to be transparent to the neo-neutron star radiation. For optical (in the absence of dust) and high X-ray energies ($\gtrsim 10$ keV), electron scattering provides the main opacity (Bahcall et al., 1970). Let’s consider a cloud of gas with mass $M_{\rm ej}$ ejected from the explosion. The cloud radius expands freely as $R=v_{\rm ej}t$ where $v_{\rm ej}=\sqrt{2E_{\rm ej}/M_{\rm ej}}$ is the characteristic velocity, $t$ is the time since ejection and $E_{\rm ej}$ is the total kinetic energy. We thus expect an homogeneous envelope to become transparent after a time

[TABLE]

which is frustratingly about the duration of the Eddington luminosity phase. Here we assume $\kappa\approx 0.1$ cm2 g*-1*, which is a reasonable value for ordinary supernova material. At lower X-ray energies ($\epsilon_{\rm x}\lesssim$ 3 keV), photoionization may delay the time at which the envelope becomes transparent by an additional factor of $\sim 5({\rm\epsilon_{\rm x}/1keV})^{-3/2}$.

Despite this, the youngest neutron star we have uncovered has an age of about 340 years (Fesen et al., 2006). In a few instances, a surviving binary companion has been detected in post-explosion deep optical imaging of extragalactic supernova (Maund & Smartt, 2009; Folatelli et al., 2014). In the case of 1993J, for example, the brightness of the transient dimmed sufficiently after about a decade so that its spectrum showed the features of a massive star superimposed on the supernova (Maund et al., 2004).

It is perhaps a stinging fact that despite the expected manifestations of neo-neutron stars, one of the main issues in the field has been that most Galactic supernova remnants have no detectable central source. Fig. 11 shows the current detections and upper limits of thermal emission in nearby neutron stars with model predictions. Observational selection effects are clearly at play when uncovering young objects yet there is the possibility that a sizable fraction of massive star collapses might produce black holes rather than neutron stars, with the clearest example being W49B (Lopez et al., 2013) and the other four examples, from Kaplan et al. (2004, 2006) all plotted in Fig. 11.

6.2 Neutron Star Mergers, Short Gamma-Ray Bursts and LIGO Events

The merger of binary neutron stars and the subsequent production of a beamed, relativistic outflow is believed to trigger short gamma-ray bursts (e.g. Abbott et al., 2017a) and expel metallic, radioactive debris referred to as a kilonova (e.g. Kasen et al., 2017). The ultimate fate of the post-merger remnant remains unclear and is dependent on the mass limit for support of a hot, differentially rotating neutron star. The merged remnant can either collapse and form a low-mass black hole (Eichler et al., 1989; Rezzolla et al., 2011; Murguia-Berthier et al., 2014) or survive as a hyper-massive neutron star (Usov, 1992; Kluźniak & Ruderman, 1998; Metzger et al., 2008, 2018). In that case, the detectability of the hyper-massive remnant will depend primarily on the orientation of the merging binary.

When the relativistic jet points in the direction of the observer, the event will likely be detected as a classical short gamma-ray burst (Nakar, 2007; Lee & Ramirez-Ruiz, 2007) and the X-ray emission emanating from the surviving remnant will be buried by the luminous afterglow emission. This can be seen in Fig. 12, where we compare the X-ray luminosity of the cooling model C shown in Fig. 9 to the on-axis X-ray afterglow luminosities of a sample of short gamma-ray bursts.

Our ability to directly uncover the emission of the newly formed, hyper-massive neutron star drastically increases when the event is off-axis, as was the case for GW170817 (Abbott et al., 2017a, b). In August 2017, Coulter et al. (2017) discovered the first optical counterpart to a gravitational wave source. In this case, the cataclysmic merger of two neutron stars. This landmark discovery initiated the field of gravitational wave astronomy and enabled an exhaustive observational campaign (Abbott et al., 2017b). In Fig. 13 we show the luminosity of the X-ray counterpart to GW170817 (Margutti et al., 2018) that is seen to be significantly dimmer than the one expected from the spin-down of a highly-magnetized, rapidly rotating remnant but only slightly brighter than the X-ray luminosity predicted for the relevant cooling model C that is plotted in Fig. 9.

The constraints imposed by the afterglow observations favor the idea that GW170817 was a typical GRB jet seen off axis (Abbott et al., 2017a, b). This interpretation assumes that our line of sight is tens of degrees from the core of the jet and thus suggests that the prospects for detecting the remnant directly might be doable for future events, in particular if they are seen further away from the axis of the jet as can be seen in Fig. 14.

As discussed in the case of core collapse supernova, one of the challenges for direct detection is that the neo-neutron star is likely to be surrounded by a thick and expanding radioactive ejecta. In the case of GW170817, the optical depth is expected to be dominated by the r-process radioactively powered, kilonova ejecta (e.g. Metzger et al., 2010; Roberts et al., 2011; Kasen et al., 2017). Given the quantities derived for the neutron star merger outflow of GW170817 (e.g. Villar et al., 2017; Kasen et al., 2017; Murguia-Berthier et al., 2017; Kilpatrick et al., 2017; Ramirez-Ruiz et al., 2019), we expect the ejecta to become transparent after a time

[TABLE]

where we have used $\kappa\approx 10$ cm2 g*-1* for the much more opaque r-process rich ejecta (Barnes & Kasen, 2013). Given the low mass ejecta, double neutron star mergers appear to be a viable system for uncovering a neo-neutron star provided that the surviving remnant is stable. Alternatively, the lack of X-ray detection of the cooling signal could be used to argue in support of a collapse to a black hole.

6.3 Neo-Neutron Stars in Accretion Induced Collapse Events

Another relevant progenitor avenue for our study is the formation of a neutron star through the collapse of oxygen-neon white dwarf stars in interacting binaries (Canal & Schatzman, 1976; Miyaji et al., 1980; Canal et al., 1990; Nomoto & Kondo, 1991; Wang, 2018a). A oxygen-neon white dwarf in a binary system might be able to augment its mass near the Chandrasekhar mass leading to accretion-induced collapse by accreting steadily or dynamically (e.g. Ruiter et al., 2019; Wang, 2018b). The formation of neutron stars from interacting oxygen-neon white dwarfs in binaries is likely to be accompanied by low mass ejecta (Woosley & Baron, 1992; Dessart et al., 2006; Metzger et al., 2009; Darbha et al., 2010), which might help direct detection.

The prospects for detection of the predicted transients appear promising (Darbha et al., 2010), yet their characterization might be difficult as they might be confused with other thermal transients predicted to occur on similar timescales ($\approx$ few days) such as failed deflagrations (Livne et al., 2005) and type .Ia supernovae (Bildsten et al., 2007). Such events are, however, not expected to be accompanied by an X-ray transient. For an ejecta mass of $M_{\rm ej}=10^{-2}\,\mathrm{M}_{\odot}$, we expect these optical transients to be uncovered by upcoming surveys to distances of a few 100 Mpc (Darbha et al., 2010), which will make the X-ray characterization of the neo-neutron star doable with current space-based facilities.

7 Summary of Results and Conclusions

We have presented a detailed study of the evolution of the outer layers of a neo-neutron star. We started just after the end of the proto-neutron star phase when the internal temperature has dropped to $\sim 2.5\times 10^{10}$ K at densities $\sim 10^{11}\,\mathrm{g\leavevmode\nobreak\ cm^{-3}}$ and above. At these temperatures the nuclei in the crust have already been formed. We developed model of the outer envelope, i.e. the region from the photosphere up to densities around $\sim 10^{5-6}\,\mathrm{g\leavevmode\nobreak\ cm^{-3}}$, at luminosities close to the Eddington luminosity, in stationary state, presented in Fig. 1. Using an extension of the neutron star cooling code NSCool (Page, 1989, 2016) we then modeled the whole neutron star, but focused on the description of the evolution of the inner envelope, at densities between $10^{5}$ to $10^{11}\,\mathrm{g\leavevmode\nobreak\ cm^{-3}}$. The evolution of the surface temperature, and hence the star’s surface thermal luminosity, during this early neo-neutron star phase in controlled by the evolution of the inner envelope and is thermally decoupled from the deeper layers on such short time scales. The initial condition of temperature $\sim 2.5\times 10^{10}$ K at high densities and surface Eddington luminosity leaves some, but not much, space for variability of the temperature and luminosity in the inner envelope as shown in Fig. 2. As a result, the surface luminosity remains close to the Eddington value, i.e., above $10^{38}$ erg s*-1*, for a few thousand seconds with effective temperatures of the order of $1.5-2\times 10^{7}$ K, as presented in Fig. 3 for a 1.4 $\mathrm{M}_{\odot}$ star. After $\sim 10^{3}$ seconds the surface temperature evolution is controlled by neutrino emission, initially by pair annihilation in the inner envelope for some $10^{5}$ seconds followed by plasmon decay until it has decreased to a few millions K and reached the “early plateau”, well known from isolated neutron star cooling theory. Models with either larger or lower surface gravity have an initially different Eddington luminosity but later follow a very similar cooling trajectory during their first year of evolution, as illustrated in Fig. 9. At ages between $10^{4}$ and $10^{7}$ seconds the luminosity drop is roughly, within a factor of a few, a power law

[TABLE]

Neutron stars in the universe could have very different origins, including core-collapse supernovae, neutron star mergers, white dwarf collapses, and the so-called electron-capture supernovae that are somewhat similar to the accretion induced collapse. In the case of birth in a core collapse supernova it is very unlikely that the neo-neutron star could be observed since it takes at least a few months till the remnant could become transparent to soft X-rays. In the case of SN 1987A, it is only after more than 30 years that a first signal of the existence of a neutron star has been possibly observed (Cigan et al., 2019). However, in the case of a born-again neo-neutron star produced by the merging of two neutron star, as is expected to be the case in short GRBs and in the GW170817 event, chances of detection of the neo-neutron star are encouraging. If the merger produces a GRB and we are strongly off-axis the neo-neutron star could be detectable once the kilonova ejecta have sufficiently expanded to become transparent to soft X-rays, which should take about a day (see Equation 26). In the case of accretion induced collapse of a white dwarf the situation is similar with little ejected material.

The direct detection of a neo-neutron star can thus be aided if the formation is followed by the ejection of low mass ejecta as in the case of neutron star mergers and accretion-induced collapses of white dwarfs. Interestingly, in both scenarios the neo-neutron star is expected to be rapidly rotating and if it has a sufficiently high magnetic field, then the spin-down luminosity might prevent the detection of the cooling signature (Rosswog et al., 2003; Price & Rosswog, 2006). Our understanding of neutron star birth has come a long way since the pioneering work by Baade & Zwicky (1934) more than 8 decades ago, but these enigmatic sources continue to remain elusive, in particular at very young ages. Neutron star mergers and accretion-induced collapses of white dwarfs provide us with an exciting opportunity to study new regimes of physics and to learn what these systems were like at the earliest epochs of formation when their luminosities are near the Eddington limit. Electromagnetic and gravitational-wave observatories over the coming years offer the potential to uncover the detailed nature of these most remarkable objects.

M.B. and D.P. are partially supported by the Consejo Nacional de Ciencia y Tecnología with a CB-2014-1 grant $\#$240512. M.B. also acknowledge support from a postdoctoral fellowship from UNAM-DGAPA. The first steps of this work we financed by a UNAM-DGAPA grant #113111. E.R.-R. acknowledge support from the Heising-Simons Foundation, the Danish National Research Foundation (DNRF132) and NSF (AST-1911206). The authors thank A. Raduta for some guidance about the nuclei specific heat and J. Schwab, R. Foley and L. Lopez for illuminating discussions. 444The version of the code used in the paper is not publicly available at the present time, but the modifications are described in Appendix C

Appendix A Physical State of Matter in the Inner Envelope

As discussed in Sect. 4 in the inner envelope we are dealing with the matter at densities $\rho_{\mathrm{b}}=10^{5}<\rho<\rho_{\mathrm{c}}=10^{11}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$ and temperatures $(1-3)\times 10^{9}<T<(1-3)\times 10^{10}$ K. At this temperatures and densities one has to take into account presence of positrons and photons. Below we briefly describe the physical ingredients employed in our model and present illustrative plots that allow to better understand our results.

A.1 Equation of State

We assume the presence of 80Ni at $\rho_{\mathrm{c}}$ (following Haensel, Zdunik, & Dobaczewski 1989) and 56Fe at $\rho_{\mathrm{b}}$, while we interpolate in both $A$ and $Z$ linearly in $\log\rho$ at intermediate densities. Pressure is obtained as the sum of radiation, free gases of electrons and positrons, a free gas of nuclei plus Coulomb interaction corrections following Potekhin & Chabrier (2010). Crystallization of ions takes place when the Coulomb coupling parameter $\Gamma\equiv\left.(Ze)^{2}\middle/(a_{\mathrm{WS}}k_{\mathrm{B}}T)\right.$ reaches 175 [$a_{\mathrm{WS}}=(\left.4\pi n_{\mathrm{i}}\middle/3\right.)^{-1/3}$ is the Wigner-Seitz cell radius, $n_{\mathrm{i}}$ being the number density of ions and $Ze$ their electric charge].

A.2 Opacity and Thermal Conductivity

The thermal conductivity $K$ is taken as the sum of the electron, $K_{\mathrm{e}}$, and photon $K_{\mathrm{ph}}$ conductivities. In the inner envelope the plasma is fully ionized. Thus, there is no need to take into account the effects of partial ionization on opacity. The radiative opacity $\kappa_{\mathrm{rad}}$ consists of two terms: free-free absorption and electron scattering. The former was calculated based on the fits of Schatz et al. (1999). The latter is based on the modern fit of Poutanen (2017), which takes into account electron degeneracy and pair production (the fit handles both Thompson and Compton scattering). A correction factor of Potekhin & Yakovlev (2001) was used for adding free-free and electron-scattering opacities. The electron thermal conductivity is taken from Yakovlev & Urpin (1980) when ions are in a liquid phase and from Potekhin et al. (1999) in the solid phase.

The resulting thermal conductivity $K$ is illustrated in Fig. 15 as well as the corresponding total opacity $\kappa$ defined by

[TABLE]

This $\kappa$, which include contributions from photons and electrons, should not be confused with the more restricted (Rosseland mean) radiative opacity.

The different shape of the contour lines in panel (a) clearly exhibit different regimes that are indicated by the boxed labels:

• •

$\gamma\!:\!ee$ – conductivity dominated by photons and controlled by Thomson/Compton scattering on electrons and positrons;

• •

$\gamma\!:\!e$ – conductivity dominated by photons and controlled by Thomson scattering on electrons;

• •

$e\!:\!l$ – conductivity dominated by electrons and controlled by scattering on ions in the liquid phase;

• •

$e\!:\!cs$ – conductivity dominated by electrons and controlled by scattering on ions in a classical Coulomb solid;

• •

$e\!:\!qs$ – conductivity dominated by electrons and controlled by scattering on ions in a quantum Coulomb solid;

In the density regime considered here for the inner envelope, photon opacity is dominated by free-free absorption only in a very narrow region at the transition between photon dominated to electron dominated transport. We have a discontinuity in the electron conductivity along the melting curve which may be fictitious as argued by Baiko et al. (1998) but which is small enough and occurring at densities high enough that it has a negligible effect on our results. Notice the dramatic effect of pairs in limiting $K$ at high temperatures simply due to the strong increase of the number of scatterers, electrons and positions, with temperature in this regime.

A.3 Specific Heat

The specific heat is computed as described in Potekhin & Chabrier (2010) to which we added the contribution of radiation and pairs when present. As we are interested in temperatures up to about $1-2$ MeV we have moreover added the contribution of nuclear excitations as discussed in Appendix B.

Contour plots of the total specific heat are shown on Fig. 17. The different regimes, as indicated by the boxed labels are:

• •

$\gamma e^{+}e^{-}$ – photons and electron-positron pairs;

• •

$\mathit{electrons}$ – electrons;

• •

$\mathit{ions:l}$ – ions in a Coulomb liquid;

• •

$\mathit{ions:cs}$ – ions in a classical Coulomb solid;

• •

$\mathit{ions:qs}$ – ions in a quantum Coulomb solid.

Nuclear excitation never dominate but make a significant contribution at the highest densities ($\geq 10^{10}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$) and temperatures ($\gg 10^{9}$ K).

A.4 Neutrino Emission

In the density range of our inner envelope, $10^{5}-10^{11}\mathrm{g\leavevmode\nobreak\ cm^{-3}}$, neutrino emission is dominated by three process, in order of decreasing temperature importance: $e^{+}-e^{-}$-pair annihilation, plasmon decay and e-ion bremsstrahlung. For the first two we follow Itoh et al. (1996) and Kaminker et al. (1999) for the third one. We present in Fig. 17 contour plots of the total neutrino emissivity. Notice the dramatic change in temperature dependence when crossing the (dotted white) line from pair annihilation to plasmon decay dominance. The very strong temperature dependence of the pair annihilation process when approaching this line is due to the exponential suppression of pairs when electrons become degenerate. Similarly, when shifting from plasmon decay to electron-ion bremsstrahlung the temperature dependence of the plasmon process increase rapidly due to the exponential suppression of plasmons below the plasmon temperature.

A.5 Time-Scales

Besides the micro-physics ingredients, $\kappa$, $C_{V}$, and $Q_{\nu}$, the two evolutionary time-scales dictated by them are also very illustrative: the neutrino cooling time scale $\tau_{\nu}\equiv C_{V}/Q_{\nu}$ and the heat diffusion time scale $\tau_{\mathrm{h}}\equiv C_{V}/K$. We display both of them in Fig. 18 as they are very helpfull to understand our results.

Appendix B Nuclei specific heat

In principle, calculation of the specific heat is straightforward as it can be directly derived from the system’s partition function $Z=Z(T)$ (see, e.g., Landau & Lifshitz 1993) as

[TABLE]

where primes denote derivative with respect to the temperature $T$ and $k_{\mathrm{B}}$ is the Boltzman constant. We describe below how do we proceed to calculate the partition function.

The nucleus has a discrete excitation energy spectrum but only low lying energy levels (up to a few MeV) are known reliably from experiments. Moreover, at higher energies the density of states grows so rapidly that it is more convenient to approximate with a continuous distribution. If $\rho_{\mathrm{LD}}(E,J)$ is the density of energy levels of angular momentum $J$ at energy $E$, the “observable level density” is $\rho_{\mathrm{LD}}(E)=\sum_{J}\rho_{\mathrm{LD}}(E,J)$ while the density of sates or “true level density”, that takes into account the spin degeneracy, is $\Omega_{\mathrm{DS}}(E)=\sum_{J}(2J+1)\,\rho_{\mathrm{LD}}(E,J)$ (Gilbert & Cameron, 1965; Huizenga & Moretto, 1972). We will follow the commonly used Back-Shifted Fermi Gas (BSFG) approximation (see, e.g., von Egidy & Bucurescu 2005), an extension of the non-interactive Fermi gas model of Bethe (1936), in which

[TABLE]

where $E_{1}$ is the energy back-shift, $a$ the level density parameter, and $\sigma$ the spin cutoff, whose values are obtained by fitting experimental data. Within this approximation we can calculate $Z(T)$ as

[TABLE]

where $E_{\mathrm{cf}}=E_{i_{\mathrm{cf}}}$ is some arbitrary energy cutoff level to switch from discrete to continuous regime; $g_{i}$ is the spin degeneracy factor, which is obtained experimentally together with the energy levels $E_{i}$. In the continuous spectrum range the spin degeneracy is in principle taken into account in $\Omega_{\mathrm{DS}}(E)$. This spin degeneracy is not experimentally determined in the high excitation (continuous) regime and there are large uncertainties in its value (see, e.g., von Egidy & Bucurescu 2009) and for this reason many authors prefer to use $\rho_{\mathrm{LD}}(E)$ instead of $\Omega_{\mathrm{DS}}(E)$ in the evaluation of $Z(T)$ in Eq. (B3).

We used the procedure described above to calculate the contributions of 56Fe and 60Ni nuclei to the heat capacity. For 56Fe we used experimental data on the energy levels and spin degeneracy factors from Junde et al. (2011) and for 60Ni – from Browne & Tuli (2013).555All these data are available online at http://www.nndc.bnl.gov/chart/ The values of $E_{1}$ and $a$ were taken from the recent fits of Bucurescu & von Egidy (2015). The choice of the cutoff energy is important, so we tried two approaches: cutoff at 5 MeV and cutoff at the first energy level for which experimental value of $g$ is not known. In the latter case the cutoff energy was $\sim 3.8$ MeV for 56Fe and $\sim 3.3$ MeV for 60Ni. The heat capacity calculated using Eqs. (B1), (B2) and (B3) is per nucleus. Thus, one have to multiply it by the ion number density $n_{\mathrm{i}}=\left.\rho\middle/(Am_{u})\right.$, where $A$ is the nucleus atomic mass number and $m_{u}$ is atomic mass unit. The results are presented on Fig. 19, which shows the heat capacity, of both the nuclei and the total, as a function of temperature at two fixed densities $\rho=10^{10}$ and $10^{11}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$. Solid black curve corresponds to 56Fe, dashed black curve – to 60Ni, dot-dashed black curve – to total heat capacity including contribution from 56Fe nuclei, dot-dashed gray curve – to total heat capacity excluding nuclei contribution. Solid gray curve demonstrates 60Ni nuclei heat capacity for the energy cutoff of 5 MeV instead of $\sim 3.3$ MeV (see previous paragraph). Dashed gray curve shows 60Ni nuclei heat capacity calculated employing the level density $\rho_{\mathrm{LD}}$ instead of the density of states $\Omega_{\mathrm{DS}}$ in Eq. (B3).

From Fig. 19 one can make several conclusions. First, at low temperatures ($T\sim 10^{9}$ K) 56Fe and 60Ni heat capacities are considerably different. This is not a surprise because at these temperatures the heat capacity is governed by a first few low lying energy levels, which can differ rather noticeably even for similar nuclei. On the other hand, this does not matter much as at $T\sim 10^{9}$ K nuclei heat capacity is much less than the total heat capacity and, thus, can be neglected. Second, the maximum contribution of the nuclei heat capacity to the total heat capacity is achieved at $T\sim 10^{9.8}$ K and can be around 50% of the total heat capacity at $\rho=10^{11}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$. Nuclei heat capacity is directly proportional to the density of matter (see paragraph before previous), thus, its contribution at lower densities is lower and at sufficiently low densities ($\rho\lesssim 10^{9}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$) can be neglected at any temperature. Third, at temperatures $T\sim 10^{9.6}-10^{10}$ K nuclei heat capacity is sensitive to the particular nuclear species and to the choice of the energy cutoff. The difference can be up to $\sim 2$ times. At higher densities one enters the neutron drip regime where the heat capacity is dominated by the dripped neutron liquid (see e.g., Page & Reddy 2012) and the contribution from the nuclear excitations can again be neglected.

Appendix C Numerical Method

We base our calculations on the code NSCool (Page, 1989, 2016) with important adjustments to solve for hydrostatic equilibrium in the inner envelope in the conditions where radiation and pairs pressures are important.

The structure equations are initially solved from the center of the star down to $\rho_{\mathrm{c}}=10^{11}\leavevmode\nobreak\ \mathrm{g\leavevmode\nobreak\ cm^{-3}}$ employing the zero temperature EOS and this interior structure is not modified afterward. At densities between $\rho_{\mathrm{c}}$ and $\rho_{\mathrm{b}}$ the structure equations are solved at every time step. The thermal evolution equations (1) and (2) are solved in the whole star, i.e., from the center down to $\rho_{\mathrm{b}}$, at every time step.

So, in the inner envelope structure and thermal evolution equations have to be solved at each time step. There are several ways to do it and we had tried some of them until we have found a suitable one. The most considerable difficulty lies in the fact that the outer parts of the inner envelope are dominated by photons and electron-positron pairs. Thus, the adiabatic index is close to 4/3 and the system is close to being unstable.

In the standard long term cooling calculation scheme (see, e.g., Page 1989; Gnedin et al. 2001) thermal evolution equations are usually solved fully implicitly employing Newton-Raphson method (Henyey scheme Henyey et al. 1959). The easiest way to modify this scheme to handle neo-neutron stars is to solve structure equations separately from thermal equations at each Newton-Raphson iteration for the thermal equations. Unfortunately, this idea does not work. The thermal equation (2) [which is basically the energy conservation law] and the hydrostatic equilibrium equation [Eq. (4) of Potekhin & Chabrier (2018)] have a tendency to create oscillations in pressure, radius and temperature. This is easy to understand: if we solve them separately, some decrease in the radius will cause an increase in the temperature due to the injection of contraction energy [Eq. (2)], which will increase the pressure and cause an increase in the radius due to the hydrostatic equilibrium equation. This will, in turn, cause the temperature and pressure to drop and a decrease in the radius. Clearly, this method is prone to instability and should not be used. We implemented it and found out that it indeed resulted in diverging iterations and in oscillations.

So, to deal with this tendency to oscillate one has to solve structure and thermal equations together in a single Newton-Raphson iteration scheme. In this case the changes in the pressure, radius and temperature are coordinated with each other at each iteration and consistent solution can be obtained. As it turns out there is now need to solve all six equations together in a single Newton-Raphson scheme. Actually, it is sufficient to solve only four equations [Eqs. (1), (2) and (1), (4) from Potekhin & Chabrier (2018)] together and equations (2) and (3) from Potekhin & Chabrier (2018) can be solved separately as they do not produce any oscillations. We implemented this approach and it worked. However, the iterations converged slowly and not for all initial conditions.

So, we improved our solver further and implemented “globally convergent” Newton scheme of Press et al. (2007), which employed backtracking line searches. It improved the situation. The iterations converged faster. But still the number of iterations for early timesteps was 5–10 times bigger than in the standard long term cooling. Probably, some even more elaborate solver can improve the situation yet it looks like that not much more. The alternative is to switch from a sequence of hydrostatic equilibriums to a hydrodynamic calculation in full GR, but this is far beyond the scope of current paper.

Another numerical complication is the necessity to match the initial and boundary conditions. For the standard long term cooling it is possible to start with constant (red-shifted) temperature profile (thus, zero luminosity inside, which is inconsistent with non-zero surface luminosity, i.e. $L_{\mathrm{s}}\neq L_{\mathrm{b}}$) and the matching will occur automatically at the first time step. In our neo-neutron stars models the Henyey method would not converge at the first time step if such inconsistent surface luminosity is employed. So, we have developed a special matching procedure for the luminosity to start with the consistent initial and boundary conditions: T_{t=0}(\rho)=F\big{(}\rho,\{p_{1},p_{2},\ldots\},p_{\mathrm{match}}\big{)}, where $p_{1},p_{2},\ldots,p_{\mathrm{match}}$ are free parameters of the parametrization of an arbitrary initial temperature profile. The procedure is as follows: we fix the values of $p_{1},p_{2},\ldots$ and use Newton-Raphson method to search for the value of $p_{\mathrm{match}}$ until the initial profile satisfies the boundary condition to the desired precision [i.e, we stop when $L_{\mathrm{s}}\{T_{\mathrm{s}}(T_{\mathrm{b}})\}=L_{\mathrm{b}}$]. Typically, this takes 5-6 Newton-Raphson iterations. Employment of such a procedure means that our initial temperature profile is no longer completely arbitrary.

In particular, as an initial temperature we take a uniform value $T_{0}$ at densities above $\rho_{\mathrm{c}}$ and in the inner envelope we choose

[TABLE]

where $\Delta T=T_{\mathrm{c},0}-T_{\mathrm{b}}$, and $\gamma>0$ is a power law index. So, if $\gamma$ is 1, then $T_{0}(\rho)$ is just linear in $\log{\rho}$. We usually fixed the value of $T_{\mathrm{c},0}$ to be $2.5\times 10^{10}$ K and for various values of $\gamma$ we solved for $\Delta T$ to match the initial and boundary conditions.

Unfortunately, with the parametrization (C1) the matching occurs only at super-Eddington surface luminosities for any tested value of $\gamma$ 666Of course, with the surface boundary condition (7) we cannot have super-Eddington surface luminosity; thus, we had to extrapolate $T_{\mathrm{b}}-T_{\mathrm{s}}$ relations of Sect. 3.1 to higher temperatures to obtain the matching.. As we do not consider mass loss and stellar winds in the current work we had two options: change the initial temperature parametrization or explicitly set the initial luminosity in the inner envelope. We decided to do the latter. Setting the initial luminosity directly requires a separate step in the algorithm to solve for the initial temperature given the initial luminosity. We incorporated matching of the initial and boundary conditions in this step. In such scheme we have lost direct control over the initial temperature, but, as we show in Sect. 5, direct control over luminosity might be more useful for studying neo-neutron stars. Besides, we can still control $T_{\mathrm{b,0}}$ via $T_{\mathrm{b}}-T_{\mathrm{s}}$ relations and the fact that $L_{\mathrm{s}}(T_{\mathrm{s}}(T_{\mathrm{b}}))=L_{\mathrm{b}}$. We can also control $T_{\mathrm{c},0}$ by adjusting the initial luminosity (see details in Sect. 5). We kept the parametrization (C1) to demonstrate how a relatively small change in the initial temperature profile can considerably affect the cooling during the first $\sim 10^{4}$ s.