Idealised simulations of the deep atmosphere of hot Jupiters

Context. The anomalously large radii of hot Jupiters has long been a mystery. However, by combining both theoretical arguments and 2D models, a recent study has suggested that the vertical advection of potential temperature leads to a hotter adiabatic temperature proﬁle in the deep atmosphere than the proﬁle obtained with standard 1D models. Aims. In order to conﬁrm the viability of that scenario, we extend this investigation to 3D, time-dependent models. Methods. We use a 3D general circulation model DYNAMICO to perform a series of calculations designed to explore the formation and structure of the driving atmospheric circulations, and detail how it responds to changes in both the upper and deep atmospheric forcing. Results. In agreement with the previous, 2D study, we ﬁnd that a hot adiabat is the natural outcome of the long-term evolution of the deep atmosphere. Integration times of the order of 1500 yr are needed for that adiabat to emerge from an isothermal atmosphere, explaining why it has not been found in previous hot Jupiter studies. Models initialised from a hotter deep atmosphere tend to evolve faster toward the same ﬁnal state. We also ﬁnd that the deep adiabat is stable against low-levels of deep heating and cooling, as long as the Newtonian cooling timescale is longer than ∼ 3000 yr at 200 bar. Conclusions. We conclude that steady-state vertical advection of potential temperature by deep atmospheric circulations constitutes a robust mechanism to explain the inﬂated radii of hot Jupiters. We suggest that future models of hot Jupiters be evolved for a longer time than currently done, and when possible that models initialised with a hot deep adiabat be included. We stress that this mechanism stems from the advection of entropy by irradiation-induced mass ﬂows and does not require a (ﬁnely tuned) dissipative process, in contrast with most previously suggested scenarios.


Introduction
The anomalously large radii of highly irradiated Jupiter-like exoplanets known as hot Jupiters remain one of the key unresolved issues in our understanding of extrasolar planetary atmospheres. The observed correlation between the stellar irradiation of a hot Jupiter and its observed inflation (for examples, see Demory & Seager 2011;Laughlin et al. 2011;Lopez & Fortney 2016;Sestovic et al. 2018) suggests that it is linked to the amount of energy deposited in the upper atmosphere. Several mechanisms have been suggested as possible explanations (see Baraffe et al. 2009Baraffe et al. , 2014Fortney & Nettelmann 2010, for a review). These solutions include tidal heating and physical (i.e. not for stabilisation reasons) dissipation (Leconte et al. 2010;Arras & Socrates 2010;Lee 2019), ohmic dissipation of electrical energy (Batygin & Stevenson 2010;Perna et al. 2010;Batygin et al. 2011;Rauscher & Menou 2012), deep deposition of kinetic energy , enhanced opacities which inhibit cooling (Burrows et al. 2007) or ongoing layered convection that reduces the efficiency of heat transport (Chabrier & Baraffe 2007). At the present time however, there is no consensus across the community on a given scenario because the majority of these solutions require finely tuned physical environments which either deposit additional energy deep within the atmosphere or affect the efficiency of vertical heat transport.
Recently, Tremblin et al. (2017, hereafter PT17) suggested a mechanism that naturally arises from first physical principles. Their argument goes as follows: consider the equation for the evolution of the potential temperature Θ, which is equivalent to entropy in this case: where D/Dt stands for the Lagrangian derivative in spherical coordinates, H is the local heating or cooling rate, c p is the heat capacity at constant pressure, and Θ is defined as a function of the temperature T and pressure, P: where P 0 is a reference pressure and γ = C p /C v is the adiabatic index. In a steady state, Eq. (1) reduces to where u is the velocity. In the deep atmosphere, radiative heating and cooling both tend to zero (i.e. H → 0) because of large atmospheric opacities. In this case (with H → 0), and if the winds do not vanish (i.e. |u| 0, see Sect. 3.2), the potential temperature Θ must remain constant for Eq.
(3) to be valid. In other words, the temperature-pressure profile must be adiabatic and satisfy the scaling: We emphasise that this adiabatic solution is an equilibrium that does not require any physical dissipation. There is an internal energy transfer to the deep atmosphere through an enthalpy flux, but there is no dissipation from kinetic, magnetic, or radiative energy reservoirs to the internal energy reservoir. Dissipative processes D dis would act as a source term with u · ∇Θ ∝ D dis and would drive the profile away from the adiabat. Physically, as discussed by PT17, this constant potential temperature profile in the deep atmosphere is driven by the vertical advection of potential temperature from the outer and highly irradiated atmosphere to the deep atmosphere by large-scale dynamical motions where it is almost completely homogenised by the residual global circulations (which themselves can be linked to the conservation of mass and momentum, and the large mass/momentum flux the super-rotating jet drives in the outer atmosphere). The key point is that it causes the temperaturepressure profile to converge to an adiabat at lower pressures than those at which the atmosphere becomes unstable to convection. As a result, the outer atmosphere connects to a hotter internal adiabat than would be obtained through a standard, "radiativeconvective" single column model. This potentially leads to a larger radius compared with the predictions born out of these 1D models.
Whilst PT17 was able to confirm this hypothesis through the use of a 2D stationary circulation model, there are still a number of limitations to their work. Perhaps most importantly, the models they used only considered the formation of the deep adiabat within a 2D equatorial slice. The steady-state temperature-pressure profiles at other latitudes remain unknown, as well as the nature of the global circulations at these high pressures in the equilibrated state. Strong ansatzes were also made about the nature of the meridional (i.e. vertical and latitudinal) wind at the equator, with the 2D models of PT17 prescribing the ratio of latitudinal to vertical mass fluxes, a presciption which could potentially affect the proposed scenario. The purpose of this paper is to reduce and constrain these assumptions and limitations and to demonstrate the viability of a deep adiabat at equilibrium. This is done by means of a series of idealised 3D general circulation model (GCM) calculations designed so as to allow us to fully explore the structure of the deep atmospheric circulations in the equilibrated atmospheres of hot Jupiters, as well as investigate the time-evolution of the deep adiabat. As we demonstrate in this work, the adiabatic profile predicted by PT17 naturally emerges from such calculations and appears to be robust against changes in the deep atmosphere radiative properties. This is the core result of this work.
The structure of the work is as follows. Our simulation properties are described in Sect. 2, where we introduce the GCM DYNAMICO used throughout this study. We then demonstrate that when using DYNAMICO, not only are we able to recover standard features observed in previous short-timescale studies of the atmospheres of hot Jupiters (Sect. 3.1), but also, when the simulations are extended to long-enough time-scales, we can show that an adiabatic profile develops within the deep atmosphere (Sect. 3.2). We then explore the robustness of our results by presenting a series of sensitivity tests, including changes in the thermal forcing in the outer and deep atmosphere (Sect. 3.3). Finally, in Sect. 4, we provide concluding remarks, including suggestions for future computational studies of the atmospheres of hot Jupiter and a discussion about implications for the evolution of highly irradiated gas giants.

Method
DYNAMICO is a highly computationally efficient GCM that solves the primitive equation of meteorology (see Vallis 2006 for a review and Dubos & Voitus 2014 for a more detailed discussion of the approach taken in DYNAMICO) on a sphere (Dubos et al. 2015). It is being developed as the next state-of-the-art dynamical core for Earth and planetary climate studies at the Laboratoire de Météorologie Dynamique and is publicly available 1 and has recently been used to model the atmosphere of Saturn at high resolution (Spiga et al. 2020). Here, we present some specificities of DYNAMICO (Sect. 2.1) as well as the required modifications we implemented to model the atmospheres of hot Jupiters (Sect. 2.2).

DYNAMICOs numerical scheme
Briefly, DYNAMICO takes an energy-conserving Hamiltonian approach to solving the primitive equations. This has been shown to be suitable for modelling the atmospheres of hot Jupiters (Showman et al. 2008;Rauscher & Menou 2012), although it may not be valid in other planetary atmospheres . Rather than the traditional latitude-longitude horizontal grid (which presents numerical issues near the poles due to singularities in the coordinate system -see the review of Williamson 2007 for more details), DYNAMICO uses a staggered horizontal-icosahedral grid (see Thuburn et al. 2014 for a discussion of the relative numerical accuracy for this type of grids) for which the number of horizontal cells N is defined by the number of subdivisions d of each edge of the main spherical icosahedral 2 : As for the vertical grid, DYNAMICO uses a pressure coordinate system whose levels can be defined by the user at runtime. Finally, the boundaries of our simulations are closed and stressfree with zero-energy transfer (i.e. the only means of energy injection and removal are the Newtonian cooling profile and the horizontal numerical dissipation). We note that unlike some other GCM models of hot Jupiters (e.g. Schneider & Liu 2009;Liu & Showman 2013;Showman et al. 2019), we do not include an additional frictional (i.e. Rayleigh) drag scheme at the bottom of our simulation domain, instead relying on the hyperviscosity and impermeable bottom boundary to stabilise the system. As a consequence of the finite difference scheme used in DYNAMICO, artificial numerical dissipation must be introduced in order to stabilise the system against the accumulation of grid-scale numerical noise. This numerical dissipation takes the form of a horizontal hyper-diffusion filter with a fixed hyperviscosity and a dissipation timescale at the grid scale, labelled τ dissip , which serves to adjust the strength of the filtering (the longer the dissipation time, the weaker the dissipation). Technically, DYNAMICO includes three dissipation timescales, each of which either diffuses scalar, vorticity, or divergence independently. However, for our models, we set all three timescales to the same value. It is important to point out that the hyperviscosity is not a direct equivalent of the physical viscosity of the planetary atmosphere, but can be viewed as a form of increased artificial dissipation that both enhances the stability of the code and accounts for motions, flows, and turbulences which are unresolved at typical grid-scale resolutions. This is known as the large eddy approximation and has long been standard practice in the stellar (e.g. Miesch 2005) and planetary (e.g Cullen & Brown 2009) atmospheric modelling communities. Because it acts at the grid-cell level, the strength of the dissipation is resolutiondependent at a fixed τ dissip (this can be seen in our results in Fig. 7).
In a series of benchmark cases, Heng et al. (2011, hereafter H11) have shown that both spectral and finite-difference-based dynamical cores which implement horizontal hyper-diffusion filters can produce differences of the order of tens of percent in the temperature and velocity fields when varying the dissipation strength. We also found such a similar sensitivity in our models: for example, the maximum super-rotating jet speed varies between 3000 and 4500 ms −1 as the dissipation strength is varied. The dissipation strength must thus be carefully calibrated. In the absence of significant constraints on hot Jupiter zonal wind velocities, this was done empirically by minimising unwanted small-scale numerical noise as well as replicating published benchmark results (an alternative, which is especially useful in scenarios where direct or indirect data comparisons are unavailable, is to plot the spectral decomposition of the energy profile and adjust the diffusion such that the energy accumulation on the smallest scales is insignificant). We found that setting τ dissip = 2500 s in our low-resolution runs leads to benchmark cases in good agreement with the results of, for example, Mayne et al. (2014a), whilst also exhibiting minimal small-scale numerical noise. This is in reasonable agreement with other studies, with our models including a hyper-diffusion of the same order of magnitude as for example H11. We note that due to differences in the dynamics between those of Saturn and that observed in hot Jupiters, and in particular due to the presence of the strong superrotating jet, we must use a significantly stronger dissipation to counter grid-scale noise than that used in previous atmospheric studies calculated using DYNAMICO (Spiga et al. 2020).

Newtonian cooling
In our simulations of hot-Jupiter atmospheres using DYNAM-ICO, we do not directly model either the incident thermal radiation on the day-side or the thermal emission on the night-side of the exoplanet. This would be prohibitively computationally expensive for the long simulations we perform in the present work. Instead we use a simple thermal relaxation scheme to model those effects, with a spatially varying equilibrium temperature profile T eq and a relaxation timescale τ that increases with pressure throughout the outer atmosphere. Specifically, this is done by adding a source term to the temperature evolution equation that takes the form This method, known as Newtonian Cooling, has long been applied within the 3D GCM exoplanetary community (i.e. , Showman et al. 2008, Showman & Polvani 2011, Mayne et al. 2014a,b, or Guerlet et al. 2014), although it is gradually being replaced by coupling with simplified but more computationally expensive radiative transfer schemes (e.g. Showman et al. 2009;Rauscher & Menou 2012, or Amundsen et al. 2016 due to its limitations (e.g. it is difficult to use to probe individual emission or absorption features such as non-equilibrium atmospheric chemistry or stellar activity). The forcing temperature and cooling timescale we use within our models have their basis in the profiles Iro et al. (2005) calculated via a series of 1D radiative-transfer models. These models were then parametrised by H11, who created simplified day-side and night-side profiles. The parametrisation used here is based on this work, albeit modified in the deep atmosphere since this is the focus of our analysis. As a result, it somewhat resembles a parametrised version of the cooling profile considered by Liu & Showman (2013).
For P > 10 bar, we set T eq = T night = T day = 1800 K.
For P > 10 bar, we consider a series of models that lie between two extremes: at one extreme we set log(τ 220 ), which we define as the decimal logarithm of the cooling timescale τ at the bottom of our model atmospheres (i.e. at P = 220 bar), to infinity, which implies that the deep atmosphere is radiatively inert, with no heating or cooling. As for the other extreme, this involves setting log(τ 220 ) = 7.5, which implies that radiative effects do not diminish below 10 bar. In Sect. 3 we explore results at the first extreme, with no deep radiative dynamics. Subsequently, in Sect. 3.3.3 we explore the sensitivity of our results to varying this prescription.

Results
The default parameters used with our models are outlined in Table 1, with the resultant models, as well as the simulation specific parameters, detailed in Table 2. In Sect. 3.2, we use the results of models A and B to demonstrate the validity of the work of PT17 in the time-dependent, three-dimensional regime. We next explore the robustness and sensitivity of our results to numerical and external effects in Sect. 3.3. We note that throughout this paper, all times are either given in seconds or in Earth years -specifically one Earth year is exactly 365 days.

Validation of the hot Jupiter model
We start by exploring the early evolution of model A, testing how well it agrees with the benchmark calculations of H11. The model is run for an initial period of 30 yr in order to reach an evolved state before we take averages over the next five years of data. We note that this model was also used to calibrate the horizontal dissipation (τ dissip ). In Fig. 1, we show zonally and temporally averaged plots of the zonal wind and the temperature as a function of both latitude and pressure.
We find that the temperature (left panel) is qualitatively similar to that reported by both H11 and Mayne et al. (2014a). The temperature range we find (∼750 K → ∼2150 K) matches the results of these latter authors (∼700 → 2000 K) to within a 10% margin of uncertainty. This is satisfactory given the differences between the various setups and numerical implementations of the GCMs, as well as the variations that occur when adjusting the length of the temporal averaging window.
The zonal wind displays a prominent, eastward, superrotating equatorial jet that extends from the top of the atmosphere down to approximately 10 bar (we note that as we continue to run this model for more time, the vertical extent of the jet increases, eventually reaching significantly deeper than 100 bar after 1700 yr). It exhibits a peak wind velocity of ≈3500 m s −1 , depending on the averaging window considered, in good agreement with the work of both H11 and Mayne et al. (2014a) who found peak jet speeds on the order of 3500 → 4000 m s −1 . In the upper atmosphere, it is balanced by counter-rotating (westward) flows at extratropical and polar latitudes. The zonal wind is also directed westwards at all latitudes below ∼50 bar, with this wind also contributing to the flows balancing the large mass and momentum transport of the super-rotating jet.
Initial temperature 1800 Table 2. Models discussed in this work.

A
The base low-resolution model, in which the deep atmosphere is isothermally initialised

J & K
Highly evolved versions of model A which have reached a deep adiabat, and then had their outer atmospheric Newtonian cooling modified to reflect a different surface temperature: 1200 K in model J and 2200 K in model K The differences we find between our models and the reference models are not unexpected. As discussed by H11, the jet speed and temperature profile are indeed highly sensitive not only to the numerical scheme adopted by the GCM (i.e. spectral vs. finite difference -see their Fig. 12) but also to the form and magnitude of horizontal dissipation and Newtonian cooling used. In our models, unlike H11, we explicitly set our deep (P > 10 bar) cooling to zero, which may explain the enhanced deep temperatures observed in our models, most likely an early manifestation of the deep adiabat we expect to eventually develop.
As noted by other works (e.g. Menou & Rauscher 2009;Mayne et al. 2014a), it takes a long time for the deep atmosphere to reach equilibrium, and the above simulation is by no means an exception: the eastward equatorial jet extends deeper and deeper as time increases, with no sign of stopping by the end of the simulated duration. This long-timescale evolution is explored in detail in the following section.
Equatorially averaged (i.e. the zonal-mean at the equator) T -P profiles (orange) for two evolved models that were either isothermally (a) or adiabatically (b) initialised. In both cases there is no forcing below 10 bars (i.e. when P > 10 bar), and the forcing above this point is plotted in dark grey. In both cases, the models have been run long enough such that their T -P profiles have fully evolved from their initial states, either isothermal (a) or adiabatic (b) for P > 10 bar, as shown by the light grey dashed line, to the same steady state, a deep adiabat that corresponds to T surface = ∼1900 K -which is ∼100 K hotter than the equilibrium temperature at 10 bar.

The formation of a deep adiabat
As discussed by PT17, and in Sect. 1, an adiabatic profile in the deep atmosphere (i.e. P > ∼1 → ∼10 bar) should be a good representation of the steady state atmosphere. In order to confirm this, we performed a series of calculations with a radiatively inert deep atmosphere (i.e. no deep heating or cooling, as required by the theory of PT17). We explore this using two models, A and B, which only differ in both their initial condition and their duration. In model A, the atmosphere, including the deep atmosphere, is initially isothermal with T =1800 K and is evolved for more than 1500 Earth years in order to reach a steady state in its T -P profile (as shown in Fig. 2a). As a consequence of the long timescales required for the model to reach equilibrium and the computational cost of such an endeavour, model A (and B) is run at relatively low resolution 3 . We investigate the sensitivity of our results to spatial resolution in Sect. 3.3.1. Model B is identical to model A except in the deep atmosphere, where it is initialised with an adiabatic T -P profile for P > 10 bar. As a result of this model being initialised close to the expected equilibrium solution, model B was then run for only 100 yr in order to confirm the stability of the steady state. In both cases, we find that the simulation time considered is long enough such that the thermodynamic structure of the atmosphere has not changed for multiple advective turnover times t adv ∼ 2πR HJ /u φ . 3 However this does not mean that our models have problems conserving angular momentum; they maintain 97.44% of the initial angular momentum after over 1700 yr of simulation time (which compares well to other GCMs: Polichtchouk et al. 2014).
The first three panels show snapshots of the profiles at different pressures (P = 0.72 mbar in a, P = 455 mbar in b, and P = 12.7 bar in c), whilst the last (d) shows the time-averaged profiles at the same pressure as panel c in order to illustrate the variable nature of the deep atmosphere. We note that the period of the time average was approximately seven and a half years. Figure 2 shows that both models have evolved to the same steady state: an outer atmosphere whose T −P profile is dictated by the Newtonian cooling profile, and a deep adiabat which is slightly hotter (∼1900 K) than the cooling profile at P = 10 bar (1800 K). This is reinforced by the latitudinal and longitudinal temperature profile throughout the simulation domain. In Fig. 3 we plot the zonal wind and temperature profile at three different heights (pressures). Here we can see that in the outer atmosphere (panels a and b) the profile is dominated by Newtonian cooling, with horizontal advection (and the resulting offset hotspot) starting to become significant when moving towards middle pressures. As for the deep atmosphere (snapshot in panel c and time average in panel d), here we start to see evidence of both the heating and near-homogenisation of the deep atmosphere. We note that we refer to the atmosphere as nearly homogenised because the temperature fluctuations at for example P = 10 bar are less than 1% of the mean temperature. Importantly, this convergence to a deep adiabat not only occurs in the absence of vertical convective mixing (an effect which is absent from our models, which contain no convective driving), but also at a significantly lower pressure (P = 10 bar) than the pressure (∼40 bar for HD 209458b; Chabrier et al. 2004) at which we would expect the atmosphere to become unstable to convection (and so, in the traditional sense, prone to an adiabatic profile).
Therefore, the characteristic entropy profile of the planet is warmer than the entropy profiles calculated from standard 1D irradiated models. We discuss the implications of this result for the evolution of a highly irradiated gas giant in Sect. 4.
In model A, the steady state described above is very slow to emerge from an initially isothermal atmosphere. This is illustrated in Fig. 4  >1500 yr of evolution is required to reach the same steady state as model B.
As further discussed in Sect. 4, this slow evolution of the deep adiabat is probably one of the main reasons why this result has not been reported by prior studies of hot Jupiter atmospheres.
The mechanism advocated by PT17 relies on the existence of vertical and latitudinal motions that efficiently redistribute potential temperature. In order to determine their spatial structure, we plot in Fig. 5 the zonally and temporally averaged meridional mass-flux stream function and zonal wind velocity for model A.
Starting with the zonal wind profile (grey lines) we can see evidence for a super-rotating jet that extends deep into the atmosphere, with balancing counter flows at the poles and near the bottom of the simulation domain. In the deep atmosphere, this jet has evolved with the deep adiabat, extending towards higher pressures as the developing adiabat (almost) homogenises (and hence barotropises) the atmosphere. This barotropisation on long timescales seems similar to the drag-free simulation started from a barotropic zonal wind in Liu & Showman (2013).
The meridional mass-flux stream function is defined according to We find that the meridional (latitudinal and vertical) circulation profile is dominated by four vertically aligned cells extending from the bottom of our simulation atmospheres to well within the thermally and radiatively active region located in the upper atmosphere. These circulation cells lead to the formation of a strong, deep downflow at the equator (which can be linked to the high equatorial temperatures in the upper atmosphere), weaker upper-atmosphere downflows near the poles, and a mass-conserving pair of upflows at mid latitudes (θ = 20 • → 30 • ). The meridional circulation not only leads to the vertical transport of potential temperature (as fluid parcels from the outer atmosphere with high potential temperature are mixed with their "cooler" deep-atmosphere counterparts), but also to the almost complete latitudinal homogenisation of the deep atmosphere (with only small temperature variations remaining).
In a fully radiative model, these circulations would also mix the outer atmosphere, leading to the equilibrium temperature profiles we instead impose via Newtonian cooling (see, e.g. Drummond et al. 2018a,b for more details about the 3D mixing in radiative atmospheres). We note that the vertical extent of the zonal wind, and the structure of the lowest cells in the mass-flux stream function, appear to be affected by the bottom boundary, suggesting that they extend deeper into the atmosphere. Whilst this is interesting and important, it should not affect the final state reached by our P-T profiles. However it does does suggest that models of hot Jupiters should be run to higher pressures to fully capture the irradiation-driven deep flow dynamics.
The primary driver of the latitudinal homogenisation is fluctuations in the meridional circulation profile, which are visible within individual profile snapshots but are averaged out when we take a temporal average. This includes contributions from spatially small-scale velocity fluctuations at the interface of the large-scale meridional cells. Evidence for these effects can be seen in snapshots of the zonal and meridional flows, in an RMS analysis of the zonal velocity, and of course in the deep temperature profile that these advective motions drive. The first reveals complex dynamics, such as zonally asymmetric and temporally variable flows, which are hidden when looking at the temporal average, but which mask the net flows when looking at a snapshot of the circulation. The second reveals spatial and temporal fluctuations on the order of 5 → 10 m s −1 in the deep atmosphere. Finally, the third (as plotted in panels c and d of Fig. 3, which show snapshots of the time average of the zonal wind and temperature profile, respectively) reveals small-scale temperature and wind fluctuations, likely associated with the deep atmosphere mixing, that are lost when looking at the average "steady" state.
However, a more detailed analysis of the dynamics of this homogenisation, as well as the exact nature of the driving flows A114, page 7 of 13 A&A 632, A114 (2019) and dynamics, is beyond the scope of this paper. Although interesting in its own right, the mechanism by which the circulation is set up in the deep atmosphere of our isothermally initialised simulations might not be relevant to the actual physical mechanism happening in hot Jupiters with hot, deep atmospheres. As a consequence of both the meridional circulations described above, and the zonal flows that form as a response to the strong dayside-night-side temperature differential, the deep atmosphere T -P profile is independent of both longitude (Fig. 6a) and latitude (Fig. 6b). Only in the upper atmosphere (P < 10 bar) do the temperature profiles start to deviate from one another, reflecting the zonally and latitudinally varying Newtonian forcing. Taken together, the two panels of Fig. 6 confirm that the latitudinal and vertical steady-state circulation, the super-rotating eastward jet, and any zonally asymmetric flows act to advect potential temperature throughout the deep atmosphere, leading at depth to the formation of a hot adiabat without the need for any convective motions.

Robustness of the results
Having confirmed that a deep adiabatic temperature profile connecting with the outer atmospheric temperature profile at P = 10 bar is a good representation of the steady state within our hot-Jupiter model atmospheres, we now explore the robustness of this result.

Sensitivity to changes in the horizontal resolution
We start our exploration of the robustness of our results by confirming that the eventual convergence of the deep atmosphere on to a deep adiabat is independent of resolution. Figure 7 shows the T -P profiles obtained for three models at the same time (t ≈ 1800 yr but with different resolutions (our "base"-resolution model, A, a "mid"-resolution model, C, and a high-resolution model, D). The mid-resolution model (C) has almost reached the exact same equilibrium adiabatic profile as the base-or lowresolution case (A): comparing this with the time-evolution of model A (Fig. 4) confirms that they are both on the path to the same equilibrium state, and that a significant amount of computational time would be required to reach it. This becomes even clearer when we look at a high-resolution model (D). Here we find that despite the long timescale of the computation, the deep atmosphere still exhibits a temperature inversion, suggesting, in comparison to Fig. 4, that the model has a long way to go until it reaches the same, deep adiabat equilibrium.
In general, we find that the higher the resolution, the more slowly the atmosphere temperature profile evolves towards the adiabatic steady-state solution. This most likely stems from the fact that horizontal numerical dissipation, on a fixed dissipation timescale, decreases with increasing resolution. We note that we kept the horizontal dissipation timescale constant due to both the computational expense of the parameter study required to set the correct dissipation at each resolution, and the numerical dissipation independence of the steady-state in the deep atmosphere.
Evidence for the impact of the smallscale flows on this slow evolution can be seen in the temporal and spatial RMS profiles of the zonal flows, which reveal that as we increase the resolution by a factor of two, the magnitude of the smallscale velocity fluctuations decreases by roughly the same factor. These results are in agreement with the effect of changing the numerical dissipation timescale (τ dissip ) at a fixed resolution, where longer timescales also slow down the circulation, thereby increasing the time required to reach a steady T -P profile in the deep atmosphere (not shown). Despite these numerical limitations, it remains clear that the presence and strength of any numerical dissipation do not affect the steady state solutions of the simulation, which remains as an adiabatic P-T profile in the deep atmosphere.

Sensitivity to changes in the forcing function of the upper atmosphere
We next explore how the deep adiabat responds to changes in the outer atmosphere irradiation and thermal emission (via the imposed Newtonian cooling). The aim is not only to test the robustness of the deep adiabat, but also to explore the response of the adiabat to changes in the atmospheric state. As part of this study, the two scenarios we consider were initialised using the evolved adiabatic profile obtained in model A, but with a modified outer atmosphere cooling profile such that T night = 1200 K (model J) or T night = 2200 K (model K). Figure 8 shows the equilibrium T -P profiles (solid lines) as well as snapshots of the T -P profiles after only 200 yr of "modified" evolution (dashed lines); a plot of model A is also included to aid comparison. Model J evolves in less than 200 yr towards a new steadystate profile that corresponds to the modified cooling profile. The deep adiabat reconnects with the outer atmospheric profile at P = 10 bar and ∼1250 K (in agreement with the relative offset found in our 1800 K models, A and B). The meridional mass circulation (not shown) displays evidence for the same qualitative flows driving the vertical advection of potential temperature as models A and B. However, it also shows signs that it is still evolving, suggesting that the steady-state meridional circulation takes longer to establish than the vertical temperature profile.
In model K, we find that, 200 yr after modifying the cooling profile of the outer atmosphere, the deep atmosphere has not yet reached a steady state. Indeed, it takes approximately 1000 yr of evolution for it to reach equilibrium, which we show as a solid line in Fig. 8. This confirms that model K, although slow to evolve relative to the cooling case (model J), does eventually settle onto a deep equilibrium adiabat. Additionally, this adjustment occurs significantly faster than the equivalent evolution of a deep adiabat from an isothermal start.
Based on the results of this section, we conclude that it takes less time for the deep atmosphere to cool than to warm when it evolves toward its adiabatic temperature profile. In order to understand this timescale ordering, we have to note that the only way for the simulation to inject or extract energy is through the fast Newtonian forcing of the upper atmosphere and also that the thermal heat content of the deep atmosphere is significantly larger than that of the outer layers. The deep ( d ) and upper ( u ) atmospheres are connected by the advection of potential temperature that we rewrite in a conservative form as an enthalpy flux: ρc p T u , and we simplify the process to two steps between the two reservoirs (assuming they have similar volumes): injection or extraction by enthalpy flux and Newtonian forcing in the upper atmosphere. . Equatorially averaged T -P profiles for three models: A (green), J (yellow) and K (orange). The orange (K) and yellow (J) models have had their outer atmosphere cooling modified such that T eq = 2200 K or 1200 K, respectively. The solid lines represent the equilibrium T -P profiles whilst the dashed lines represent the T -P profiles 200 yr after the outer atmosphere forcing was adjusted (shown in dark grey for each model). We note that after 200 yr of "modified" evolution, only the 2200 K model has not reached equilibrium.
-In the case of cooling, the deep atmosphere contains too much energy and needs to evacuate it. It will setup a circulation to evacuate this extra energy to the upper layers with an enthalpy flux that would lead to an upper energy content set by ρ u c v T u ∼ ρ u c v T u,init + ρ d c v (T d,init − T d,eq ) if we ignore first Newtonian cooling. Here, T u would then be very large essentially because of the density difference between the upper and lower atmosphere. The Newtonian forcing term proportional to −(T u − T u,eq )/τ is then very large and can efficiently remove the energy from the system.
-In the case of heating, the deep atmosphere does not contain enough energy and needs an injection from the upper layers. This injection is coming from the Newtonian forcing and can at log(τ 220 )=7.5 (E) log(τ 220 )=11 (F) log(τ 220 )=15 (G) log(τ 220 )=20 (H) log(τ 220 )=22.5 (I) Fig. 9. Newtonian cooling relaxation timescale profiles used in the models shown in Fig. 10. We note that a smaller value of τ means more rapid forcing towards the imposed cooling profile (which in all cases is isothermal in the deep atmosphere, where P > 10 bar), and that the relaxation profiles are identical for P < 10 bar (grey line).
first only inject ρ u c v (T u,eq − T u,init ) in the system. The enthalpy flux will then lead to an energy content in the deep atmosphere given by we assume that all the extra energy is pumped by the deep atmosphere. Because of the density difference and the limited variations in the temperature caused by the forcing, the temperature change in the deep atmosphere is small and will require more injection from the upper layers to reach equilibrium. However, even in the most favourable scenario in which all the extra energy is transferred, the Newtonian forcing cannot exceed −(T u,init − T u,eq )/τ which explains it takes a much longer time to heat the deep atmosphere than to cool it.

Sensitivity to the addition of Newtonian cooling to the deep atmosphere
It is unlikely that the atmosphere will suddenly turn thermally inert at pressures greater than 10 bar. Rather, we expect the thermal timescale to gradually increase with increasing pressure. In this section, we examine the sensitivity of the deep atmospheric flows, circulations, and thermal structure to varying levels of Newtonian cooling. Additionally, we are motivated to quantify the maximum amount of Newtonian cooling under which the deep atmosphere is still able to maintain a deep adiabat.
To explore this, we consider five models each with different cooling timescales at the bottom of the atmosphere (i.e. five different values of log (τ 220 )). From this, we can then linearly interpolate the relaxation timescale in log (P) space between 10 and 220 bar. The resultant profiles are plotted in Fig. 9, and can be split into three distinct groups: (1) The relaxation profile with log(τ 220 ) = 7.5 (model E) represents a case with rapid Newtonian cooling that does not decrease with increasing pressure.
(2) The case log (τ 220 ) = 11 (model F) is a simple linear continuation of the relaxation profile we use between P = 1 bar and 10 bar. It is the simplest possible extrapolation of the upper atmosphere thermal timescale profile, and likely represents the strongest realistic forcing in the deep atmosphere. Pressure (bar) log(τ 220 ) = 7.5 (E) log(τ 220 ) = 11 (F) log(τ 220 ) = 15 (G) log(τ 220 ) = 20 (H) log(τ 220 ) = 22.5 (I) T forced T init Fig. 10. Snapshots of the T -P profile for five initially adiabatic simulations (coloured lines -based on model B, and with the same outer atmosphere cooling profile (dark grey)) which are then forced to a deep isothermal profile (grey dashed line) with varying log (τ 220 ) (Eq. (10)).
accordance with expectations born out from 1D atmospheric models of hot Jupiter atmospheres (see, e.g. Iro et al. 2005).
The results we obtained are summarised by the T -P profiles we plot in Fig. 10. For low levels of heating and cooling in the deep atmosphere (models G, H and I), the results are almost indistinguishable from models A and B, with only a decrease in the outer atmosphere connection temperature of a few Kelvin in model G. We find a more significant reduction in the temperature of the T -P when we investigate model F, in which we set log (τ 220 ) = 11. In particular, there is a deepening of the connection point between the outer atmosphere and the deep adiabat, which only becomes apparent for P > 20 bar in this model. This result suggest that model F falls near the pivot point between models in which the deep atmosphere is adiabatic and those that relax toward the imposed temperature profile. This is confirmed by model E, in which τ 220 = 7.5, where we find that the deep adiabat has been rapidly destroyed (in <30 yr), such that the deep T -P profile corresponds to the imposed cooling profile throughout the atmosphere. This occurs because the Newtonian timescale has become smaller than the advective timescale, which means that the imposed temperature profiles dominate over any dynamical effects.
Before closing this section, let us briefly comment on the meridional circulation profiles obtained in those models that converge onto a similar deep adiabatic temperature profile (models G, H and I). For all of them, we recover the same qualitative structure we found for model A, characterised by meridional cells of alternating direction that extend from the deep atmosphere to the outer regions. However, the finer details of the circulations differ from the ones seen in model A. This is illustrated in Fig. 11 which displays the meridional circulation and zonal flow profiles for models G (Fig. 11a) and H (Fig. 11b). As the Newtonian cooling becomes faster in the deep atmosphere, the number of meridional cells increases (see also Fig. 5), to the point that in model E, no deep meridional circulation cells exists and the deep circulation profile is essentially unstructured. Despite these differences in the shape of the meridional circulation, the steady-state profile obtained in these simulations in the deep atmosphere is again an adiabatic PT profile provided the Newtonian cooling is not (unphysically) strong.

Conclusions based on the simulation results
By carrying out a series of 3D GCM simulations of irradiated atmospheres, we have shown in the present paper that if the deep atmosphere is initialised on an adiabatic PT profile, it remains, as a steady state, on this profile. We have also demonstrated the fact that if the deep atmosphere is initialised at an overly hot state, it rapidly cools down to the same steady state adiabatic profile, whereas if the deep atmosphere is initialised at an overly cold state, it slowly evolves towards the steady state adiabatic profile. Furthermore, in all the above cases, the deep adiabat forms at lower pressures than those at which we would expect the atmosphere to be convectively unstable based on information from 1D models. We have also shown that this steady-state adiabatic profile is unaffected by changes in the deep Newtonian cooling and is independent of the details of the flow structures, provided that the velocities are not completely negligible. The hot adiabatic deep atmosphere is the natural final outcome of the simulations for various resolutions, even though the timescale to reach steady state is longer at higher resolution when starting from an overly cold initial state.
When the simulations are initialised on an overly cold profile, the time needed to reach the steady state is of the order of t ∼ 1000 yr, explaining why the formation of a deep adiabat has not been seen in previous GCM studies: this timescale is far greater than the time taken for the outer atmosphere to reach an equilibrium state (t < ∼1 yr for P < 1 bar). As a result, the vast majority of published GCM models only contain a partially evolved deep atmosphere, the structure of which is directly comparable to the early outputs of our isothermally initialised calculation. Examples of this early evolution of the deep atmosphere towards a deep adiabat (as seen in the early outputs plotted in Fig. 4) include Fig. 6 of , where the deep temperature profile shows signs of heating from its initial isothermal state, albeit only on the irradiated side of the planet, Fig. 7 of Amundsen et al. (2016), where we see a clear shift from their initially isothermal deep atmosphere towards a deep adiabat, and Fig. 8 of Kataria et al. (2015), where we again  (Mayne et al. 2014b), and including a robust two-stream radiation scheme (Amundsen et al. 2014). Here we show snapshots of the T-P profile at 0.25 (purple), 2.5 (green), and 25 (orange) Earth years, along with two example adiabats (grey dotted and dashed lines) designed to show how the deep atmosphere gets warmer and connects to steadily warmer adiabats as the simulation progresses. We note that at the end of the simulated time this progression is ongoing towards a deep, hot adiabat, albeit at an increasingly slow rate. see a temperature inversion and a push towards a deep adiabat for Wasp-43b. It is tempting to think that if these simulations were run longer, they would evolve to a similar, deep adiabatic structure (with a corresponding increase in the exoplanetary radius). In order to investigate this possibility, we extended the model of Amundsen et al. 2016, run with the Unified Model of the Met Office (which includes a robust two-stream radiation scheme that replaces the Newtonian Cooling in our models), for a total of ≈25 Earth years. The results are shown in Fig. 12, where we plot the pressure-temperature profile at three different times, A114, page 11 of 13 A&A 632, A114 (2019) along with examples of the approximate deep adiabat that best matches each snapshot. We see that the deep atmosphere rapidly converges towards a deep adiabat with further vertical advection of potential temperature warming up this adiabat as the simulation goes on. Since this process continues during the simulation, the result not only reinforces our conclusions but suggests that our primary Newtonian cooling profile represents a reasonable approximation of the incident irradiation and radiative loss.
The results obtained in the present simulations suggest that future hot Jupiter atmosphere studies should be initialised with a hot, deep adiabat starting at the bottom of the surface irradiation zone (P ∼ 10 bar for HD 209458b). Furthermore, in a situation where the equilibrium profile in the deep atmosphere is uncertain, we suggest that this profile should be initialised with a hotter adiabat than expected, rather than a cooler one. The simulation should then be run long enough for the deep atmosphere to reach equilibrium. This is in agreement with the results of Amundsen et al. (2016), who also suggested that future GCM models should be initialised with hotter profiles than currently considered. For instance, recent 3D simulations of HD 209458b were initialised with a hotter interior T-P profile (e.g. one of the models of the aforementioned authors was initialised with an isotherm that is 800 K hotter than typically used in GCM studies, thus bringing the deep atmosphere closer towards its deep adiabat equilibrium temperature), and show important differences, on the timescales considered, between the internal dynamics obtained with this setup, and the ones obtained with a cooler, more standard, deep atmospheric profile (see, Lines et al. 2018aLines et al. ,b, 2019. Using the aforementioned, more accurate atmosphere initial profiles should not only bring these models towards a more physical hot Jupiter parameter regime (with a correct inflated radius), but should also provide a wealth of information on how the deep adiabat responds to changes in parameter and computational regime.

Evolution of highly irradiated gas giants
The results obtained in the present GCM simulations have strong implications for our understanding of the evolution of highly irradiated gas giants. As mentioned above, we first emphasise that simulations initialised from an overly cold state are not relevant for the evolution of inflated hot Jupiters (although they could be of some interest for re-inflation, but this is beyond the scope of this paper). Indeed, inflated hot Jupiters are primarily in a hot initial state, and as far as the evolution is concerned only the steady state of the atmosphere is important; the shorter timescales needed to reach this steady state are irrelevant for the evolution (with a typical Kelvin Helmholtz timescale of ∼1 Myr).
As shown in the present simulations, provided they are run long enough, hot Jupiter atmospheres converge at depth, that is in the optically thick domain, to a hot adiabatic steady-state profile without the need to invoke any dissipation mechanism such as ohmic or kinetic-energy dissipation. These 3D dynamical calculations confirm the 2D steady-state calculations of PT17. Importantly enough, the transition to an adiabatic atmospheric profile occurs at lower pressures than the ones at which the medium is expected to become convectively unstable (and therefore adiabatic according to the Schwarzchild criterion). This means that the planet lies on a hotter internal entropy profile than suggested by 1D irradiation models, yielding a larger radius. The mechanism of potential temperature advection in the atmosphere of irradiated planets thus provides a robust solution to the radius inflation problem.
As mentioned previously, almost all scenarios suggested so far to resolve the anomalously inflated planet problem rely on the (uncomfortable) necessity to introduce finely tuned parameters. This is true in particular for all the different dissipation mechanisms, whether they involve kinetic energy or ohmic and tidal dissipation. This is in stark contrast with the present mechanism, in which entropy (potential temperature) is advected from the top to the bottom of the atmosphere. High-entropy fluid parcels are moved from the upper to the deep atmosphere and toward high latitude while low-entropy fluid parcels come from the deep atmosphere and are deposited in the upper atmosphere. This gradually changes the entropy profile until a steady-state situation is obtained. Although an enthalpy (and mass and momentum) flux is associated with this process down to the bottom of the atmosphere (characterised by some specific heat reservoir), this does not require a dissipative process (from kinetic, magnetic, or radiative energy reservoirs into the internal energy reservoir).
In order to both characterise this deep heating flux and to confirm that our hot, deep adiabat would not be unstable due to high-temperature radiative losses, we also explored the vertical enthalpy flux in our model and compared it to the radiative flux, as calculated for a deep adiabat using ATMO (PT17). This analysis reveals that the vertical enthalpy flux dominates the radiative flux at all P > 1 bar: for example, averaging over a pressure surface at P = 10 bar, we find a net vertical enthalpy flux (ρc p T u z) of −1.04 × 10 8 ergs −1 cm −2 compared to an outgoing radiative flux of 7.68 × 10 6 ergs −1 cm −2 , suggesting that any deep radiative losses are well compensated by energy (enthalpy) transport from the highly irradiated outer atmosphere. This result is reinforced by the UM calculation shown in Fig. 12, which intrinsically includes this deep radiative loss and shows no evidence of cooling due to deep radiative effects.
This (lack of a requirement for additional dissipative processes) is of prime importance when trying to understand the evolution of irradiated planets. Whereas dissipative processes imply an extra energy source in the evolution ( M˙ dm, wherė is the energy dissipation rate, to be finely tuned), to slow down the contraction of the planet, there is no need for such a term in the present process. Indeed, as an isolated substellar object (i.e. without nuclear energy source) cools down, its gravitational potential energy is converted into radiation at the surface, with a flux σT 4 eff . Let us now suppose that the same object is immersed into an isotropic medium characterised by a pressure of P = 220 bars and a temperature of T ∼ 4000 K, typical conditions in the deep atmosphere of 51Peg-B-like hot Jupiters. Once the original inner adiabat of the object (after its birth) has cooled down to 4000 K at 220 bars, the thermal gradient between the external and internal media will be null, which essentially reduces the local convective flux and the local optically thick radiative flux to zero. Thus, the cooling flux will be reduced to almost zero. At this point, the core can no longer significantly cool and is simply in thermal equilibrium with the surrounding medium. Both the contraction and the cooling flux are essentially insignificant: dR/dt ≈ 0, σT 4 int ≈ 0, in which we define σT 4 int as the radiative and convective cooling flux at the interior-atmosphere boundary. Indeed, convection will become inefficient in transporting energy and any remaining radiative loss in the optically thick deep core will be compensated by downward energy transport from the hot outer atmosphere. For a highly irradiated gas giant, the irradiation flux is not isotropic, but the combination of irradiation and atmospheric circulation will lead to a similar situation, with a deep atmosphere adiabatic profile of ∼4000 K at 220 bars for all latitudes and longitudes. Therefore, there is no further significant cooling of the planet's interior and we also have σT 4 int ≈ 0. The evolution of the planet is stopped (dS /dt ≈ 0), or in other words its cooling time is now prohibitively long, and the planet lies on a constant adiabat determined by the equilibrium between the inner and atmospheric ones at the interior-atmosphere boundary. The situation will last as long as the planet-star characteristics remain the same, illustrating the robustness of the potential temperature advection mechanism to explain the anomalous inflation of these bodies.
Therefore, irradiation-induced advection of potential temperature appears to be the most natural and robust process to resolve the radius-inflation puzzle. We note that this does not exclude other processes (e.g. dissipative ones) from operating within the atmospheres of hot Jupiters, but they are unlikely to be the dominant mechanisms responsible for the radius inflation.