Evaporation of planetary atmospheres due to XUV illumination by quasars

TL;DR

This paper investigates how XUV radiation from growing supermassive black holes in galactic centers can cause significant atmospheric mass loss on planets, with potential impacts varying by galaxy type and mass.

Contribution

It introduces a model estimating planetary atmosphere loss due to black hole XUV irradiation across different galaxy types and cosmic epochs.

Findings

Approximately 50% of planets may lose a Martian atmosphere.

About 10% could lose an Earth's atmosphere.

0.2% might lose Earth's ocean mass.

Abstract

Planetary atmospheres are subject to mass loss through a variety of mechanisms including irradiation by XUV photons from their host star. Here we explore the consequences of XUV irradiation by supermassive black holes as they grow by the accretion of gas in galactic nuclei. Based on the mass distribution of stars in galactic bulges and disks and the luminosity history of individual black holes, we estimate the probability distribution function of XUV fluences as a function of galaxy halo mass, redshift, and stellar component. We find that about 50% of all planets in the universe may lose the equivalent of a Martian atmosphere, 10% may lose an Earth's atmosphere, and 0.2% may lose the mass of Earth's oceans. The fractions are appreciably higher in the spheroidal components of galaxies, and depend strongly on galaxy mass, but only weakly on redshift.

Tables

Table 1. Table 1: Examples of fluxes, fluences, mass loss rates, and accumulated mass lost. In the top half of the Table, we use a timescale appropriate for the growth of supermassive black holes (40 Myr), while for the bottom half we use a timescale appropriate for planets around older stars (5 Gyr). Bold entries indicate the starting point from which the other three entries are derived for a given row.

	${\log }_{10}{F}_{\mathrm{XUV}}\left(\frac{\mathrm{erg}}{s{\mathrm{cm}}^2}\right)$	${\log }_{10}{\Psi }_X\left(\frac{\mathrm{erg}}{{\mathrm{cm}}^2}\right)$	${\log }_{10}{\dot{M}}_{\mathrm{loss}}$ $\left(\frac{g}{s}\right)$	${\log }_{10}{M}_{\mathrm{lost}}$ (g)
Earth’s Oceans	3.7	18.8	9.0	24.15
Earth’s Atmosphere	1.3	16.4	6.6	21.7
Mars’s Atmosphere	-1.0	14.1	4.3	19.4
Flux at Earth from Sgr A* if ${\eta }_X=0.1$	-2.18	12.9	-	-
Current XUV Flux at Proxima b	2.47	19.7	7.8	25.0
Current XUV Flux at Earth	0.67	17.9	6.0	23.2
Current Earth Mass Loss Rate	-1.8	15.4	3.5	20.7

Equations

$$\dot{M}=ϵ\frac{F_{\mathrm{XUV}}\pi R_p^2}{G{M}_p/R_p}$$

$$M_{\mathrm{lost}}\approx 3ϵ(4\pi G{\rho }_P{)}^{-1}\int F_{X}dt$$

$$M_{BH}=\int \frac{(1-{ϵ}_{\mathrm{BH}})L_{BH}}{{ϵ}_{\mathrm{BH}}{c}^2}dt=\frac{4\pi r^2}{c^2}\int \frac{(1-{ϵ}_{\mathrm{BH}})F_X}{{ϵ}_{\mathrm{BH}}{\eta }_X}dt$$

$${\Psi }_X\approx \frac{{ϵ}_{\mathrm{BH}}}{1-{ϵ}_{\mathrm{BH}}}\frac{{\eta }_X{M}_{\mathrm{BH}}{c}^2}{4\pi r^2}$$

$$r=1.8\mathrm{kpc}{\left(\frac{M_{\mathrm{BH}}}{10^9{M}_{\odot }}\right)}^{1/2}{\left(\frac{{\Psi }_X}{10^{17}\mathrm{erg}{\mathrm{cm}}^{-2}}\right)}^{-1/2}$$

$$\frac{d{M}_{\ast }}{d{\Psi }_X}=\frac{d{M}_{\ast }}{dr}\left|\frac{dr}{d{\Psi }_X}\right|$$

$$\Sigma ={\Sigma }_0\exp \left(-b(r/r_E{)}^{1/n}\right)$$

$$\frac{dP}{d{\log }_{10}{\Psi }_X}=\frac{\ln (10)}{\Gamma (2n)}\frac{E_X{b}^{2n}}{8\pi r_E^{2}n{\Psi }_X}\exp \left(-b{\left(\frac{E_X}{4\pi r_E^2{\Psi }_X}\right)}^{1/2n}\right)$$

$$\frac{dP}{d{\log }_{10}{\Psi }_X}=\frac{\ln (10)E_X/(4\pi {\Psi }_X{a}^2)}{(\sqrt{E_X/(4\pi a^2{\Psi }_X)}+1{)}^3}$$

$$\frac{d{\rho }_{\ast }}{d{\log }_{10}\Psi }$$

$$BT\sim \{\begin{array}{l}{\mathcal{N}}_{0.01-0.99}\left(0.2+0.2{\log }_{10}\left(\frac{M_{\ast }}{10^{10}{M}_{\odot }}\right),0.33\right)\mathrm{if}\mathrm{SF}\\ {\mathcal{N}}_{0.01-0.99}\left(0.55+0.1{\log }_{10}\left(\frac{M_{\ast }}{10^{10}{M}_{\odot }}\right),0.33\right)\mathrm{if}\mathrm{not}\mathrm{SF},\end{array}$$

$${\log }_{10}{M}_{\mathrm{BH}}\sim \mathcal{N}(8.46+1.05{\log }_{10}\left(\frac{M_b}{10^{11}{M}_{\odot }}\right),0.34)$$

$${\log }_{10}{r}_{E,b}\sim \{\begin{array}{l}\mathcal{N}\left(-0.097+0.097{\log }_{10}\left(\frac{M_{\ast }}{10^{10}{M}_{\odot }}\right),0.38\right)\mathrm{if}\mathrm{SF}\\ \mathcal{N}\left(-0.35+0.25{\log }_{10}\left(\frac{M_{\ast }}{10^{10}{M}_{\odot }}\right),0.38\right)\mathrm{not}\mathrm{if}\mathrm{SF}\end{array}$$

$${\log }_{10}{r}_{E,d}\sim \{\begin{array}{l}\mathcal{N}\left(0.48+0.125{\log }_{10}\left(\frac{M_{\ast }}{10^{10}{M}_{\odot }}\right),0.378\right)\mathrm{if}\mathrm{SF}\\ \mathcal{N}\left(0.28+0.42{\log }_{10}\left(\frac{M_{\ast }}{10^{10}{M}_{\odot }}\right),0.4\right)\mathrm{if}\mathrm{not}\mathrm{SF}\end{array}$$

$$p(v|r)\equiv \frac{p(v,r)}{p(r)}=\frac{4\pi r^{2}4\pi v^{2}f}{4\pi r^2\rho }$$

$$\frac{dr}{dt}=v_r=\pm \sqrt{2(E-\varphi (r))-\frac{L^2}{r^2}}.$$

$$M_{\mathrm{lost}}=7.7\times 10^{25}\mathrm{g}{\left(\frac{{\eta }_X}{0.1}\frac{F_{\mathrm{bol}}}{F_E}\right)}^{1/2}\frac{t_0}{1\mathrm{Gyr}}$$

Full text

1 Introduction

The recent discoveries of the TRAPPIST1 planetary system (Gillon et al., 2017), LHS 1140b (Dittmann et al., 2017), and Proxima Centauri b (Anglada-Escudé et al., 2016), each with planets in the putative habitable zone of their low-mass host stars, has generated intense interest in the susceptibility of their atmospheres to erosion and more generally the viability of life on their surfaces (Ribas et al., 2016; Garraffo et al., 2016; Airapetian et al., 2017; Dong et al., 2017; Ranjan et al., 2017). Despite their low masses, M-dwarfs are capable of generating substantial high-energy radiation through frequent flares (e.g. Barnes et al., 2016; Meadows et al., 2016; Davenport et al., 2016; Vida et al., 2017; Ribas et al., 2017; Miles & Shkolnik, 2017). The large flux of high-energy radiation to which the close-in planets are subjected is capable of ionizing their upper atmospheres and driving hydrodynamic winds, the same mechanism responsible for occasionally dramatic mass loss among hot Jupiters (Murray-Clay et al., 2009).

A given planet may be subject to sources of high energy radiation besides its host star. There are many such sources in the universe, but we will argue that the most important, aside from the host star of the planet, is radiation from the growth of supermassive black holes in galactic nuclei. In their most extreme phase of growth, when the black hole is potentially visible as a quasar, we will argue that the radiation from the black hole is sufficient to substantially affect planets on scales comparable to the entire galaxy.

In this paper we will apply results from the study of stellar XUV radiation on planetary atmospheres to the problem of quasars illuminating the same atmospheres. This will result in a straightforward relationship between the planet’s distance from the center of the galaxy, the mass of the supermassive black hole at the galaxy’s center, and the mass lost by the planet’s atmosphere. Based on this relationship, we can estimate the probability distribution function of mass loss, i.e. the probability a random planet in the universe will be subject to a given amount of mass loss, as a function of redshift, mass, and the stellar component in which the planet’s host star resides. We will estimate the effect introduced by the fact that many stars in the universe have non-circular orbits, and hence will occasionally pass much closer to the central SMBH. Finally, we discuss the implications for habitability - on Earth, and in galaxies in general.

2 The relationship between XUV flux from black holes and mass loss from planetary atmospheres

Planetary atmospheres subjected to the high-energy radiation of Active Galactic Nuclei (AGN) may experience a number of effects harmful to life on the planet. The primary effect we consider here is the loss of substantial mass from the atmosphere due to heating by photons near or exceeding the ionization potential of hydrogen.

In principle the mass loss from planetary atmospheres can be quite complicated and difficult to predict. In the solar system, the Earth loses mass via Jeans mass loss, namely the loss of individual particles on the high-velocity tail of the Maxwell-Boltzmann distribution, while Mars has lost much of its atmosphere to the solar wind via sputtering as a result of a minimal protective magnetic field and low surface gravity (Jakosky et al., 2017). In addition to mass lost one particle at a time by these methods, planets subject to heavy irradiation from high-energy photons can experience hydrodynamic winds, wherein the upper atmosphere is heated and flows outward in the fluid regime (e.g. Lammer et al., 2003; Erkaev et al., 2007; Murray-Clay et al., 2009; Lopez et al., 2012; Owen & Jackson, 2012; Owen & Alvarez, 2016; Chen & Rogers, 2016; Zahnle & Catling, 2017).

As is standard in the literature, we parameterize the mass loss by the incident flux of XUV radiation (meaning radiation with energies encompassing extreme ultraviolet radiation starting at 10 eV up to and including soft X-rays),

[TABLE]

Here $F_{\mathrm{XUV}}$ is the XUV flux, and $R_{p}$ and $M_{p}$ are the plaent’s radius and mass respectively. Physically, this equation assumes that some fraction $\epsilon$ of the incident power available to heat the atmosphere, $\pi R_{p}^{2}F_{\mathrm{XUV}}$, is used in driving the atmosphere out of the planet’s gravitational potential. Crucially $\dot{M}$ depends on the planet’s properties (excluding composition and other higher-order complications) only via the bulk density $\rho_{p}\propto M_{p}/R_{p}^{3}$, which does not vary immensely, particularly for rocky planets of the most interest for habitability.

The value of $\epsilon$ can be estimated theoretically (Murray-Clay et al., 2009; Owen & Alvarez, 2016) to be of order $10\%$ over a wide range of fluxes (Bolmont et al., 2017). Below a critical value of $\sim 0.1\ \mathrm{erg}\ \mathrm{s}^{-1}\mathrm{cm}^{-2}$, the flux is insufficient to drive a hydrodynamic wind and $\epsilon$ falls precipitously. Above about $1000\ \mathrm{erg}\ \mathrm{s}^{-1}\mathrm{cm}^{-2}$, $\epsilon$ begins to fall slowly as the outflow solution begins to lose substantial energy to radiative cooling. In this regime $\dot{M}\propto F_{\mathrm{XUV}}^{1/2}$ (Murray-Clay et al., 2009). The Appendix discusses an example where the outflow is subject to high radiative losses, namely planets in the habitable zone around host objects with high XUV luminosities. In the regime where $\epsilon\approx 10\%\approx\mathrm{const.}$, one only needs a timescale over which a given XUV flux acts to estimate the total mass lost by hydrodynamic winds.

The total mass lost due to hydrodynamical winds in the regime where $\epsilon\approx\mathrm{const}$ is simply the time integral of equation (1),

[TABLE]

where the integral is evaluated from the time a planet develops a steady atmosphere to the epoch under consideration. We have assumed that the planet’s properties are roughly constant in time. This may not be true in the early stages of planetary evolution, as even terrestrial planets may retain substantial gas accreted directly from the protoplanetary disk. However, we expect that for planets that may host life, the properties of the planet itself should change relatively little over long time periods.

Equation (2) implies that the critical remaining unknown quantity is the fluence, $\Psi_{X}\equiv\int F_{X}dt$, integrated over the relevant timescale. In order to proceed, one path is to specify a source luminosity $L_{X}$ and lifetime $t_{0}$, which translates into a fluence given a distance $r$, i.e. $\Psi_{X}=t_{0}L_{X}/(4\pi r^{2})$. In the case of a supermassive black hole, we take $L_{X}=\eta_{X}\eta_{\mathrm{Edd}}L_{\mathrm{Edd}}$, where $\eta_{X}$ is the fraction of the spectrum of sufficiently high energy to aid in driving hydrodynamic winds, $\eta_{\mathrm{Edd}}$ is the Eddington ratio, namely the ratio of the bolometric luminosity to the Eddington luminosity, and $L_{\mathrm{Edd}}\approx 1.3\times 10^{38}(M_{\mathrm{BH}}/M_{\odot})\ \mathrm{erg}\ \mathrm{s}^{-1}$ is the Eddington luminosity itself for a black hole of mass $M_{\mathrm{BH}}$.

The value of $t_{0}$ has been extensively studied over many decades as the quasar duty cycle or lifetime. It can be estimated by matching the present-day mass function of supermassive black holes to the luminosity function of quasars at higher redshift (Richstone et al., 1998; Haehnelt et al., 1998; Salucci et al., 1999; Wyithe & Loeb, 2002; Aird et al., 2017) the clustering properties of quasars (Martini & Weinberg, 2001), or via hydrodynamic simulations (Hopkins et al., 2005). These models are generally consistent with a picture in which the black hole radiates at or near the Eddington luminosity during the quasar phase for a time of order the Salpeter time, i.e. the mass-independent time it takes to double the black hole’s mass during Eddington-limited accretion with radiative efficiency $\epsilon_{\mathrm{BH}}\sim 10\%$, roughly $t_{0}=40\ \mathrm{Myr}$. In this case, $\Psi_{X}=t_{0}\eta_{X}L_{\mathrm{edd}}/(4\pi r^{2})$.

At the same time, we note that the current mass of the black hole may be written directly in terms of the fluence as follows.

[TABLE]

The remaining fraction, $1-\epsilon_{\mathrm{BH}}$, of this energy is gained by the black hole thereby increasing its mass. Canonically $\epsilon_{\mathrm{BH}}\sim 0.1$. If we assume that both $\epsilon_{\mathrm{BH}}$ and $\eta_{X}$ are constant over the lifetime of the black hole, we can rearrange the previous equation to obtain

[TABLE]

So long as $\eta_{X}$ and $\epsilon_{\mathrm{BH}}$ are roughly constant, this expression will hold regardless of the growth history of the black hole. This is a nice feature when compared to our previous estimate of $\Psi_{X}$, which requires explicitly assuming that the black hole accretes at the Eddington-limited rate for a fixed time $t_{0}$. However, the two estimates are quite closely related, since some methods for estimating $t_{0}$ rely on equation (3). Moreover, the two estimates depend on $r$, $M_{\mathrm{BH}}$, and $\eta_{X}$ in exactly the same way, so we can find a $t_{0}$ for which the two are identical for any given value of $\epsilon_{\mathrm{BH}}$. For the canonical value of $\epsilon_{\mathrm{BH}}=0.1$, the fluences are equal when $t_{0}=50$ Myr, remarkably close to our fiducial choice above.

Because these estimates are so congruous, for convenience and definitiveness we will adopt $t_{0}=40\ \mathrm{Myr}$. This will allow us to make a more sophisticated estimate of $\eta_{X}$, since we will know in what state the black hole is accreting. This will also allow us to easily check which fluences correspond to fluxes above or below the critical values for $\epsilon$ discussed above. A corollary is that fluxes, fluences, and total mass lost are directly proportional to each other over several orders of magnitude in each.

In Table 1, we list various masses whose loss may be important for a planet’s habitability, and the corresponding high-energy fluences and fluxes that would be required to remove that much mass assuming $\Psi_{X}=F_{X}t_{0}$. One can also go in the other direction, starting with a known flux, such as the high-energy flux at Proxima Centauri (Ribas et al., 2017), or the current high energy flux at the earth from the sun (Ribas et al., 2005), or from Sgr A* at Eddington luminosity. To translate between fluxes and fluences, we usually assume our fiducial value of $t_{0}=40\ \mathrm{Myr}$, but for cases in the table where the high-energy flux originates from a star, we adopt a substantially longer timescale of $5\ \mathrm{Gyr}$ to reflect the fact that the planets in question have experienced that flux for a much longer period of time. For the purposes of these simple estimates, we neglect the fact that younger stars have substantially higher fractions of their luminosity emitted in the XUV (Ribas et al., 2005) – for the sun this increases the XUV fluence by about an order of magnitude (Zahnle & Catling, 2017). Values for the current XUV fluxes at the Earth from the Sun, and at Proxima b from Proxima Centauri are taken from Ribas et al. (2005) and Ribas et al. (2017) respectively.

We shall use the first three entries in the Table as references in the figures throughout the rest of the paper. In addition to being a useful reference, this Table also immediately raises an important issue. The current mass loss rate from the Earth’s atmosphere is a factor of 300 lower than predicted by Equation (1), despite the fact that the current XUV flux at the Earth is well within the range for which Bolmont et al. (2017) predict that $\epsilon\sim 0.1$. This is understood to be a feature of the current structure of Earth’s atmosphere, wherein hydrogen and helium only escape by diffusion. With warmer, still formally habitable, surface temperatures, substantially more water vapor could make it to the upper atmosphere where it would be subject to photodissociation and the loss of hydrogen (Kasting et al., 2015). This discrepancy emphasizes the importance of other factors in the structure of the atmosphere which we do not consider. It is therefore prudent to expect that our results will only hold for the subset of planetary atmospheres to which Equation (1) applies.

3 The distribution of xuv fluences from AGN

In the previous section we showed that the mass lost from a given planet’s atmosphere satisfies $M_{\mathrm{lost}}\propto\Psi_{X}=t_{0}\eta_{X}\eta_{\mathrm{Edd}}L_{\mathrm{Edd}}/(4\pi r^{2})$, where $t_{0}$, $\eta_{X}$, and $\eta_{\mathrm{Edd}}$ can each be estimated. This relation can be re-arranged to obtain the three-dimensional distance $r$ between the planet and the black hole at which the planet is subject to a given fluence for a particular choice of $E_{X}\equiv t_{0}\eta_{X}\eta_{\mathrm{Edd}}L_{\mathrm{Edd}}$. Scaling to relevant physical values,

[TABLE]

Clearly this set of values suggests that this effect can be important, at least in the cores of galaxies.

As a first estimate of how relevant this effect is, we adapt a plot from Brodie et al. (2011), namely masses of spherical stellar systems plotted against their effective radii. Lines of constant $\Psi_{X}$ are superimposed, suggesting which systems will contain large fractions of planets strongly effected by XUV radiation from their supermassive black hole. Figure 1 shows the result, where a few systems of known black hole mass and half light radii are highlighted in color. The smaller background points show a large compilation of data spanning globular clusters to elliptical galaxies111http://sages.ucolick.org/spectral_database.html retrieved May 10, 2017 (including data from Barmby et al., 2007; Chilingarian et al., 2009; Harris, 2010; Brodie et al., 2011; Strader et al., 2011; Misgeld & Hilker, 2011; McConnachie, 2012; Usher et al., 2013; Jennings et al., 2014, and many others), where their known V-band luminosities have been translated into black hole masses by assuming that the correlation between spheroid luminosity and black hole mass (here we use Häring & Rix, 2004) established at high masses extends to arbitrarily small masses. This is a highly speculative assumption (see e.g. Kormendy & Ho, 2013; Graham & Scott, 2015; Pihajoki, 2017, for steeper versions of this relation).

We highlight a few systems in Figure 1 that do not require translation through this empirical correlation because their black hole masses are estimated based on direct modeling of the system in question. For the Milky Way, we assume a half mass radius of $0.75\ \mathrm{kpc}$ for the bulge (e.g. Li, 2016), and a black hole mass of $4\times 10^{6}M_{\odot}$ (e.g. Boehle et al., 2016). For M87, we take the black hole mass from Walsh et al. (2013) and the effective radius from Agnello et al. (2014). We also plot the recently-reported intermediate mass black hole from in 47 Tuc Kızıltan et al. (2017) and two black holes in the ultra-compact dwarf galaxies reported in Ahn et al. (2017).

Based purely on Figure 1, focusing largely on the regime $M_{\mathrm{BH}}\ga 10^{6}M_{\odot}$ where the relationship between luminosity and black hole mass is reasonably well-callibrated, we can already see a few important results. The half-mass radii of these stellar systems are generally large enough that many of the planets they presumably contain will be subject to only moderate quantities of XUV fluence – enough to damage the present-day atmosphere of Mars, but not much more. There are several interesting outliers, namely M87, and the ultra-compact dwarfs. These are each peculiar systems in their own way, respectively one of the most massive galaxies in the local universe, and galaxies with remarkably small radii. This suggests that regardless of the fate of the typical planet in the universe, there will be some unlucky enough to reside in unusual galaxies such as these. Globular clusters also stray towards the $\log_{10}\Psi_{X}=17$ line, suggesting that globular clusters may be less hospitable to life than the typical star. Of course it is not clear that any of them actually contain the intermediate mass black holes suggested by the extrapolation, but interestingly the reported result of Kızıltan et al. (2017) is very much in line with the extrapolation.

Figure 1 is a simple first look at how planets in nearby stellar systems might fare under the influence of XUV radiation from their black holes. To take it a step further, we will model the full stellar distributions (rather than simply looking at half-light radii), in order to compute the probability density that a random planet at a particular epoch in the universe will be subject to a fluence $\Psi_{X}$. The first step is to compute the distribution $dP/d\log_{10}\Psi_{X}$ for a single stellar component. With these distributions, we can then construct the full distribution of expected fluences for a random planet in the universe (as opposed to a random planet in a given stellar system) by summing up the individual distributions with appropriate weights.

The fluence distribution for a given stellar system is straightforward to compute. Given that there is a relatively simple relationship between $r$ and $\Psi_{X}$, namely Equation (5), we can change variables so that

[TABLE]

where $M_{*}$ is the cumulative mass in stars, either exposed to a given $\Psi_{X}$ or within radius $r$. The radial mass distribution $dM_{*}/dr$ has a number of well-known forms, and $dr/d\Psi_{X}$ follows immediately from equation (5). We should note that in every plot in this paper, we are actually showing either $dM_{*}/d\log_{10}\Psi_{X}=\ln(10)\Psi_{X}dM_{*}/d\Psi_{X}$, or the probability density per $\log_{10}\Psi_{X}$, namely $dP/d\log_{10}\Psi_{X}=(1/M_{\mathrm{tot}})dM_{*}/d\log_{10}\Psi_{X}$, since the $x$-axis of each plot will be logarithmically scaled. Here $M_{\mathrm{tot}}$ is the total stellar mass in a given system.

For a thin axisymmetric disk with surface density $\Sigma$, we have $dM_{*}/dr=2\pi r\Sigma$, while for a spherical distribution with density $\rho$, $dM_{*}/dr=4\pi r^{2}\rho$. A standard choice for $\Sigma$ is the Sersic density profile (Sérsic, 1963),

[TABLE]

where $\Sigma_{0}$ is a normalization constant, $r$ is the distance from the center of the galaxy, $r_{E}$ is the radius containing half the mass, $n$ is the Sersic index, and $b$ is a function of $n$ set such that $r_{E}$ does in fact contain half the mass. Plugging this profile into equation (6) and adding the appropriate factors to change variables to log-scale, we find

[TABLE]

Here $\Gamma(x)=\int_{0}^{\infty}t^{x-1}e^{-}tdt$ is the gamma function. This equation is quite straightforward to evaluate, and essentially reduces to a power law with an exponential cutoff in the quantity $(E_{X}/r_{E}^{2}\Psi_{X})^{(1/2n)}$. This immediately says that changes in $E_{X}$ and $r_{E}^{2}$ affect the distribution in the same way but with opposite sign – increasing $E_{X}$ or decreasing $r_{E}^{2}$ shift the distribution to higher $\Psi_{X}$. Increasing $n$, meanwhile, increases the width of the distribution. Each of these effects is shown in turn in Figure 2.

Figure 2 also shows the difference between the infinitely thin Sersic profile with $n=4$, and the three-dimensional Hernquist (1990) profile, which when viewed in projection yields an excellent approximation to the $n=4$ Sersic profile. Since the two profiles are nearly the same in projection, clearly each particle in the 3D distribution will be strictly farther away from the center than the corresponding particle in the flattened distribution. It is therefore unsurprising that the flat distribution extends to higher $\Psi_{X}$, while the 3D distribution contains contributions at much lower $\Psi_{X}$. Following the same procedure as for the Sersic profile, but with $dM_{*}/dr=4\pi r^{2}\rho$, and $\rho=M_{\mathrm{tot}}a/(2\pi r(r+a)^{3})$, we find that the fluence distribution is

[TABLE]

Here $a$ is the characteristic scale of the potential-density pair, which is related to the effective radius (i.e. the projected half-mass radius) via $r_{E}\approx 1.8153a$ (Hernquist, 1990). Just as in equation (8), the dimensionless quantity determining the shape of the distribution is $E_{X}/(4\pi r_{E}^{2}\Psi_{X})$, which is the ratio of the fluence in question $\Psi_{X}$ to the fluence at the half-mass radius $r_{E}$.

While these distributions show the fluence expected for a single component of a single galaxy, a random planet in the universe will be located in a multi-component galaxy with a distribution of values of $r_{E}$, $n$ and $E_{X}$. To estimate the fluence distribution for all planets, we proceed via the Monte Carlo method in the spirit of Guillochon & Loeb (2015), namely drawing values of each of these quantities from an appropriate distribution and summing the results.

Explicitly, we proceed using the following notation: $M_{h}$ is the mass of the dark matter halo, $M_{*}$ is the stellar mass, $M_{b}$ is the bulge mass, $M_{d}$ is the disk mass, $SF$ is either true or false, true meaning that the galaxy is star-forming and false meaning that the galaxy is quenched, the bulge to total ratio is $BT=M_{b}/(M_{b}+M_{d})$, and $r_{E,d}$ and $r_{E,b}$ are the effective radii of the disk and bulge respectively. Our calculation proceeds as follows:

Set up a grid of halo masses, equally spaced in $\log_{10}M_{h}$ with spacing $\Delta\log_{10}M_{h}$. For each $M_{h}$… 2. 2.

Record the density of halos per dex in halo mass, $\phi_{10}(M_{h})$ 3. 3.

Draw from the conditional probability distribution $p(M_{*}|M_{h},z)$ 4. 4.

Draw from $P(SF|M_{*},z)$ 5. 5.

Draw from $p(BT|M_{*},SF)$, implying $M_{d}$ and $M_{b}$. 6. 6.

Draw from $p(M_{\mathrm{BH}}|M_{b})$, implying $E_{X}$. 7. 7.

Draw from $p(r_{E,d}|M_{d},SF,z)$ 8. 8.

Draw from $p(r_{E,b}|M_{b},SF,z)$ 9. 9.

Use Equations (8) and (9) to compute $dP/d\log_{10}\Psi$ for both the disk and bulge component respectively.

Each of these distributions is then summed up with the appropriate weight to yield the density of stellar mass in the universe subject to a given fluence, $d\rho_{*}/d\log_{10}\Psi$, which has units of solar masses per cubic Mpc per dex in fluence,

[TABLE]

This quantity is proportional to the number of planets per volume in the universe per dex in fluence, assuming that planet formation has no dependence on additional parameters. The validity of this assumption and the consequences of relaxing it will be discussed in section 5.

In order to proceed, we need to specify the various probability density functions in the list above. For the most part these will be normal distributions in the logarithm of the variables, reflecting the fact that observationally many of these relations are consistent with having a fixed lognormal scatter.

First, we have the halo mass function $\phi_{10}(M_{h})$, where the subscript is a reminder that this quantity is $dN/d\log_{10}M_{h}$, not $dN/dM_{h}$. For these we use the online tool HMFcalc (Murray et al., 2013) with the default settings, namely a Planck-SMT cosmology (Planck Collaboration et al., 2016) cosmology and a Sheth et al. (2001) fitting function, but we expand the tabulated range of masses to $10^{8}$ to $10^{17}M_{\odot}/h$, decrease the tabulated mass intervals, and obtain values at redshifts between $z=0$ and $z=10$ in steps of $\Delta z=1$. The value of the mass function here simply serves as a weight in the sum given by equation (3).

For the distribution of stellar masses given a halo mass, $p(M_{*}|M_{h})$, we employ the redshift-dependent relation from Moster et al. (2013), with a log-normal scatter of $0.15\mathrm{dex}$ independent of $M_{h}$ or $z$. The scatter in this relation has been constrained with a variety of techniques at $z=0$ (Reddick et al., 2013; Zu & Mandelbaum, 2015), but its redshift evolution is not well-constrained. The relationship between scatter and the low-mass slope of the stellar mass halo mass relation is also degenerate - that is, the data are consistent with a family of solutions where the slope and scatter are increased simultaneously (Garrison-Kimmel et al., 2017). However, as we will see, the low-mass end of the relation is not incredibly important to our final result.

Next, we draw a random number to determine if the galaxy is star-forming or quenched. To do so we employ the stellar mass function fits provided by Leja et al. (2015). In particular the stellar mass function is split into star-forming and quenched galaxies, so the probability that a galaxy is star-forming or quenched as a function of stellar mass is just $P_{\mathrm{SF}}=\phi_{SF}(M_{*})/(\phi_{SF}(M_{*})+\phi_{\mathrm{quenched}}(M_{*}))$.

Once the star-forming status of the galaxy is known, we can draw from the observed relationship between bulge mass and stellar mass from Lang et al. (2014), which shows little dependence on redshift. For star-forming galaxies,

[TABLE]

meaning that $BT$ is drawn from a clipped normal distribution, which is a normal distribution with the mean given by the first argument and standard deviation given by the second. The allowed range is given by the subscript, and if a value is drawn outside the allowed range, it is replaced with the appropriate bound, e.g. $BT=0.01$ if the draw from the normal distribution is below $0.01$.

With the mass partitioned, we rely on the observed correlation between bulge mass and black hole mass. In particular we adopt the relation reported in McConnell & Ma (2013), namely

[TABLE]

where $\mathcal{N}$ with no subscript denotes a normal distribution, meaning that $M_{\mathrm{BH}}$ itself is drawn from a log-normal distribution. Note that the scatter in the relation is quite uncertain, with values quoted between 0.2 and 0.6 depending on the statistical method, and with some indication that the scatter is larger at smaller masses. With the mass of the black hole, we can now estimate $E_{X}=t_{0}\eta_{X}\eta_{\mathrm{Edd}}L_{\mathrm{Edd}}$. As discussed in the previous section, we take $t_{0}=40\ \mathrm{Myr}$, $\eta_{\mathrm{Edd}}=1$, and $L_{\mathrm{Edd}}=1.3\times 10^{38}(M_{\mathrm{BH}}/M_{\odot})\ \mathrm{erg}\ \mathrm{s}^{-1}$. For $\eta_{X}$, rather than assuming a constant value, we adopt a model for the spectrum of AGN from Thomas et al. (2016), assuming in their model that the proportion of flux in a non-thermal component is $p_{NT}=0.1$, the photon power-las index of the non-thermal emission is $\Gamma=1.8$, and the coronal radius in units of gravitational radii is $\log_{10}r_{\mathrm{cor}}=1.5$.

The final ingredients entering our model for the cosmological fluence distribution are the effective radii of each component. Once again we rely on the results of Lang et al. (2014)

[TABLE]

for the spheroidal component, or

[TABLE]

for the disk. The medians of all four of these distributions are assumed to evolve as $r_{E}\propto(1+z)^{-1}$ (van der Wel et al., 2014; Allen et al., 2017), where equations (13) and (14) are evaluated at $z=1$.

We carry out this procedure in a number of halo mass ranges, drawing 30,000 samples in evenly spaced intervals of $\log_{10}M_{h}$. The results are presented in Figure 3, which shows $d\rho_{*}/d\log_{10}\Psi_{X}$. As in Figure 2, we include three vertical reference lines corresponding to some relevant masses from our solar system, namely the current mass of the Martian atmosphere, the current mass of Earth’s atmosphere, and the current mass of Earth’s oceans. We show lines for $z=0$, 2, 4, 6, 8, and 10 with increasingly redder lines, and we also separate out the spheroidal components. As the redshift decreases, $d\rho_{*}/d\log_{10}\Psi_{X}$ generally increases at all values of $\Psi_{X}$ as the mass of stars in the universe increases. The contribution to the high-fluence end of the distribution is dominated by spheroids, and non-negligible parts of the distribution lie at values of $\Psi_{X}$ above the reference lines. Based on the drop-off in $\epsilon$ for XUV fluxes below $0.1\ \mathrm{erg}\ \mathrm{s}^{-1}\ \mathrm{cm}^{-2}$, fluences below the line for the Martian atmosphere are unlikely to correspond any longer to the values of $M_{\mathrm{lost}}$ on the lower $x-$axis.

The effect can also be broken down by mass and redshift. Figure 4 shows the fraction of stars at a given halo mass (in a given component at a given redshift) for each of the three reference quantities of mass loss. Graphically one can imagine that these lines are where the cumulative distribution functions (CDFs) computed from Figure 3 cross the vertical reference lines of that figure. Although the effect is slightly different at different redshifts, the trend from Figure 4 is quite clear: planets in higher-mass galaxies are more likely to be subject to high XUV fluences. This is a consequence of the black hole masses in these galaxies growing super-linearly with stellar mass (since the bulge-to-total ratio increases), but the effective radii growing at a much more modest rate.

The redshift dependence, shown in Figure 5, is much less dramatic. In each mass bin, shown by the non-solid lines, there tends to be a moderate decrease in the influence of XUV radiation at lower redshifts. This is the result of the increase in effective radius for galaxies at low redshift. However, over the whole mass range considered (the solid lines), the trend is reversed. Except at the very lowest redshifts, the fraction of stars in the universe subject to high XUV radiation increases with time. This apparent paradox is merely the result of the increasing characteristic mass of the halo mass function – in other words, at low $z$, high-mass galaxies become more common, and these are the galaxies where the largest fraction of stars are subject to high XUV fluences.

4 The influence of non-circular orbits

Implicit in the calculations presented in the previous section is that the stars follow roughly circular orbits, so that their current radii are representative of their typical radii. In other words, we assumed that when computing the fluence $\Psi_{X}=\int L_{X}/(4\pi r^{2})dt$, the $r^{2}$ factor could be taken out of the integral. This is likely to be a good approximation in the disk components of galaxies, modulo the issue of stellar migration (e.g. Sellwood & Binney, 2002), but in spheroidal systems where the velocity distribution is often closer to isotropic, stars can easily traverse substantially different radii over the course of a single orbit.

To get an idea for the order of magnitude of this effect, we examine orbits in the Hernquist potential we adopted in the previous section for the density distribution of spheroidal systems. For each orbit, we compute the ratio of $\int_{0}^{t_{0}}r^{-2}dt$ to $r_{0}^{-2}t_{0}$, where $r_{0}$ is the particle’s initial radius. The latter quantity is proportional to the fluence if the particle is on a circular orbit at radius $r_{0}$, while the former is the analogous quantity taking into account the fact that particle’s radius changes over time. The ratio is therefore equal to $\Psi_{X}/\Psi_{X,\mathrm{circ}}$, the actual fluence over the fluence assuming a circular orbit.

We compute this distribution using a Monte Carlo method. In particular, we draw a large number of particles at random from the 6D phase space density $f$ of the system, and integrate them forward in time for $40\mathrm{Myr}$. During the integrations, in addition to the particle’s state, we integrate $r^{-2}dt$ to obtain the numerator of the ratio defined above. We draw particles from the phase space density in three steps. First, the particle’s radius is drawn from the probability density $4\pi r^{2}\rho/M_{\mathrm{tot}}$. Since $\rho$ has a particularly simple form for the Hernquist profile, this equation can be easily integrated and inverted to solve for $r$ as a function of the CDF, making it trivial to draw a random value of $r$ with the inverse CDF method. Once $r=r_{0}$ is known, we draw the magnitude of the velocity from the conditional probability density

[TABLE]

Once again we have made use of the spherical symmetry to integrate out the spatial angles, and similarly we assume the velocity is isotropic to integrate out the velocity angles. Note that $f$ is normalized such that $fd^{3}\vec{x}d^{3}\vec{v}$ is the mass, as opposed to the probability, contained in the differential 6D volume about a given point. For a given $r$ and $v$, this expression can be evaluated analytically using the expression from Hernquist (1990) for $f$ as a function of the energy per unit mass $E=(1/2)v^{2}+\phi(r)$ where $\phi(r)=GM/(r+a)$ is the gravitational potential. Note that $E$ is strictly a function of $r$ and $v$. We sample from $p(v|r)$ using rejection sampling from a uniform distribution, making use of the fact that $p(v|r)$ is a broad distribution, and there is a finite range of allowed velocities, namely $0<v<\sqrt{2GM_{\mathrm{tot}}/(r+a)}$, for which a star is bound. Finally, the star can be assigned a velocity vector $\vec{v}$ by drawing a random point from the surface of the unit sphere, again by the assumption of isotropy.

Once the star’s initial radius and velocity vector are chosen, we immediately know the star’s specific energy (as in the previous paragraph) and initial specific angular momentum $L=r(v_{x}+v_{y})$ where, because of the spherical symmetry, without loss of generality we have assumed that the star lies on the $z$-axis. On timescales similar to or shorter than the dynamical time, as is the case for our integration, $E$ and $L$ are constants. We can therefore compute $r$ as a function of time by integrating

[TABLE]

The two solutions correspond to the particle moving inward and outward along its orbit. In order to choose the correct root, we flag the initial sign of $v_{r}$ as determined by the initial draw of the velocity, and we change the sign when the argument of the square root becomes negative. In this unphysical regime, we evaluate the square root using the absolute value of the argument, and the particle quickly re-enters the physically allowed region. We verify that the distribution of $r$ for a large sample of these particles agrees with the analytic distribution both in the initial conditions and after the $40$ Myr evolution.

The result of this calculation is shown in Figure 6 for masses between $10^{8}M_{\odot}$ and $10^{12}M_{\odot}$, and $a$’s set according to the median of the distribution for quenched galaxies in equation (13). Although there is a narrow peak at $\Psi_{X}/\Psi_{X,\mathrm{circ}}=1$ containing about 5% of the stars, the full distribution is much broader, with a distribution reasonably fit by a log-normal with median $\log_{10}\Psi_{X}/\Psi_{X,\mathrm{circ}}\approx 0.18$ and width $0.43$ dex. The distribution is also roughly independent of mass. Qualitatively, the distribution of $\Psi_{X}$ itself is so broad that the additional width introduced by the effect of non-radial orbits is not likely to be important.

5 Implications for Habitability

The habitability of any given planet is complex and time-dependent. In some sense, one must understand the planet’s climate over Gyr timescales where there are many more unknowns than in terrestrial climate science – the planet’s mass, spin, inclination, period, surface configuration, atmospheric composition and pressure, and so forth (see section 4.1 of Shields et al., 2016, and references therein). In this section we shall nonetheless comment on a few aspects of habitability implied by our results.

5.1 Galactic Habitable Zones

For many readers, our results may be reminiscent of the concept of the galactic habitable zone (Gonzalez et al., 2001; Lineweaver et al., 2004; Forgan et al., 2017; Spitoni et al., 2017). Just as each star has a circumstellar habitable zone, a set of planetary orbits that are likely to admit stable climates with liquid water at the surface (Kopparapu et al., 2013), it has been suggested that each galaxy has a particular range of radii where life is more likely to develop. In most formulations, it is assumed that the metallicity of the gas from which a planetary system forms must be in a certain range near the solar value. This supposition is based on the fact that terrestrial planets form from heavy elements, so planets are more likely to form when the metallicity is higher. If the metallicity is too high, planet formation may make gas giants more common, leading to the ejection of many terrestrial planets (e.g. Barclay et al., 2017). Suppositions about habitability based on metallicity are quite speculative, since observationally there is no indication of any dependence of terrestrial planet frequency on metallicity (Buchhave et al., 2012; Schlaufman, 2015). There is a well-known correlation between metallicity and gas giant frequency (Fischer & Valenti, 2005), but the consequences for habitability are unclear.

Another key component of galactic habitable zones is the proximity of the planet to supernovae and gamma-ray bursts. Most stellar explosions occur in proportion to the star formation rate averaged over at least 40 Myr timescales, so the centers of galaxies, particularly near $z=2$ are supposed to be particularly dangerous (Lineweaver et al., 2004) though the absolute normalization of the effect, i.e. the exact probability that a supernova at a given distance will render a planet uninhabitable, is quite uncertain. A useful starting point is the computation by Gehrels et al. (2003), which examines the effects on Earth’s atmosphere by running a global circulation model coupled to a photochemical network. The primary effect of a supernova according to this result is the production of odd nitrogen molecules (namely a single nitrogen paired with some number of oxygen atoms), which catalyze the destruction of ozone. The magnitude of the effect is such that a supernova at $\sim 8\ \mathrm{pc}$ would cause order unity destruction of the Earth’s ozone layer over a multi-year period. Although this would undoubtedly cause damage to Earth’s current biosphere, it is less clear that it would render the Earth uninhabitable, much less a generic terrestrial planet, whose atmospheric chemistry and thickness may be substantially different.

The mechanism we discuss in this work, namely the loss of substantial atmospheric mass due to winds driven by XUV radiation from AGN, is arguably a much stronger and more robust effect than the ozone destruction by supernovae typically considered. However, it is unclear whether the levels of mass loss we find in this paper are in fact destructive for life. A useful comparison is Proxima Centauri b, whose host star is prone to X-ray flares, and is therefore likely to have experienced a large fluence of high-energy radiation over the course of its life. In spite of this substantial mass loss, there are many scenarios (Ribas et al., 2016; Barnes et al., 2016) where the planet retains substantial volatiles. On the other hand, it has recently been shown (Jakosky et al., 2017) that mass loss from the Martian atmosphere (through a different mechanism) was likely responsible for a dramatic change in character of the planet, from temperate with substantial liquid water on the surface, to the current cold dusty planet where nearly all of the water near the surface is frozen. It is likely that two different planets subject to the same atmospheric mass loss will experience dramatically different climate and habitability effects, ranging from negligible to dramatic. It is also possible that atmospheric mass loss may increase the habitability of a planet, if the planet initially forms with a substantial hydrogen atmosphere (Luger et al., 2015).

5.2 Implications for Earth

As shown in Table 1, the high-energy flux at Earth’s current distance from the Milky Way’s central supermassive black hole is small, even if Sgr A∗ were emitting an unobscured Eddington luminosity. It follows that even though there are multiple lines of evidence pointing to recent accretion events (Su et al., 2010; Bland-Hawthorn et al., 2013; Churazov et al., 2017), planets in the solar system itself have probably never been subject to this particular mechanism of mass loss. There is at least one caveat to this conclusion, however. While stars in disk galaxies typically follow circular orbits, they may be perturbed when they pass through resonances with spiral arms or bars (Sellwood & Binney, 2002). Simulations show that the metallicity distribution function of the Milky Way has properties (in particular skewness) consistent with a substantial fraction of stars in the outer disk having experienced migration (Loebman et al., 2016). In fact, the high measured metallicity of Alpha Centauri has led Barnes et al. (2016) to suggest that the system has migrated outwards by at least 3.5 kpc. The Sun may have experienced such a process, but it is also quite consistent with having been born near its current galactocentric radius (Allende Prieto et al., 2004). Given the mass of Sgr A∗, and using $\Psi_{X}\sim 10^{17}\mathrm{erg}\ \mathrm{cm}^{-2}$ in Equation (5), planets in the Milky Way are only likely to be affected if their stars orbit within $\sim 0.1$ kpc of the galactic center. It follows that even planets in the Alpha Centauri system, where there is evidence for substantial past radial migration, are unlikely to have been affected by mass loss due to XUV irradiation by Sgr A∗.

Another possibility arises during the course of a merger between two galaxies. Each one initially contains a black hole at its center. Over time, the black holes will be driven to the center of the post-merger galaxy by dynamical friction (Begelman et al., 1980; Kelley et al., 2017) on a timescale comparable to the dynamical time. During this phase, if the black holes are actively accreting, they could essentially extend the reach of their XUV radiation far beyond the center of the galaxy. There are numerous examples of pairs of AGN on kpc scales (Fu et al., 2012; McGurk et al., 2015; Comerford et al., 2015). If the black holes are no more likely than usual to be active as they settle to the center of the galaxy, then the small quasar duty cycle would suggest that this is a minor effect. However, mergers are exactly when the black holes may be most likely to be active. Even so, the effect of these off-center AGN is not likely to be larger than the central AGN simply because the stellar density at the center of the galaxy is larger than in its outskirts.

Another aspect of mergers is the rearrangement of stellar orbits. Although the Milky Way likely has not had a substantial merger in the past 10 Gyr (Stewart et al., 2008), it is destined to merge with Andromeda in a few Gyr (Cox & Loeb, 2008; Sohn et al., 2012; van der Marel et al., 2012). Simulations with orbital parameters consistent with constraints on the relative motion of the Milky Way and Andromeda predict that the merger remnant will be a large elliptical galaxy. The simulations also predict the distribution of final orbits for stars at similar present-day galactocentric radius to the sun. These stars are predicted to have substantially larger semimajor axes than $\sim 8$ kpc (van der Marel et al., 2012), suggesting once again that solar system planets are unlikely to ever face enough XUV radiation from a supermassive black hole to harm their atmospheres.

6 Conclusion

In this paper we have brought together three different astrophysical subjects: the accretion history of supermassive black holes, the spatial distribution of stars in galaxies, and the erosion of planetary atmospheres. We assumed that most galaxies in the universe have a supermassive black hole which has grown largely through the direct accretion of gas, largely over a small number of instances where the accretion luminosity reaches the Eddington limit for a short time, about 40 Myr. Next, we followed observational relationships between stellar and halo mass, bulge to total ratio, the spatial sizes of the bulge and disk, and black hole mass. From these two sets of assumptions, we can compute the PDF of XUV fluences to which random stars in the universe are subject by adding up the analytically-determined PDFs associated with individual stellar components. The scatter in each observational relation is included by a Monte Carlo method, in which the sum is carried out over many draws from each distribution.

The atmospheric mass loss PDF can be computed at a variety of epochs, and the contribution to the PDF from bulges vs. disks, and from different ranges of halo mass can be tracked as well. We have found that the largest fluences are contributed by the highest-mass halos, which have the largest black holes, but not dramatically larger spatial sizes. The PDF changes relatively little across cosmic time because, although galaxies at a given mass are growing in spatial size, the increasing contribution of higher halo masses as the characteristic halo mass increases largely cancels out the effect. In the end the XUV fluence from supermassive black holes is appreciably higher in the spheroidal components of galaxies, and depends strongly on galaxy mass, but only weakly on redshift.

The XUV fluence can be directly related to the mass of atmosphere lost by a planet of a given density over a range of about 4 orders of magnitude in fluence. Using this translation, we find that a substantial fraction of planets in the universe may have lost potentially non-negligible quantities of atmosphere (see Figure 5). In particular, about 50% of all planets in the universe may lose the equivalent of a Martian atmosphere, 10% may lose an Earth’s atmosphere, and 0.2% may lose the mass of Earth’s oceans. Although ultimately the effect on habitability is uncertain, we note that the damage inflicted to planetary atmospheres by XUV irradiation from black holes may be substantially greater than what some past works have assumed would render the Earth uninhabitable.

Acknowledgements

The authors would like to thank Lars Hernquist, Carl Rodriguez, James Guillochon, James Owen, Aaron Romanowsky, and Rory Barnes for helpful discussions. This work has benefitted from numerous freely available resources, including the arXiv, NASA ADS, numpy, scipy, and matplotlib; we thank those responsible for making these resources available. JCF is supported by an ITC Fellowship.

Appendix A Habitable zones around compact objects

The first exoplanets ever discovered were found orbiting a millisecond pulsar (Wolszczan & Frail, 1992), and the case for searching for planetary transits around X-ray binaries has recently been made by Imara & Di Stefano (2017). The influence of XUV radiation of such sources on nearby planets is therefore an interesting question. In particular, how much damage is inflicted on the atmosphere of a planet in the habitable zone of an object emitting high energy radiation?

Since the location of the habitable zone is, to a first approximation, just set by the bolometric flux the planet receives, the corresponding XUV flux is just the fraction of the bolometric luminosity emitted in high-energy radiation, $\eta_{X}$, times the flux corresponding to habitability. For reference, the bolometric flux at Earth is $F_{E}=1.4\times 10^{7}\mathrm{erg}\ \mathrm{cm}^{-2}\ \mathrm{s}^{-1}$. For $\eta_{X}\ga 10^{-3}$, we are well within the regime where the mass loss scales as $F_{X}^{1/2}$ instead of directly with $F_{X}$. Following the results of Bolmont et al. (2017), we will reduce $\epsilon$ by $\sqrt{F_{X}/10^{3}\ \mathrm{erg}/\mathrm{cm}^{2}/\mathrm{s}}$, in which case the mass lost over a period of 1 Gyr is

[TABLE]

It follows that for planets around such objects, water may be liquid on their surfaces, but it is unlikely to remain there for very long.