# Advection of Potential Temperature in the Atmosphere of Irradiated Exoplanets: A Robust Mechanism to Explain Radius Inflation

Pascal Tremblin, G. Chabrier, J. Mayne, S. Amundsen, I. Baraffe, F. Debras, B. Drummond, J. Manners, S. Fromang

► To cite this version:

Pascal Tremblin, G. Chabrier, J. Mayne, S. Amundsen, I. Baraffe, et al. Advection of Potential Temperature in the Atmosphere of Irradiated Exoplanets: A Robust Mechanism to Explain Radius Inflation. The Astrophysical Journal, 2017, 841 (1), 10.3847/1538-4357/aa6e57. hal-01834142

HAL Id: hal-01834142

https://hal.science/hal-01834142v1

Submitted on 28 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Advection of Potential Temperature in the Atmosphere of Irradiated Exoplanets: A Robust Mechanism to Explain Radius Inflation

P. Tremblin<fnref n="1" />,<fnref n="2" />, G. Chabrier<fnref n="2" />,<fnref n="3" />, N. J. Mayne<fnref n="2" />, D. S. Amundsen<fnref n="4" />,<fnref n="5" />, I. Baraffe<fnref n="2" />,<fnref n="3" />, F. Debras<fnref n="2" />,<fnref n="3" />, B. Drummond<fnref n="2" />, J. Manners<fnref n="2" />,<fnref n="6" />, and S. Fromang<fnref n="7" />

Received 2017 February 20; revised 2017 April 18; accepted 2017 April 18; published 2017 May 19

## Abstract

The anomalously large radii of strongly irradiated exoplanets have remained a major puzzle in astronomy. Based on a two-dimensional steady-state atmospheric circulation model, the validity of which is assessed by comparison to three-dimensional calculations, we reveal a new mechanism, namely the advection of the potential temperature due to mass and longitudinal momentum conservation, a process occurring in the Earth's atmosphere or oceans. In the deep atmosphere, the vanishing heating flux forces the atmospheric structure to converge to a hotter adiabat than the one obtained with 1D calculations, implying a larger radius for the planet. Not only do the calculations reproduce the observed radius of HD 209458b, but also reproduce the observed correlation between radius inflation and irradiation for transiting planets. Vertical advection of potential temperature induced by non-uniform atmospheric heating thus provides a robust mechanism to explain the inflated radii of irradiated hot Jupiters.

**Key words:** atmospheric effects – methods: numerical – planets and satellites: general

## 1. Introduction

The anomalously large radii of irradiated hot Jupiters are one of the most intriguing problems in our understanding of extrasolar giant planets. While various physical mechanisms have been proposed to resolve this puzzle (see Baraffe et al. 2010, 2014; Fortney & Nettelmann 2010, for reviews), none of them thus far provide a satisfactory explanation. Generally speaking, they either lack a description of a robust mechanism (e.g., to explain downward transport of kinetic energy deep enough to reach the internal adiabat (Showman & Guillot 2002), or they must invoke fine-tuned conditions for resistivity (ohmic resistivity, Batygin et al. 2011; enhanced opacities, Burrows et al. 2007; ongoing layered convection, Chabrier & Baraffe 2007). On the other hand, there is now growing observational evidence to suggest a correlation between inflated radii and incident stellar flux (Demory & Seager 2011; Laughlin et al. 2011; Miller & Fortney 2011; Weiss et al. 2013; Figueira et al. 2014) and the recently observed re-inflated hot Jupiter EPIC 211351816.01 (Grunblatt et al. 2016) also supports this correlation (even though other explanations for this object might be possible). In the meantime, much effort in the community has been dedicated to the development of general circulation models (GCMs) (e.g., Showman et al. 2009; Dobbs-Dixon et al. 2012; Mayne et al. 2014, see Showman et al. 2010 for a review), but so far none of these models have been able to reveal the long-sought mechanism responsible for the inflated radii of hot Jupiters.

In this paper, we propose a new approach to the problem of the atmospheric circulation of tidally locked exoplanets. Starting from our 1D radiative/convective atmosphere code ATMO (Amundsen et al. 2014, 2017; Tremblin et al. 2015, 2016), recently applied to irradiated giant planets (Drummond et al. 2016), we construct a stationary two-dimensional (2D) circulation model of the equatorial region of the planet (see Sections 2 and 3). The

originality (and key) of this approach is to look for a stationary solution, while all existing 2D/three-dimensional (3D) atmospheric studies of irradiated hot Jupiters have only considered time-dependent simulations that cannot be run long enough to reach a steady-state in the deep atmosphere, where the thermal timescale becomes prohibitively long. In Section 4, we show that the resulting pressure/temperature and wind structures from the 2D solution compare very well with the equatorial solution from the 3D simulation. In the deep atmosphere, the 2D solution naturally converges to an adiabatic structure before becoming unstable to convection and thus connects to a hotter adiabat than that obtained using standard 1D models, which, in turn, corresponds to a significantly larger radius for the planet. We derive a reduced model to explain this behavior, providing a consistent solution for the inflated radii of hot Jupiters: the circulation induces a vertical mass flux in the deep atmosphere, advecting the potential temperature in a region with a low heating rate, thus imposing a hot adiabatic interior for the planet (see Section 5). In Section 6, we show that the models reproduce the observed trend of increasing radii with increasing irradiation, and we provide conclusions and perspectives in Section 7.

## 2. 2D Stationary Circulation Model

To describe the atmosphere, we look for stationary solutions of the Euler equations in spherical coordinates, *r* (radius), *φ* (latitude), and *λ* (longitude) (Mayne et al. 2014), where the velocity components are denoted by *u*<sub>*r*</sub>, *u*<sub>*φ*</sub>, and *u*<sub>*λ*</sub>, density is *ρ*, pressure is *P*, *Ω* is the angular velocity, and *g*(*r*) is the gravitational acceleration. A heating rate *H*<sub>rad</sub> is included, consistently calculated from the radiative transfer equation (Amundsen et al. 2016; see the Appendix). We then look for a solution with *u*<sub>*φ*</sub> = 0 (because of the north-south hemispheric symmetry) but with ∂*u*<sub>*φ*</sub>/∂*φ* ≠ 0 and *u*<sub>*λ*</sub>, *u*<sub>*r*</sub> ≠ 0 at the equator *φ* = 0.

<footnote n="1">¹ Maison de la Simulation, CEA-CNRS-INRIA-UPS-UVSQ, USR 3441, CEA Paris-Saclay, F-91191 Gif-Sur-Yvette, France; pascal.tremblin@cea.fr</footnote>

<footnote n="2">²Astrophysics Group, University of Exeter, EX4 4QL Exeter, UK</footnote>

<footnote n="3">³Ecole Normale Supérieure de Lyon, CRAL, UMR CNRS 5574, F-69364 Lyon Cedex 07, France</footnote>

<footnote n="4">⁴Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10025, USA</footnote>

<footnote n="5">⁵NASA Goddard Institute for Space Studies, 2880 Broadway, New York, NY 10025, USA</footnote>

<footnote n="6">⁶Met Office, Exeter, EX1 3PB, UK</footnote>

<footnote n="7">⁷Laboratoire AIM, CEA/DSM-CNRS-Université Paris 7, Irfu/Service d'Astrophysique, CEA Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>
The meridional momentum conservation condition trivially vanishes at the equator, so we lack an equation to find the longitudinal (zonal), vertical velocities, and the derivative of the latitudinal (meridional) velocity at the equator. The wind is essentially driven by the longitudinal pressure gradients, which implies gradients of longitudinal velocities through the longitudinal momentum conservation condition. In a steady-state, mass conservation implies that the resulting longitudinal mass flow is balanced by a combination of the vertical and meridional ones. We assume that the geometry of the wind is defined by a constant form factor $\alpha$, which corresponds to the ratio of the meridional over vertical mass fluxes. This yields the extra equation:

$$ \frac{1}{r^2} \frac{\partial r^2 \rho u_r}{\partial r} = \frac{1}{r\alpha} \frac{\partial \rho u_\phi}{\partial \phi}. \quad (1) $$

The physical significance of $\alpha$ can be highlighted by considering the following two limits, namely $\alpha \to \infty$ for a purely longitudinal and meridional wind, and $\alpha \to 0$ for a purely longitudinal and vertical wind, respectively. We use the 3D GCM simulations of Amundsen et al. (2016) to calibrate the profile of $\alpha$. Note that this factor is likely to depend on the gravity of the object. By injecting the mass conservation into the internal energy equation (see the Appendix) we get, after some algebraic manipulation, the equation of entropy conservation written as a function of the advection of the potential temperature (defined as $\theta = T(P_0/P)^{(\gamma-1)/\gamma}$, where $P_0$ is the reference pressure, and $\gamma = C_p/C_v$ is the ratio of the specific heats at constant pressure and volume, respectively):

$$ -\frac{\gamma}{\gamma-1} \mathbf{P} u \cdot \nabla (\ln \theta) = H_{\text{rad}}. \quad (2) $$

The full system of equations can be rewritten in a more compact form (see the Appendix). When there is no flow (i.e., zero velocity), the system of equations reduces to the equation of hydrostatic balance $\partial P/\partial r = -\rho g$ and a constant flux transport $F_{\text{rad}} = \text{cst}$, corresponding to our 1D scheme (in ATMO) in the absence of convection. We thus extend the same solver to include the three equations needed to compute $u_\lambda$, $u_r$, and $\partial u_\phi/\partial \phi$ at the equator, and we replace the constant flux transport equation of the 1D scheme by energy conservation including a heating rate (i.e., the derivative of the flux transport equation; see the Appendix for details).

### 3. Boundary Conditions and Input Parameters

The system of equations is solved on a fixed height grid between $R_{\text{min}}$ and $R_{\text{max}}$. The boundary conditions at the bottom of the atmosphere are (i) a linear extrapolation of the velocities from inside the computational domain (i.e., we keep a constant gradient at the boundary), (ii) an imposed maximum pressure, and (iii) an imposed radiative flux $F_{\text{rad}} = F_{\text{int}}$, where $F_{\text{int}}$ is the same as in the 3D calculation. At the top of the atmosphere, we use a linear extrapolation of all the variables ($\log(P)$ for the pressure). For an axisymmetric irradiation (or absence of irradiation), an axisymmetric solution exists with $\partial u_\phi/\partial \phi = 0$, $u_r = 0$, $u_\lambda = \text{cst}$, and any 1D solution that obeys hydrostatic balance and constant flux transport conditions. Thus, the system as proposed is not able to fully specify the wind solution, since any constant zonal wind is a solution with isotropic irradiation. This is not surprising, as in 3D models the magnitude of the zonal wind is known to depend on the full system and on the dissipation processes (physical and

<figure><img src="image_3.png" /><figcaption>Figure 1. Input profiles for the 2D stationary circulation model. The first input is the zonal wind profile at the substellar point ($\lambda = 0$), and the second input is the $\alpha$ profile, defined as the ratio of the meridional to vertical mass flux at the equator.</figcaption></figure>

numerical). Calculation of the magnitude of the zonal wind is a very active field of research with 3D GCMs, and the determination of the zonal wind strength is still a matter of debate (e.g., Fromang et al. 2016).

In our 2D approach, we opt not to try to constrain the zonal wind magnitude but to parameterize it, a different but complementary approach to those trying to explain the equatorial superrotation (e.g., Vallis 2006; Showman & Polvani 2011). We impose a zonal wind vertical profile at a given longitude of the equator (we choose the substellar point) and the 2D circulation model computes the full equatorial solution starting from this imposed profile. Therefore, the 2D circulation model depends on two inputs: a zonal wind profile at the substellar point and an $\alpha$ profile (the ratio of meridional to vertical mass flux) for the entire equatorial region. We choose input parameters that are derived from the 3D models to test our 2D solution against the equatorial solution of a 3D simulation; the profiles are displayed in Figure 1. In order to determine $\alpha$, we decided to interpolate at the equator the values of the quantities of Equation (1) from the 3D GCM run of HD 209458b from Amundsen et al. (2016). Then, at each vertical layer, we chose $\alpha$ such that it minimized the residual of Equation (1) with a least squares method on the longitude. The zonally averaged residual is smaller than 40% of the zonally average vertical mass flux; hence we conclude that this approximation is relatively good. The data from the 3D GCM were obtained after 1600 earth days of simulation, where the atmosphere has only converged above 1 bar. The values of $\alpha$ are nearly constant with time above that limit (we checked from 1000 to 1600 days), but are still evolving at deeper levels. In the region where the $\alpha$ value has not converged (for pressures larger than ~20 bars), we decided to choose a constant value. We have tested $\alpha = 1, 30, 100$ and found that the result is independent of this value. The other usual input parameters include the gravity at the surface, the internal heat flux, the elemental composition, and the stellar and orbital properties, which are all given in Table 1.

### 4. The 2D Equatorial Model versus the 3D Equator Calculation

In this section we compare the result from our 2D circulation model with those from the equatorial region of our 3D GCM simulations of HD 209458b (Amundsen et al. 2016). The results
<table><caption>**Table 1**
List of the Input Parameters Used in the 2D Circulation Model</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>Vertical resolution</td><td>50</td></tr><tr><td>Longitudinal res.</td><td>10</td></tr><tr><td>P<sub>max</sub> at R<sub>min</sub></td><td>200 bar</td></tr><tr><td>R<sub>min</sub></td><td>1.287 R<sub>jup</sub></td></tr><tr><td>R<sub>max</sub></td><td>1.48 R<sub>jup</sub></td></tr><tr><td>T<sub>int</sub></td><td>100 K</td></tr><tr><td>log g</td><td>2.97</td></tr><tr><td>γ</td><td>1.377</td></tr><tr><td>μ</td><td>2.34</td></tr><tr><td>Ω</td><td>2 × 10<sup>-5</sup> s<sup>-1</sup></td></tr><tr><td>F<sub>*</sub> Star spectrum</td><td>HD 209458 (Kurucz)</td></tr><tr><td>R<sub>*</sub></td><td>1.118 R<sub>☉</sub></td></tr><tr><td>a</td><td>0.0475 au</td></tr></tbody></table>

are portrayed in Figures 2(A) and (B), respectively. A number of characteristics are well-reproduced by the 2D solution, namely:

1. The day-night temperature gradient.

2. The shift of the hot spot in the upper atmosphere.

3. The pressure at which the horizontal temperature gradients vanish.

4. The variations of the longitudinal wind velocity.

5. The global structure and magnitude of the vertical wind.

The slight differences in the uppermost part of the atmosphere ($P \le 10^3$ bar) are likely caused by the different boundary treatments between the dynamical model and the 2D stationary model. The main differences between the 2D and 3D results are (i) the pressure/temperature structure of the deep atmosphere ($P \ge 10$ bar), and (ii) the differences between the vertical velocities in the deep atmosphere (mainly downward between 1 and 10 bar). These differences stem from the fact that the 3D time-dependent solution has not yet converged toward a steady-state solution. The 3D system tries to heat the deep interior from the initial 1D profile. This leads to artificial high vertical velocities in the deep atmosphere. Even these velocities, however, are unable to heat up the deep layers because of the extremely long thermal relaxation timescale at such depths. The system then tries to increase its energy and the internal heat flow is too small to efficiently heat up the deep layers on the timescale of the 3D simulation. This is an artificial problem, due to the fact that the initial profile is too cold. In reality, these objects lose energy by cooling, so it should be much easier to reach a steady-state in the 3D runs by initializing with a hotter PT profile and allowing the system to lose the extra energy (Amundsen et al. 2016). Even though the vertical velocities in the 2D deep atmosphere are small, however, the vertical mass flux remains significant, as shown in Figure 3, and this behavior has important consequences, as explained in the next section.

## 5. Solution to the Inflated Radii Problem

As illustrated by the dashed cyan line in Figure 2(B), the 2D solution converges toward an adiabatic profile at significantly smaller pressures than the deeper inner convectively unstable regions (convection is not included in this 2D solution), which starts at $P \sim 40$ bars for HD 209458b with the proper $1.4 R_{jup}$ radius (see Figure 1 of Chabrier et al. 2004). At these levels ($\sim 10$ bars), the atmospheric and interior adiabats connect naturally, leading to a hotter characteristic entropy profile for

the planet than the one found in standard irradiated models. This in turn implies a larger radius for a given mass, $\ge 1.4 R_{jup}$ for HD 209458b. The radius of $1.4 R_{jup}$ corresponds to the adiabat illustrated by the pink dashed line in the P/T plots, which is the one obtained in the evolutionary calculations, including an extra energy equal to 1% of the stellar luminosity in the planet thermal budget (Chabrier et al. 2004) and the dashed blue line corresponds to a radius of $1.6 R_{jup}$. This solves the radius inflation problem: 2D atmospheric circulation leads to an adiabatic structure in the deep atmosphere before convection sets in, leading to a hotter internal adiabat. This behavior can be explained with the simple following model. In deep enough layers (deeper than $\sim 1$ bar), when the gradients of longitudinal velocities are negligible, the system of equations reduces to

$$ \frac{\partial r^2 \rho u_r}{\partial r} = 0 = \frac{\partial p}{\partial \lambda} = 0 = \frac{\partial p}{\partial r} \approx -\rho g \\ -u_r P \frac{\gamma}{\gamma - 1} \frac{\partial \ln(\theta)}{\partial r} = H_{\text{rad}}, \quad (3) $$

which can be rewritten as:

$$ \begin{gathered} r^2 \rho u_r = R_0^2 \rho_0 u_{r,0} = \dot{M}_r \\ \frac{\partial \ln \theta}{\partial \ln P} = \frac{\gamma - 1}{\gamma} \frac{H_{\text{rad}} r^2}{\dot{M}_r g}. \end{gathered} \quad (4) $$

Therefore, in the deep layers, as soon as the heating rate becomes sufficiently small, $H_{\text{rad}} \ll \dot{M}_r g \gamma / (r^2(\gamma - 1))$, $\theta = \text{cst}$, and thus $P \propto T^{\gamma/(\gamma-1)}$, which implies that the profile becomes adiabatic, a direct consequence of the vertical mass flux and the vertical transport of the potential temperature, even in the absence of convection. In turn, layers at pressure $P$ will be at temperature $T_0(P/P_0)^{(\gamma-1)/\gamma}$. From the 2D models, we get an estimate of $H_{\text{rad}} r^2 (\gamma - 1) / \dot{M}_r g \gamma \approx 10^{-2}$ between 10 and 100 bars. We do not show the deepest layers in Figure 2 in order to focus on the comparison with the 3D run; however, our results confirm that the profiles follow the adiabat down to $\sim 2$ kbars. Furthermore, we show in the Figure 3 that the vertical mass flux follows the $1/r^2$ behavior for $P \ge 1$ bar below the hot spot ($\lambda = 36^\circ$ and $72^\circ$), as expected from mass conservation, confirming the validity of the reduced model at these longitudes. For other longitudes, the PT profile is maintained on the same adiabat because of the longitudinal advection of the potential temperature (see Equation (2)). In the absence of radiatively driven vertical velocities, convection sets in at $P \ge 40$ bars (Chabrier et al. 2004), but in the presence of circulation, the transition between the vertical global circulation and the convective zone might take place at deeper levels. In any case, our atmospheric adiabatic profile reconnects with a hotter internal adiabat than the conventional one. At this stage, it is worth stressing that this process differs entirely from the one proposed by Showman & Guillot (2002). Indeed, it is not related to the vertical transport of kinetic energy, a process hard to realize in a strongly gravitationally stratified medium, but is instead related to the vertical advection of the potential temperature (i.e., of entropy) due to a vertical mass flux and a small heating rate. The adiabatic adjustment is similar to the one proposed in Youdin & Mitchell (2010) with turbulent
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Temperature, zonal wind, and vertical wind profiles as a function of pressure for different longitudes $\lambda$. (A): values extracted from a 3D run of the UM (21). (B): 2D stationary circulation model. The pink dashed line in the *P/T* plots portrays the thermal profile obtained in the evolutionary models with extra energy of 1% of the stellar luminosity, giving a radius of 1.4 $R_{\text{jup}}$ (Chabrier et al. 2004), while the dashed cyan line in the *P/T* 2D plot corresponds to the deep atmosphere adiabat that all profiles converge to ($\sim$1.6 $R_{\text{jup}}$).</figcaption></figure>

mixing, but in our 2D solutions the adjustment is simply induced by the vertical velocities of the global circulation itself. This process is very similar to the one taking place in the Earth's stratified atmosphere or in deep oceanic layers in cases of strong surface irradiation. The large-scale non-uniform heating induces longitudinal motions (by conservation of longitudinal momentum), which in turn produce large air (in the atmosphere) or water (beneath the sea surface) mass fluxes that transport adiabatically the potential temperature to

deep layers and thereby warm them up (see chapter 16 in Vallis 2006 and chapter 9.15 in Gill 1982).

We have verified that this process and the reduced system are independent of the value of $\alpha$ in the deep layers. Tests with $\alpha = \infty$, deeper than 1 bar, and $\alpha = 2$, in the upper layers, yields the same conclusions. Hence, our results are independent of the degree of convergence of the deep atmosphere in the 3D model. In the upper layers, the magnitude of the vertical wind depends directly on the value
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Vertical mass fluxes as a function of pressure for different longitudes at the equator in the 2D steady-state model. The dashed magenta line represents the scaling of the mass fluxes over inverse squared radius when vertical motions are dominant in the mass conservation equation.</figcaption></figure>

of $\alpha$ used at these levels. Hence $\alpha$ is also directly constrained by the high-altitude vertical wind of the 3D model, which is well-converged for $P \le 10^{-1}$ bar. Only when $\alpha$ is high everywhere in the atmosphere do the vertical velocities vanish, and we get back to the standard flat isothermal profile in the deep atmosphere (see the Appendix).

The 2D model applies only to the equator of the planet. Since the longitudinal advection of the potential temperature maintains the 2D PT structure on the adiabat reached under the hot spot, we can expect that the meridional advection of potential temperature (see Equation (2)) maintains the 3D PT structure on the same adiabat. Other latitudes of the full 3D solution in the UM simulation indeed seem to converge to the same adiabat as the 2D model, most likely because of the deep circulation identified in previous simulations (Heng et al. 2011; Mayne et al. 2014). Since the vertical mass flux is mainly downwelling in our equatorial solution, we expect the full 3D steady-state to produce upwelling flux at higher latitudes to globally conserve mass. We illustrate this possibility in Figure 4 by displaying the upward vertical velocities at the high latitudes around the hot spot and by computing the deep vertical velocities with the vertical mass conservation condition. Current computational capabilities, however, prevent the possibility for 3D models to probe the long timescales needed for the circulation to reach this steady-state down to the convection zone. The combination of both 3D runs and our 2D stationary circulation model thus enables us to better understand the short-term and long-term dynamics of hot jupiters. It demonstrates that irradiation-induced atmospheric circulation explains the inflated radii of hot jupiters.

Since the 2D solution reaches high temperatures in the deep atmosphere, the electric conductivity is high and ohmic dissipation might affect our results. We have estimated the efficiency of ohmic dissipation with Equation (17) of Batygin et al. (2011) and we get an efficiency of the order of $10^{-5}$ below ~8 bars. The reason is that even though the electric conductivity is high, the wind solution from the 3D simulation has almost zero velocity for $P \ge 8$ bars, and hence no kinetic energy to be dissipated. As a consequence, ohmic dissipation should not impact our result, as also found in other studies (Rogers & Komacek 2014; Rogers & Showman 2014), even though future (lacking) constraints on the magnitude of the

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** (A) For pressures less than 1 bar, the vertical wind profiles are extracted at different latitudes between 0° and 72° of the 3D UM simulation (at the longitude of the hot spot $\lambda = 72°$). At higher pressures the steady-state vertical velocities are approximated by integration of the vertical mass conservation. (B) shows the corresponding vertical mass fluxes.</figcaption></figure>

deep circulation in the atmosphere will help to precisely quantify this effect.

## 6. Correlation of the Inflated Radii with Increasing Irradiation Flux

There is now growing observational evidence suggesting a correlation between inflated radii and incident stellar flux (Demory & Seager 2011; Laughlin et al. 2011; Miller & Fortney 2011; Weiss et al. 2013; Figueira et al. 2014). Laughlin et al. (2011) have found a best-fit dependence of $R_{\text{obs}} - R_{\text{pred}} \propto T_{\text{eff}}^{\beta}$ with $\beta = 1.4 \pm 0.6$ ($T_{\text{eff}}$ is the equilibrium temperature of the planet defined as $T_{\text{star}} \times (R_{\text{star}}/a)^{1/2}$). In Figure 5, we show the PT profiles obtained with the 1D and 2D models as a function of irradiation flux. All the other parameters are kept constant. We keep the same wind profile used for HD 209458b (see Figure 1), although the deep adiabat does not depend much on the wind magnitude (see Figure 7 in the Appendix). It is clear that the stronger the irradiation, the hotter the internal adiabat the averaged 2D profiles connect with. We can easily estimate the magnitude of the temperature change at these levels (~200 bars) compared with the 1D models. For HD 209458b (brown profiles) the profile becomes isothermal at about 1800 K, ~1 bar in 1D, whereas in 2D it becomes adiabatic at ~5 bars, with an adiabatic gradient of $(\gamma - 1)/\gamma \approx 0.21$ ($\gamma \approx 1.27$). The temperature at 200 bars is thus now $1800 \times (200/5)^{0.21} \approx 4000$ K. Using evolutionary
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Average PT profiles of the 2D circulation model as a function of irradiation flux (solid lines); the standard irradiation for HD209 is $9.41 \times 10^8$ erg s$^{-1}$ cm$^{-2}$. For comparison, we show the PT profiles obtained in 1D with dashed lines.</figcaption></figure>

models (Chabrier et al. 2004), we can calculate the interior profiles and the corresponding radii that connect to the different atmospheric structures at the interior-atmosphere boundary. Figure 6 displays the observed correlation between these inflated radii and the irradiation flux (red lines), compared with the ones obtained in 1D (brown). The systems displayed in the plot are taken from www.exoplanet.eu. Note that we focus here on irradiated hot jupiters, i.e., objects with a dominant gaseous envelope with a lower-limit radius of 1 $R_{\text{jup}}$ in the absence of irradiation. Planets with smaller radii are less massive, e.g., hot Neptunes, and have a smaller gaseous atmosphere. Radius inflation for such objects is not an issue and can be easily explained by properly taking into account the irradiating flux (e.g., Baraffe et al. 2008, 2010). For the sake of simplicity and to illustrate our purpose, we keep the same gravity in the atmospheric models, log (g) = 2.97, as a fiducial value, while varying the incident flux. In the evolutionary models, we use a mass of 0.7 $M_{\text{jup}}$, the mass of HD 209458b, except for the most irradiated object, for which we take a mass of 2 $M_{\text{jup}}$, so all models are consistent with log (g) = 2.97 ± 0.15. As shown in the figure (and in Table 2) the 1D models hardly exceed 1.05 $R_{\text{jup}}$. In contrast, the 2D models, all based on the same physical setup, reproduce the whole trend of increasing radius with increasing irradiation, including for the most inflated objects, with radii ~2 $R_{\text{jup}}$. Overall, the radii obtained in 2D seem to be a bit larger than the observed ones, a result for which we see different possible explanations:

1. We used the same input profiles as the one obtained in 3D for HD 209458b. Variations in the input parameters might yield slightly smaller radii.

2. The steady-state 3D solutions might converge to intermediate profiles between the 1D and 2D ones. This is indeed suggested by our 3D simulation, which seems to converge toward the 1.4 $R_{\text{jup}}$ (the observed value) interior adiabat, whereas the 2D solution predicts a slightly larger value, 1.6 $R_{\text{jup}}$.

3. Dissipation processes in the atmosphere (e.g., shear, ohmic dissipation that could reduce the magnitude of the vertical velocities) might lead to a slightly cooler adiabat and thus a smaller radius.

Nevertheless, the 2D steady-state circulation models well reproduce the overall observed trend of enhanced radii of hot

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Evolution of the radius as a function of irradiation for the 2D models (red) and for the 1D irradiated models (brown). All the atmospheric models are calculated with log (g) = 2.97. All the evolutionary models are calculated for an object mass of 0.7 $M_{\text{jup}}$, except for the most irradiated one, for which we use a mass of 2 $M_{\text{jup}}$ in order to keep log (g) = 2.97 ± 0.15.</figcaption></figure>

<table><caption>**Table 2**
1D and 2D Radii for a Planet with log (g) = 2.97 for Different Irradiation Fluxes</caption><thead><tr><th>Flux<br>(erg s<sup>-1</sup> cm<sup>-2</sup>)</th><th>1D Radius<br>(R<sub>jup</sub>)</th><th>2D Radius<br>(R<sub>jup</sub>)</th></tr></thead><tbody><tr><td>1.47 × 10<sup>7</sup></td><td>1.05</td><td>1.15</td></tr><tr><td>5.88 × 10<sup>7</sup></td><td>1.05</td><td>1.28</td></tr><tr><td>2.35 × 10<sup>8</sup></td><td>1.05</td><td>1.36</td></tr><tr><td>9.41 × 10<sup>8</sup></td><td>1.05</td><td>1.60</td></tr><tr><td>3.77 × 10<sup>9</sup></td><td>1.1</td><td>1.82</td></tr></tbody></table>

jupiters as a function of irradiation, based on a robust mechanism, namely atmospheric circulation induced by irradiation. The relation we obtain is $R_{2D}-R_{1D} \propto T_{\text{eff}}^{\beta}$ with $\beta = 1.43 \pm 0.14$, in good agreement with the aforementioned one inferred from observations.

## 7. Discussion and Conclusions

In this paper, we have presented a new approach to the problem of atmospheric circulation in tidally locked, strongly irradiated exoplanets. We have constructed a 2D stationary circulation model that is able to predict the pressure/temperature and wind structure of the atmosphere at the equator, given a zonal wind profile at a given longitude and an estimate of the meridional to vertical mass flux ratio. Since the 2D steady-state solution agrees well with the 3D simulation in the upper atmosphere, observational evidence of this mechanism remains relatively elusive at this stage. Nevertheless, it is worth mentioning that the steady-state solution does not have the isothermal plateau characteristic of the 1D solution, nor does it present cold traps as in Parmentier et al. (2016). The absence of these features might provide indirect observational constrains.

With the increasing computational power of massively parallel architectures, a full 3D steady-state might be reached sufficiently deep to better characterize the nature of the flow in the deep atmosphere. Nonetheless, the steady-state nature of our 2D circulation model provides a very useful complementary tool to current 3D simulations to understand the long-term physical processes at play in the atmospheres of hot jupiters. With this model we have demonstrated that atmospheric circulation induced
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Pressure–temperature profiles of the 2D circulation model with constant profiles $\alpha = 2$ and $u_{\lambda} = 6$ km s$^{-1}$ (A), $\alpha = 100$ and $u_{\lambda} = 6$ km s$^{-1}$ (B), $\alpha = 2$ and $u_{\lambda} = 2$ km s$^{-1}$ (C), $\alpha = 2$ and $u_{\lambda} = 8$ km s$^{-1}$ (D). The magenta dashed line is the same as in Figure 2.</figcaption></figure>

by irradiation can explain the inflated radii of these planets. Indeed, the induced vertical mass flux and the advection of the potential temperature naturally constrain the atmosphere pressure/temperature profile to become adiabatic, even in the absence of convection, when the heating rate is small, as is the case in the deep atmosphere. Since this structure becomes adiabatic at smaller pressures than the one at which the atmosphere becomes unstable to convection, the 2D steady-state atmospheric profile reconnects with a hotter internal adiabat. This in turn implies a larger radius for a given mass, explaining the long-standing problem of the anomalously inflated radii of these objects.

N.M. acknowledges funding from the Leverhulme Trust via a Research Project Grant. James Manners acknowledges the support of a Met Office Academic Partnership secondment. This work is partly funded by the ERC grants No. 320478-TOFU and No. 247060-PEPS. Some of the calculations used the STFC DIRAC HPC service and the new University of Exeter Supercomputer, ISCA. D.S.A. acknowledges support from the NASA Astrobiology Program through the Nexus for Exoplanet System Science.

## Appendix

The stationary Euler equations in spherical coordinates $r$ (radius), $\phi$ (latitude), and $\lambda$ (longitude), supplemented by the energy equation for an ideal gas (whose internal energy $e$ is

given by the equation of state $P = (1 - \gamma)e$) read,

$$ \frac{1}{r^2} \frac{\partial^2 \rho u_r}{\partial r} + \frac{1}{r} \frac{\partial \cos \phi \rho u_\phi}{\partial \phi} + \frac{1}{r \cos \phi} \frac{\partial \rho u_\lambda}{\partial \lambda} = 0 $$

$$ u_r \frac{\partial u_r}{\partial r} + \frac{u_\phi}{r} \frac{\partial u_r}{\partial \phi} + \frac{u_\lambda}{r \cos \phi} \frac{\partial u_r}{\partial \lambda} - \frac{u_\phi^2 + u_\lambda^2}{r} + \frac{1}{\rho} \frac{\partial p}{\partial r} + g(r) - 2\Omega u_\lambda \cos \phi = 0 $$

$$ u_r \frac{\partial u_\phi}{\partial r} + \frac{u_\phi}{r} \frac{\partial u_\phi}{\partial \phi} + \frac{u_\lambda}{r \cos \phi} \frac{\partial u_\phi}{\partial \lambda} + \frac{u_r u_\phi}{r} + \frac{u_\lambda^2 \tan \phi}{r} + \frac{1}{\rho r} \frac{\partial p}{\partial \phi} + 2\Omega u_\lambda \sin \phi = 0 $$

$$ u_r \frac{\partial u_\lambda}{\partial r} + \frac{u_\phi}{r} \frac{\partial u_\lambda}{\partial \phi} + \frac{u_\lambda}{r \cos \phi} \frac{\partial u_\lambda}{\partial \lambda} + \frac{u_r u_\lambda}{r} - \frac{u_\lambda u_\phi \tan \phi}{r} + \frac{1}{\rho r \cos \phi} \frac{\partial p}{\partial \lambda} - 2\Omega u_\phi \sin \phi + 2\Omega u_r \cos \phi = 0 $$

$$ \frac{u_r}{\gamma - 1} \frac{\partial p}{\partial r} + \frac{u_\phi}{\gamma - 1} \frac{1}{r \cos \phi} \frac{\partial p}{\partial \phi} + \frac{u_\lambda}{\gamma - 1} \frac{1}{r \cos \phi} \frac{\partial p}{\partial \lambda} + \frac{\gamma p}{\gamma - 1} \left( \frac{1}{r^2} \frac{\partial^2 u_r}{\partial r} + \frac{1}{r \cos \phi} \frac{\partial \cos \phi u_\phi}{\partial \phi} + \frac{1}{r \cos \phi} \frac{\partial u_\lambda}{\partial \lambda} \right) = H_{\text{rad}} $$

(5)
where the velocity components are denoted by $u_r$, $u_\phi$, and $u_\lambda$, density is $\rho$, pressure is $P$, $\Omega$ is the angular velocity, and $g(r)$ is the gravitational acceleration. Since our radiative solver is plane-parallel, the heating rate at the equator is approximated as $H_{\text{rad}}(r, \cos\phi, \cos\lambda) = -(\partial r^2 F_{\text{rad}}/\partial r)/r^2 \approx -(\partial F_{\text{rad}}(r, \cos\lambda))/\partial r$ for the conservation of energy. At the equator $\phi = 0$, and assuming a constant residual meridional to vertical mass flux ratio $\alpha$:

$$ \frac{1}{r^2} \frac{\partial r^2 \rho u_r}{\partial r} = \frac{1}{r\alpha} \frac{\partial \rho u_\phi}{\partial \phi}. \qquad (6) $$

We get the reduced system:

$$ \begin{aligned} & \frac{1}{r} \frac{\partial^2 u_r}{\partial r} + r \frac{u_r}{\rho} \frac{\partial \rho}{\partial r} + \frac{\partial u_\phi}{\partial \phi} + \frac{\partial u_\lambda}{\partial \lambda} + \frac{u_\lambda}{\rho} \frac{\partial \rho}{\partial \lambda} = 0 \\ & \frac{1}{r} \frac{\partial^2 u_r}{\partial r} + r \frac{u_r}{\rho} \frac{\partial \rho}{\partial r} - \frac{1}{\alpha} \frac{\partial u_\phi}{\partial \phi} = 0 \\ & ru_r \frac{\partial u_\lambda}{\partial r} + u_\lambda \frac{\partial u_\lambda}{\partial \lambda} + u_r u_\lambda + 2\Omega r u_r + \frac{1}{\rho} \frac{\partial p}{\partial \lambda} = 0 \\ & u_r \frac{\partial u_r}{\partial r} + \frac{u_\lambda}{r} \frac{\partial u_r}{\partial \lambda} - \frac{u_\lambda^2}{r} - 2\Omega u_\lambda + \frac{1}{\rho} \frac{\partial p}{\partial r} + g(r) = 0 \\ & \frac{ru_r}{\gamma - 1} \frac{\partial p}{\partial r} + \frac{u_\lambda}{\gamma - 1} \frac{\partial p}{\partial \lambda} \\ & \qquad - \frac{\gamma p}{\gamma - 1} \left( r \frac{u_r}{\rho} \frac{\partial \rho}{\partial r} + \frac{u_\lambda}{\rho} \frac{\partial \rho}{\partial \lambda} \right) + r \frac{\partial F_{\text{rad}}}{\partial r} = 0. \end{aligned} \qquad (7) $$

All the details of our radiative scheme with irradiation and isotropic scattering are given in Drummond et al. (2016). For the advection terms of the type $u_x \partial(\rho)/\partial x$, we used an upwind scheme. For stability reasons, we also needed to use an upwind (relative to the radiative flux transport, i.e., from bottom to top in the atmosphere) heating rate of the type $H_{\text{idepth}} = -(F_{\text{idepth+1}} - F_{\text{idepth}})/(r_{\text{idepth+1}} - r_{\text{idepth}})$. For axisymmetric irradiation (or absence of irradiation), an axisymmetric solution exists with $\partial u_\phi/\partial \phi = 0$, $u_r = 0$, $u_\lambda$ constant, and any 1D solution that fulfills the hydrostatic balance and the constant flux transport. Then, the system as proposed is not able to fully specify the wind solution, since any constant zonal wind is a solution. Therefore, we impose a zonal wind profile at a given longitude of the equator (the substellar point), which is equivalent to imposing the total momentum in the zonal wind.

In Figure 7, we show the 2D PT structures with constant profiles at $\alpha = 2$ and $u_\lambda = 6 \text{ km s}^{-1}$ (A), $\alpha = 100$ and $u_\lambda = 6 \text{ km s}^{-1}$ (B), $\alpha = 2$ and $u_\lambda = 2 \text{ km s}^{-1}$ (C), $\alpha = 2$ and $u_\lambda = 8 \text{ km s}^{-1}$ (D). With $\alpha = 100$ throughout the whole

atmosphere, we recover the usual isothermal plateau (at $T \approx 1800$ K), similar to the 1D solutions (e.g., Drummond et al. 2016). The occurrence of an isothermal plateau is indeed expected, as $\alpha = 100$ corresponds to a case for which the mass flux would be essentially purely meridional, with a negligible vertical contribution, a rather unlikely solution. Indeed, as seen in Figure 1, the 3D run has already converged to an adiabatic structure at higher temperatures ($T \approx 2000$ K), for a pressure $P \le 10$ bars. This suggests that the full 3D solution converges to the same adiabat as the 2D steady-state circulation model. The PT structures with different wind velocities at the substellar point show that the deep adiabat is only weakly dependent on this parameter.

## References

* Amundsen, D. S., Baraffe, I., Tremblin, P., et al. 2014, A&A, 564, 59

* Amundsen, D. S., Mayne, N. J., Baraffe, I., et al. 2016, A&A, 595, A36

* Amundsen, D. S., Tremblin, P., Manners, J., Baraffe, I., & Mayne, N. J. 2017, A&A, 598, A97

* Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315

* Baraffe, I., Chabrier, G., & Barman, T. 2010, RPh, 73, 016901

* Baraffe, I., Chabrier, G., Fortney, J., & Sotin, C. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 763

* Batygin, K., Stevenson, D. J., & Bodenheimer, P. H. 2011, ApJ, 738, 1

* Burrows, A., Hubeny, I., Budaj, J., & Hubbard, W. B. 2007, ApJ, 661, 502

* Chabrier, G., & Baraffe, I. 2007, ApJL, 661, L81

* Chabrier, G., Barman, T., Baraffe, I., Allard, F., & Hauschildt, P. H. 2004, ApJL, 603, L53

* Demory, B.-O., & Seager, S. 2011, ApJS, 197, 12

* Dobbs-Dixon, I., Agol, E., & Burrows, A. 2012, ApJ, 751, 87

* Drummond, B., Tremblin, P., Baraffe, I., et al. 2016, A&A, 594, A69

* Figueira, P., Oshagh, M., Adibekyan, V. Z., & Santos, N. C. 2014, A&A, 572, A51

* Fortney, J. J., & Nettelmann, N. 2010, SSRv, 152, 423

* Fromang, S., Leconte, J., & Heng, K. 2016, A&A, 591, A144

* Gill, A. E. 1982, JFM, 157, 532

* Grunblatt, S. K., Gaidos, E., Huber, D., & Howard, A. 2016, AAS, 228, 201.04

* Heng, K., Menou, K., & Phillipps, P. J. 2011, MNRAS, 413, 2380

* Laughlin, G., Crismani, M., & Adams, F. C. 2011, ApJL, 729, L7

* Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014, A&A, 561, 1

* Miller, N., & Fortney, J. J. 2011, ApJL, 736, L29

* Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, ApJ, 828, 22

* Rogers, T. M., & Komacek, T. D. 2014, ApJ, 794, 132

* Rogers, T. M., & Showman, A. P. 2014, ApJL, 782, L4

* Showman, A. P., Cho, J. Y. K., & Menou, K. 2010, in Exoplanets, ed. S. Seager (Tucson, AZ: Univ. Arizona Press), 471

* Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564

* Showman, A. P., & Guillot, T. 2002, A&A, 385, 166

* Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71

* Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJL, 817, L19

* Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJL, 804, L17

* Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics (Cambridge: Cambridge Univ. Press)

* Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14

* Youdin, A. N., & Mitchell, J. L. 2010, ApJ, 721, 1113