Dusty outflows in planetary atmospheres: Understanding "super-puffs" and transmission spectra of sub-Neptunes

TL;DR

This paper proposes that dusty, non-static atmospheric outflows in super-puff exoplanets explain their large apparent radii and flat transmission spectra, challenging static atmospheric models and impacting exoplanet characterization.

Contribution

It introduces a novel scenario of dusty outflows carrying small dust grains to high altitudes, explaining observations of super-puffs and their flat spectra, which static models cannot.

Findings

Dusty outflows can inflate transit radii.

High-altitude dust flattens transmission spectra.

Scenario explains super-puff observations.

Abstract

`Super-puffs' are planets with anomalously low mean densities ($\lesssim 10^{-1}~{\rm g\ cm}^{-3}$). With a low surface gravity, the extended atmosphere is susceptible to extreme hydrodynamic mass loss (`boil off') on a timescale much shorter than the system's age. Even more puzzling, super-puffs are estimated to have a scale height of $\sim 3000~{\rm km}$, yet recent observations revealed completely flat transmission spectra for Kepler 51b and 51d. We investigate a new scenario that explains both observations: non-static outflowing ($\dot{M}\gtrsim 10^{-10}~M_\oplus~{\rm yr}^{-1}$) atmospheres that carry very small dust grains ($\sim 10~{\rm A}$ in size, $\sim 10^{-2}$ in mass fraction) to high altitudes ($\lesssim 10^{-6}~{\rm bar}$). Dust at high altitudes inflates the observed transit radius of the planet while flattens the transmission spectra.Previous static atmospheric models struggles to achieve cloud elevation and production of photochemical haze at such high altitudes. We propose to test this scenario by extending the wavelength coverage of transmission spectra. If true, dusty atmospheric outflows may affect many young ($\lesssim 10^9~{\rm yr}$), low mass ($\lesssim 10~M_\oplus$) exoplanets, thereby limit our ability to study the atmospheric composition in transmission, and inflate the observed transit radius of a planet hence obscure the underlying mass-radius relationship.

Figures

Tables

Table 1. Table 1: Properties of the representative models

Model	${\tau }_{\mathrm{kh}}$	${\rho }_{\mathrm{ini}}\left(r_{\mathrm{in}}\right)$	$m_d/m_{\mathrm{gas}}$	${\dot{M}}_{10}^{†}$
	(${10}^9\mathrm{yr}$)	(${10}^{-8}g{\mathrm{cm}}^{-3}$)	(${10}^{-2}$)
0 (Parker wind)	2.3	18.9	1.7	4.0
1 (Optical & IR)	6.2	0.37	2.6	5.4
2 (UV & X-ray)			4.8	6.5

Equations

$$\begin{array}{rl}& p=p_{\mathrm{p}}\exp \left[{\beta }_{\mathrm{p}}\left(\frac{r_{\mathrm{p}}}{r}-1\right)\right]\ge p_{\infty }=p_{\mathrm{p}}{\mathrm{e}}^{-{\beta }_{\mathrm{p}}};\\ & \rho =\frac{p\mu }{k_{\mathrm{B}}{T}_{\mathrm{eq}}};{\beta }_{\mathrm{p}}\equiv \frac{G{M}_{\mathrm{p}}\mu }{r_{\mathrm{p}}{k}_{\mathrm{B}}{T}_{\mathrm{eq}}},\end{array}$$

$$\begin{array}{rl}& {\dot{M}}_{\mathrm{Parker}}\sim 4\pi r_{\mathrm{s}}^2{c}_{\mathrm{s}}{\rho }_{\mathrm{p}}\exp \left(\frac{3}{2}-\frac{2r_s}{r_{\mathrm{p}}}\right),\\ & {\dot{M}}_{\mathrm{rad}}\sim \left(\frac{L_{\ast }}{4\pi a^2}\right)\pi r_{\mathrm{p}}^2{\left(\frac{c_s^2}{2}\right)}^{-1}.\end{array}$$

$$\begin{array}{rl}{v}_{\mathrm{term}}& \sim 3\mathrm{m}{\mathrm{s}}^{-1}\times {\left(\frac{\mu }{m_{\mathrm{H}}}\right)}^{-1/2}\left(\frac{M_{\mathrm{p}}}{M_{\oplus }}\right){\left(\frac{r}{R_{\oplus }}\right)}^{-2}\\ & \times \left(\frac{r_{\mathrm{d}}}{10\mathring{\mathrm{A}}}\right)\left(\frac{\rho /m_{\mathrm{H}}}{10^{12}{\mathrm{cm}}^{-3}}\right){\left(\frac{T}{10^3\mathrm{K}}\right)}^{-1/2};\\ {\tau }_{\mathrm{prec},\mathrm{p}}& \equiv \frac{r_{\mathrm{p}}}{v_{\mathrm{term},\mathrm{p}}}\sim 10^{-1}\mathrm{yr}{\left(\frac{r_{\mathrm{d}}}{10\mathring{\mathrm{A}}}\right)}^{-1}\left(\frac{{\rho }_{\mathrm{p}}/m_{\mathrm{H}}}{10^{12}{\mathrm{cm}}^{-3}}\right).\end{array}$$

$$\begin{array}{rl}{\dot{M}}_{\mathrm{crit}}& \equiv 4\pi r^2\rho v_{\mathrm{term}}\simeq 2\times 10^{-11}{M}_{\oplus }{\mathrm{yr}}^{-1}\left(\frac{M_c}{M_{\oplus }}\right)\\ & \times \left(\frac{r_{\mathrm{d}}}{10\mathring{\mathrm{A}}}\right){\left(\frac{T}{10^3\mathrm{K}}\right)}^{-1/2}{\left(\frac{\mu }{m_{\mathrm{H}}}\right)}^{1/2}.\end{array}$$

$$⟨r_{\mathrm{eff}}⟩\simeq {\left\{\frac{1}{\pi }{\int }_0^{R_{\mathrm{H}}}\mathrm{d}b2\pi b\left[1-{\mathrm{e}}^{-\tau (b)}\right]\right\}}^{1/2},$$

$${\sigma }_{\mathrm{d},\mathrm{ext}}\simeq {\sigma }_{-16}\times 10^{-16}{\mathrm{cm}}^2{\left(\frac{r_{\mathrm{d}}}{10\mathring{\mathrm{A}}}\right)}^3{\left(\frac{\lambda }{\mu \mathrm{m}}\right)}^{-\delta },$$

$$\frac{m_{\mathrm{d}}}{m_{\mathrm{gas}}}\simeq X_{\mathrm{d}}\left(\frac{m_{\mathrm{C}}{N}_{\mathrm{C},\mathrm{dust}}}{m_{\mathrm{H}}}\right)\simeq 0.56{\left(\frac{r_{\mathrm{d}}}{10\mathring{\mathrm{A}}}\right)}^3\left(\frac{X_{\mathrm{d}}}{10^{-4}}\right),$$

$$\frac{\mathrm{d}{M}_{\mathrm{a}}}{\mathrm{d}r}=4\pi r^2\rho ,\frac{\mathrm{d}p}{\mathrm{d}r}=-\frac{G(M_{\mathrm{a}}+M_c)\rho }{r^2},p=\kappa {\rho }^{\gamma },$$

$$\exp \left(-\frac{{\mathcal{M}}^2}{2}\right)=\varrho \exp \left(\frac{3}{2}-\frac{2}{\eta }\right),\varrho =\frac{1}{\mathcal{M}{\eta }^2},$$

$$\sum _{h\nu }F(h\nu ){\sigma }_{\mathrm{d},\mathrm{ext}}(h\nu )=4\pi r_{\mathrm{d}}^2{\sigma }_{\text{textsc}sb}{\tilde{T}}_{\mathrm{d}}^{4}q({\tilde{T}}_{\mathrm{d}}).$$

Full text

planets and satellites: atmospheres — planets and satellites: composition — planets and satellites: formation — planets and satellites: physical evolution — method: numerical

11footnotetext: Center for Computational Astrophysics, Flatiron Institute,

New York, NY 10010; lwang@flatironinstitute.org22footnotetext: Princeton University Observatory, Princeton, NJ 0854433footnotetext: Department of Physics and Kavli Institute for Astrophysics

and Space Research, Massachusetts Institute of Technology,

Cambridge, MA 02139

1 Introduction

“Super-puffs” are planets that have sub-Neptune masses ($\lesssim 5~{}M_{\oplus}$) but gas-giant transit radii ($\gtrsim 5~{}R_{\oplus}$), and thus extremely low mean densities ($<10^{-1}~{}{\rm\,g}~{}{\rm\,cm}^{-3}$) and large scale heights ($\sim 3000~{}{\rm\,km}$). A prime example is Kepler 51b, which has a $\sim 7R_{\oplus}$ transit radius but a mass of only $\sim 2.1~{}M_{\oplus}$ (consolidated by independent transit timing variation analyses of several groups e.g. Roberts et al. in prep; Masuda, 2014, M14 hereafter). The ensemble of discovered super-puffs include Kepler 51c, 51d; Kepler 79d, 79e (Jontof-Hutter et al., 2014); and Kepler 87c (Ofir et al., 2014). In this letter we concentrate our discussions on the well-studied Kepler 51b unless specially noted.

Recent works (Owen & Wu 2017; Wang & Dai 2018, WD18 hereafter) suggest that hydrodynamic and photoevaporative loss of atmospheres might be a ubiquitous effect responsible of the observed bimodal radius distribution of close-in sub-Neptune planets (Fulton et al., 2017). Given their low surface gravity, super-puffs are expected to have excessive hydrodynamic mass-loss even without stellar high energy radiation (“boil-off”, see also Owen & Wu, 2016), and should disperse on a timescale of $\sim 10^{3}~{}\mathrm{yr}$ (§2.1), much shorter than the system’s age ($\sim 0.3~{}{\rm Gyr}$ for Kepler 51 from gyrochronology; M14). Similarly, Lammer et al. (2016) noted that CoRoT-24b must also have high-altitude aerosols to increase the apparent transit radius, thereby lowering the implied mass loss rate. However, they were agnostic of how aerosols could form or be lifted to such hight altitudes.

Super-puffs, with their large scale heights, are considered ideal targets for transmission spectroscopy. However, the HST WFC3 observation of Kepler 51b and 51d yield flat transmission spectra in the near-infrared (Roberts et al. in prep). This is reminiscent of the flat spectrum of GJ1214b (Kreidberg et al., 2014). If cloud/haze are invoked to mute the absorption features, they have to be advected to or produced at such a high altitude that current models would struggle with (§2.2). We hereby consider a non-static atmosphere characterized by a slow hydrodynamic outflow ($\gtrsim 10^{-10}~{}M_{\oplus}~{}\mathrm{yr}^{-1}$), producing a relatively small mass-loss over the age of Kepler 51. Dust grains can be carried to much higher altitude in this outflow, increasing the observed transit radius to $\sim 7~{}R_{\oplus}$ while muting signatures of other species in the atmosphere.

2 Basic ideas

2.1 Isothermal atmosphere: Inevitable escape

Generally a planetary atmosphere can be divided into a convective isentropic interior and a radiation-dominated, approximately isothermal exterior (Rafikov, 2006; Owen & Wu, 2016; Ginzburg et al., 2016). In the isothermal layer, hydrostatic density and pressure profiles are given by,

[TABLE]

where the subscripts “p” and “$\infty$” denote the quantities at the planetary radius and infinite radius respectively, $G$ is the gravitational constant, $k_{\mathrm{B}}$ is the Boltzmann constant, $M_{\mathrm{p}}$ is the planetary mass (core and atmosphere combined), $\mu$ is the (dimensional) mean molecular mass, and $T_{\mathrm{eq}}\simeq 886~{}{\rm\,K}\ \left({L_{*}}/{L_{\odot}}\right)^{1/4}\left({a}/{0.1~{}{\rm\,au}}\right)^{-1/2}$ is the equilibrium temperature at planetary orbit radius $a$ and host star luminosity $L_{*}$. The dimensionless parameter $\beta_{\mathrm{p}}$ is also called the “restricted Jeans parameter” (e.g. Fossati et al., 2017; Cubillos et al., 2017). We also remind the reader that $p_{\infty}$ in eq. 1 serves as a confining term preventing the isothermal atmosphere from a spontaneous outflow. If one naively assumes a clear atmosphere (free of cloud/haze) of solar abundance, $p_{\mathrm{p}}\sim 20-100~{}{\rm\,mbar}$ is required at the observed transit radius (e.g. Lopez & Fortney, 2014; Lammer et al., 2016). For Kepler 51b, this leads to $\beta_{\mathrm{p}}\simeq 9.8$ and $p_{\infty}\sim 10^{-6}~{}{\rm\,bar}$ using eq. (1) Such $p_{\infty}$ is a few orders of magnitude greater than any plausible sun-like stellar wind total pressure (Murray-Clay et al., 2009). Unconfined atmospheres hydrodynamically lose mass at $\dot{M}\sim\mathrm{min}\{\dot{M}_{\rm Parker},\ \dot{M}_{\rm rad}\}$, where (e.g. Parker, 1958),

[TABLE]

Here $c_{s}=(k_{\mathrm{B}}T_{\mathrm{eq}}/\mu)^{1/2}$ is the isothermal sound speed and $r_{s}=GM_{\mathrm{p}}/(2c_{\mathrm{s}}^{2})$ is the sonic radius. We find $\dot{M}\sim 10^{-3}~{}M_{\oplus}~{}\mathrm{yr}^{-1}$ with $p_{\mathrm{p}}\sim 10~{}{\rm\,mbar}$ for Kepler 51b, dispersing the atmosphere in $\lesssim 10^{3}~{}\mathrm{yr}$—much shorter than the estimated age of the system ($\sim 0.3~{}{\rm Gyr}$), which in turn questions the earlier assumption of “clear” atmosphere.

2.2 Dusts in the Atmospheres

Aerosols, which could consist of dusts and liquid droplets, could dramatically increase the opacity of gas. The enhanced opacity lowers the required pressure at the apparent planet radius $p_{\mathrm{p}}$ by several orders of magnitude, giving rise to a much slower outflow. However, maintaining aerosol particles at a radius as high as $7~{}R_{\oplus}$ over Kepler 51b is difficult in a static atmosphere. In-situ formation of dusts (for clouds/haze) demands rather high gas density; photochemical calculations reveal that dust formation is very inefficient below $p\sim 10^{-7}-10^{-6}~{}{\rm\,bar}$ (Morley et al., 2012, 2013; Fortney et al., 2013; Kawashima & Ikoma, 2018). Aerosols are also subject to planetary gravity; dust grains with radius $r_{\mathrm{d}}$ precipitate at terminal velocity $v_{\rm term}$ and timescale $\tau_{\rm prec}$ (Baines et al., 1965; Draine, 2011),

[TABLE]

The eddy diffusion coeffient required to lift $r_{\mathrm{d}}=10~{}\mathrm{\AA}$ dusts to $\sim 7~{}R_{\oplus}$ is at least $K_{zz}\sim 10^{11}-10^{12}~{}{\rm\,cm}^{2}~{}\mathrm{s}^{-1}$, which is significantly greater than the values observed on the Earth (Pilinski & Crowley, 2015) and modeled on exoplanets (Morley et al., 2013). Even if dust formation at high altitudes were sufficient to compensate dust precipitation, in a static atmosphere with the eq. (1) density profile, heavy elements in this layer are rapidly depleted at timescale $\lesssim(m_{\rm metal}/m_{\mathrm{d}})\tau_{\rm prec}$ [here $(m_{\rm metal}/m_{\mathrm{d}})$ is the atmospheric mass ratio of metal elements to dusts].

We thus consider non-static atmospheres in which aerosols are co-moving with outflows. The critical mass-loss rate, at which $v_{r}=v_{\rm term}$ (note that this Equation does not depend on $r$; see also WD18)

[TABLE]

Whenever $\dot{M}\gg\dot{M}_{\mathrm{crit}}$, dusts experience neglibible precipitation, and can be considered as co-moving with gas. $\dot{M}$ must also satisfy $\dot{M}<\dot{M}_{\mathrm{max}}\sim(M_{\rm atm}/\tau_{\mathrm{p}})$, where $M_{\rm atm}$ is the total mass of atmosphere and $\tau_{\mathrm{p}}$ is the planet’s age (approximated by the host star’s age $\tau_{*}$; for Kepler 51b, $\dot{M}_{\mathrm{max}}\sim 10^{-9}~{}M_{\oplus}~{}\mathrm{yr}^{-1}$). Dusts of $\sim 10~{}\mathrm{\AA}$ sizes should be abundantly produced by geological activities, while laboratory experiments (Zhao et al., 2018) show that gas-phase formation of tiny graphites and polycyclic aromatic hydrocarbon (PAH) can also be very efficient even at relatively low temperatures and UV intensities. Meanwhile, the temperature throughout most of the internal atmosphere (§3.1) is higher than dust sublimation temperature ($\sim 1500~{}{\rm\,K}$), preventing tiny grains from fast coagulating: larger grains fall back to the internal atmosphere and are broken into gaseous species.

2.3 Effective transit radii

High-altitude aerosols lead to extra extinction on stellar light from the observer’s view, thus effectively increases the planet trasiting radii. To ease later discussion, we define the effective transit radius:

[TABLE]

where $\tau(b)$ is the optical depth along the line-of-sight (LoS) at impact parameter $b$ relative to the planet geometric center. The upper limit of the integral is $R_{\mathrm{H}}$ (the planet’s Hill radius) where the assumption of excluding host star gravitation likely breaks down. We estimate the optical depth by $\tau(b)\simeq\Sigma(b)X_{\mathrm{d}}\sigma_{\rm d,ext}$, where $\Sigma(b)$ is the column density along the LoS, $X_{\mathrm{d}}\simeq n_{\mathrm{d}}/(\rho/m_{\mathrm{H}})$ is the number fraction of dust particles relative to hydrogen nuclei, and $\sigma_{\rm d,ext}$ is the extinction cross section of a single dust particle. At optical and infrared (IR) wavelengths $0.2\lesssim(\lambda/{\rm\,\mu m})\lesssim 2$, the extinction cross section of very small grains is well approximately given by a smooth power-law function,

[TABLE]

where $(\sigma_{-16},\delta)\simeq(0.92,1.55)$ for graphites, and $(0.11,0.93)$ for silicates (Draine & Malhotra, 1993). PAH grains at $r_{\mathrm{d}}\sim 10~{}\mathrm{\AA}$ have an absorption edge at $\lambda\sim 1~{}{\rm\,\mu m}$, and are optically similar to graphites at shorter wavelengths (Li & Draine, 2001). For simplicity we assume that all aerosols consist of graphite dusts. The dust-to-gas mass ratio corresponding to number ratio $X_{\mathrm{d}}$ is, assuming hydrogen atmosphere,

[TABLE]

where $N_{\rm C,dust}\simeq 470~{}(r_{\mathrm{d}}/10~{}\mathrm{\AA})^{3}$ is the number of carbon atoms per dust grain.

3 Detailed Modeling

3.1 Isentropic interior

Although all interesting atmospheric dynamics take place in the radiative exterior, hydrodynamic structures of the convective interior should still be consistently calculated by solving,

[TABLE]

where $M_{\rm a}$ is the mass of atmosphere enclosed by radius $r$, $M_{c}$ is the mass of the solid planet core, $\kappa$ is the specific entropy parameter, and $\gamma$ is the adiabatic index (we take $\gamma=1.4$ for the molecular atmospheres in this letter).

The gravitation in the radiative exterior of atmosphere depends on both $M_{c}$ and the total mass of isentropic atmosphere $M_{\rm atm}$, while the self-gravity of the gas in that layer is usually negligible. In practice, we first pick an $M_{c}$ and an $M_{\rm atm}$ and obtain a model of the external radiative atmosphere. Then, we solve eq. (8) as a boundary value problem such that (1) $M_{\rm a}(R_{c})=0$ ($R_{c}$ is the planet core radius), and (2) $p$ and $\rho$ match the external atmosphere profiles at the radiative-convective boundary $r_{\mathrm{rcb}}$, which is adjusted so that $M_{\rm a}(r_{\mathrm{rcb}})=M_{\rm atm}$. The isentropic atmosphere is characterized by its Kelvin-Helmholtz timescale $\tau_{\rm kh}$ (e.g. Owen & Wu, 2017).

3.2 Dusty outflowing exterior

The model planet orbits the host star (for simplicity, we round off to $M_{*}=M_{\odot}$, $L_{*}=0.88L_{\odot}$ from M14) on a $a=0.25~{}{\rm\,au}$ circular orbit ($T_{\mathrm{eq}}=543~{}{\rm\,K}$). The planet combines an $M_{c}=1.7~{}M_{\oplus}$, $R_{c}=1.14~{}R_{\oplus}$ solid core and an $M_{\rm atm}=0.4~{}M_{\oplus}$ convective atmosphere.

3.2.1 Model 0: Isothermal Parker wind

The first model (Model 0) that we consider is constructed analytically. If we assume an isothermal $T=T_{\mathrm{eq}}$, the well-known Parker wind solution satisfies (Parker, 1958),

[TABLE]

where $\mathcal{M}\equiv v_{r}/c_{s}$ is the radial Mach number, $\eta\equiv r/r_{s}$ is the dimensionless radius normalized by the sonic radius $r_{s}$, and $\varrho\equiv\rho/\rho_{s}$ is the dimensionless density normalized by $\rho_{s}$ (the density at sonic radius).

3.2.2 Consistent thermochemical simulations

Models 1 and 2 involve full hydrodynamic simulations that incorporate radiation and thermochemistry described in WD18. The axisymmetric 2.5-dimensional spherical-polar mesh centers at the planet, whose polar axis points to the host star. It spans $(r,\theta)\in[3~{}R_{\oplus},400~{}R_{\oplus}]\otimes[0,\pi]$ at resolution $256\times 128$ (radial zones are spaced logarithmically and latitudinal zones evenly), to guarantee that all relevant physical processes are included in the simulation domain. The initial conditions obey the isothermal hydrostatics at $T_{\mathrm{eq}}$ in eq. (1), where $\rho_{\rm ini}(r_{\rm in})$ (the initial mass density at the inner boundary $r_{\rm in}=3~{}R_{\oplus}$) is the variable parameter. Initial abundances of chemical species are uniform across the simulation domain; they are identical to WD18, except for the dusts. We adjust $\rho_{\rm ini}(r_{\rm in})$ and the dust-to-gas mass ratio $(m_{\mathrm{d}}/m_{\rm gas})$ for each simulation so that $\langle r_{\mathrm{eff}}\rangle=r_{\mathrm{p}}=7~{}R_{\oplus}$ and $\dot{M}_{\mathrm{max}}>\dot{M}\gg\dot{M}_{\mathrm{crit}}$ in steady states.

Both models include the host star luminosity $L(2~{}{\rm\,eV})=0.88~{}L_{\odot}$, representing infrared and optical radiation. Model 2 also involves high-energy photons represented by four photon energy bins ($h\nu=7~{}{\rm\,eV}$ for soft FUV, $12~{}{\rm\,eV}$ for Lyman-Werner band FUV, $25~{}{\rm\,eV}$ for EUV, and $3~{}{\rm\,keV}$ for the X-ray) at luminosities111These high-energy luminosities are estimated with the recipes in Owen & Wu (2017) and WD18, adopting Ribas et al. (2005) for $L(t<10^{8}~{}\mathrm{yr})$ and assuming $L/L(t<10^{8}~{}\mathrm{yr})=\mathrm{min}\{1,(\tau_{*}/10^{8}~{}\mathrm{yr})^{-1.5}\}$.: $L(7~{}{\rm\,eV})=L(25~{}{\rm\,eV})=L(3~{}{\rm\,keV})=8\times 10^{28}~{}{\rm\,erg}~{}\mathrm{s}^{-1}$, and $L(12~{}{\rm\,eV})=8\times 10^{27}~{}{\rm\,erg}~{}\mathrm{s}^{-1}$. Rays are parallel to the polar axis, entering the simulation domain at the outer radial boundary with fluxes $F(h\nu)=L(h\nu)/(4\pi a^{2})$.

We include $r_{\mathrm{d}}=10~{}\mathrm{\AA}$ graphites in these two models as a proxy of dusts of all sizes and components. Dust temperature is estimated by the dual-temperature profile $T_{\mathrm{d}}=\mathrm{max}\{T_{\mathrm{eq}},\tilde{T}_{\mathrm{d}}\}$ (similar to Chiang & Goldreich 1997), where $\tilde{T}_{\mathrm{d}}$ is obtained by solving

[TABLE]

Here $\sigma_{\textsc{sb}}$ is the Stefan-Boltzmann constant, and $q(T_{\mathrm{d}})$ is the dust emissivity.

3.3 Results

3.3.1 Model profiles

Table 1 summarizes the key properties and results of our models. All models demand $\sim 10^{-2}$ of atmospheric mass in dusts to achieve $\langle r_{\mathrm{eff}}\rangle=7~{}R_{\oplus}$ with $M_{\mathrm{p}}=2.1~{}M_{\oplus}$. The gas pressure required at $r_{\mathrm{p}}$ is merely $\sim 10^{-8}-10^{-9}~{}{\rm\,bar}$, while the $p=100~{}{\rm\,mbar}$ radii is much lower (compared to §2.1): $r_{100~{}{\rm\,mbar}}\simeq 2.5~{}R_{\oplus}$ (Model 0) or $2.2~{}R_{\oplus}$ (Models 1 and 2). Density, temperature and radial velocity profiles along the radial ray at $\theta=\pi/2$ (i.e. perpendicular to the direction to the host star) of all models are presented by Figure 1. Figure 2 illustrates the meridional plots of density, temperature and velocity profiels for Model 2 in steady state, which are similar to the EUV photoevaporation models discussed in WD18: a hot ($T>10^{4}~{}{\rm\,K}$), anisotropic EUV-dominated outflow, a warm ($T\lesssim 10^{3}~{}{\rm\,K}$) intermediate layer, and a “tail” behind the night hemisphere.

Curiously, there are also day-night meridional motions in Models 1 and 2. This is the consequence of dust temperature excess: in regions accessible by $h\nu=2~{}{\rm\,eV}$ photons, dust temperature $T_{\mathrm{d}}\simeq 1050~{}{\rm\,K}\gg T_{\mathrm{eq}}$ due to $q(T_{\mathrm{d}})\ll 1$ (eq. 10), causing gas temperature $T\gg T_{\mathrm{eq}}$ via dust-gas thermal accommodation. Figure 2 illustrate such meridional motion, which never leaves the planetary gravity potential, but still satisfies $|v|\gg v_{\rm term}$ hence can keep the dusts aloft. We nonetheless choose not to over-interpret this result: atmospheric circulation requires proper treatment of radiative transfer, dimensionality and planet spin to model, which are postponed to future works.

3.3.2 Transit light curves and model consistency

Figure 3 illustrates the synthetic transit light curves (limb darkening profile adopted from M14), plus a simple “hard sphere” for reference. All models have a extended but gentler ingress/egress than the hard-sphere. Model 2 has a relatively sharper ingress/egress, because EUV photons carve a cliff in density and temperature by launching a photoevaporative wind. The synthetic light curve is symmetric about the mid-transit, as $\langle r_{\mathrm{eff}}\rangle=7~{}R_{\oplus}$ is still deep in the planet’s potential well. To analyze the detectability of the difference in the light curves, we re-sample systhetic light curves with 1-minute cadence and add a white noise component of $1500~{}{\rm ppm}$ to mimic the Kepler observation of Kepler 51b. The resultant light curves were analyzed with a conventional Mandel & Agol (2002) transit model similar to that employed by M14. We found that more extended and gentler ingress/egress of the synthetic light curves can be accommodated by a combination of higher impact parameter $b$ and slightly different limb darkening coefficients than those reported by M14. A future observation of the system with higher photometric precision is required to distinguish Models 0 through 2 which differs by only $\sim 200~{}{\rm ppm}$.

4 Discussion and Summary

In this letter, we showed that a dusty outflow of a planetary atmosphere could enhance the opacity at high altitudes, therefore successfully explains the puffy Kepler 51b, and flat transmission spectrum of super-puff exoplanets. The dusty outflow scenario relies on the mass-loss rate $\dot{M}$, which should stay in a proper range ($\dot{M}_{\mathrm{crit}}\ll\dot{M}\ll\dot{M}_{\mathrm{max}}$; see §2.2), favoring the class of young, low-mass sub-Neptunes. Cubillos et al. (2017) suggests that $\sim 15\%$ of sub-Neptunes are too puffy and may be currently experiencing mass loss. The mechanism is maximized when the atmospheric dispersal timescale is similar to the age of the system [e.g. $\sim 0.3~{}{\rm Gyr}$ for Kepler 51 (M14), and $\lesssim 1~{}{\rm Gyr}$ for Kepler 79 (Walkowicz & Basri, 2013)].

Dusty outflows have several implications. First, extinction cross sections of small grains are smooth function of wavelengths in optical and near-infrared (see also Draine & Lee, 1984; Draine & Malhotra, 1993; Li & Draine, 2001). Dusts therefore obscure the signatures of some other chemical species in planetary atmospheres, limiting the ability of transmission spectroscopy. Figure 4 plots the strength of water features against planet mass for sub-Neptune planets (Crossfield & Kreidberg, 2017). We note a possible dichotomy that only low-mass ($\lesssim 10~{}M_{\oplus}$) planets tend to have muted absorption features. One explanation is that planets more massive than 10 M have gravitational wells too strong to allow adequate atmospheric loss, as seen in numerical explorations of WD18. Meanwhile, due to the large optical depths in Ly$\alpha~{}$(Draine, 2011) and the metastable helium line (Oklopcic & Hirata, 2018), a simple calculation show that both lines should still be observable by transmission spectra for planets undergoing dust outflows. Second, the observed $r_{\mathrm{p}}$ may differ significantly from the predicted radius assuming a clear atmosphere (§3.3). A key objective of the TESS mission is to accurately measure the masses and radii of $>50$ sub-Neptunes, followed by ensemble analyses of their compositions, which may be significantly biased if leaving dusty outflows unaccounted for. Third, as $\sigma_{\rm d,ext}$ increases at shorter wavelengths, $\langle r_{\mathrm{eff}}\rangle$ in optical bands should be greater than infrared. The transiting radii yielded by eq. (6) at $\lambda=0.5~{}{\rm\,\mu m}$ are $\sim 10-20\%$ greater than $\lambda=1~{}{\rm\,\mu m}$. Such phenomenon has been observed for a few exoplanets (e.g. Ehrenreich et al., 2014). Extending wavelength coverage of transmission spectra (e.g. Spitzer) should also be able to detect more dust-specific signatures.

This work is supported by the Center for Computational Astrophysics of Flatiron Institute, and the Department of Astrophysical Sciences of Princeton University. We thank our colleagues (alphabetical order): Xue-Ning Bai, Adam Burrows, Jeremy Goodman, Xiao Hu and Kento Masuda, for helpful discussions and comments.