© ESO, 2017

1. Introduction

Neutron stars are created in supernova explosions and begin their lives surrounded by a fallback disk with a mass of order 0.1 to 0.2M_⊙ (Fryer & Heger 2000), giving a disk-to-primary mass ratio similar to protoplanetary disks close to gravitational instability. These disks have a metallicity equal to or above the Galactic average, and are thus dust-rich. The presence of a dusty disk with 10 M_⊕ of material has been proposed to explain a mid-infrared excess around a young isolated neutron star (Wang et al. 2006) and a mid-IR counterpart has been also detected around the magnetar 1E 2259+286 (Kaplan et al. 2009). These discoveries are consistent with dusty fallback disks around neutron stars ~1 Myr after formation.

The neutron star family broadly consists of four categories, according to the main mechanism that powers their emission:

• young radio pulsars (young PSRs, about ≈ 2200 objects known), powered by their rotational energy;

• millisecond radio pulsars (MSPs, ≈ 400 objects), which have accreted matter from a companion star;

• thermally emitting dim isolated neutron stars (DINSs, 7 objects known), powered by their thermal cooling and accreting from the interstellar medium;

• accreting neutron stars (ANSs, a few hundred known objects), powered by accretion of gas from a companion star.

These comprise nearly all known neutron stars, with rare exceptions like rotating radio transients, magnetars and central compact objects. These distinctions are important because different physical processes shape the environment around different neutron stars, leading to a range of effects on their planets.

The first exoplanets discovered were the three low-mass objects found around the millisecond radio pulsar B1257+12 (Wolszczan & Frail 1992). One of these is a tiny object with ~0.02 M_⊕ whereas the other two are Super-Earths of ≈4 M_⊕ (Wolszczan 1994). After B1257+12, other neutron stars have been found to host sub-stellar companions. The “diamond-planet” system PSR J1719–1438 is a millisecond pulsar surrounded by a Jupiter-mass companion thought to have formed via ablation of its donor star (Bailes et al. 2011). The system PSR B1620–26 instead is a millisecond pulsar-white dwarf binary surrounded by a Jupiter-mass planet in a 40 yr orbit (Thorsett et al. 1993). The latter system is in the globular cluster NGC 6121 (M4) and its formation is probably related to dynamical interactions that occurred in the cluster (Sigurdsson et al. 2003). It has also been recently proposed that stochastic timing variations observed around the millisecond pulsar PSR J1937+21 might be related to the presence of an asteroid belt (Shannon et al. 2013), and a similar idea has been suggested for the young PSR J0738–4042 (Brook et al. 2014). There is even a gamma ray burst (GRB 101225A) which has been proposed to be the consequence of a tidal disruption of a minor body around a neutron star (Campana et al. 2011) and some authors have explained the enigmatic phenomenon of fast radio bursts as asteroid/comet collisions with neutron stars (Geng & Huang 2015).

Neutron star planets can be first-, second- or third-generation. First generation planets would be formed in the usual manner, as a by-product of the star formation process, and would likely be ablated or unbound during stellar death. Second generation objects would form in the supernova fallback disk around a freshly-formed neutron star. Third generation planets would form from a disk consisting of a disrupted binary companion (possibly previously overflowing its Roche lobe), thought to be essential for producing millisecond pulsars such as B1257+12. The supernova explosion, the accretion from a companion for millions up to billion years that MSPs undergo, and the emission of high energy X-ray/γ-ray radiation and MeV–TeV particles (the pulsar wind) are all disruptive processes that might destroy planets or disrupt their orbits. Although, based on the fact that very few planets have been found to date, Martin et al. (2016) have recently suggested that the formation of planets around pulsars is an inherently rare phenomenon, it is clear that planet formation can happen around neutron stars. Furthermore, existing timing measurements of MSPs allow for low-mass planets on larger orbits, and the presence of planets around most non-MSP neutron stars is still essentially unconstrained.

Even if the fraction of neutron star systems that form planets is as small as the current detection fraction, the large number of neutron stars in the galaxy (~ 10⁹; Bisnovatyi-Kogan 1992) still guarantees that a relatively large number, perhaps ≳10⁷, of such planetary systems exist. We move one step forward and ask how different pulsar planets are with respect to those found around main sequence stars, whether they can retain an atmosphere, and whether they may be habitable.

Planetary systems around neutron stars need not be similar to the thousands planets that have been found around main sequence stars (e.g. Borucki et al. 2010; Batalha et al. 2013; Burke et al. 2014). For example, the habitability of a planet is generally defined in terms of equilibrium surface temperature set by the radiant energy it receives from its host star. Such stellar radiation is in first approximation blackbody radiation peaking at near-IR, optical or UV and typical habitable zones are set at a fraction of, and up to a few, astronomical units. A habitable zone of a drastically smaller size than around main sequence stars was derived for white dwarfs by Monteiro (2010), Agol (2011), Barnes & Heller (2013). They found that white dwarfs with a surface temperature of less than ≈ 10⁴K may host a habitable planet for ≈ 3 Gyr if it lies in a narrow band ~ 0.005–0.02 au from the white dwarf. In the case of neutron stars the blackbody radiation emitted might peak in X-rays, copious amounts of ionizing high energy particles might be present and almost no near-IR/optical/UV radiation is emitted and therefore different mechanisms need to be considered.

In this paper we address two problems related to neutron star planetary systems. The first is to verify whether there is any further evidence for the presence of a cloud of debris/gas around the planet-hosting PSR B1257+12, as suggested by Pavlov et al. (2007). In Sects. 2 and 3 we thus present our analysis on archival 2007 Chandra observations of this system and discuss the limits we can set on the presence of a cloud of absorbing material. The second goal is to verify whether neutron star planets can harbor an atmosphere and what are the conditions on their surfaces, which are bathed in ionizing radiation and energetic particles.

We discuss the physical conditions found in the environment around neutron stars in Sect. 4. We then discuss in Sect. 5 how such harsh conditions can affect a planetary atmosphere. Then we apply our calculations to the observed population of pulsar planets (Sect. 6) and infer in which conditions they are likely to harbor an atmosphere, which timescales are required for their evaporation and finally we touch the question of a neutron star habitable zone. In Sect. 7 we discuss potential differences in composition of neutron star planets, related to the different type of environment in which they are formed.

2. X-ray observations

We used archival Chandra observations (ObsID: 7577) carried on 2007 May 3 (start time 01:03:41 UTC) for a total of 18.31 ks. The data were collected with the Advanced CCD Imaging Spectrometer (ACIS) in very faint mode. The data were processed using the CIAO software (v.4.9) along with the calibration database (CALDB; v.4.7.3). The X-ray photons were identified with the SAOImage ds9 tool (v7.4) by selecting a circular source region of 1.5′′ in radius centered around the brightest pixel close to the coordinates of the pulsar (RA:13^h 00^m 03.0810^s, Dec:12° 40′55.875′′; see Fig. 1).

Fig. 1 PSR B1257+12  as seen by Chandra on 2007, May 3. The green circle is the 1.5′′ radius extraction region used to find the source counts. The region is centered around the white pixel, which is the brightest one in the image. The color scale refers to number of counts.

The background was obtained from an annular region centered around the source with inner radius of 3′′ and outer radius of 30′′. The photons were then extracted with the tool specextr which calculates also response and ancillary files for X-ray spectral analysis. The X-ray spectra were analyzed with the XSPEC software (v12.9.1). We also use the archival Chandra observation (ObsID 5518) taken on May 22, 2005 (20.05 ks long) to reproduce the results discussed in Pavlov et al. (2007).

3. Data analysis and results

We begin with the 2005 observations to try to reproduce the results of Pavlov et al. (2007). We find 25 photons in the 0.3–8 keV band and we fit the spectrum by using two models, a blackbody with normalized area (bbodyrad in XSPEC) and a simple power law. Both models are fitted with absorption where we keep the hydrogen column density N_H fixed at a value of 3 × 10²⁰ cm^-2 as inferred from the pulsar timing (via the dispersion measure). We find best fit parameters compatible within the 1σ error with those reported by the authors.

The 2007 observation shows a total of 21 source photons in the 0.3–8 keV energy range. Given the small number of photons, we do not group our data and use the C-statistics (Cash 1979). The results of our spectral fit are reported in Table 1.

Both a power-law and a blackbody model give results which are compatible, within the statistical uncertainties, to those obtained for the 2005 Chandra data. The 0.3–8.0 keV luminosity of the power-law model corresponds to an X-ray luminosity of $\mathrm{3.}{\mathrm{1}}_{-1.3}^{\mathrm{+}\mathrm{2.8}}\mathrm{\times }{\mathrm{10}}^{\mathrm{29}}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1}$ for an assumed source distance of 700 pc, whereas the blackbody model gives a very small emission area of the order of 2300 m², analogous to the previous findings of 2005¹. Since the projected area is very small, we explore the degeneracy of the blackbody temperature and its normalization K by plotting the confidence contours for the 68%, 90% and 99% level (see Fig. 2). Even when using the most extreme values for the normalization (at the boundary of the 99% confidence contour, K ≈ 2.2) the size of the emitting hot spot is A_{⊥ ,a} ≈ 3 × 10⁴ m², more than two orders of magnitude smaller than the predicted theoretical values of ~ 10⁷ m².

Fig. 2 Contour plot for temperature (units of keV) vs. normalization K. The contours refer to the 68% (red line), 90% (green line) and 99% (blue line) confidence levels. Degeneracy between temperature and normalization is observed, but it is still not sufficient to reconcile the small emitting area with theoretical models of millisecond pulsars.

Since we see no difference in the spectral parameters between the 2005 and 2007 observations, we try to fit the data simultaneously to increase the signal-to-noise ratio, which is valid under the assumption that the underlying spectral shape has remained truly the same in the two observations. The results of our fit are reported in Table 2.

We also try to leave the N _H value free in the fit, but the strong degeneracy with the blackbody normalization K (see Fig. 3) does allow to place only a 90% confidence upper limit of N_H < 4 × 10²¹cm^-2. To better constrain the surface of the emitting area we also try to fit a neutron star atmosphere model (nsatmos), where we fix the neutron star mass, radius, distance, and column density at 1.4 M_⊙, 10 km, 700 pc and 3 × 10²⁰ cm^-2, respectively. We obtain a good fit for a hot cap with a temperature of 1.1 × 10⁶ K and a 90% confidence level upper limit on the emitting area of less than 4 × 10⁶ m². The results of this model are shown in Table 2.

Fig. 3 Contour plot for the joint 2005–2007 fit of the column density NH vs. blackbody normalization. The contours and symbols are the same as in Fig. 2. The degeneracy between column density and blackbody normalization is evident.

Following Pavlov et al. (2007, see their Eq. (2)), we can set an upper limit on the putative cloud of absorbing material around the pulsar from our constraints on the hydrogen column density. We assume that the cloud has a spherical shape and that the volume is truncated at around the location of the outermost planet (≈ 0.5 au). This gives an upper limit on the mass of the cloud of: ${\mathit{M}}_{\mathrm{cl}}\mathit{<}\mathrm{1.5}\mathrm{\times }{\mathrm{10}}^{-4}{\left(\frac{\mathrm{ℒ}}{\mathrm{0.5}\hspace{0.17em}\mathrm{au}}\right)}^{\mathrm{2}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\oplus }}\mathit{,}$(1)at the 90% confidence level, where ℒ is the size of the cloud in units of au. For sizes of order 1 au and above, this is similar to dusty debris disk masses around main-sequence stars (Wyatt 2008). Thus, if there is a dust cloud around B1257+12, it is comfortably in the mass regime of debris disks.

4. Neutron star family

4.1. Radio pulsars

Young PSRs are relatively young (0–1 Gyr) neutron stars formed after a supernova explosion. Their spin period is usually in the 0.01–10 s range and they possess a strong magnetic field of the order of 10¹⁰–10¹³ G. The rotation of their strong dipolar magnetic field causes the emission of low frequency radiation and a relativistic particle wind with a loss of energy from the system with a consequent spin down. A pulsar loses energy at a rate: ${\mathit{L}}_{\mathrm{sd}}\mathrm{=}\mathit{I\omega }\mathit{\omega ̇}\mathrm{=}\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}\mathit{I}\frac{\mathit{Ṗ}}{{\mathit{P}}^{\mathrm{3}}}\mathit{,}$(2)where L_sd is the spin-down power, I is the moment of inertia of the neutron star, ω is the neutron star angular frequency, its spin down rate and P and Ṗ are the spin period and the spin period derivative.

MSPs are old neutron stars that are formed in binaries and have been spun up to millisecond spin periods by accretion processes (see reviews Bhattacharya & van den Heuvel 1991; Patruno & Watts 2012). These objects can then either stop accreting and become binary MSPs or ablate completely their companion and become isolated MSPs. Some formation scenarios for the pulsar planets around PSR B1257+12  discuss indeed their formation after the millisecond pulsar ablated its companion leaving a post-ablation disk of ~ 0.1 M_⊙ (Stevens et al. 1992; see also Podsiadlowski 1993; Veras 2016 for reviews of different models proposed). These pulsars are among the most precise clocks in the Universe (e.g., Taylor 1991) and timing their pulse arrival times allows the measurement of their spin and orbital variations with exquisite precision. Thanks to this stability, small perturbations that emerge in the timing analysis of MSPs can highlight the presence of planets or other small bodies which would be impossible to detect otherwise.

A large fraction of MSPs (~ 80%) are found in binaries with a companion star that can be another neutron star, a main sequence star, a brown dwarf or a white dwarf. Beside their spin period, MSPs differ from young PSRs also because their dipolar magnetic field is of the order of 10⁸–10⁹ G, whereas the general physics of the pulsar wind and radiative emission is believed to be similar. Their weak magnetic field means also that they are spinning down at a significantly slower pace when compared to the young pulsars, and indeed most MSPs are expected to emit radio pulsations for at least a few more billion years. Most of the energy lost is carried away by an intense flux of low radio frequency radiation and a relativistic particle wind. A small part of the energy budget (~ 0.01–1%; Pavlov et al. 2007) is instead converted into high energy radiation consisting mostly of X-rays.

A large fraction of the spin-down power of pulsars is quickly converted into an energetic flow of relativistic particles (pulsar wind), composed mostly by electron-positron pairs. The magnetic field of the pulsar is instead dipolar (scaling with distance as r^-3) up to the so-called light cylinder, which is the location where the magnetic field lines of the pulsar have an angular velocity equal to the speed of light. After this point the field lines open and wind up creating a toroidal field whose strength scales as r^-1. The pulsar wind exerts a ram pressure p_w that is a function of the total spin-down power of the pulsar (e.g., Harding & Gaisser 1990): ${\mathit{p}}_{\mathit{w}}\mathrm{=}\mathit{\xi }\hspace{0.17em}{\mathit{L}}_{\mathrm{sd}}\mathit{/}\mathrm{(}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\mathit{c}\mathrm{)}\mathit{,}$(3)where c is the speed of light and ξ is an efficiency factor that accounts for the fraction of spin-down power that is effectively transformed into relativistic wind power. Such parameter is different than 1 because a fraction 0.01–0.8 of the pulsar energy output can be converted into high energy gamma ray photons, a discovery recently made thanks to the Fermi gamma-ray telescope (see e.g., Abdo et al. 2013). A small fraction of power is instead carried away in form of X-rays and an even smaller fraction is transformed into coherent radio waves detectable as radio pulses from Earth, comprising only a tiny percentage of the energy budget. Therefore realistic values for ξ are in the 0.2–0.99 range. Such emission processes proceed until the pulsar crosses the so-called “death-line”, i.e., a region in the P−Ṗ diagram where pulsars turn off since their spin-down power is insufficient to sustain the pulsar wind and the production of electron-positron pairs.

4.2. Isolated neutron stars

DINSs are also known as thermally emitting neutron stars since their primary photon energy output comes from the thermal emission from their surface. DINSs however, emit energy also in form of neutrinos, with their cooling rate that can be described as: $\frac{\mathrm{d}\mathit{U}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}{\mathit{L}}_{\mathit{\gamma }}\mathrm{-}{\mathit{L}}_{\mathit{\nu }}\mathit{,}$(4)where U is the internal energy of the neutron star and L_γ and L_ν are the photon and neutrino power. After a neutron star is formed, the neutrino energy loss is the dominant mechanism, whereas photon cooling becomes dominant at a later stage. To understand the typical timescale of this transition one needs to know the exact interior composition of neutron stars since the neutrino production strongly depends on the exact particle interactions that occur in the neutron star core. For a number of plausible models, the neutrino cooling dominates in the first 100 000 yr (e.g., Page et al. 2006).

A good approximation for L_γ is given by blackbody emission: ${\mathit{L}}_{\mathit{\gamma }}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{R}}^{\mathrm{2}}{\mathit{\sigma }}_{\hspace{0.17em}\mathrm{SB}}\hspace{0.17em}{\mathit{T}}_{\mathrm{e}}^{\mathrm{4}}\mathit{,}$(5)where T_e is the temperature of the neutron star envelope and σ_SB the Stefan-Boltzmann constant. If the neutron star is in the photon cooling era, then L_ν ≈ 0 and the integration of Eq. (4) gives (Tsuruta 1979): $\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{=}\mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{3}}{\mathit{\alpha }}^{\mathrm{2}}{\left(\frac{{\mathit{M}}_{\mathrm{NS}}}{{\mathit{M}}_{\mathrm{\odot }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{T}}_{\mathrm{e}\mathit{,}\mathrm{7}}^{-2}\mathrm{\left(}\mathit{f}\mathrm{\right)}\left[\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{T}}_{\mathrm{e}\mathit{,}\mathrm{7}}\mathrm{\left(}\mathit{f}\mathrm{\right)}}{{\mathit{T}}_{\mathrm{e}\mathit{,}\mathrm{7}}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\right)}^{\mathrm{2}}\right]\mathit{,}$(6)where the indices f and i indicate the final and initial temperatures and M_NS is the neutron star mass (in the ≈ 1.2–2.0 M_⊙ range). The parameter α ≈ 0.1–1 relates the internal (core) temperature to the envelope temperature. The T ∝ t^{− 1/2} dependence of the temperature on time means that the surface temperature halves as the age of the neutron star quadruples. As an example, we start with a neutron star that has just finished to cool mainly via neutrino emission, and we assume it has an age of 0.1 Myr. If the surface temperature of this object is of the order of 10⁵ K (as observed in some isolated neutron stars of similar age, Kaplan 2008) then its temperature would have dropped to ≈ 6000 K when its age is ~ 1 Gyr. The thermal luminosity of the 1-Gyr old neutron star would then be L ≈ 10²⁴ erg s^-1.

The minimum distance that a planet can have from the neutron star is the so-called tidal radius, which sets the region where an object would be disrupted by the tidal forces generated by the neutron star. This radius is: ${\mathit{R}}_{\mathrm{t}}\mathrm{=}{\left(\frac{{\mathit{M}}_{\mathrm{NS}}}{{\mathit{M}}_{\hspace{0.17em}\mathrm{p}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}\mathit{,}$(7)where M_p and R_p and the mass and radius of the planet. The fraction of thermal radiation incident onto such planet at the tidal radius distance is: $\mathit{f}\mathrm{=}\frac{\mathit{\pi }\hspace{0.17em}{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{R}}_{\mathrm{t}}^{\mathrm{2}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}{\left(\frac{{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}}{{\mathit{R}}_{\mathrm{t}}}\right)}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}{\left(\frac{{\mathit{M}}_{\hspace{0.17em}\mathrm{p}}}{{\mathit{M}}_{\mathrm{NS}}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{·}$(8)For an Earth-like planet f ≈ 10^-5 and thus the amount of energy incident on the planet would be similar to that received by the Earth from the Sun only if the neutron star temperature is slightly larger than 10⁶ K. However, for a temperature of the neutron star of 6000 K the amount of power received by the planet is equivalent to what Earth would receive if it were at a distance of more than 100 au from the Sun. Under these particular conditions planets will be frozen worlds for most of their lifetime.

4.3. Extra sources of heat and X-ray radiation

Isolated neutron stars are directly exposed to the interstellar medium and it is expected that all of them would accrete some of this material. Such accretion process generates extra power due to the conversion of the accreted gas rest mass into energy, with a typical efficiency of the order of 10–20%. This so-called Bondi-Hoyle accretion process should be continuous and might be the main source of power for these type of systems. Toropina et al. (2012) performed MHD calculations to estimate the Bondi-Hoyle accretion luminosity of an isolated neutron star with a magnetic field of 10¹²G (typical for PSRs) and found that ${\mathit{L}}_{\mathrm{acc}}\mathrm{=}\frac{\mathit{G}\hspace{0.17em}{\mathit{M}}_{\mathrm{NS}}\mathit{Ṁ}}{{\mathit{R}}_{\hspace{0.17em}\mathrm{NS}}}\mathrm{\approx }\mathrm{6}\mathrm{\times }{\mathrm{10}}^{\mathrm{29}}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1}$ (although we caution that different calcualtions show a much lower X-ray luminosity for the same process, see e.g. Popov et al. 2015). If the neutron star has a magnetosphere, as expected from most neutron stars, then this power might be mainly emitted as X-ray photons pulsating at the rotational rate of the neutron star.

In the case of pulsars, the X-ray radiation might also be formed as a consequence of non-thermal processes originating in the pulsar magnetosphere that heat up the polar caps (and which are unrelated to accretion processes). To date few tens of MSPs have been observed by the Chandra and XMM-Newton observatories to emit 0.1–10 keV X-ray radiation with a power of the order of 10²⁹–10³¹ erg s^-1 (Pavlov et al. 2007). A few outliers with a power of 10³²–10³³ erg s^-1 are also observed, but they constitute a small minority of the sample. Even if the X-ray radiation represents a negligible fraction of the total power emitted by pulsars, its effect on a planetary atmosphere might be severe given that a substantial amount of energy can be deposited into the outer layers of the atmosphere. Indeed the X-rays interact with the electron K-shell of atoms, producing ions and energetic electrons that drive the heating of the atmosphere (Cecchi-Pestellini et al. 2006). The fast electrons thus generated induce secondary ionizations that contribute substantially to the atmospheric chemistry and can provide a substantial contribution to the atmospheric heating.

5. Effect of irradiation on neutron star planets

As we have discussed above, neutron stars can irradiate a planet via two main channels: relativistic winds (when the neutron star is active as a pulsar) and/or via high energy X-ray radiation produced either during the pulsar phase or via accretion from the interstellar medium for DINSs.

5.1. Effect of pulsar winds

For simplicity we begin by assuming that a planet around a young pulsar possesses an isothermal atmosphere dominated by gas pressure and no magnetosphere. The gas pressure is: ${\mathit{p}}_{\mathrm{g}}\mathrm{=}\mathit{\rho }\hspace{0.17em}\mathit{k}\hspace{0.17em}\mathit{T}\mathit{/}\mathit{m,}$(9)where T is the temperature of the atmosphere, k is the Boltzmann constant, ρ is the density of the atmosphere which can be assumed to have an exponential profile ρ = ρ₀ exp(−z/h) and h = kT/mg is the scale-height (with g being the gravitational acceleration on the planet). The mean molecular mass is m = μm_H (where μ is the mean molecular weight and m_H is the hydrogen mass). The pulsar wind and atmospheric pressures are equal at a height of about: $\mathit{ẑ}\mathrm{=}\mathit{h}\hspace{0.17em}\hspace{0.17em}\ln \hspace{0.17em}\left(\frac{{\mathit{\rho }}_{\mathit{0}}\hspace{0.17em}\mathit{k}\hspace{0.17em}\mathit{T}\hspace{0.17em}\mathit{4}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{r}}^{\mathit{2}}\hspace{0.17em}\mathit{c}}{\mathit{\xi }\hspace{0.17em}{\mathit{L}}_{\mathrm{sd}}\hspace{0.17em}\mathit{m}}\right)\mathrm{·}$(10)This is the location where the relativistic particles that compose the wind will form a shock. The density after the shock can be found by equating Eqs. (3) and (9): $\mathit{\rho }\mathrm{=}\frac{\mathit{\xi }\hspace{0.17em}{\mathit{L}}_{\mathrm{sd}}\mathit{m}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\mathit{c}\hspace{0.17em}\mathit{k}\hspace{0.17em}\mathit{T}}\mathrm{·}$(11)The temperature of the shocked gas can be calculated from the perfect gas-law by knowing the pulsar-wind velocity v_pw: ${\mathit{T}}_{\mathrm{sh}}\mathrm{=}\mathit{\psi }\frac{\mathit{mp}}{\mathit{k}\hspace{0.17em}\mathit{\rho }}\mathrm{=}\mathit{\psi }\frac{\mathrm{3}}{\mathrm{16}}\frac{\mathit{m}}{\mathit{k}}{\mathit{v}}_{\mathrm{pw}}^{\mathrm{2}}\mathit{,}$(12)where the parameter ψ is an efficiency parameter which is equal to one when the whole spin down power is carried by the relativistic wind and the whole energy carried by the wind particles is transferred into the atmospheric gas. In the remainder of the paper we assume for simplicity that ψ ~ 1. Since the pulsar-wind is ultra-relativistic, we can assume that v_pw ≈ c and thus the shocked gas will have an enormous temperature T_sh ~ 10¹³ K. The post-shock gas will be optically thin, thus the emerging radiation will be characterized by a very small Compton y-parameter and the main cooling mechanism for the shocked gas will be thermal Bremsstrahlung. At these temperatures, the typical photon energy is ≈ 1 GeV (i.e., ≈ 1.6 × 10^-3 erg). Such high energy photons will penetrate deep into the atmosphere and deposit their energy in a layer whose depth depends on the photon energy (see Sect. 5.2).

Since the pulsar is injecting a fraction f of its spin down energy into the planet, then the mass loss Ṁ_pw can reach its maximum allowed value if we assume that all the kinetic energy of the relativistic wind is transferred via collisions to the atmospheric particles (and thus we neglect the photon energy deposition discussed above): $\begin{array}{ccc}\mathit{Ṁ}\mathrm{pw}& \mathrm{=}& \mathit{\psi }\frac{\mathit{f}\hspace{0.17em}{\mathit{L}}_{\mathrm{sd}}}{{\mathit{U}}_{\mathrm{b}}}\mathit{,}\\ & \mathrm{=}& \mathrm{5}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}\mathit{\psi }\left(\frac{{\mathit{L}}_{\mathrm{sd}}}{{\mathrm{10}}^{\mathrm{34}}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1}}\right){\left(\frac{{\mathit{R}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{\oplus }}}\right)}^{\mathrm{3}}{\left(\frac{\mathit{D}}{\mathrm{1}\hspace{0.17em}\mathrm{au}}\right)}^{-2}{\left(\frac{{\mathit{M}}_{\mathrm{p}}}{\mathrm{5}\hspace{0.17em}{\mathit{M}}_{\mathrm{\oplus }}}\right)}^{-1}\hspace{0.17em}\mathrm{g}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}\end{array}$(13)where U_b is the specific gravitational binding energy, $\mathit{f}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}{\left(\frac{{\mathit{R}}_{\mathrm{p}}}{\mathit{D}}\right)}^{\mathrm{2}}$ depends on the planet radius and distance (D) from the pulsar. The total binding energy per unit mass on a planet is: ${\mathit{U}}_{\mathrm{b}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{p}}}\mathrm{·}$(14)By taking as a typical value the Earth’s radius and M_p = 5 M_⊕, then the binding energy is U_b ≈ 3 × 10¹² erg g^-1. The high energy photons are much more energetic than the typical mean molecular energy U_bm ≈ 5 × 10^-10 erg and thus might potentially cause a large mass outflow from the atmosphere.

5.2. Effect of high energy radiation

Smith et al. (2004) and Cecchi-Pestellini et al. (2009) have shown that the X-ray irradiation of the planetary atmospheres can produce a substantial heating. Usually the high energy X-rays and gamma-rays can penetrate much deeper in the atmosphere than the UV and soft X-rays. However, for sufficiently thick atmospheres even the hardest X-ray/gamma photons (>10 keV) might not reach the surface of the planet. Those planets orbiting around isolated neutron stars (which are not emitting pulsar winds, Sect. 4.1) might therefore undergo a similar evolution as the planets exposed to intense X-ray radiation from young main sequence stars. On Earth, the flow of X-rays is quickly blocked by the upper atmosphere (thermosphere) which is a very optically thin layer of gas that becomes ionized when interacting with X-rays and UV radiation. Such layer has a relatively large temperature of hundreds up to thousands degrees but it is very inefficient at conducting heat since it is very rarefied. The penetration depth of high energy radiation depends on three quantities, the initial intensity I_ν,0 of the radiation field, the energy E of the high energy photon and the composition of the atmosphere. In general the attenuation of X-ray/gamma radiation is characterized by the so-called linear attenuation coefficient β and can be described as: ${\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathit{\nu }\hspace{0.17em}\mathrm{0}}{\mathit{e}}^{\mathrm{-}\mathit{\beta }\hspace{0.17em}\mathit{x}}\mathit{,}$(15)where x is the distance traveled by the X-ray photon in the atmosphere. The linear attenuation coefficient is also dependent on the density of the medium being crossed by the high energy radiation. Since the atmospheric density scales exponentially, the linear attenuation coefficient can be written: $\mathit{\beta }\mathrm{=}{\mathit{\beta }}_{\mathit{m}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}\mathit{z}\mathit{/}\mathit{h}\mathrm{)}\mathit{,}$(16)where β_m is the mass attenuation coefficient and can be found in tables². To find the height z_abs (measured from the planet surface) where the specific intensity will be reduced by a factor e, we define the total atmospheric thickness as S such that x = S−z and combine Eqs. (15) and (16): $\hspace{0.17em}\ln \hspace{0.17em}\left(\frac{{\mathit{I}}_{\mathit{\nu }}}{{\mathit{I}}_{\mathit{\nu ,}\mathit{0}}}\right)\mathrm{=}\mathrm{-}{\mathit{\beta }}_{\mathit{m}}{\mathit{\rho }}_{\mathit{0}}\mathrm{(}\mathit{S}\mathrm{-}\mathit{z}\mathrm{)}{\mathit{e}}^{\mathrm{-}\mathit{z}\mathit{/}\mathit{h}}\mathit{.}$(17)The solution for this equation can be written as: ${\mathit{z}}_{\mathrm{abs}}\mathrm{=}\mathit{S}\mathrm{-}\mathit{h}\hspace{0.17em}\mathit{W}\left(\frac{{\mathit{e}}^{\mathit{S}\mathit{/}\mathit{h}}}{{\mathit{\rho }}_{\mathrm{0}}{\mathit{\beta }}_{\mathit{m}}\hspace{0.17em}\mathit{h}}\right)\mathit{,}$(18)where the function W is the Lambert W function. A heating of the atmosphere by high energy radiation can penetrate in deeper atmospheric layers than UV and/or soft X-ray radiation and thus generate a hydrodynamic escape of the atmosphere, beside the thermal jeans escape due to direct collisions with atoms/molecules.

Fig. 4 Absorption height of high energy radiation for Earths and Super-Earths with total atmospheric thickness of S = 0.1 Rp. The solid curves show the penetration depth of radiation of different energy (hard X-rays/γ-rays βm = 0.01 → dark red; hard X-rays βm = 0.1 → orange; soft X-rays βm = 100 → yellow). The dashed gray line identifies the threshold height where the number density reaches 1015 cm-3. Below this threshold height, collisions play an important role in transferring heat to the lower atmosphere.

The typical density of the atmospheric layer where the X-ray photons are absorbed will depend on a number of factors like the physical state and chemical composition of the atmosphere. If we choose as an illustrative example an Earth-like atmosphere, then the height z_abs will occur at around 50–70 km from the surface, where we have used β_m = 0.1–0.01 which is suitable for photons with energies larger than ~ 10 keV. The number density can be calculated as: $\mathit{n}\mathrm{=}\frac{\mathit{\rho }}{\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}}\mathit{,}$(19)and gives a value of the order of 10¹⁷ cm^-3. At these densities the metals in the atmosphere will undergo frequent collisions and will be de-excited while being in turbulent or thermally induced motion and transfer their energy to other atoms/molecules. An illustrative example obtained for Super-Earths of different mass is shown in Fig. 4. The density of the atmosphere at the surface of the planet is calculated by imposing hydrostatic equilibrium. This implies that the atmospheres of larger planets have a strong density gradient, whereas smaller planets have more homogeneous atmospheres. Since Earth-like planets can retain only a tiny fraction of atmospheric mass we have used a smoothing function to connect the properties of Earth-like planets and Super-Earths. The smoothing function is defined as: $\mathit{f}\mathrm{\left(}{\mathit{M}}_{\mathrm{p}}\mathrm{\right)}\mathrm{=}\hspace{0.17em}\tan \mathrm{h}\mathrm{(}\hspace{0.17em}{\mathit{M}}_{\mathrm{p}}\mathrm{-}{\mathit{M}}_{\mathrm{o}}\mathrm{)}\mathrm{+}\mathit{ϵ}$(20)where M₀ = 1M_⊕  and ϵ = 10^-6 is the fraction of planetary mass retained in the atmosphere of a 1 M_⊕  planet. Since there might be also a strong compositional difference of the atmospheric elements, we used a similar smoothing function for the molecular weight (that we have chosen to go from μ = 29 for a 1 M_⊕  planet down to 3 for a Super-Earth). Finally, the radius of the planet is also increased gradually from 1 R_⊕ (for an Earth like planet) to 2 R_⊕ for Super-Earths.

6. Application to pulsar planetary systems

Fig. 5 Neutron star planets habitable zone. The diagram shows the luminosity produced by the neutron star (as pulsar wind or X-ray radiation), which is relevant to heat the planet atmosphere as a function of the planetary distance. The habitable zone is the orange area between the oblique dashed gray lines. Open symbols are upper limits.

6.1. Atmospheric heating from high energy radiation

There are several definitions in the literature about what constitutes the boundaries of an habitable zone around main sequence stars. In this work we follow the definition given by Kasting et al. (1993): the habitable zone is the region around a star where a Earth-like planet (with a CO₂/ H₂O/N₂ atmosphere) can have large amounts of liquid water on its surface.

According to Selsis et al. (2007) a necessary (but not sufficient) condition for habitability of a planet is that its equilibrium temperature T_{e q} is below 270 K. Such equilibrium temperature can be defined as: ${\mathit{T}}_{\hspace{0.17em}\mathrm{eq}}\mathrm{=}{\left[\frac{\mathrm{\Sigma }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{a}\mathrm{)}}{\mathit{\phi \sigma }}\right]}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}$(21)were Σ is the incident stellar energy flux, a is the Bond albedo of the planet, σ is the Stefan-Boltzmann constant and φ is a geometric factor of order unity that accounts for the redistribution of the heat on the surface. This equilibrium temperature does not relate necessarily to the surface temperature of the planet (see Selsis et al. 2007, for a discussion). In this paper we follow Kaltenegger & Sasselov (2011) and set the inner and outer planet habitable zone boundaries by requiring that the equilibrium temperature lies in the range 175 K < T_eq < 270 K.

We can thus rewrite the above equilibrium temperature expression in terms of the fraction η of X-rays and gamma-rays which are able to penetrate the atmosphere in a layer which is dense enough to allow the thermalization of the released energy (rather than the thermal escape of the gas). If we rewrite Eq. (21) as a function of the input neutron star luminosity (L_NS) and distance of the planet from the pulsar then we obtain: ${\mathit{T}}_{\hspace{0.17em}\mathrm{eq}}\mathrm{=}\mathit{\eta }{\left[\frac{{\mathit{L}}_{\mathrm{NS}}}{\mathrm{16}\mathit{\pi \sigma }\hspace{0.17em}{\mathit{D}}^{\mathrm{2}}}\right]}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{·}$(22)The biggest uncertainty in this equation lies in the parameter η which depends on the energy of the photon and the physical parameters of the atmosphere.

Fig. 6 Evaporation timescale of planets irradiated by a source of X-rays with LX = 1029 erg s-1 (top row), pulsar wind Lsd = 4 × 1032 erg s (middle row) and Lsd = 1034 erg s-1 (bottom row). All figures on the left show the case for a very thick atmosphere holding 30% of the planet total mass. The right panels show the case for a thinner atmosphere of 1% of the total planet mass (which is a reasonable assumption for the atmospheric mass of Super-Earths; see Lopez & Fortney 2013; Elkins-Tanton & Seager 2008). The gray circles represent the mass and distance of four of the five pulsar planets with the exclusion of PSR B1257+12  a, which is too small to hold any atmosphere. The dots represent (from left to right): PSR B1257+12 b, PSR B1257+12 c, PSR J1719–1438 b, PSR B1620–26 b. The dotted black lines represent isochrones of 0.1, 1 and 10 Gyr (from bottom to top, respectively). The radius of all planets has been fixed to Rp = 2 R⊕ and XM = 10-3. The evaporation timescales have been calculated with Eq. (24) (top figures) and 13 (middle and bottom figures).

To have an equilibrium temperature between 175 and 270 K, a planet needs to lie between the lower (D_l) and upper (D_u) habitable zone boundaries. We start as an initial example with the planets in the PSR B1257+12  system where the total input neutron star luminosity needs to be between ηL_X ≈ η 5.5 × 10²⁹ erg s^-1 (see Pavlov et al. 2007; Yan et al. 2013 and Sect. 3 of this work) and L_sd ≈ 1.5 × 10³⁴ erg s^-1. If the main source of power comes from the X-ray luminosity alone, then even if we take η = 1 then D_l = 0.02 au and D_u = 0.06 au and all three pulsar planets in PSR B1257+12 would be too cold. However, if sufficient gamma radiation is produced in the upper atmosphere due to the presence of the pulsar wind shock, then the boundaries of the habitable zone for planets b and c in PSR B1257+12  shift outwards between 2 and 5 au (see Table 3). In Fig. 5 we show the neutron star habitable zone (calculated with Eq. (22)) defined as the location where a pulsar planet surface temperature is between 175 and 270 K. The blue circles refer to the measured X-ray luminosity and distance of the pulsar planets in PSR B1257+12 and PSR B1620–26. The red squares (the open square is an upper limit) refer to the total spin-down luminosity of the same pulsar planets with the addition of the “diamond planet” pulsar PSR J1719–1438. A plausible irradiation luminosity for the planet lies between the red squares and the blue circle if the X-ray luminosity is isotropic and the pulsar wind (partially) hits the planet and the gamma-ray luminosity is absorbed by the planet atmosphere. The planet around PSR B1620–26 is too cold even in the most optimistic case, whereas the 175–270 K band lies in between the values for the pulsar planets in PSR B1257+12. No X-ray information is currently available for PSR J1719–1438 although the X-ray luminosity required to fall in the temperate 175–270 K band needs to be smaller than any other known X-ray pulsar.

Furthermore, most isolated neutron stars with Bondi-Hoyle accretion have an X-ray luminosity larger than PSR B1257+12  (Pfahl & Rappaport 2001) which is on the lower end of the X-ray luminosity spectrum rather than a typical system. For rocky planets similar to the Earth, such habitable zone would exist for a very brief amount of time due to the atmospheric evaporation whereas for Super-Earths with dense atmospheres this phase could be potentially very long. To evaluate the amount of matter lost by the planet’s atmosphere as a consequence of high energy irradiation we adapt the calculations of London et al. (1981) and Ruderman et al. (1989), who considered the evaporation rate (Ṁ_w) of a light stellar companion due to X-ray heating from a pulsar. Differently from Ruderman et al. (1989), we do not use the mass-radius relation for low mass stars, but we leave the radius and mass dependencies of the planet explicitly in the equations. We assume that the velocity of the emerging particle wind due to X-ray evaporation is equal to the escape velocity from the planet: ${\mathit{v}}_{\mathrm{e}}\mathrm{=}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{p}}}}\mathrm{=}\mathrm{2.6}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-4}{\left(\frac{\mathit{M}}{\mathit{R}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathrm{cm}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{.}$(23)For an Earth-like planet as the one under consideration here, the total mass loss rate is therefore: $\mathit{Ṁ}\mathrm{w}\mathrm{=}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{-14}\mathit{\chi }\left(\frac{{\mathit{X}}_{\mathrm{M}}}{{\mathrm{10}}^{-3}}\right){\left(\frac{{\mathit{M}}_{\mathrm{p}}}{{\mathit{M}}_{\mathrm{\oplus }}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{{\mathit{R}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{\oplus }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{{\mathit{L}}_{\mathit{x}}}{\mathrm{(}\mathrm{erg}\mathit{/}\mathrm{s}\mathrm{)}}\frac{{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{4}{\mathit{D}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{g}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}$(24)where X_M is the fractional metal abundance in the atmosphere, χ is the ratio between the soft (0.2–1 keV) X-ray intensity to the total incident X-ray intensity, relative to the same ratio for the X-ray binary pulsar Her X-1 (which was used by Ruderman et al. 1989 as a reference source for the X-ray emission). If we assume an atmosphere completely dominated by heavy elements (as is the case for the Earth) and we use χ = 1 and X_M = 1 then the expression above reduces to: $\mathit{Ṁ}\mathrm{w}\mathrm{\simeq }\mathrm{3}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}\left(\frac{{\mathit{L}}_{\mathrm{x}}}{{\mathrm{10}}^{\mathrm{30}}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1}}\right){\left(\frac{\mathit{D}}{\mathrm{1}\hspace{0.17em}\mathrm{au}}\right)}^{-2}\left(\frac{{\mathit{X}}_{\mathrm{M}}}{{\mathrm{10}}^{-3}}\right)\hspace{0.17em}\mathrm{g}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{.}$(25)We note that these estimates are compatible with mass loss rates observed in exoplanets exposed to intense levels of radiation (e.g., Ehrenreich & Désert 2011). The evaporation timescale to lose an Earth-like atmosphere mass (M_atm) would therefore be τ_evp = M_atm/Ṁ_w ~ 5 × 10²¹ g /Ṁ_w. Typical values for τ_evp lie in the range ~ 10⁵–10⁷ yr for Earth-like planets at distances of 1–10 au and up to trillion years for the most extreme Super-Earths with thick atmospheres (see Fig. 6).

6.2. Habitable zone parameter range

We can now ask the question about the physical parameters that a planet needs to have in order to fall within the pulsar habitable zone (as discussed in Sect. 6.1). We begin with a pulsar that has a spin-down luminosity among the smallest possible values (see Fig. 7 for the distribution of spin-down luminosity for all measured PSRs and MSPs), i.e. L_sd ≈ 10²⁹ erg s^-1 and we consider the effect of the wind on its atmosphere.

Fig. 7 Histogram of the spin-down luminosity of radio pulsars. The median value is ≈ 4 × 1032 erg s-1.

If we assume that the planet is similar to the Earth (in mass, atmospheric composition, density and distance) then the shock will form at an height ẑ of approximately 200 km (see Eq. (10)) which would correspond to the thermosphere of the planet. We now require that the planet equilibrium temperature lies within the range 175–275 K and we use Eq. (21) to find a range of allowed distances: 0.02–0.06 au. Using Eq. (13) for a typical Earth mass/radius planet, gives a mass loss rate of the order of 10¹² g s^-1 which means that an Earth-like atmosphere would be potentially consumed in a time as short as few hundred years (for spin down luminosities of 10³⁴ erg s^-1). The most optimistic case that will give the longest possible survival time of the atmosphere is the case of a Super-Earth which, in the most extreme cases, can have an atmosphere of about 30% the planet’s mass (Elkins-Tanton & Seager 2008). The binding energy U_b for such planet would be similar to that of an Earth-like planet and the corresponding mass loss would be Ṁ_pw ≈ 10¹² g s^-1. The atmospheric mass of such planets can reach up to ~ 10²⁸ g, which would give an evaporation timescale of the order of 0.1–1 Gyr (assuming that all pulsar wind power is used to eject the atmospheric gas).

It is important also to stress that the pulsar wind is not a process that will continue indefinitely. Indeed once the pulsar reaches a sufficiently slow spin period it will cross the so-called “death-line” in the P−Ṗ diagram, meaning that the pulsar wind will turn off. For young pulsars, this occurs on a timescale of the order of million years, whereas for MSPs it is of the order of billion years. However, this would also turn off the energy source of the planet and therefore the temperature will drop dramatically thus removing any possibility to define an habitable zone, unless a Bondi-Hoyle accretion process generates a sufficiently large amount of X-ray radiation as discussed in Sect. 4.3 or other effects like radiogenic heating (Sect. 7) or tidal heating play a dominant role.

6.3. Presence of a magnetosphere

If a magnetosphere is present around pulsar planets, then its shielding effect can be fundamental to deflect and/or confine the incoming relativistic pulsar wind particles. If we assume that the planetary magnetic field B_p(r) is dipolar, then when the magnetic pressure balances the pulsar wind ram pressure, a shock will form at a distance r_s from the planet: $\frac{{\mathit{B}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{\left(}{\mathit{r}}_{\mathrm{s}}\mathrm{\right)}}{\mathrm{8}\mathit{\pi }}\mathrm{=}\mathit{\xi }\hspace{0.17em}{\mathit{L}}_{\mathrm{sd}}\mathit{/}\mathrm{(}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}\mathit{c}\mathrm{)}\mathit{.}$(26)If B_p is the value of the magnetic field on the planet’s surface, then r_s can be found imposing the condition that it is larger than the radius of the planet R_p. In Table 4 we calculate the minimum surface magnetic field B_p,min required for each of the known pulsar planets to withstand the ram pressure of the pulsar wind.

6.4. Are pulsar winds really hitting the planets?

In the preceding discussion we have assumed that the pulsar wind is emitted isotropically. The mechanism to produce a radio pulsar beam and particle wind is still not completely understood, but the current physical picture requires the production of electron-positron pairs. When pair creation happens, the negative current flows along the poles whereas the positive charge/current flows in a sheet on the pulsar equator (Chen & Beloborodov 2014). In the open filed zone there is a jump in the toroidal field component B_φ that requires charge density to be negative (Lyubarskii 1990). There is a matter flow of ions in the equatorial plane with energy density which is about twice the magnetic field energy density. In the aligned rotator geometry (i.e., when spin and magnetic axes are nearly aligned), the current sheet of ions in the equator extends indefinitely if nothing stops this particle flux. If the misalignment is instead substantial, the flow of particles oscillates (Philippov et al. 2015; Tchekhovskoy et al. 2016). It is thus possible to speculate that the wind could completely miss the planets orbiting around the pulsar because of a geometric misalignment of the orbital plane with respect to the ion current sheet or negative charge flow. In this case the planet will be only irradiated by the X-ray emission of the neutron star.

For PSR B1257+12  we can calculate with Eq. (24) what would be the total mass lost by the planets in this case. By using the current characteristic age of the pulsar (≈ 850 Myr) we expect that the two outermost planets should have lost a total of ≈ 5 × 10^-4M_⊕. We caution that given the many uncertainties on the quantities involved, it is not possible to conclude that this is the correct scenario for PSR B1257+12 . However, we note that the value calculated above is close to the estimate made in Eq. (1) based on X-ray observations.

7. Formation and composition of neutron star planets

The low-mass planetary system around PSR B1257+12 is unique in its architecture among the ~ 400 millisecond pulsars³. This has been interpreted as evidence for a very low probability formation channel. However, low-mass planetary systems at larger orbital separations are still unconstrained around MSPs and other neutron stars. Therefore, the PSR B1257+12 planets may represent a rare case of either forming at, or migrating to, sub-au distances, while planet formation around neutron stars could still be a general phenomenon.

Neutron star planets may originate from one of three epochs: a protoplanetary disk during the star formation process (1st generation); a post-supernova fallback disk (2nd generation); or an accretion disk formed from material stripped by the neutron star from a binary companion (3rd generation). Any 1st generation planets, if they exist, are likely destroyed or have their orbits disrupted during the supernova. Millisecond pulsar planets cannot easily be explained by 2nd generation formation, as spin-up of the pulsar by a companion star which migrated to within 1 au would disrupt their orbits. This is justified by the fact that, since PSR B1257+12  contains a fully recycled pulsar, its companion star must have been a low-mass star which has undergone Roche lobe overflow during the main sequence phase. This requires that the orbital separation had been less than about ~ 1 au. Thus, 3rd generation formation, from the remnants of the companion star, is preferred (Podsiadlowski 1993).

The gas-to-dust ratio, Δ_g/d, i.e. the mass ratio of gas-phase volatiles (H, He, N, etc.) to solids (Fe, Mg, Si, Al etc.) is a major controlling parameter of planet formation. Models of the formation of the PSR B1257+12 planets have found that protosolar nebula like conditions with Δ_g/d = 100 are preferred over a supernova fallback disk with Δ_g/d ≈ 4 (e.g. Currie & Hansen 2007). Indeed, the architecture of the PSR B1257+12 system resembles the inner solar system, and it may have formed from a radially confined planetesimal belt similar to that proposed for the solar system (Hansen et al. 2009).

Based on main sequence and supernova nucleosynthesis models for progenitor masses ≤ 25 M_⊙, the gas-to-dust ratio of post-supernova material is ≳ 20 (Nomoto et al. 2006). Assuming a sufficiently massive (M_dust ≳ 10^-4M_⊙) fallback disk remains around the newly-formed neutron star, the high gas-to-dust ratio would speed up dust growth and favour the occurrence of the streaming instability (Δ_g/d ≳ 1), potentially leading to a rapid build-up of planetesimals (Johansen et al. 2007). Such post-supernova material would also have a peculiar composition. The mass ratio of O to Si, Fe and Mg in the supernova ejecta increases from roughly solar (≈unity) for a 13 M_⊙ progenitor to ~ 3 for 25 M_⊙ (Nomoto et al. 2006). The fraction of O left over from rock formation and available for H₂O formation is different in these cases, so 2nd generation neutron star planets could potentially be very water-rich.

A major difference between neutron stars and main sequence ones is the shape of their radiation field. Radio waves would not participate significantly in disk heating or ionization, and the ultraviolet and visible luminosity is negligible for evolved systems although it may play a role when the disk is accreting onto the neutron star (e.g. Martin et al. 2016). The high flux of high-energy photons and relativistic wind particles may lead to a high ionization fraction which could mean there is no dead zone in the disk. Martin et al. (2016) suggest that this would mean a low formation efficiency of rocky planets, as dead zones are thought to aid growth (Johansen et al. 2014). As discussed earlier, the shape of the relativistic particle wind is currently not known, and its role in the disk ionization is thus not well constrained, so we leave these issues open.

7.1. Evolution of internal temperature: radiogenic heating

The internal energy flux of Earth is (47 ± 2) × 10¹⁹ erg s^-1 or 7.87 × 10^-8 erg s^-1 g^-1 (Davies & Davies 2010). This energy comes from radioactive isotopes, currently totalling about 5 × 10^-8 erg s^-1 g^-1, of which 2.9 × 10^-8 erg s^-1 g^-1 is from ⁴⁰K and the rest largely ²³⁸U and ²³²Th, and from the heat of formation. For Earth, the total internal energy production is only 0.027% of solar insolation and does not directly control the surface temperature. It is, however, important for habitability because it drives plate tectonics and the geodynamo which gives rise to the planetary magnetic field. This may offer protection against atmospheric ablation by the pulsar wind (Sect. 6.3).

Second-generation planets around a neutron star would be composed of the ejecta of the progenitor supernova, which can produce a wide range of radioisotope abundances. Notably, ²³⁸U and ²³²Th may be less relevant in SNII ejecta than in average galactic material, because they are thought to be produced in neutron star mergers rather than Type II supernovae, although this is still being debated (Eichler et al. 1989; Freiburghaus et al. 1999; Cowan et al. 2005; Tanvir et al. 2013). The production of the long-lived planetary energy source ⁴⁰K in SNII is well established (e.g. Rauscher et al. 2002; Kobayashi et al. 2006). Table 5 shows the energy production per gram of planetary core (refractory) mass for Earth and for planets composed of ejecta from SNII with different progenitor masses.

8. Conclusions

We have shown that the harsh environment around neutron stars can still accommodate planets with warm atmospheres, provided they are Super-Earths. We addressed the effect of radiation produced by neutron stars on planetary atmospheres and defined a habitable zone for the particular case of low UV/optical radiation, moderate X-ray, and high gamma-ray and relativistic particle flux. We find that, if part of the pulsar power is injected in the atmosphere, all of the three PSR B1257+12 planets may lie in the habitable zone. In particular, the two Super-Earths may have retained their atmosphere for at least a hundred million years provided they contain a large atmospheric fraction of the total planet mass, with the atmosphere possibly still being present to these days. We also find that if a moderately strong planetary magnetosphere is present, the atmospheres can survive the strong pulsar winds and reach survival timescales of several billion years. The same argument applies to possible pulsar planets around more powerful objects than PSR B1257+12. Alternatively, a similar result can be achieved if a non-isotropic pulsar wind is present in the system. Furthermore, planets which lie within a band of ≈ 0.01–1 au can be in their habitable zone provided that the neutron stars transfer their energy to the planetary atmosphere via X-rays alone (e.g., through Bondi-Hoyle accretion). We have also briefly discussed the formation of neutron star planets, their potential high water content, and we have highlighted the fact that radiogenic heating could play an important role in these systems.

Acknowledgments

We would like to thanks O. Shorttle, S. Popov, Y. Cavecchi, M. Kenworthy, J. Hessels, R. Wijnands, S. Portegies Zwart, B. Czerny and T. Stolker for useful discussions. A.P. acknowledges support from a NWO Vidi Fellowship.

References

All Tables

All Figures

	Fig. 1 PSR B1257+12  as seen by Chandra on 2007, May 3. The green circle is the 1.5′′ radius extraction region used to find the source counts. The region is centered around the white pixel, which is the brightest one in the image. The color scale refers to number of counts.
In the text

	Fig. 2 Contour plot for temperature (units of keV) vs. normalization K. The contours refer to the 68% (red line), 90% (green line) and 99% (blue line) confidence levels. Degeneracy between temperature and normalization is observed, but it is still not sufficient to reconcile the small emitting area with theoretical models of millisecond pulsars.
In the text

	Fig. 3 Contour plot for the joint 2005–2007 fit of the column density NH vs. blackbody normalization. The contours and symbols are the same as in Fig. 2. The degeneracy between column density and blackbody normalization is evident.
In the text

Fig. 4 Absorption height of high energy radiation for Earths and Super-Earths with total atmospheric thickness of S = 0.1 Rp. The solid curves show the penetration depth of radiation of different energy (hard X-rays/γ-rays βm = 0.01 → dark red; hard X-rays βm = 0.1 → orange; soft X-rays βm = 100 → yellow). The dashed gray line identifies the threshold height where the number density reaches 1015 cm-3. Below this threshold height, collisions play an important role in transferring heat to the lower atmosphere.

In the text

	Fig. 5 Neutron star planets habitable zone. The diagram shows the luminosity produced by the neutron star (as pulsar wind or X-ray radiation), which is relevant to heat the planet atmosphere as a function of the planetary distance. The habitable zone is the orange area between the oblique dashed gray lines. Open symbols are upper limits.
In the text

Fig. 6 Evaporation timescale of planets irradiated by a source of X-rays with LX = 1029 erg s-1 (top row), pulsar wind Lsd = 4 × 1032 erg s (middle row) and Lsd = 1034 erg s-1 (bottom row). All figures on the left show the case for a very thick atmosphere holding 30% of the planet total mass. The right panels show the case for a thinner atmosphere of 1% of the total planet mass (which is a reasonable assumption for the atmospheric mass of Super-Earths; see Lopez & Fortney 2013; Elkins-Tanton & Seager 2008). The gray circles represent the mass and distance of four of the five pulsar planets with the exclusion of PSR B1257+12  a, which is too small to hold any atmosphere. The dots represent (from left to right): PSR B1257+12 b, PSR B1257+12 c, PSR J1719–1438 b, PSR B1620–26 b. The dotted black lines represent isochrones of 0.1, 1 and 10 Gyr (from bottom to top, respectively). The radius of all planets has been fixed to Rp = 2 R⊕ and XM = 10-3. The evaporation timescales have been calculated with Eq. (24) (top figures) and 13 (middle and bottom figures).

In the text

	Fig. 7 Histogram of the spin-down luminosity of radio pulsars. The median value is ≈ 4 × 1032 erg s-1.
In the text