A heavy molecular weight atmosphere for the super-Earth $\pi$ Men c

We improved the published transit ephemeris (Gandolfi et al., 2018) using TESS data from Sectors 1, 4, 8, 11–13 and the code pyaneti (Barragán et al., 2019), which allows for parameter estimation from posterior distributions.

Figure 4 shows the spectrum obtained during the second HST observation (top), marking the stellar features at which we looked for absorption in the $-$70 to $+$10 km/s velocity range, and the resulting signal-to-noise ratio per spectral bin (middle). The bottom panel compares the in- and out-of-transit spectra.

Having detected absorption of the C ii 1335 Å feature, we looked for a similar signal at the C ii 1334 Å line, without success. Next, we show that the absorption signal at 1334 Å is hidden by ISM contamination (Figure 4, bottom), which affects the S/N. We first fitted the C ii 1335 Å stellar feature using a Gaussian profile and the out-of-transit spectrum (black dashed line) and then employed a further Gaussian profile to fit the position and strength of the planetary absorption on the in-transit spectrum (yellow solid line). We obtained that the Gaussian profile simulating the planetary absorption lies at a velocity of about $-$48.5 km/s with respect to the position of the main C ii 1335 Å feature, has a normalized amplitude of $\sim$0.18, and a width of $\sim$0.11 Å. These fits were performed considering the unbinned spectra. We then derived the strength of the C ii 1334 Å line, prior to ISM absorption by scaling the Gaussian fit to the C ii 1335 Å feature by the ratio of the oscillator strength times the statistical weight of the two lines (about 1.8; blue dashed line). We further simulated the C ii ISM absorption profile at the position of the C ii 1334 Å resonance line employing a Voigt profile (purple solid line), in which we set the position of the line equal to that obtained from the reconstruction of the stellar Ly$\alpha$ line (García Muñoz et al., 2020) and a C ii ISM column density equal to that of hydrogen scaled to the expected ISM carbon-to-hydrogen abundance ratio and ISM C ionization fraction (Frisch & Slavin, 2003) (green solid line). Figure 4 (bottom) indicates that the simulated profile is a good match to the out-of-transit spectrum, particularly accounting for the uncertainties involved in inferring the C ii ISM column density. Finally, we added to this line the planetary absorption obtained from fitting the C ii 1335 Å feature, rescaled by 1.8 (orange solid line). Lastly, Figure 4 (bottom) shows that the difference between the simulated C ii 1334 Å line profiles before and after adding the planetary absorption is significantly smaller than the observational uncertainties and hence undetectable in the data. In summary, reduced S/N due to ISM contamination hides the planetary absorption signal at 1334 Å.

To estimate the likelihood that our C ii signal is due to intrinsic stellar variability, we did a Monte Carlo simulation where we assumed the population of intrinsic stellar variability between the mean in-transit and mean out-of-transit spectrum is represented by the measured in-transit absorption from each line (Table 1). We excluded from this representative population the C ii lines because of the putative planetary absorption and the O i triplet because the noise properties and intrinsic variability of these airglow-contaminated lines are different. For the remaining eight emission lines, we calculate a weighted mean of 2.67$\pm$1.27% for the in-transit absorption. We made 10⁶ realizations of nine emission lines (the eight ’unaffected’ lines of the representative population plus C ii 1335.7 Å, the line with the 3.4-sigma in-transit absorption) each with in-transit absorption randomly drawn from a normal distribution with mean 2.67% and standard deviation 1.27%. We find that the likelihood of 1 or more emission lines having an absorption $\geq$6.76% due to intrinsic stellar variability as measured by our COS spectra is 2%.

We inspected the O i lines for evidence of planetary absorption. The O i triplet at 1302–1306 Å is comprised of three emission lines that share an upper energy level. O i 1302.168 Å is the resonance line, and is significantly affected by ISM absorption. O i 1304.858 and 1306.029 Å each have similar oscillator strengths to the resonance line. COS’s wide aperture leads to significant contamination of each line of the O i triplet by geocoronal airglow. We corrected each orbit’s pipeline-reduced x1d spectrum’s O i emission lines for geocoronal airglow by performing a least-squares fit of the airglow templates downloaded from MAST (Bourrier et al., 2018a) to the spectrum (Figure 5, top). We created mean in-transit and out-of-transit spectra in the same way as the C ii lines, and measured absorptions in each of the three lines (shortest wavelength to longest wavelength) of $-$3.65$\pm$40.65 %, $-$8.33$\pm$6.32%, and 5.17$\pm$3.66% over the same velocity range as C ii (Figure 5, bottom). No significant variation was observed in the O i triplet. The O i line at 1355 Å similarly does not show any significant variation between in-transit and out-of-transit. Lastly, we produced time series for the flux of some stellar lines (Figure 6). The series reveals no obvious variability that could cause a false transit detection.

Our nominal transit depth for C ii of 6.8% translates into a projected area equivalent to a disk of radius $R_{\rm{C\;II}}$/$R_{\rm{p}}$=15.1. This is about the extent of $\pi$ Men c’s Roche lobe in the substellar direction ($R_{\rm{L_{1}}}$/$R_{\rm{p}}$=13.3). The spectroscopic velocity of $+$10 km/s for the C ii ions is consistent with our photochemical-hydrodynamic predictions, and likely traces absorption in the vicinity of the planet’s dayside as the gas accelerates toward the star. Comparable velocities are predicted by 3D models on the planet’s dayside (Shaikhislamov et al., 2020b) before the escaping gas interacts with the impinging stellar wind. Negative velocities of $-$70 km/s (and probably faster, as the stellar line becomes weak and the S/N poor at the corresponding wavelengths) suggest that the C ii ions are accelerated away from the star by other mechanisms such as tidal forces and magnetohydrodynamic interactions with the stellar wind. For reference, the velocity of the solar wind at the distance of $\pi$ Men c is on the order of 250 km/s (Shaikhislamov et al., 2020b). There is observational evidence for HD 209458 b, GJ 436 b and GJ 3470 b that their escaping atmospheres also result in preferential blue absorption (particularly the latter two) (Vidal-Madjar et al., 2003; Ehrenreich et al., 2015; Bourrier et al., 2018b). Models considering the 3D geometry of the interacting stellar and planet winds also favor blue absorption, especially when the stellar wind is stronger (Shaikhislamov et al., 2020a, b) and arranges the escaping atoms into a tail trailing the planet. We consider a C ii tail to be realistic scenario for $\pi$ Men c.

To gain intuition, we have built the phenomenological model of $\pi$ Men c’s tail sketched in Figure 7 (top). The C ii ions are injected into a cylindrical tail of fixed radius 15$R_{\rm{p}}$ and are subject to a prescribed velocity $U$=$-$($U_{0}$ $-$ $x$/$t_{\rm{acc}}$) km/s that varies in the tail direction. This geometry surely simplifies the true morphology of the escaping gas, which is likely to resemble an opening cone tilted with respect to the star-planet line (Shaikhislamov et al., 2020b). This crude description allows us at the very least to obtain analytical expressions for some of the relevant diagnostics. Here, $U_{0}$=10 km/s (a typical value from the photochemical-hydrodynamic simulations; Figure 2) and $t_{\rm{acc}}$ is an acceleration time scale. (Note: $x$$<$0 in the tail, and the ions are permanently accelerating.) Related accelerations have been predicted by physically-motivated 3D models (Ehrenreich et al., 2015; Shaikhislamov et al., 2020a, b), under the combined effect of gravitational, inertial and radiative forces. Magnetic interactions with the stellar wind might additionally affect ion acceleration.

The C ii ions photoionize further into C iii with a time scale $t_{\rm{C\;II}}$$\sim$20 hours (calculated for unattenuated radiation from $\pi$ Men at wavelengths shorter than the 508-Å threshold and the corresponding C ii cross section (Verner & Yakovlev, 1995)). The collision of stellar wind particles with the planetary C ii ions might ionize further the latter (ionization potential of 24 eV), especially at the mixing layer between the two flows (Tremblin & Chiang, 2013), but it remains unclear whether collisional ionization can compete on the full-tail scale with photoionization. This should be assessed in future work. The continuity equations that govern the model are:

As the flow accelerates through the tail, the total density decays and the C ii ions are converted into C iii. The solution for the C ii ion number density is:

This highly simplified exercise aims to estimate reasonable values for the free parameters [C ii]₀ and $t_{\rm{acc}}$ that reproduce the transit depth and Doppler velocities from the COS data. To produce the wavelength-dependent opacity, we integrate [C ii]($x$) along the line of sight keeping track of the Doppler shifts introduced by the varying $U$($x$). We represent the C ii cross section at rest for the 1335.7 Å line by a Voigt lineshape with thermal ($T$=6000 K, a reference temperature in the tail; our findings do not depend sensitively on this temperature, as the absorption over a broad range of velocities is caused by the bulk velocity of the ions rather than by their thermal broadening) and natural (Einstein coefficient $A_{\rm{ul}}$=2.88$\times$10⁸ s^-1) broadening components, and a fractional population of the substate from which the line arises based on statistical weights and equal to 2/3.

Figure 7 (bottom) shows how [C ii]₀ and $t_{\rm{acc}}$/$t_{\rm{C\;II}}$ affect the absorption of the stellar line. Absorptions consistent with the measurements are generally found for [C ii]₀$\geq$10⁴ cm^-3 and $t_{\rm{acc}}$/$t_{\rm{C\;II}}$$\leq$1, and result in C ii tails longer than $\sim$50$R_{\rm{p}}$ and large amounts of C ii moving at $-$70 km/s. A key outcome of the model is that for $t_{\rm{acc}}$/$t_{\rm{C\;II}}$$>$1 the C ii ion photoionizes too quickly to produce significant absorption at the faster velocities. This could in principle be compensated with increasing amounts of C ii ions entering the tail, but the photochemical-hydrodynamic model indicates that it is not possible to go beyond [C ii]₀$\sim$10⁶ cm^-3 without obtaining unrealistically large escape rates. For [C ii]₀$\sim$10⁴$-$10⁶ cm^-3, the mass fluxes of carbon atoms through the tail range from $\sim$10⁸ to $\sim$10¹⁰ g s^-1. These mass fluxes are comparable to the loss rates of carbon nuclei predicted by our photochemical-hydrodynamic model for atmospheres that have carbon abundances larger than solar (Table 2). This sets a weak constraint on the carbon abundance that prevents us for the time being from assessing whether carbon is a major constituent of $\pi$ Men c’s atmosphere.

This insight into $\pi$ Men c’s carbon abundance, even though of limited diagnostic value, is consistent with the interpretation of in-transit absorption by C ii at the hot Jupiter HD 209458 b (Vidal-Madjar et al., 2004; Ben-Jaffel & Sona Hosseini, 2010; Linsky et al., 2010). Both HD 209458 b and $\pi$ Men c transit Sun-like stars and exhibit similar transit depths. 3D models of HD 209458 b (Shaikhislamov et al., 2020a) show that its C ii absorption signal is explained by carbon in solar abundance. Because the predicted mass loss rate of HD 209458 b is an order of magnitude higher than for $\pi$ Men c (a result mainly from its lower density), a reasonable guess is that a $\times$10 solar enrichment for $\pi$ Men c will result in comparable transit depths. Future 3D modeling of the interaction of the C ii ions with the stellar wind and the planet’s magnetic field lines should help refine these conclusions.

We investigated $\pi$ Men c’s extended atmosphere with a 1D photochemical-hydrodynamic model that solves the gas equations at pressures $\leq$1 $\mu$bar (García Muñoz et al., 2020). Heating occurs from absorption of stellar photons by the atmospheric neutrals. Cooling is parameterized as described in our previous work, and includes emission from H${}_{3}^{+}$ in the infrared, Ly$\alpha$ in the FUV, atomic oxygen at 63 and 147 $\mu$m, and rotational cooling from H₂O, OH and CO also in the infrared. We adopted a NLTE formulation of H${}_{3}^{+}$ cooling that captures the departure of the ion’s population from equilibrium at low H₂ densities. We included a parameterization of CO₂ cooling at 15 $\mu$m (Gordiets et al., 1982). The model considers 26 chemicals (H₂, H, He, CO₂, CO, C, H₂O, OH, O, O₂, H${}_{2}^{+}$, H⁺, H${}_{3}^{+}$, He⁺, HeH⁺, CO${}_{2}^{+}$, CO⁺, C⁺, HCO${}_{2}^{+}$, HCO⁺, H₂O⁺, H₃O⁺, OH⁺, O⁺, O${}_{2}^{+}$ and electrons) that participate in 154 chemical processes. It does not include hydrocarbon chemistry, although this omission is not important as hydrocarbons are rapidly lost at low pressures in favor of other carbon-bearing species (Moses et al., 2013).

To explore a broad range of compositions, we adopted nuclei fractions $Y_{\rm{C}^{*}}$ and $Y_{\rm{O}^{*}}$ such that $Y_{\rm{C}^{*}}$+$Y_{\rm{O}^{*}}$ goes from a few times 10^-6 to $\sim$0.4, and $Y_{\rm{C}^{*}}$/$Y_{\rm{O}^{*}}$ ratios from 0.1 to 10. We imposed $Y_{\rm{He}^{*}}$=0.1$Y_{\rm{H}^{*}}$. By definition $Y_{\rm{H}^{*}}$+$Y_{\rm{He}^{*}}$+$Y_{\rm{C}^{*}}$+$Y_{\rm{O}^{*}}$=1, and thus the composition is specified by only two nuclei fractions. With the above information, we estimated the corresponding vmrs from Figure 7 of Hu & Seager (2014) and assigned them as boundary conditions of our extended atmosphere model. For the other gases, we adopted zero vmrs. Hydrocarbons (e.g. C₂H₂, CH₄) become abundant in the bulk atmosphere for $Y_{\rm{C}^{*}}$/$Y_{\rm{O}^{*}}$$\geq$2. Because our model does not currently include hydrocarbons, we transferred all the C nuclei at the base of the extended atmosphere from C₂H₂ and CH₄ into CO and C.

Table 2 summarizes 30 cases, each for a different bulk atmospheric composition. For cases 05-06 and 11-12 we first assumed for their vmrs in the bulk atmosphere:

• •

  (05) H₂: 7.66$\times$10^-1; He: 1.66$\times$10^-1; CO: 6.90$\times$10^-3; H₂O: 6.21$\times$10^-2.

• •

  (06) He: 1.46$\times$10^-1; CO₂: 9.76$\times$10^-2; H₂O: 7.32$\times$10^-1; O₂: 2.44$\times$10^-2.

• •

  (11) H₂: 7.88$\times$10^-1; He: 1.63$\times$10^-1; CO: 2.48$\times$10^-2; H₂O: 2.48$\times$10^-2.

• •

  (12) H₂: 4.47$\times$10^-1; He: 1.18$\times$10^-1; CO₂: 1.45$\times$10^-1; CO: 1.45$\times$10^-1; H₂O: 1.45$\times$10^-1.

However, our numerical experiments showed that H₂O (but also H₂ and O₂) were unstable for these cases, and their number densities dropped by orders of magnitude in a few elements of the spatial grid. This is evidence that these molecules chemically transform before reaching the $\mu$bar pressure level. To avoid numerical difficulties, in these four cases we prescribed the vmrs at the $\mu$bar pressure level by replacing the unstable molecules by their atomic constituents while preserving the original nuclei fractions. They are indicated with a $\dagger$ in Table 2, which lists the adopted vmrs.

The mass fraction of e.g. H₂ and He is given by the ratio vmr(H₂)$\;\mu_{\rm{H}_{2}}$/($\sum_{s}$vmr_s$\;\mu_{s}$) and vmr(He)$\;\mu_{\rm{He}}$/($\sum_{s}$vmr_s$\;\mu_{s}$), with the summation extending over all species, respectively. The mass fraction of heavy volatiles in the bulk atmosphere $Z$ is calculated as one minus the added mass fractions of H₂ and He.

The right hand side of Table 2 summarizes the model outputs. $\dot{m}_{\rm{H}}$($@$Mach=1) quotes the mass flux of neutral H i atoms at the sonic point. It departs from the mass loss rate of H nuclei ($\dot{m}_{\rm{H}^{*}}$) if at the sonic point hydrogen is ionized. $\dot{m}_{\rm{He}^{*}}$, $\dot{m}_{\rm{C}^{*}}$, $\dot{m}_{\rm{O}^{*}}$ and $\dot{m}$ quote the mass loss rates of the specified nuclei and the net mass loss rate of the atmosphere. All mass loss rates are calculated over a 4$\pi$ solid angle.

Table 2 shows that $\dot{m}_{\rm{H}}$($@$Mach=1) is about 5$\times$10⁹ g s^-1 for a H₂/He atmosphere. 3D simulations of $\pi$ Men c’s extended atmosphere (Shaikhislamov et al., 2020b) show that the Ly$\alpha$ absorption signal varies monotonically with the mass flux of protons in the stellar wind. According to their numerical experiments (their Figure 4), increases in the stellar wind flux by a factor of a few result in deeper transits by a similar factor of a few. Because ENA generation depends on both the flux of protons in the stellar wind and the flux of neutral hydrogen in the planet wind, we estimate that a factor of a few (we take $\times$1/4) decrease in $\dot{m}_{\rm{H}}$($@$Mach=1) with respect to the case of an H₂/He atmosphere suffices to bring ENA generation to undetectable levels for HST/STIS. We take $\dot{m}_{\rm{H}}$($@$Mach=1)=1.25$\times$10⁹ g s^-1 as our approximate threshold for non-detection of Ly$\alpha$ absorption. A refined estimate of this threshold calls for 3D simulations over a variety of bulk atmospheric compositions.

As the adopted intrinsic temperature affects the predicted atmospheric mass (because of its impact on the scale height), we constrain its possible values through thermal evolution calculations. For an adiabatic planetary envelope that cools efficiently by convection under the moderating effects of atmospheric opacity, a pattern emerges: The lower the present-day $T_{\rm int}$ and the more massive the atmosphere is, the longer it takes to cool down to that state from an initial hot state ($T_{\rm int}$ much larger than 100 K) after formation (time $t$=0).

Figure 8 (top left) shows evolution tracks for different present-day $T_{\rm int}$ and a heavy volatile mass fraction $Z$=0.85. The largest $T_{\rm int}$ that yields a cooling time in agreement with the system’s age is 52 K, although $T_{\rm int}$=60 K might yield a solution where the radius is matched within the measurement uncertainties. Assuming $T_{\rm int}$=100 K requires an extra heating source that maintains such a high heat flux at present times. For inflated hot Jupiters, extra heating of $\sim$0.1% to a few percent of the incident stellar irradiation is required to explain their large radius, which is consistent with Ohmic heating (Thorngren & Fortney, 2018). For $\pi$ Men c, less extra heating is required ($\sim$0.01$\%$; gray dashed curve). However, the mechanism that may provide this extra heating to the smaller $\pi$ Men c is not obvious. While Ohmic heating may occur at sub-Neptunes (Pu & Valencia, 2017), the predicted amount falls short by two orders of magnitude of what is needed for $\pi$ Men c (Figure 8, top right). Another option is tidal heating provided that the planet is on an eccentric orbit, although the orbital eccentricity remains poorly constrained (Gandolfi et al., 2018; Huang et al., 2018; Damasso et al., 2020). For an eccentricity $e$=$0.001$ and tidal quality factor $Q_{\rm{p}}$$\geq$10³ we find that tidal heating affects negligibly the planetary cooling. Tidal heating however becomes effective at extending the cooling time of an atmosphere with $T_{\rm int}$=100 K up to the current system’s age if, for example, $e$=0.02 and $Q_{\rm{p}}$$\sim$10³ (about three times the Earth’s tidal quality factor), or if $e$=0.1 and $Q_{\rm{p}}$$\sim$5$\times$10⁴ (about the Saturn/Uranus/Neptune value). Further, we estimated the circularization time for these two configurations (Jackson et al., 2008) to be $\tau_{\rm{circ}}$=0.6 and 30 Gyr respectively. Even the shortest of them is on the order of our prescribed survival time. It is thus reasonable to expect that the recent history of $\pi$ Men c’s atmosphere may have occurred while the planet followed an eccentric orbit that could have sustained $T_{\rm int}$$\sim$100 K. It is also conceivable that the outer companion in the $\pi$ Men planetary system may endow a non-negligible eccentricity to the innermost planet’s orbit (De Rosa et al., 2020; Kunovac Hodžić et al., 2020; Xuan & Wyatt, 2020), as is found in some close-in sub-Neptune-plus-cold Jupiter systems. Importantly, to the effects of planet mass loss, higher $T_{\rm int}$ values that imply lower $M_{\rm atm}$ are less likely to survive over time scales of Gigayears, which sets a limit to how high $T_{\rm int}$ can be given that $\pi$ Men c still hosts an atmosphere. Finally, we adopt $T_{\rm int}$=100 K as a reasonable upper limit for $\pi$ Men c.

Figure 8 (bottom left) shows evolution tracks for $Z$ from 0.5 to 0.88 at present-day $T_{\rm int}$=52 K. The lower the adopted $Z$, the less massive the atmosphere is and the quicker it cools down and contracts. Eventually, the planet adopts a state of (nearly) equilibrium evolution with the incident flux, where contraction slows down and further cooling progresses on very long timescales. In this state, the planet radiates ($\sim$(100/1150)⁴$\sim$6$\times$10^-5) only 0.006% more than if it was in true equilibrium with the stellar irradiation. Mass fractions $Z$$<$0.5 may be possible for present $T_{\rm int}$$<$44 K. For comparison, Saturn has $T_{\rm int}$$\sim$77 K and Neptune $\sim$53 K, and thus we do not expect that $T_{\rm int}$ at $\pi$ Men c will be much higher than for them given that its mass is much lower. Also, a non-zero eccentricity may keep $T_{\rm int}$ above such values. Mass loss, fixed in the preparation of Figure 8 (bottom left) to 2$\times$10¹⁰ g/s, prolongs the cooling time immediately after formation but speeds up the contraction well into the future as the planet loses its envelope.

We have simulated the future evolution of $\pi$ Men c for some of the compositions deemed realistic for present-day $\pi$ Men c with the goal of exploring whether the planet will ever cross the radius valley. As an example, Figure 8 (bottom right) shows that for $Z$(CO₂)=0.50 and $T_{\rm{int}}$=44 K, $\pi$ Men c will turn into a bare rocky core in about 0.5 Gyr. The same conclusion is found for other $Z$–$T_{\rm{int}}$ combinations.

The interior structure model calculates the atmospheric pressure-temperature ($p$–$T$) profile using an analytical formulation (Guillot, 2010) that depends on the ratio $\gamma$=$\kappa_{\rm vis}$/$\kappa_{\rm IR}$ of constant visible and IR opacities in the $T$–optical depth relation and on the local opacity. For the IR opacities, we take the Rosseland mean $\kappa_{\rm R}$ for solar composition (Freedman et al., 2008) as our baseline. In our opacity 1 setting, we adopt $\gamma$=0.123 and $\kappa_{\rm IR}$=1$\times$$<$$\kappa_{\rm R}$$>$ (bracket $<$$>$ denotes average over wavelength) which was confirmed to reproduce well the $p$–$T$ profiles published for the H₂/He/H₂O atmosphere of K2-18 (Scheucher et al., 2020). As $\pi$ Men c’s atmosphere may also contain large abundances of CO₂, our opacity 2 setting adopts $\gamma$=0.500 and $\kappa_{\rm IR}$=2$\times$$<$$\kappa_{\rm R}$$>$, which is appropriate for a CO₂-dominated atmosphere with some admixture of H₂/He. Using two sets of opacities for each interior model calculation is a pragmatic way of bracketing the real opacity of $\pi$ Men c’s atmosphere. In a conservative spirit, our constraints on the planet’s bulk composition from the atmospheric stability argument utilize the opacity setting that results in the lowest $Z$.

Caption for Table 2.

Summary of photochemical-hydrodynamic model runs. Columns 2–4 quote the assumed nuclei fractions for H, C and O in the bulk atmosphere. Columns 6–13 quote the adopted vmrs at the $\mu$bar pressure level for H₂, H, He, CO₂, CO, C, H₂O and O, and 14 quotes the corresponding heavy mass fraction. Column 15 quotes the mass flux of H i atoms at the location of the sonic point, which occurs for all cases within the range $r_{{Mach}=1}$/$R_{\rm{p}}$=3.5–4.6. Columns 16–20 quote the mass loss rates of H, He, C and O nuclei, and 21 quotes the net mass loss rate. $\dagger$: For these cases, we specified the vmrs at 1 $\mu$bar considering that molecules such as H₂, H₂O and O₂ become photochemically unstable before reaching the $\mu$bar pressure level. We also assume that the 1-$\mu$bar pressure level occurs at the radial distance $\sim$2.06$R_{\earth}$ for the TESS radius of the planet. This neglects the vertical extent of the region between a few mbars and 1 $\mu$bar, a reasonable approximation for moderate atmospheric temperatures. We impose that at the 1-$\mu$bar level the temperature coincides with the planet’s $T_{\rm{eq}}$ (1150 K).

Caption for Figure 8.

Top Left. Temporal evolution of planet radius $R_{\rm{p}}$($t$) scaled by the measured $R_{\rm{p}}$=2.06$R_{\earth}$ for models with $Z$(CO₂)=0.85 and present-day $T_{\rm{int}}$ values as labeled. Both $R_{\rm{p}}$ and $T_{\rm{int}}$ decrease with time as the planet cools. The gray box indicates the uncertainty in present radius ($\pm$0.03$R_{\earth}$) and stellar age (5.2$\pm$1.1 Gyr). For the gray dashed curve of $T_{\rm{int}}=100$ K, the evolution calculations include an extra energy source that converts 0.01% of the incident energy flux into heating of the interior. Top Right. Cooling time, defined as the time to reach the measured planet radius $R_{\rm{p}}$. Symbols are placed at 5.1 Gyr (13 Gyr) if contraction is halted, i.e. if $R_{p}$ becomes independent of time and agrees with (remains larger than) the measured $R_{p}$. All runs are for interior models with $T_{\rm{int}}$=100 K, $Z$=0.85, and no mass loss. Bottom Left. Same as top left panel but for models with a lower $T_{\rm int}$=52 K and $Z$ values as labeled. The moderate $Z$=0.50 model seems to fall short in cooling time. However, evolving that model further to $T_{\rm int}$=44 K leads to equilibrium evolution and $R_{\rm{p}}$(t) still in agreement with the measured value. Bottom Right. Planet size evolution considering a fiducial mass loss of 2$\times$10¹⁰ g s^-1 (blue) and omitting mass loss (red) for $Z$(CO₂)=0.5 and $T_{\rm{int}}$=44 K. Equilibrium is reached when mass loss is omitted. When mass loss is taken into account, the planet continues shrinking as it sheds its atmosphere until it turns into a bare rocky core, which is predicted to happen within 0.5 Gyr from now. The squares indicate potential present states of $\pi$ Men c, from where the evolution is calculated backward and forward in time.