On the effect of pulsar evaporation on the cooling of white dwarfs

To produce a He WD with single star evolution requires a time longer than Hubble time, so it is widely accepted that He WDs result from binary interaction through stable mass-transfer or common-envelope evolution. The short spin periods and low magnetic fileds of the NSs in LMBPs imply that LMBPs likely formed through the former process (Bhattacharya & van den Heuvel, 1991).

There is mounting evidence that mass transfer in LMXBs is likely non-conservative (e.g. Jacoby et al. 2005; Antoniadis et al. 2012). So we construct a model for the evolution of LMXBs assuming that 30 percent of the transferred matter from the low-mass donor star is accreted by the NS, and the remaining 70 percent of the material is ejected from the vicinity of the NS in the form of isotropic winds, taking away the specific angular momentum of the NS. Meanwhile, the mass accretion rate of the NS is also limited by the Eddington accretion rate. Moreover, angular momentum loss caused by magnetic braking and gravitational wave radiation is included (see Paxtion et al. 2015 for details).

When the companion star detaches from its Roche-lobe and the mass transfer terminates, we assume that the NS behaves as a MSP, and evaporation commences (van den Heuvel & van Paradijs 1988; Ruderman et al. 1989). The mass loss rate ($\dot{M}_{\rm 2,evap}$) of the companion star caused by evaporation is given as follows ( van den Heuvel & van Paradijs 1988; Stevens et al. 1992):

where $f(<1)$ measures the fraction of the incident radiation energy that is invested in overcoming the companion’s gravity. In principle it should be calculated by solving the wind’s structure considering the interplay between thermal and dynamical processes, although it is usually left as a free parameter. In addition, $v_{\rm 2,esc}$ is the escape velocity at the surface of the companion star, $R_{2}$ is the radius of the companion star, and $a$ is the orbital separation of the system. The spin-down luminosity of the radio pulsar is $L_{\rm p}=4\pi^{2}I\dot{P}_{\rm s}/P_{\rm s}^{3}$, where $I$ (taken to be $10^{45}\rm g\,cm^{2}$) is the moment of inertia of the pulsar, $P_{\rm s}$ the spin period of the pulsar and $\dot{P}_{\rm s}$ the derivative of the spin period (e.g. Liu & Li 2017):

with $B$ being the surface magnetic field of the pulsar. Substitute Eq. (2) into Eq. (1), we obtain

with

which combines the evaporation efficiency and the pulsar’s parameters into a single parameter. Here $B_{8}$ and $P_{\rm s,ms}$ are the magnetic field in units of $10^{8}$ G and the spin period in units of ms, respectively. We stop the evaporation process either when the envelope of the He core has been completely stripped or a Hubble time is reached. Meanwhile, we neglect any accretion onto the NS during the evaporation phase, and we will give a short discussion about this assumption in Sect. 4.

The evolutionary tracks presented in this paper are calculated using the stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA, version 10398; see Paxton et al. 2011, 2013, 2015, 2018). Cooling of WDs depends on the equations of state (EOS) and opacity of the atmosphere of the WDs. In MESA, different EOS and opacity tables are used on different regions in the density-temperature ($\rho-T$) diagram. In the region of $2.5\lesssim\log T(\rm K)\lesssim 7.5$ and $-10\lesssim\log\rho(\rm g/cm^{3})\lesssim 2$, OPAL EOS tables and SCVH tables, named as MESA EOS tables, are employed (Rogers & Nayfonov 2002; Saumon et al. 1995). Outside the MESA EOS region, HELM tables (Timmes & Swesty 2000) and PC tables (Potekhin & Chabrier 2010) are used (see Fig. 1 of Paxton et al. 2011). The opacity tables are constructed from the work by Cassisi et al. (2007), which combines radiation opacities with the electron conduction opacities. Radiation opacities are taken from Ferguson et al. (2005) for $2.7\leq{\rm log}T(\rm K)\leq 4.5$, and OPAL (Iglesias & Rogers 1993, 1996) for $3.75\leq{\rm log}T(\rm K)\leq 8.5$. The energy-loss rates and their derivatives from neutrinos are taken from Itoh et al. (1996). In addition, convection is treated by standard mixing length theory by default (Cox & Giuli 1968).

In our calculations, the initial masses of the neutron stars ($M_{\rm NS}^{\rm i}$) and the companions ($M_{\rm 2}^{\rm i}$) are set to be 1.4 $M_{\rm\odot}$ and 1.3 $M_{\rm\odot}$, respectively. We take the initial orbital period $P_{\rm orb}^{\rm i}$ to be longer than the bifurcation periods to guarantee the formation of He WDs for the evolution of LMXBs (Pylyser & Savonije 1988, 1989). The ratio $\alpha$ of the mixing length to the pressure scale height is set to be 2.0 (e.g. Shao & Li 2012; Wang et al. 2014).

We calculate the binary evolution with three metallicities: $Z=0.02\,({\rm Pop.\,I}),\,0.001({\rm Pop.\,II}),\,{\rm and}\,0.0001$ (Pop. III) considering the fact that some LMBPs may originate from a relatively low metallicity environment (Rivera-Sandoval et al. 2015; Cadelano et al. 2015; Cromartie et al. 2020).

One of the most important parameters in our work is $\Re$, because it determines the intensity of the evaporation power for a specific binary system. We constrain the possible range of $\Re$ from the observed parameters of radio pulsars. The $P_{\rm s}-\Re$ diagram is plotted in Fig. 1 for currently observed MSPs with $P_{\rm s}\leq 10\,\rm ms$, with $f$ set to be 1. We can roughly estimate the mean value $\overline{\Re}\sim 0.137$ and the maximum value $\Re_{\rm max}\sim 6$. In our calculations we adopt four values for $\Re=0$, 0.037, 0.37, and 1.0.

It is widely known that there is a $P_{\rm orb}\,-\,M_{\rm WD}$ relation for the final outcomes of LMXB evolution (e.g. Tauris & Savonije 1999; Lin et al. 2011). Our simulated $P_{\rm orb}\,-\,M_{\rm WD}$ relations with different values of $\Re$ and $Z$ are shown in Fig. 2, and the LMBPs in which the He WD temperatures have been derived from observation are plotted with stars. Most of the observed sources follow the relations for Pops. I and II metallicities. However, we note that three sources, PSRs J0740+6620, J1641+3627F and J1816+4510, seem to follow the relation of Pop. III, implying that they may originate from extremely low metallicity stars. In addition, the $P_{\rm orb}-M_{\rm WD}$ relation becomes slightly scattered when different values of $\Re$ are adopted, especially for systems with short orbital periods, reflecting the decrease in the WD mass (or the increase in the orbital period) due to evaporation-induced mass loss. Therefore, the evaporation process can considerably alter the evolutionary tracks of the He WDs.

In Fig. 3, the final mass ($M_{\rm Env}$) of the envelope on the He WDs is shown as a function of the final WD mass. The panels from left to right show the cases of $Z=0.02,\,\,0.001$, and $0.0001$, respectively. Beyond the left boundaries the calculation becomes non-convergent because of hitting the $\log Q$-limit (except for the case of $\Re=0$)¹¹1$\log Q=\log\rho-2\log T+12$, where $Q$ and $T$ are the density and temperature of a specific zone in the star model, respectively.. The reason for the non-convergence may be that strong evaporation dramatically changes the structure of the WD’s atmosphere. Fig. 3 shows the following features:

(1) The values of $M_{\rm Env}$ generally decrease with increasing $\Re$ for the same $Z$. This is because the mass loss rate is higher for larger value of $\Re$. However, this tendency becomes less significant with decreasing $Z$.

(2) There is more mass lost in narrower orbits (and with less massive WDs), but the amount of mass loss does not change monotonically with $M_{\rm WD}$ for the models with different values of $\Re$ (see below for the explanation).

The values of $M_{\rm Env}$ decrease monotonically with increasing WD mass in the case of $\Re=0$. This can be understood as follows. If there is no evaporation, the factors that can affect the final hydrogen content in the envelope are the initial hydrogen content in the proto-WDs and the residual hydrogen burning process. The amount of the residual hydrogen envelope at the end of the mass-transfer decreases with increasing proto-WD mass (see Fig. 10 of Istrate et al. 2016). Moreover, more massive WDs (with masses larger than $M_{\rm th}$) experience more flash cycles, which further consume the residual hydrogen (Althus et al. 2000, 2001, 2001; Benvenuto & Vito 2004, 2005; Istrate et al. 2016). Thus, more massive WDs have a less massive envelope and cool faster compared with less massive ones.

However, the amount of the decrease in $M_{\rm Env}$ does not change monotonically with $M_{\rm WD}$ for the models with $\Re=0.037$, 0.37, and 1.0. This is actually the result of the competition between the hydrogen burning at the bottom layer of the envelope and mass loss from the outer layer stripped by evaporation. The pressure at the bottom layer of the envelope reduces due to mass loss, which results in insufficient hydrogen burning. The initially massive envelope combined with insufficient hydrogen burning may lead to a more massive envelope for less massive WDs. Beyond a threshold mass, the mass loss induced by evaporation becomes small, so evaporation can hardly affect the evolution of the envelope. The reasons are as follows. Firstly, a more massive WD has a smaller radius and a larger escape velocity. Secondly, a more massive WD corresponds to a wider orbit. These factors give rise to a smaller wind mass loss rate induced by evaporation. Therefore, the evaporation process hardly affects the WD cooling in this case.

Fig. 4 shows the distribution of the WDs in the $M_{\rm WD}-T_{\rm eff}$ diagram. Same as Fig. 3, the panels from left to right show the cases with $Z=0.02$, 0.001, and 0.0001, respectively. The upper boundaries represent the effective temperature ($T_{\rm eff,max}$) at the beginning of the WD cooling process, and the lower boundaries represent the effective temperature at Hubble time. The closed left boundaries show the full cooling tracks.

It can be seen from Fig 4 that, for the models with the same $Z$, the WDs tend to cool down to lower $T_{\rm eff}$ with increasing $\Re$ for the same $M_{\rm WD}$. This is easy to understand because larger $\Re$ causes more envelope mass to be stripped and a smaller hydrogen envelope mass²²2The envelope is actually the H/He mixing atmosphere. left behind, as shown in Fig. 3. This clearly indicates that $M_{\rm Env}$ dramatically decreases with increasing $\Re$ for the same WD mass. For example, the ratio of $M_{\rm Env}$ in the cases of $\Re_{\rm max}(=1.0)$ and $\Re_{\rm min}(=0)$ can be as large as an order of magnitude for short orbit period systems, which can lead to a significant difference (about several $10^{3}\,\rm K$) in $T_{\rm eff}$. In the cases of $\Re\neq 0$, the lowest value of $T_{\rm eff}$ decreases with increasing $\Re$, same as $M_{\rm Env}$, and the decrease in $T_{\rm eff}$ also depends on $M_{\rm WD}$.

Furthermore, comparing the models with different $Z$ reveals two remarkable features. (1) The WDs with lower $Z$ can cool down to lower values of $T_{\rm eff}$ with the same $\Re$ and $M_{\rm WD}$. For example, in the case of $\Re=1$, a 0.26$M_{\rm\odot}$ WD can cool down to an effective temperature $\,\sim$ 4000 K and 3000 K for Pop. II and Pop. III metallicities respectively, but > 4000 K for Pop. I metallicities. (2) From Figs. 3 and 4, it seems that smaller $M_{\rm Env}$ does not correspond to lower $T_{\rm eff}$ when we compare the models with different metallicities. Although the WDs with lower $Z$ have a slightly more massive envelope, their temperatures are lower than those with higher $Z$ but the same WD mass. For example, for a 0.26 $M_{\rm\odot}$ WD, the final values of $M_{\rm Env}$ increase with decreasing metallicity, i.e., about $10^{-3},\,10^{-2.6}$ and $10^{-2.5}\,M_{\rm\odot}$ for Pops. I, II and III metallicities, respectively. However, the final values of $T_{\rm eff}$ decrease with decreasing metallicity.

The reasons for the above two features are as follows. Firstly, from the $P_{\rm orb}-M_{\rm WD}$ relation (see Fig. 2), the orbital separation is smaller for lower metallicity for the same WD mass, which leads to a larger wind mass-loss rate. However, the envelope at the end of the mass transfer phase is more massive in the lower metallicity models (see Fig. 10 of Istrate et al. 2016). This explains why slightly more massive envelope is left behind in the lower metallicity models, though the systems have larger wind mass-loss rate. In addition, the opacity is smaller for lower metallicity. Thus, the WDs with lower metallicity cool faster than with higher metallicity.

We compare our results with the observed He WDs in LMBPs. In Fig. 4, we can see that for the case without evaporation (red dots), only part of the observed sources with relatively high temperatures can be accounted for by the calculated results. With increasing $\Re$, more observed sources, especially those with low temperatures, are covered by the calculated results. In addition, there are several sources with extremely low temperatures which can only be explained by Pops. II or III models. We summarize that although the models without evaporation considered likely reproduce the masses and temperatures of WDs with $T_{\rm eff}>7000\,\rm K$, evaporation is probably required to explain those with $T_{\rm eff}\lesssim 5000\,\rm K$.