The Sub-Neptune Desert and Its Dependence on Stellar Type: Controlled by Lifetime X-Ray Irradiation

With the exception of the Sun, EUV (10–124 nm) observations of stars are scarce. While studies have used models derived from solar observations to estimate XUV fluxes for old, Sun-like stars (Lecavelier Des Etangs 2007; Sanz-Forcada et al. 2011), such extrapolations are not possible for other stellar types. Thus, we focus specifically on the evolution of X-rays. X-rays are the primary driver for hydrodynamic escape (Owen & Jackson 2012). Furthermore, the XUV emission for Sun-type stars (Sanz-Forcada et al. 2011) and the trends of the near and far-ultraviolet for M dwarfs (Shkolnik & Barman 2014) are all suggestive of the XUV evolution qualitatively following the same saturation and decay behavior as the X-rays.

To track the evolution of X-ray flux as a function of time for different stellar spectral types, we use the observationally derived relations in Jackson et al. (2012) and Shkolnik & Barman (2014). Jackson et al. (2012) used X-ray survey observations of open clusters to report the ratio of the X-ray to bolometric luminosity across stellar ages of 10 Myr to 4.5 Gyr, with bins spanning stellar spectral types from K5 to F0. Shkolnik & Barman (2014) used ROSAT (Röntgensatellit) observations to determine the X-ray to J-band (centered at 1.25 μm) flux ratio for M4 to K4 spectral types with stellar ages of 10 Myr to 5 Gyr. Combining these two data sets, the rough range of stellar masses for which the lifetime X-ray evolution has been measured is 0.3–1.4 M_⊙. These data measure the steady-state X-ray emission and do not account for additional X-ray emissions from large flares or coronal mass ejections.

The same trend of a saturated phase of X-ray emission early in the lifetime of the star, tens to hundreds of millions of years, followed by a decay is found across all stellar spectral types (Figure 1). The main difference between spectral types is twofold: as stellar mass decreases, stars saturate with a greater proportion of their luminosity in the X-rays, while the saturation phase also increases in duration (Jackson et al. 2012; Shkolnik & Barman 2014). Both trends are visible in Figure 1. Although this behavior is not monotonic across spectral types as a result of the uncertainties in the observations, the overall trend is apparent. We note that the binning by spectral type of Jackson et al. (2012) was carried using B − V color—we convert the main-sequence B − V color bins to stellar mass bins solely for reference in Figure 1. The conversion from B − V color to stellar mass varies as a function of stellar age, and the full treatment for our calculation of lifetime-integrated X-ray fluxes is more complex and described in detail in Section 2.1.3.

Zoom In Zoom Out Reset image size

Both Jackson et al. (2012) and Shkolnik & Barman (2014) parameterized the X-ray to bolometric luminosity ratio (L_X/L_bol) and X-ray to J-band luminosity ratio (L_X/L_J), respectively, as a function of stellar age (for a given spectral type) by fitting a broken power law. The saturated phase is fit by a flat line, and the subsequent decay is fit with a power law. The uncertainties in the X-ray measurements of open clusters and field stars, which are used to define individual data points, can span orders of magnitude for the least precise measurements, with the precision varying as a function of stellar age and spectral type—both parameters are independent variables in our analysis. In order to accurately account for these errors in our calculations, we carry out a Markov-chain Monte Carlo (hereafter MCMC) fit to all data points in each of the stellar bins defined by Jackson et al. (2012) and Shkolnik & Barman (2014). In carrying out the MCMC, we take the least-squares fit of a broken power law (using a piecewise function to fit to all data points simultaneously) to the data and their errors in L_X/L_bol to be the maximum likelihood result.

The fit itself is carried out in the log(age(yr)) and log(L_X/L_bol) parameter space, such that we are fitting to the data with a piecewise function of two lines (as in Figure 1), one with a slope of zero representing the saturation phase and one a line of negative slope representing the decay in X-ray luminosity as a function of time. The parameters that we fit for are the log(L_X/L_bol) intercept for the line of negative slope, the logarithm of the time at which the saturation phase ends (log τ_sat), and the slope of the line representing the X-ray luminosity decay—under the condition that the piecewise function is continuous. These fits, which are carried out in logarithmic space, are then converted to the following: the ratio of the X-ray to the bolometric, or to the J-band for Shkolnik & Barman (2014), luminosity during the saturation phase (L_X/L_bol); the time at which the saturation phase ends τ_sat; and the power-law index α for the decay in L_X/L_bol after τ_sat. This approach returned fit parameters in many cases identical to those reported in Jackson et al. (2012) and Shkolnik & Barman (2014). So as to ensure agreement of our maximum likelihood fits with those reported in Jackson et al. (2012), we carried out the "lower weighting" of data points referenced with arrows described in Figure 2 of Jackson et al. (2012)—we found that multiplying the errors for these data points by a factor of 5 produced fit parameters in closest agreement with theirs.

With our maximum likelihood fit parameters in hand, we sampled the parameter space with an MCMC of 100 samplers over 1200 steps using the emcee routine of Foreman-Mackey et al. (2012). Discarding the first 200 steps during which the samplers have not yet sampled the full parameter space, we are left with 100,000 fits, comprising three posterior probability distributions, one for each fit parameter. We randomly sample the fit parameters for 1000 of the fits, taking note that this number of subsamples still accurately represents the shape of the original posterior probability distributions. These 1000 fits are used to account for the effects that the errors in the X-ray measurements have on our calculation of the lifetime-integrated X-ray luminosities and fluxes.

We ultimately seek to calculate estimates of the lifetime-integrated X-ray flux for all planets in the Kepler sample. However, the calculation of the lifetime-integrated X-ray flux for each planet, as well as a Monte Carlo sampling of each planet's errors, is precluded by the time-intensive nature of the calculation. Because the resolution of the stellar masses in the isochrones that we use as inputs as well as the ages at which lifetime-integrated X-ray fluxes have been measured are known, we are aware of the maximum resolution at which lifetime-integrated X-ray luminosities can be accurately calculated as a function of stellar mass and age. For this reason, we precalculate a lookup table—an array of lifetime-integrated X-ray luminosities for given stellar masses and ages, which can be quickly queried for use in estimating the lifetime-integrated X-ray flux for each Kepler planet.

As discussed in Section 2.1.2, the data on X-ray evolution are reported in terms of L_X/L_bol or L_X/L_J. As such, it becomes necessary to model the temporal evolution of the bolometric and J-band luminosities for the stellar types of interest, in order to determine the temporal evolution of the X-ray luminosity itself. Stellar isochrone models are well suited for this, and we use the PARSEC isochrones obtained from the CMD 2.8 web interface (Bressan et al. 2012; Chen et al. 2014). We utilize the PARSEC isochrones for their fine mass resolution (∼0.03) over our 0.3 < M_⊙ < 1.4 range of interest, as well as the availability of tracks for stars with nonsolar metallicities.

In generating our lookup table, we define 40 linearly spaced mass bins over the solar-mass range 0.3 < M_⊙ < 1.4, as well as 50 different stellar ages up to which we integrate. Twenty ages are logarithmically spaced between 200 Myr (limited on the low end by the lowest ages in the PARSEC isochrones) and 1 Gyr, and the final 30 are linearly spaced from 1 to 10 Gyr so as to finely sample the age range in which most Kepler planets fall (Johnson et al. 2017; also see Figure 3). Finally, five bins in metallicity are specified at Z = 0.0025, 0.0075, 0.0125, 0.0175, and 0.03.

For each mass, age, and metallicity bin, we track the star's luminosity and B − V color at each time step in the isochrone (spaced as log(time(yr)) = 0.03). At each time step, we use the B − V color to place the star in the appropriate bin of X-ray evolution as described in Section 2.1.2. We then multiply the isochrone luminosity by one of the L_X/L_bol ratios among the 1000 fits generated for this bin in Section 2.1.2. This leaves us with L_X as a function of time, which we integrate over the lifetime of the star by quadrature (a numerical integration using polynomial interpolation functions) to get a lifetime-integrated X-ray luminosity. For stars of M_* > 0.9 M_⊙ that are older than 1 Gyr, we monitor the evolution of the star's stellar effective temperature to determine whether it leaves the main sequence during its lifetime. We consider a greater than 10% change in stellar effective temperature versus log(time(yr)) to be a signature that the star has left the main sequence. This criterion accurately captures the point in time before dramatic changes to the star's other physical parameters, e.g., R_*, also occur. We stop the integration of lifetime-integrated X-ray flux at the time the star leaves the main sequence, due to the lack of comprehensive observations of X-ray emissions from post-main-sequence stars. This process of calculating lifetime-integrated X-ray luminosity is repeated for all mass, age, and metallicity bins, and then repeated for each one of the 1000 sets of X-ray evolution fit parameters from Section 2.1.2. This effectively provides 1000 realizations of the lookup table that sample the extent of the errors in the X-ray measurements.

In Figure 2, we show our lookup table for our closest-to-solar metallicity bin (Z = 0.0175). Figure 2(a) shows the mean lifetime-integrated X-ray luminosity of the 1000 realizations in X-ray evolution fit parameters, while Figure 2(b) shows the relative error in the lifetime-integrated X-ray luminosities. The physical implications of the results in the lookup table will be discussed in Section 3.1. Over the majority of the parameter space, relative errors are around 30%, although for stars of M_* ∼ 0.7 M_⊙ older than 3 Gyr, relative errors are greater than 100%. This is due to poor constraints on the post-saturation phase X-ray evolution for this particular bin of 0.935 <B − V < 1.275 (see Jackson et al. 2012 Figure 2(f)). For reference and to facilitate comparisons with the results of evolution models, we also show a lookup table of instantaneous X-ray luminosity (Figure 2(c)). This uses the same methodology used to generate the lifetime-integrated X-ray luminosity lookup table, but without an integration over time.

Zoom In Zoom Out Reset image size

Although we utilize and report results using the PARSEC isochrones, due to their mass and metallicity resolution, as well as their direct outputting of B − V color, we have verified our results by repeating the above analysis using the Baraffe et al. (2015) stellar isochrones, combined with estimates of B − V color from model stellar effective temperature using the observations of Pecaut & Mamajek (2013). We find that our results are largely independent of the isochrones used.

In this section, we calculate the lifetime-integrated X-ray flux for stars in the Kepler sample. We utilize two sets of Kepler data. For the purposes of defining the photoevaporation desert with the greatest possible precision, we use data from the CKS (Johnson et al. 2017; Petigura et al. 2017). The CKS used follow-up spectroscopic observations to improve the precision on the stellar parameters of 1305 Kepler stars, and in turn the measurements of 2025 Kepler planets. The CKS measurements have resulted in a factor of ∼3 improvement in the precision in the measurement of stellar masses, as well as a factor of ∼2 improvement in the stellar radius precision that has resulted in a factor of ∼2 improvement in precision for the planetary radii (using the planet-to-star radius ratios from the latest Q1-17 DR 25 Kepler catalog; Johnson et al. 2017). Furthermore, the CKS is the first comprehensive survey of Kepler stellar ages (see Figure 3), increasing the sample size of stars with ages by almost two orders of magnitude over previous studies. This enables, for the first time, estimates of parameters, such as the lifetime-integrated X-ray flux, to be calculated for a majority of planets in the Kepler data set.

Zoom In Zoom Out Reset image size

Because the CKS is a magnitude-limited survey (less than magnitude 14.2 in the Kepler bandpass; Petigura et al. 2017), it excludes most low-mass stars from the Kepler sample and consists largely of Sun-like stars. To explore differences in the photoevaporation desert and valley as a function of stellar type, we use the latest release of the general Kepler data set: the Q1-17 DR 25 catalog (Thompson et al. 2018). The recent calculation of uniform false-positive probabilities for all Kepler Objects of Interest in the Q1-17 DR 25 catalog using the methodology of Morton et al. (2016) has resulted in a large increase in the number of confirmed Kepler planets, totaling 3182 planets.  The stellar parameters utilize those corrected from the original DR 25 catalog, i.e., the DR 25 Supplemental Stellar data.

The relative errors on the properties of the Kepler planets and host stars can be large. Even with the improved precision of the CKS, the typical 5% errors on stellar masses, 11% uncertainty factor on stellar radius, and factor of 2 uncertainty in ages can result in significant variations in estimates of the lifetime-integrated X-ray flux for a given planet. The effect of these errors on defining a region free of planets (the photoevaporation desert) must be taken into account, especially when dealing with samples limited to several dozen to a couple hundred planets (the reason for the smaller samples of planets will be discussed in Section 2.2.3). The consideration of errors becomes even more important when using the regular Kepler sample to compare properties dependent on host-star stellar type, in which typical errors on planet radii are ∼40% (Johnson et al. 2017).

In consideration of these errors, as well as the fact that our calculations of lifetime-integrated X-ray luminosities rely on the numerical integration of binned observational data (making it impractical to propagate errors), we take a Monte Carlo approach to modeling the errors on planets. We replace each planet in the specific data set being considered (CKS or Q1-17 DR 25) with 1000 "fractional planets," each weighted to be 1/1000th of a planet. The physical parameters for each of the 1000 planet fractions and their respective host star are determined by taking posteriors directly from the relevant table: the Q1-17 DR 25 stellar catalog (Mathur et al. 2016), the Q1-17 DR 25 planet catalog (Hoffman & Rowe 2017), or the CKS (Johnson et al. 2017; Petigura et al. 2017, posteriors were provided by Erik Petigura upon request). The lifetime-integrated X-ray flux for each fractional planet is determined by using the lookup table described in Section 2.1.3 to get a lifetime-integrated X-ray luminosity for its MC'd host star. One thousand realizations of the lifetime X-ray luminosity lookup table were previously generated to account for uncertainties in the X-ray observations. Each of the 1000 realizations of the lookup table is used for each of the 1000 planet fractions representing one planet, such that the uncertainties in the lifetime X-ray luminosity are propagated. The time-integrated X-ray luminosity for a fractional planet is then divided by the posterior for its semimajor axis to get a time-integrated X-ray flux.

For the CKS data, posteriors of the ages of each planet's host star are available and drawn from for each fractional planet. In the case of the Q1-17 DR 25 catalog, ages have not been measured. While many of the G-type host stars in the catalog have ages that have been measured by the CKS, the majority of the M&K-type stars, which we wish to compare to, have no age information. In order to treat the ages, and in turn the lifetime integration over X-ray luminosity, consistently for all stellar types in the Q1-17 DR 25 catalog, we take the approach of randomly drawing an age for the host star of each fractional planet in the catalog from the full distribution of ages of the CKS sample (whereby the errors in the CKS ages are accounted for by including 1000 age posteriors for each star in the CKS sample). Because the ages measured by the CKS form a single unimodal distribution around ∼6 Gyr (Figure 3), a random drawing of ages for each host star in the Q1-17 DR 25 catalog should produce similar average behavior over a large sample. 

Using the lifetime-integrated X-ray fluxes calculated for each of the 1000 fractional planets, we calculate planet occurrences in the parameter space of planet radius (R_p) versus lifetime-integrated X-ray flux.

Planet occurrence calculations must account for the completeness of the Kepler data set (Howard et al. 2012), i.e., whether every planet is likely to have been detected around the sample of stars of interest. We therefore focus our analysis on a subsample of the Kepler data set that is nearly complete. To obtain this sample, we follow the methodology of Wolfgang & Lopez (2015)  and perform a series of cuts in orbital period, planet radius, stellar radius, and stellar noise (combined differential photometric precision, or CDPP). Specifically, we restrict the range of orbital periods to <10 days, the planet radii to R_p > 1.2 R_⊕, the stellar radii to R_* < 1.2 R_s, and the CDPP binned to the timescale of the planets' transit <200 ppm. Our CDPP cut is less restrictive than that of Wolfgang & Lopez (2015) (200 ppm versus 100 ppm) because we are focusing on planets with shorter orbital periods (10 versus 25 days). To assess the completeness of this subsample, we compared it to estimates of completeness from Thompson et al. (2018). For the subsample of planets orbiting FGK dwarfs in the Q1-17 DR 25 catalog, with period <100 days and R_p < 6 R_⊕ (most similar to the sample in our study), Thompson et al. (2018) found that the lowest completeness was found to be 89%, with the highest being 96%. This equated to a 7% variation in completeness across the sample. Because sample incompleteness increases with orbital period, and the largest period in our sample is 10 days versus the 100 days in Thompson et al. (2018), these completeness estimates are conservative. Relative to other sources of uncertainty in our sample, the estimated errors due to pipeline incompleteness are small. The relative error due to lifetime-integrated X-ray flux is near 30% (see Figure 2(b)), and the relative errors due to binomial statistics are greater than 20% (Figure 7(b)). We therefore conclude that the effects of pipeline incompleteness will not significantly bias the main conclusions of our analysis.

The number of planets, as well as stars, left by the above cuts is tabulated in Table 1 for the two data sets we use, as well as their stellar subsamples. Around 150 planets form the sample with host stars of all spectral types for both the CKS and the Q1-17 DR 25  data sets. The extremely small sample size of planets around F-type stars (6.2) is a result of a large fraction of these stars having left the main sequence for the Kepler ages and in turn being discarded. We include our analyses on the F-type stars for reference, but focus our interpretations on the M&K- and G-type subsamples due to their larger sample sizes.

We create two grids to calculate planet occurrence: one in the parameter space of planet radius and lifetime-integrated X-ray flux, and the other in planet radius and bolometric flux space, for comparison with previous studies. We define our grids so that the parameter space defined by our bounds in bolometric flux are equivalent to the bounds defined for lifetime-integrated X-ray flux, for a solar-mass star of age 5 Gyr (defined at the low end to correspond to our 10 days period cut, at the high end based on the maximum lifetime-integrated X-ray flux calculated for the Kepler planets). To calculate the occurrence within a given grid cell, we adopt the methodology of Howard et al. (2012). We reproduce Equation (2) from Howard et al. (2012), where the average planet occurrence within a cell of R_p and lifetime-integrated X-ray flux is

$$\begin{eqnarray}&&{f}_{\mathrm{cell}}=\displaystyle \sum _{j=1}{n}_{{pl},\mathrm{cell}}\displaystyle \frac{1/{p}_{j}}{{n}_{* ,j}},\end{eqnarray} \tag{ 1 }$$

where p_j = (R_*/a)_j is the transit probability for a fractional planet j, with a being its semimajor axis. Because we sum over fractional planets, we weight each fractional planet to be 1/1000th of a full planet. n_*,j is the number of stars that pass the cuts described in Section 2.2.3 for the transit duration of the specific fractional planet under consideration.

Again, following Howard et al. (2012), we calculate the error in f_cell using binomial statistics. The standard deviation of f_cell is calculated from the probability distribution of drawing n_pl,cell planets from n_*,eff,cell = n_pl,cell/f_cell effective stars:

$$\begin{eqnarray}&&{f}_{\mathrm{cell},\mathrm{std}}=\sqrt{{n}_{* ,\mathrm{eff},\mathrm{cell}}\,{f}_{\mathrm{cell},\mathrm{std}}(1-{f}_{\mathrm{cell},\mathrm{std}})}.\end{eqnarray} \tag{ 2 }$$

In reporting occurrences, we choose to report occurrence densities (Foreman-Mackey et al. 2014), where we are dividing the occurrence in each grid cell by the area of the grid in dR_p dFlux. This way, the values that we report are not dependent on the specific grid resolution that is chosen.

We wish to search for a region of the parameter space with no planets, an occurrence desert. We expect a desert to exist for gaseous planets at high fluxes, due to atmospheric escape (e.g., Lopez & Fortney 2013; Owen & Wu 2013; Chen & Rogers 2016). Given the large uncertainties on planet properties and fluxes, no region of the occurrence grid is exactly zero. To define the desert, we take the largest region of the occurrence grid that is consistent with zero at 2σ confidence. To find this region, we adopt the following algorithm. We start by taking note of the planet occurrence and its 2σ deviation in the cell of largest planetary radii (3.9–5 R_⊕) and highest flux value (i.e., the upper right corner of the plots in Figure 7). Due to the low occurrence and large uncertainties arising from the small number of planets in this cell (typically a handful of fractional planets), its occurrence is always consistent at the 2σ level with an occurrence of zero. This cell is considered the starting point for our desert. We then consider the occurrences in each of the grid cells sharing one border with our initial desert grid cell. The cell with the lowest occurrence is added to the desert. The total occurrence in the new desert consisting of two cells is calculated, and the standard deviation of the desert occurrence is recalculated using binomial statistics for all fractional planets in the desert grid cells, using Equation (2). This process is repeated until we have drawn the largest possible region that is still consistent with an occurrence of zero at the 2σ level. We define this region to be the desert.

We initially calculate an occurrence grid with a resolution chosen such that the largest relative errors in occurrence per grid cell for the full CKS and Q1-17 DR 25  data sets are one order of magnitude; we use this grid in studies of the photoevaporation desert. In investigating the photoevaporation valley, we select a much finer resolution that lends itself to much larger relative errors in individual grid cells. However, as our interest for the valley is in occurrence summed over the radius axis, the magnitude of the relative errors for the summed occurrences still meet our criteria of <1 dex.

We facilitate the comparison of the photoevaporation desert among different parameters, present-day bolometric flux versus lifetime-integrated X-ray flux, as well as among different subsamples of stellar type and planet radius, by examining our full data set of planets in one dimension. We take all fractional planets from our data set of interest with R_p: 1.8 < R_⊕ < 4 (our definition of sub-Neptunes), and generate a cumulative distribution function (CDF) over the flux dimension, beginning at the fractional planets with the highest fluxes.

We wish to compare the CDFs for bolometric flux versus lifetime-integrated X-ray flux, which are parameters with different units. For this comparison, the bolometric flux data are renormalized to the lifetime-integrated X-ray flux data such that they have identical values at the 10th and 90th percentile of flux.

To determine whether the given CDFs in different physical parameter spaces, or for different subsamples, have statistically significant differences, we utilize the Anderson–Darling test (hereafter the AD test) for two samples. The AD test was developed as an alternative to the Kolmogorov–Smirnov test. The AD test was designed to be sensitive to differences in the tails of two CDFs, making this test appropriate for our analysis of the photoevaporation desert, which is represented in the high-flux tail. The Kolmogorov–Smirnov test is most sensitive to differences in the center of the distribution (Anderson & Darling 1954). Nevertheless, we confirm our results regarding statistical differences in the distributions with the Kolmogorov–Smirnov test and find consistent conclusions. We choose, however, to report results for the AD test, for the reasons cited above.

To compare the shape of the transition from the desert to the center of the planet distribution, we also consider the probability distribution function (PDF)—simply histograms of the planet distributions in flux space.