Probabilistic Constraints on the Mass and Composition of Proxima b

It can be shown (e.g., Ho & Turner 2011) that the probability distribution of $\sin (i)$ corresponding to an isotropic inclination distribution is

$$\begin{eqnarray}&&P(\sin (i))=\sin (i)/\sqrt{1-{\sin }^{2}(i)}.\end{eqnarray} \tag{ 1 }$$

Since this distribution peaks at $\sin (i)=1$, the mass distribution of an RV-detected planet—assuming no prior constraints on the mass—peaks at the minimum mass M₀.

In their models of the possible orbital histories of Proxima b, Barnes et al. (2016) find that galactic tides could have inflated the eccentricity of the host star's (at the time unconfirmed) orbit around the α Cen binary, leading to encounters within a few hundred au and the possible disruption of Proxima's planetary system. If so, this could affect the likely inclination of the planet in a non-isotropic way. However, Kervella et al. (2017) have presented RV measurements showing that Proxima is gravitationally bound to the α Cen system with an orbital period of 550,000 years, an eccentricity of ∼0.5, and a periapsis distance of 4200 au. At this distance, the ratio of Proxima's gravitational field to that of α Cen at the planet's orbit (∼0.05 au) is greater than 10⁸; unless Proxima's orbit was significantly more eccentric in the past, it seems unlikely that α Cen would have disrupted the system.

Mulders et al. (2015) provide up-to-date occurrence rates of planets around M dwarf stars from the Kepler mission. The sample is limited to 2 < P < 50 days, over which they find the occurrence rates to be mostly independent of the period. The binned rates and a regression curve, as well as their uncertainties, are presented in Figure 1.

Zoom In Zoom Out Reset image size

Kepler statistics for M dwarfs remain incomplete below 1 R_⊕, but complete statistics for earlier-type stars suggest a flat distribution for 0.7 < R < 1 R_⊕ (Mulders et al. 2015). Since mass–radius relationships typically find a strong dependence of mass on radius (M ∝ R^3–4; e.g., Weiss & Marcy 2014; Wolfgang et al. 2016), we assume a priori that Proxima b (M ≳ 1.3 M_⊕) is larger than 0.7 R_⊕. Therefore, for this Letter we adopt the regression curve fitted to the binned data, but set the occurrence rates to be flat for R < 1 R_⊕.

Multiple works (e.g., Weiss & Marcy 2014; Rogers 2015) have determined the existence of two distinct populations of exoplanets smaller than Neptune (R ≲ 4 R_⊕): a small radius population with densities consistent with an entirely iron and silicate composition (hereafter "rocky"), and a large radius population with lower-density planets that must have significant amounts of ice or a thick H/He atmosphere (hereafter "sub-Neptunes").

Rogers (2015) studies the abundance of planets of each composition as a function of radius. They define f_α(R) as the likelihood that a planet of radius R is dense enough to be consistent with a rocky composition and determine f_α for a sample of planets with known masses and radii. They suggest fitting the data with a two-parameter linear model:

$$\begin{eqnarray}&&{f}_{\alpha }({R}_{{\rm{P}}},{R}_{\mathrm{thresh}},{\rm{\Delta }}R)\\ =\,\left\{\begin{array}{ll}1 & {R}_{{\rm{P}}}\lt {R}_{\mathrm{thresh}}-\tfrac{1}{2}{\rm{\Delta }}R\\ 0.5+\displaystyle \frac{{R}_{\mathrm{thresh}}-{R}_{{\rm{P}}}}{{\rm{\Delta }}R} & | {R}_{{\rm{P}}}-{R}_{\mathrm{thresh}}| \lt \tfrac{1}{2}{\rm{\Delta }}R\\ 0 & {R}_{{\rm{P}}}\gt {R}_{\mathrm{thresh}}+\tfrac{1}{2}{\rm{\Delta }}R.\end{array}\right.\end{eqnarray} \tag{ 2 }$$

They find a step function to best describe the data, with ΔR fixed at zero and R_thresh ≈ 1.5 R_⊕. For the purposes of this Letter, we prefer this fit, but will also vary R_thresh and ΔR to see how they affect our results.

We stress that a planet for which f_α = 1 is only sufficiently dense to be rocky; we still cannot necessarily exclude an ice or volatile component. Here, we will assume that all planets for which f_α = 1 follow the low-radius M–R relationships given in the following section, which were empirically fitted without prior knowledge of the planets' compositions. For simplicity, we refer to these as "rocky" planets, and the other population as "sub-Neptunes," but we will revisit this distinction later on.

Since Proxima b is in the habitable zone, it receives an amount of stellar flux comparable to that received by Earth, so we should bear in mind the possibility that the volatile envelope of a sub-Neptune could be lost due to photoevaporation. Owen & Mohanty (2016) model rocky planets with thick H/He envelopes in the habitable zones of M dwarfs, finding that planets with M > 0.8 M_⊕ maintain their envelopes over Gyr timescales and are therefore uninhabitable. The 2σ lower limit on the minimum mass of Proxima b is 0.93 M_⊕, so it is unlikely that any H/He envelope on the planet would evaporate under this rule. However, we note that this study focuses on planets with a primarily rocky composition, so it may not be directly applicable to habitable zone sub-Neptunes.

Additionally, Zahnle & Catling (2013) empirically define boundaries for atmospheric evaporation as a function of stellar heating, escape velocity, and atmospheric composition. In particular, a planet receiving an Earth-like flux must have an escape velocity above ∼8 km s⁻¹ in order to maintain an H₂ atmosphere for 5 Gyr. We will revisit this requirement in Section 4.2.

Empirically derived relationships between exoplanet masses and radii rely on RV or transit-timing variation (TTV) measurements of transiting exoplanet masses. Weiss & Marcy (2014) fit a mass–radius (hereafter M–R) relationship to a sample of 65 transiting exoplanets, in which they find evidence for the two populations discussed in Section 2.3. Through least-squares regression, they find the densities of the rocky planets to increase linearly with planet radius:

$$\begin{eqnarray}&&{\rho }_{{\rm{P}}}=2.43+3.39\left(\displaystyle \frac{{R}_{{\rm{P}}}}{{R}_{E}}\right){\rm{g}}\,{\mathrm{cm}}^{-3}\end{eqnarray} \tag{ 3 }$$

while the RV-measured masses of sub-Neptunes increase nearly linearly with planet radius:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{{\rm{P}}}}{{M}_{\oplus }}=4.87{\left(\displaystyle \frac{{R}_{{\rm{P}}}}{{R}_{E}}\right)}^{0.63}.\end{eqnarray} \tag{ 4 }$$

Wolfgang et al. (2016) use an expanded version of this data set to fit power-law M–R relationships using a more statistically robust Bayesian method. For the rocky planets, they find

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{{\rm{P}}}}{{M}_{\oplus }}=1.4{\left(\displaystyle \frac{{R}_{{\rm{P}}}}{{R}_{E}}\right)}^{2.3}\end{eqnarray} \tag{ 5 }$$

and for the sub-Neptunes with RV-measured masses,

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{{\rm{P}}}}{{M}_{\oplus }}=2.7{\left(\displaystyle \frac{{R}_{{\rm{P}}}}{{R}_{E}}\right)}^{1.3}.\end{eqnarray} \tag{ 6 }$$

Due to the larger sample size and more robust fitting procedure, we adopt Equations (5) and (6) as our preferred M–R relationships, but for completeness we consider the Weiss & Marcy (2014) relationships as well. We find that the choice of M–R relationships has a minimal impact on our final results. Both sets of relationships are plotted in Figure 1.

It is important to note that the above relationships for sub-Neptunes exclude masses measured by TTV, since TTV masses have been found to be systematically lower than RV masses. This could indicate a selection bias or systematic error in the method used, but since Proxima b was detected through RV measurements, we believe it is proper to exclude the TTV masses either way.

It is also clear that there is a significant spread in the masses of the observed planets. Wolfgang et al. (2016) suggest a spread of ±1.9 M_⊕ for the sub-Neptune planets, which we adopt for our simulations. For rocky planets, the spread is noticeably smaller. There are too few planets to constrain this spread, but it should most likely increase with mass, so we arbitrarily define the spread to be 30% of the calculated mass.

The fitted occurrence rates and their uncertainties (f ± df) are given in even bins in log-space. We use them to generate a random sample of radii, where the number of radii in each bin (r₀) is selected from a normal distribution with mean value f(r₀) and standard deviation df (r₀). We find that the results converge for 1000 samples of the occurrence rates, with each sample containing ∼10⁶ radii.

To each radius, we assign a composition ("rocky" or "sub-Neptune") based on the model of Rogers (2015; Equation (2)), with R_thresh = 1.5 R_⊕ and ΔR = 0. We then assign a mass to each radius and composition from a Gaussian distribution with mean value M(R)—calculated using our chosen M–R relationships (Equations (5) and (6))—and a standard deviation dM that represents the spread. We choose a spread proportional to the calculated mass for rocky planets (dM = 0.3 × M(R)), but a constant spread for sub-Neptunes (dM = 1.9 M_⊕). We also reject negative masses, which could in principle bias the assigned masses toward higher-than-average values—however, we find that only a negligible number of masses are rejected.

Finally, we assign a line-of-sight inclination parameter $\sin (i)$ to each planet, drawn from the isotropic probability distribution discussed in Section 2.1.

The prior mass and radius distributions, P(M) and P(R), can be derived directly from the simulated sample. Factoring in the projected minimum mass M₀, the posterior distributions $P(M| {M}_{0})$ and $P(R| {M}_{0})$ can be calculated from Bayes's formula:

$$\begin{eqnarray}&&P(X| {M}_{0})=\displaystyle \frac{P({M}_{0}| X)P(X)}{P({M}_{0})},\end{eqnarray} \tag{ 7 }$$

where X represents mass or radius. Since M₀ is known, P(M₀) is just a normalizing constant. Taking ${M}_{0}={1.27}_{-0.17}^{+0.20}$ M_⊕ as the projected mass of Proxima b and the upper limit ${\sigma }_{{M}_{0}}=0.20$ M_⊕ as its standard deviation, we calculate for each simulated planet:

$$\begin{eqnarray}&&{P}_{i}({M}_{0}| X)=\exp (-{({M}_{0}-{M}_{i}{\sin }_{i}(i))}^{2}/2{\sigma }_{{M}_{0}}^{2}).\end{eqnarray} \tag{ 8 }$$

Then $P({M}_{0}| M)$ and $P({M}_{0}| R)$ are the average values of ${P}_{i}({M}_{0}| X)$ for each bin in mass or radius. The prior and posterior distributions are calculated for each sample of 10⁶ planets, and the final results are taken to be the mean result of 1000 samples.

The prior probability that a planet in a given mass bin is rocky is equal to the number of simulated rocky planets in that bin divided by the total number of planets in the same bin. Since we want to know the likelihood that Proxima b belongs to the "rocky" population, we multiply this prior composition probability distribution by the posterior mass distribution from the previous section and integrate over all masses.

The prior and posterior mass probability distributions for Proxima b are plotted in Figure 2. The shaded regions demonstrate the relative contributions of the populations at each mass. The prior distribution is valid for RV-detected planets around M dwarfs with intermediate periods (2 < P < 50 days) and radii (0.7 < R < 4 R_⊕), while the posterior distribution can be taken as the mass probability distribution for Proxima b.

Zoom In Zoom Out Reset image size

For reference, we include the posterior distribution given no prior constraints on the mass; that is, the distribution resulting from an isotropic $\sin (i)$ distribution and the measured M₀ with its uncertainty. We find that this nearly matches our result, since both $P(\sin (i))$ and P(M) are bottom heavy.

Figure 3 shows the cumulative probability that M < X for both of the considered M–R relationships (Section 2.4) as well as for the case of no prior mass distribution. We find that there is little difference between the results for each M–R relationship.

Zoom In Zoom Out Reset image size

In order to verify that sub-Neptune planets can maintain H₂ envelopes in the habitable zone, we compare the escape velocities of our simulated sub-Neptunes to the ∼8 km s⁻¹ cutoff for H₂ atmospheric escape (assuming an Earth-like stellar flux) defined by Zahnle & Catling (2013). In both the prior and posterior distributions of escape velocities, we find that fewer than 1% of the sub-Neptunes have escape velocities below this threshold, with most having v_e ≳ 15 km s⁻¹. Therefore, we do not believe that Proxima b will be subject to significant atmospheric loss if it has a sub-Neptune composition.

Table 1 lists the sets of parameters for which we run the simulation, including the mass spread dM for each composition and the central value (R_thresh) and width (ΔR, if nonzero) of the transition region defined by Equation (2). The following results for each case are given: the probability P_rocky that Proxima b belongs to the "rocky" category of planets, i.e., that its density is consistent with a fully iron and silicate composition, and the expectation values ${\langle M\rangle }_{\mathrm{rocky}}$ and ${\langle R\rangle }_{\mathrm{rocky}}$ of the mass and radius under the assumption that it belongs to this population.

Table 1.  Monte Carlo Simulation Parameters and Results

Case	Parameters				Results
	dM(rocky)	dM(sub-Neptune)	Rthresh	ΔR	Procky	${\langle M\rangle }_{\mathrm{rocky}}$	${\langle R\rangle }_{\mathrm{rocky}}$
Case A	0.3 × M	1.9 M⊕	1.5 R⊕	...	89.9%	${1.63}_{-0.72}^{+1.66}\,{M}_{\oplus }$	${1.07}_{-0.31}^{+0.38}\,{R}_{\oplus }$
Case B	0.6 × M	3.8 M⊕	1.5 R⊕	...	84.6%	${1.65}_{-0.73}^{+1.95}$ M⊕	${1.03}_{-0.36}^{+0.42}$ R⊕
Case C	0.1 × M	0.7 M⊕	1.5 R⊕	...	93.6%	${1.65}_{-0.73}^{+1.52}$ M⊕	${1.06}_{-0.24}^{+0.36}$ R⊕
Case D	0.3 × M	1.9 M⊕	1.7 R⊕	...	94.6%	${1.71}_{-0.79}^{+2.13}$ M⊕	${1.10}_{-0.33}^{+0.50}$ R⊕
Case E	0.3 × M	1.9 M⊕	1.3 R⊕	...	81.6%	${1.52}_{-0.62}^{+1.15}$ M⊕	${1.02}_{-0.27}^{+0.26}$ R⊕
Case F	0.3 × M	1.9 M⊕	1.5 R⊕	1.2 R⊕	84.8%	${1.64}_{-0.73}^{+1.99}$ M⊕	${1.06}_{-0.30}^{+0.53}$ R⊕
Case G	0.6 × M	3.8 M⊕	1.5 R⊕	1.2 R⊕	81.1%	${1.65}_{-0.75}^{+2.13}$ M⊕	${1.02}_{-0.36}^{+0.63}$ R⊕

Note. The resulting values of P_rocky, ${\langle M\rangle }_{\mathrm{rocky}}$, and ${\langle R\rangle }_{\mathrm{rocky}}$ for different mass spreads dM and compositional parameters R_thresh and ΔR. The expectation values are reported with 95% confidence intervals.

Download table as:  ASCIITypeset image

Case A is most consistent with the previous work we have cited, so we take it as our primary result. In this case, there is a ∼90% probability that Proxima b belongs to the "rocky" population, with an ∼10% likelihood that it belongs to the "sub-Neptune" population. In the case that it is rocky, the expectation values (and 95% confidence intervals) for the mass and radius are ${\langle M\rangle }_{\mathrm{rocky}}={1.63}_{-0.72}^{+1.66}$ M_⊕ and ${\langle R\rangle }_{\mathrm{rocky}}={1.07}_{-0.31}^{+0.38}$ R_⊕.

We investigate the effect of increasing (Case B) and decreasing (Case C) the mass spread for each composition, which results in lower and higher values of P_rocky, respectively. This results from low-radius (R ∼ 1.5 R_⊕), low-mass sub-Neptunes; when dM is large, they can lie significantly below the M–R relation with masses between 1 and 2 M_⊕, so that they are indistinguishable from the rocky planets in the mass domain.

In Cases D and E, we determine the effect of raising or lowering the threshold radius R_thresh at which the rocky and sub-Neptune populations are split. A 0.2 R_⊕ offset in either direction, which encompasses most of the values suggested in the literature, results in a ∼5% to 8% shift in P_rocky, where higher threshold radii allow for more rocky planets and therefore a higher probability of a rocky composition. Furthermore, allowing for a nonzero width ΔR to the cutoff region allows sub-Neptunes to exist with lower radii and masses, thereby decreasing P_rocky.

In all cases, we find P_rocky to be between 80% and 95% using the Wolfgang et al. (2016) M–R relationship, and we find similar values using the Weiss & Marcy (2014) relationship (e.g., P_rocky = 90.7% for Case A), so this result does not vary substantially over the range of reasonable values for the input parameters. This is consistent with the results of Kane et al. (2017), who determine that Proxima b is 84% likely to be terrestrial assuming an isotropic probability distribution for its inclination and a simple mass boundary between the two compositions at 5.1 M_⊕.