Is Proxima Centauri b habitable? -- A study of atmospheric loss

TL;DR

This study evaluates whether Proxima Centauri b can retain its atmosphere by modeling ion escape rates, revealing high atmospheric loss potential especially if unmagnetized, impacting its habitability prospects.

Contribution

It adapts a multi-species MHD model to assess atmospheric loss on Proxima Centauri b, providing new insights into its potential habitability based on magnetic field presence.

Findings

Ion escape rates are about 100 times higher than Earth's if unmagnetized.

Magnetic fields can reduce escape rates but remain high compared to Solar System planets.

Results are uncertain due to unknown planetary parameters.

Abstract

We address the important question of whether the newly discovered exoplanet, Proxima Centauri b (PCb), is capable of retaining an atmosphere over long periods of time. This is done by adapting a sophisticated multi-species MHD model originally developed for Venus and Mars, and computing the ion escape losses from PCb. The results suggest that the ion escape rates are about two orders of magnitude higher than the terrestrial planets of our Solar system if PCb is unmagnetized. In contrast, if the planet does have an intrinsic dipole magnetic field, the rates are lowered for certain values of the stellar wind dynamic pressure, but they are still higher than the observed values for our Solar system's terrestrial planets. These results must be interpreted with due caution, since most of the relevant parameters for PCb remain partly or wholly unknown.

Tables

Table 1. Table 1: Chemical reactions and associated rates (the reaction rates are adopted from Schunk & Nagy 2009 ).

Chemical Reaction		Rate Coefficient 222Electron impact ionization is neglected in the calculation, H+ density is from the stellar wind, the neutral hydrogen is neglected.
Primary Photolysis and Particle Impact333The photoionization frequencies are rescaled to PCb values, using the EUV estimate from Ribas et al. (2016), the latter being around 33 times that received by the Earth. in s-1
	CO2 + $h\nu$ $\to$ CO${}_{}{}^{+}{}_2{}^{}$ + $e^{-}$	5.55 $\times$ 10-5
	O + $h\nu$ $\to$ O+ + $e^{-}$	2.08 $\times$ 10-5
Ion-Neutral Chemistry in cm3 s-1
	CO${}_{}{}^{+}{}_2{}^{}$ + O $\to$ O${}_{}{}^{+}{}_2{}^{}$ + CO	$1.64\times {10}^{-10}$
	CO${}_{}{}^{+}{}_2{}^{}$ + O $\to$ O+ + CO2	$9.60\times {10}^{-11}$
	O+ + CO2 $\to$ O${}_{}{}^{+}{}_2{}^{}$ + CO	$1.1\times {10}^{-9}$ (800/Ti)0.39
	H+ + O $\to$ O+ + H 444Rate coefficient from Fox & Sung (2001).	$5.08\times {10}^{-10}$
Electron Recombination Chemistry in cm3 s-1
	O${}_{}{}^{+}{}_2{}^{}$ + $e^{-}$ $\to$ O + O	$7.38\times {10}^{-8}$ (1200/Te)0.56
	CO${}_{}{}^{+}{}_2{}^{}$ + $e^{-}$ $\to$ CO + O	$3.10\times {10}^{-7}$ (300/Te)0.5

Table 2. Table 2: The stellar wind input parameters at PCb (Garraffo et al., 2016 ) for two case studies in the PSO coordinate system. Case 1 (C1) corresponds to the maximum P dyn and P tot over one orbital period of PCb. Case 2 (C2) corresponds to minimum P dyn and P tot , but with the maximum P mag .

	nsw (cm-3)	Tsw (K)	vsw555The $y$-component of the velocity arises primarily from the orbital motion of PCb. (km/s)	IMF (nT)
C1	21400	8.42$\times$105	(-833, 150, 0)	(0, 0, -227)
C2	2460	9.53$\times$105	(-1080, 150, 0)	(0, 0, -997)

Table 3. Table 3: Ion escape rates in sec -1 .

	O+	O${}_{}{}^{+}{}_2{}^{}$	CO${}_{}{}^{+}{}_2{}^{}$	Total
PCb with 1 bar surface pressure
C1-UnM666Unmagnetized Case 1.	1.8$\times$1027	2.4$\times$1026	3.3$\times$1026	2.4$\times$1027
C2-UnM	1.1$\times$1027	9.5$\times$1025	8.2$\times$1025	1.3$\times$1027
C1-M777Magnetized Case 1.	7.3$\times$1026	5.4$\times$1026	5.8$\times$1026	1.8$\times$1027
C2-M	5.9$\times$1025	8.7$\times$1025	5.3$\times$1025	2.0$\times$1026
PCb with 93 bar surface pressure888Scale height, HPCb, is still equal to 0.85 HVenus.
C1-UnM93	3.7$\times$1027	4.1$\times$1024	1.4$\times$1023	3.7$\times$1027

Equations

$$[C{O}_2]$$

$$[O]$$

$$H_{\mathrm{PCb}}=\frac{kT}{mg}=0.85H_{\mathrm{Venus}},$$

$$[C{O}_2]$$

$$[O]$$

Full text

1 Introduction

Currently, thousands of exoplanets with diverse characteristics have been detected (Winn & Fabrycky, 2015). Arguably, one of the most important reasons for studying them is to determine whether they could be habitable. The Habitable Zone (HZ) is broadly defined as the region around a star where a planet can support liquid water for certain values of the atmospheric pressure (Lammer et al., 2009), and has been extensively investigated (e.g. Kasting et al., 1993; Kasting & Catling, 2003; Kopparapu et al., 2013). As per current estimates, it is expected that there are $\sim 10^{10}$ Earth-sized planets that lie within the HZ of their stars, in the Milky Way (Marcy et al., 2014; Winn & Fabrycky, 2015).

Most of the attention has focused on HZ studies of exoplanets around M-dwarfs for two reasons: (i) Earth-sized planets in the HZ are quite common (Dressing & Charbonneau, 2015), and (ii) detailed characterizations of their atmospheres are easier since they are situated much closer to the host star (Seager & Deming, 2010; Shields et al., 2016). Hence, there have been numerous HZ studies of M-dwarf exoplanets (Tarter et al., 2007; Scalo et al., 2007; Shields et al., 2016). This area received a major impetus recently with the potential discovery of a super-Earth (a planet slightly larger than the Earth) orbiting Proxima Centauri (Anglada-Escudé et al., 2016), which was christened Proxima Centauri b (PCb). The result is significant because Proxima Centauri is the nearest star to the Earth, and plans are already underway for exploring PCb.111https://breakthroughinitiatives.org/Initiative/3

Exoplanetary atmospheres are an important area of research since spectroscopic analyses can potentially reveal the presence of life via biosignatures (Seager & Deming, 2010). This raises an important question: at what rates are the atmospheres of close-in exoplanets (such as those in the HZ of M-dwarfs) eroded by the stellar wind? The question is of paramount importance, as high escape rates might deplete the planet of its atmosphere after a short period of time, if one excludes mechanisms such as outgassing (Kasting & Catling, 2003).

In our own Solar system, the importance of atmospheric loss mechanisms has been thoroughly documented for Mars (e.g. Lammer, 2013; Lillis et al., 2015; Dong et al., 2015b). One of the primary scientific goals of NASA’s Mars Atmosphere and Volatile EvolutioN (MAVEN) mission is to determine the escape rates, which are crucial because of their potential impact on the long-term evolution of the Martian atmosphere (e.g., loss of water). Brain et al. (2015) estimated a net ion escape rate of $\sim$2.5 $\times$ 1024 s*-1* by choosing a spherical shell at $\sim$1000 km above the planet with energies $>$25 eV over a 4-month MAVEN period. The ion loss increases by more than one order of magnitude during interplanetary corona mass ejection events (Dong et al., 2015a; Jakosky et al., 2015).

There exist several papers that have explored atmospheric losses from HZ exoplanets (e.g. Khodachenko et al., 2007; Tian, 2009; Erkaev et al., 2013), some of which relied upon hydrodynamic models. As we have learned from our Solar system, realistic estimates appear to necessitate sophisticated magnetohydrodynamical (MHD) models for studying the stellar wind-exoplanet interaction, especially for Earth-sized (or above) planets where neutral losses (e.g., thermal and photochemical escape) are not the dominant mechanism (Brain et al., 2016). Such studies have been quite sparse in exoplanetary research as they necessitate highly detailed numerical codes. However, recent publications (e.g Cohen et al., 2014; Kislyakova et al., 2014; Cohen et al., 2015) on magnetized and unmagnetized exoplanets are notable for tackling this issue.

In this Letter, we investigate atmospheric losses arising from stellar wind interaction with a planet that has PCb-like parameters. Since PCb may lack a magnetic field (Barnes et al., 2016), we examine both the magnetized and unmagnetized cases. Although our final results are outwardly simple, we note that this interaction process is actually highly complex in nature. Our model self-consistently accounts for the ionospheric chemistries and electromagnetic forces (Ma et al., 2004, 2013).

2 The simulation setup

In this Section, we describe our code, and the rationale behind our choices of the parameters.

2.1 Physical model and computational methodology

The 3-D Block Adaptive Tree Solar-wind Roe Upwind Scheme (BATS-R-US) multi-species MHD (MS-MHD) model was initially developed for Mars (Ma et al., 2004) and Venus (Ma et al., 2013). The latter is employed herein, and the neutral atmospheric profiles are based on the solar maximum conditions. The MS-MHD model solves a separate continuity equation for each ion species, and one momentum and one energy equation for the four ion fluids H+, O+, O${}_{2}^{+}$, CO${}_{2}^{+}$ (Ma et al., 2004, 2013). In contrast to most global magnetosphere models applied to Earth that start from $2$-$3$ Earth radii, the Mars/Venus MS-MHD model contains a self-consistent ionosphere, and thus the lower boundary extends down to 100 km altitude above the planetary surface. The MS-MHD model includes ionospheric chemical processes such as charge exchange and photoionization.

The chemical reactions and rate coefficients are based on Schunk & Nagy (2009), except the photoionization rates; the latter are rescaled to PCb values based on the EUV estimate in Ribas et al. (2016). The O+, O${}_{2}^{+}$, CO${}_{2}^{+}$ ion densities at the model lower boundary (i.e., on the 100 km altitude sphere) satisfy the photochemical equilibrium condition, where float boundary conditions for the velocity $\mathbf{u}$ and the magnetic field $\mathbf{B}$ have been applied. Given the high collision frequency near the inner boundary, ions, electrons and neutrals are expected to have roughly the same temperature. Thus, the plasma temperature (sum of ion and electron temperatures) is twice the neutral temperature.

As this code is capable of handling a wide array of chemical and physical processes, and has been thoroughly benchmarked and validated for Venus and Mars, we use it for studying the ion escape rates of PCb. As a consequence, we assume that the atmospheric composition of PCb is similar to that of Venus and Mars. In reality, we wish to note there are many possibilities for the atmospheric composition of PCb, provided that it even exists (Ribas et al., 2016; Turbet et al., 2016; Goldblatt, 2016). A clearer picture regarding the existence and composition of the atmosphere may emerge via the JWST (Kreidberg & Loeb, 2016).

In the model, a nonuniform, spherical grid is adopted in order to accurately capture the multi-scale physics in different regions. The radial resolution varies from 5 km ($\sim 1/2$ scale height) at the inner boundary to several thousands of km at the outer boundary. The horizontal resolution is $3.0^{\circ}$ (in both longitude and latitude). The simulation domain is defined by $-$45 R${}_{P}\leq$ X $\leq$ 15 RP and $-$30 R${}_{P}\leq$ Y, Z $\leq$ 30 RP, where RP denotes the planetary radius. The code is run in the PCb-Star-Orbital (PSO) coordinate system, where the $x$ axis points from PCb toward Proxima Centauri, the $z$ axis is perpendicular to the PCb’s orbital plane, and the y axis completes the right-hand system.

Table 1 summarizes the chemical reaction schemes and the associated rates for inelastic collisions used in the MS-MHD calculations.

2.2 The selection of the simulation parameters

We adopt stellar and planetary parameters that are consistent with PCb. However, we wish to caution the readers that our exoplanet is not necessarily representative of PCb itself, as many of these parameters are partly or wholly unknown.

The Venusian neutral atmospheric profile is used as a prototype since PCb may experience significant losses of water and $H_{2}$, which are potentially major components, quite rapidly (Ribas et al., 2016; Barnes et al., 2016). As a result, the overall composition could eventually resemble that of Venus. It may also be relatively more suited to modeling the short distance ($\sim 0.05$ AU) between Proxima Centauri and PCb (Anglada-Escudé et al., 2016). But, we note that the Venusian atmosphere is possibly denser (and colder) compared to PCb. Both of these factors will alter the scale height, which in turn affects the escape rates. The Venusian neutral atmosphere from Ma et al. (2013) is:

[TABLE]

where $z_{0}=100$ km. We reconstruct the neutral atmospheric profiles based on the following parameters. First, the study by Anglada-Escudé et al. (2016) concluded that its minimum mass was $1.27\,M_{\oplus}$, which yields a minimum radius of around $1.1\,R_{\oplus}$ (Anglada-Escudé et al., 2016) assuming a rocky composition akin to that of the Earth, i.e., adopting the empirical relationship $R/R_{\oplus}\sim(M/M_{\oplus})^{1/3.7}$ (e.g. Valencia et al., 2006). If we assume that only the gravity changes, the scale height of PCb is

[TABLE]

and we will express all quantities in terms of the Venusian values detailed in Ma et al. (2013), since the code was calibrated for Venus. The second parameter required as an input is the density at the base of the model, which requires the surface atmospheric pressure, psur. A precise value of psur remains unknown at this stage - we choose the Earth value of $1$ bar. A similar value was also considered in the simulations by Meadows et al. (2016). Thus, we obtain a density that is $0.011$ times the Venus value at the model lower boundary, the reduced value primarily arising from the much lower atmospheric pressure. PCb’s neutral atmospheric profiles are as follows:

[TABLE]

It is noteworthy that the density of the dominant neutral component at 100 km altitude on the Earth (N2) and PCb (CO2) are about 1$\times$1013 cm*-3*, nearly equal to one another. Furthermore, the scale height of N2 at 100 km is about 5.8 km which is close to the scale height of CO2 for Venus at 100 km, before rescaling it to account for the higher gravity of PCb. Therefore, the reconstructed neutral atmosphere mimics certain features of Earth’s atmosphere. The sensitivity of our results to the atmospheric pressure is evaluated by also considering the case with Venusian surface pressure ($93$ bar).

Next, the stellar wind parameters at PCb must be specified, namely the density $\left(n_{sw}\right)$, temperature $\left(T_{sw}\right)$, velocity $\left(v_{sw}\right)$ and interplanetary magnetic field ($\mathrm{B_{IMF}}$). The dynamic pressure $P_{dyn}=m_{p}n_{sw}v_{sw}^{2}$ ($m_{p}$ is proton mass), magnetic pressure $P_{mag}=\mathrm{B_{IMF}^{2}}/(2\mu_{0})$ and total pressure $P_{tot}=P_{dyn}+P_{mag}$ are also introduced. The relevant values have been listed in Table 2, and have been chosen from Garraffo et al. (2016) with either maximum or minimum total pressure (and dynamic pressure), over one orbit of PCb. The stellar wind MHD model is driven at its inner boundary by the observable stellar surface magnetic field (i.e., magnetogram).

For the magnetized case, we have chosen PCb’s dipole moment to be approximately one-third of the Earth with the same dipole polarity. A similar value has also been adopted for PCb in recent studies (Zuluaga & Bustamante, 2016). As PCb rotates slower than the Earth (Ribas et al., 2016), this choice is consistent with certain dynamo scaling laws (Christensen, 2010).

3 The Simulation Results

The most important results of our simulations, i.e., the escape rates (computed by using a sphere with the radius equal to 10 RP) are summarized in Table 3. The contour plots of the O+ ion density, the magnetic field strength $\left(|\mathbf{B}|\right)$ and the magnetic field lines for unmagnetized case 1 (C1-UnM), magnetized case 1 (C1-M) and magnetized case 2 (C2-M) are presented in Fig. 1. As seen from Fig. 1, the case with a dipole field is characterized by a relatively large barrier to the stellar wind (from $\mid$B$\mid$ color contours). In all three instances, the plasma boundaries are greatly compressed, primarily by the stellar wind dynamic pressure, Pdyn, especially for the magnetized case. The second row in Fig. 1 illustrates the plasma boundaries clearly by zooming in.

The main conclusions are as follows: 1) Ions escape along the open magnetic field lines in all three cases (top row of Fig. 1), from the forces (mostly pressure gradient and $\mathbf{J}$ $\times$ $\mathbf{B}$ terms) in the ion momentum equation. 2) Although C1-M has a global dipole magnetic field, the dipole field strength is not strong enough to fully protect the exoplanet from the stellar wind interaction due to the enormous stellar wind dynamic pressure, and the southward stellar wind Bz that enables dayside magnetic reconnection. 3) Closed dipole field lines are observed for C2-M, and not C1-M, because the value of Pdyn for case 2 ($\sim$ 5 $\times$ 104 nP) is five times smaller than that of case 1 ($\sim$ 2.5 $\times$ 105 nP). 4) In both C1-M and C2-M, the polar cap regions (i.e., the regions poleward of the auroral oval, Schunk & Nagy 2009, pp. 26) are greatly extended to lower latitudes compared to that of the Earth. Thus, the escape rates in both cases are expected to be larger than the typical values observed in our Solar system.

Interestingly, the ionospheric profiles, especially the heavy ions O${}_{2}^{+}$ and CO${}_{2}^{+}$, are not significantly affected by the stellar wind conditions $\lesssim$ 200 km (i.e., the photochemical region) regardless of the stellar wind total pressure - see Fig. 2. The slight increase in O+ density (in both panels of Fig. 2) results from the charge exchange (see Table 1) with the enhanced penetrating stellar wind protons (not shown in the plot).

From Table 3, we see that the total escape rates for the unmagnetized cases are $\sim 10^{27}$/sec. As expected from Fig. 1, the ion escape rates for C1-UnM and C1-M are similar, although the dipole fields are responsible for reducing the total ion losses to some extent. However, the total ion escape rate for the magnetized case C2-M (compared to C2-UnM) is decreased by about one order of magnitude (mainly caused by decrease in O+) and is on the order of $\sim 10^{26}$/sec. This is consistent with the notion that a global dipole magnetic field protects Venus-like exoplanets from the stellar wind to some extent, although the final rates are still higher than those observed in our Solar system. But, this conclusion is not necessarily valid for exoplanets with different sizes, neutral atmosphere profiles and stellar wind parameters.

In Fig. 2, CO${}_{2}^{+}$ is one of the major ion species in the PCb’s ionosphere as a direct consequence of the photoionization of CO2. However, CO${}_{2}^{+}$ quickly reacts with neutral O to produce mostly O${}_{2}^{+}$ which becomes the other major ion at a similar altitude as CO${}_{2}^{+}$. Due to the large scale height of neutral O (because of the lighter mass), the O+ ionospheric peak is located at a relatively high altitude, and is thus picked up first by the stellar wind. Due to momentum conservation, if the stellar wind picks up a large number of light ions, it cannot pick up too many heavy ions. In addition, the electromagnetic shielding caused by the high-altitude ionized O+ ions (e.g, as a conductor layer) also helps in preventing the stellar wind and IMF from penetrating deeply into the planetary ionosphere and picking up the heavy ions. This helps explain why there are higher O+ but lower heavy ion losses in C1-UnM93, as compared to C1-UnM. C1-UnM93 is the unmagnetized case with the Venusian surface pressure of $93$ bar instead of the Earth value of $1$ bar with the scale height chosen in accordance with (2). Overall, the total ion escape rate does not change too much, indicating that the surface pressure does not significantly affect the ion escape rates under extreme stellar wind conditions.

In contrast to our preceding discussion, it is noteworthy that the terrestrial planets in our Solar system exhibit loss rates that are typically $\sim$1025/sec (Lammer, 2013; Brain et al., 2016).

4 Discussion and Conclusion

It was argued in Seki et al. (2001) that approximately $2\%$ of the atmospheric oxygen had been lost over 3 billion years on Earth. As the escape losses for Proxima Centauri b (PCb) in the unmagnetized case are about two orders of magnitude higher,999Since our choices of PCb’s radius and surface atmospheric pressure are close to that of the Earth, the total amount of atmosphere will also be similar to Earth. all of the atmosphere could be depleted much faster - possibly in a span of $\mathcal{O}\left(10^{8}\right)$ years. In turn, this has very important ramifications for surface-based life-as-we-know-it, given the importance of elements such as oxygen (Raymond & Segrè, 2006; Decker & Van Holde, 2010).

If gases such as oxygen are depleted on these short timescales, sufficient time may not exist for complex life to evolve. Our simulations indicate that the escape rates for PCb in the magnetized case are higher than that of the Earth, implying that some of the above conclusions for the unmagnetized case are also valid here. But, it is equally important to recognize that the magnetized case is quite sensitive to the values of the stellar wind parameters (see Table 3). The atmosphere depletion could occur over $\mathcal{O}\left(10^{8}\right)$ and $\mathcal{O}\left(10^{9}\right)$ years for C1-M and C2-M, respectively.

Thus, our results appear to indicate that PCb, and other similar M-dwarf exoplanets, are not generally capable of supporting an atmosphere over long (Gyr) timescales in both the unmagnetized or magnetized cases when the total stellar wind pressure is high. In contrast, for lower values of the pressure, the magnetized case can potentially sustain an atmosphere over Gyr timescales. It has been widely established that the HZ of M-dwarfs evolves over time (Shields et al., 2016). If the stellar wind pressure was higher at earlier epochs, this would have led to increased escape rates. Hence, even with a dipole field, the atmosphere may be eroded on sub-Gyr timescales.

We have also demonstrated that the ionospheric profiles (especially the heavy ions O${}_{2}^{+}$ and CO${}_{2}^{+}$) are mostly unaffected by the stellar wind conditions at $\lesssim$ 200 km. This is important because the stellar wind pressure is highly variable, as noted in Garraffo et al. (2016). Hence, the fact that the stellar wind variability does not greatly influence the lower regions lends some credence to the hypothesis that the surface biology may not be significantly affected (Grießmeier et al., 2016).

There are some important caveats, with regards to the above conclusions, that must be reiterated at this stage. Most of the parameters for PCb and Proxima Centauri are currently unknown. The composition and thickness of the atmosphere could be different, enabling the retention of the atmosphere over longer periods of time. The existence of processes, such as outgassing which occurs on Titan (Tobie et al., 2006), may counteract the atmospheric losses by serving as sources. Thus, a complete understanding of this problem necessitates a knowledge of both source and loss mechanisms. Our model does not include kinetic processes, which contribute to the ion loss mechanisms (Strangeway et al., 2005).

Finally, as observed in Sec. 1, there are $\sim 10^{10}$ exoplanets in the HZ around M-dwarfs. Thus, even if our parameter choices are invalidated by future observations of PCb, these values would, in all likelihood, be valid for other M-dwarf HZ exoplanets whose parameters are similar to the ones considered herein.

To summarize, this Letter sheds light on atmospheric losses from M-dwarf exoplanets in the HZ such as PCb, and clearly delineates the role of the planetary magnetic field and the stellar wind parameters. We have shown that Venus-like exoplanets (with PCb-like parameters) are characterized by high escape rates in the unmagnetized limit, but these values are reduced to an extent when a dipole magnetic field exists. Hence, PCb may undergo significant atmospheric erosion over Gyr timescales in both the magnetized and unmagnetized cases, but this statement must be taken with due caution since there are many inherent uncertainties. Future missions such as JWST will be essential in placing constraints on stellar winds and exoplanet atmospheres, thereby paving the way for more accurate estimations of atmospheric losses.

The authors thank Adam Burrows, Zhenguang Huang, Abraham Loeb, Janet Luhmann and Edwin Turner for the helpful and encouraging comments. This research was supported by the NASA Living With a Star Jack Eddy Postdoctoral Fellowship Program, administered by the University Corporation for Atmospheric Research. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. The Space Weather Modeling Framework that contains the BATS-R-US code used in this study is publicly available from http://csem.engin.umich.edu/tools/swmf. For distribution of the model results used in this study, please contact the corresponding author.