© F. Sainsbury-Martinez et al. 2021

1 Introduction

The anomalously large radii of highly irradiated gaseous exoplanets remains one of the key unresolved areas in our understanding of extrasolar giant planetary atmospheres. In recent years, however, this state of affairs has started to change, at least for highly irradiated Jupiter-like planets (i.e. hot Jupiters). In particular, a large number of different mechanisms have been suggested as possible explanations for the observed inflation (see Baraffe et al. 2014 and Fortney et al. 2021, for recent reviews of this topic), as well as the correlation between said inflation and the planetary irradiation (for examples, see Demory & Seager 2011; Laughlin et al. 2011; Lopez & Fortney 2016; Sestovic et al. 2018; Thorngren & Fortney 2018). However, much less work has gone into trying to understand the inflated radii of lower-mass brown dwarfs. At first glance, it seems reasonable to suggest that the same mechanism responsible for radius inflation in hot Jupiters should also operate within these brown dwarf atmospheres given how qualitatively similar they are. However, this is not guaranteed as not only must the inflation mechanism operate in an atmosphere in which the influence of gravity is significantly enhanced, it must also be able to explain the observed correlation between inflation and irradiation for brown dwarfs (Casewell et al. 2020b). This correlation is not simply a linear increase in radius with irradiation, but, rather, a combination of two trends: an increase in the observed radius with irradiation for brown dwarfs orbiting main sequence stars and a lack of inflation for very highly irradiated brown dwarfs orbiting hot white dwarfs (see Sect. 2). We note that a third trend, namely, that of decreasing inflation with increasing giant planetary mass, has also been observed; however, Thorngren et al. (2021) suggests that this effect ranks rather low in importance when compared to the link between irradiation and inflation.

Recently, Tremblin et al. (2017, hereafter PT17) and Sainsbury-Martinez et al. (2019, hereafter FS19) explored in steady-state 2D or long-timescale 3D models, respectively, a hot Jupiter inflation mechanism that naturally arises from first physical principles, namely: the formation of a hot, deep, and non-convective adiabat via the vertical transport of potential temperature.

Rather than developing via convective flows, the formation of this adiabatic region (i.e. region of constant potential temperature) in the deep atmosphere is instead driven by the vertical advection of potential temperature from the highly irradiated outer atmosphere to the deep atmosphere by large-scale dynamical motions, such as the equatorial super-rotating jet and its associated mass balancing flows. Here, it is almost completely homogenised by deep zonal and meridional flows and circulations.Critically this results in the ‘deep’ temperature-pressure profile converging onto an adiabat at lower pressures than that at which the atmosphere might be expected to be unstable to convection. As a result, the outer atmosphere connects to a hotter internal adiabat than would be obtained from a standard 1D ‘radiative-convective’ model potentially leading to a highly irradiated exoplanet having a larger radius than predicted. PT17 showed that this mechanism is able to reproduce the observed correlation between radius inflation and irradiation in the case of hot Jupiters.

The purpose of this paper is to extend the work of PT17 and FS19 into the brown dwarf regime and explore whether the presence of a deep, hot, and non-convective adiabat can also explain the radius inflation trend observed for these objects. To that end, we performed a series of somewhat idealised, long-timescale, 3D global circulation model (GCM) calculations of three brown dwarfs: Kepler-13Ab, KELT-1b, and SDSS1411B, the latter of which shows no signs of radius inflation in observations and is in a very short orbit around a hot white dwarf. These simulations will allow us to not only explore if a stable, inflating, advective adiabat is able to form, but also to investigate the flows, circulations, and fluxes associated with said development, or a lack thereof.

The structure of this work is as follows: We start, in Sect. 2, by detailing some key observations of brown dwarfs. This includes an overview of the mass-radius-irradiation relation for known transiting hot brown dwarfs as well as key observational parameters for the three brown dwarfs (and their host-stars) under consideration here: Kepler-13Ab, KELT-1b, and SDSS1411B. Next, in Sect. 3, we introduce the model used here to explore this inflation mechanism. This includes an overview of both the numerical scheme used in our 3D GCM (DYNAMICO), as well as the 1D models that we use to define and parametrise said models radiative dynamics. Then, in Sect. 4, we detail our results. We start by exploring the differences between models of each of our three browndwarfs, after which we introduce a series of test cases designed to explore at what pressure (referred to as the horizontal convergence pressure) each brown dwarf model atmosphere might be expected to develop a deep, hot, advective adiabat (Sect. 4.1). Having calculated a horizontal convergence pressure range for each brown dwarf, we then explore how robust each brown dwarfs advective adiabat is to the presence of deep radiative forcing on a similar (or slightly weaker) timescale to that seen in 1D ‘radiative-equilibrium’ models (Sect. 4.2). Finally, in Sect. 4.3, we explore the flows, circulations, and fluxes which lead to the development, or lack thereof, of a radiatively stable deep adiabat for each of our brown dwarfs. We conclude, in Sect. 5, by providing concluding remarks, discussing the implications of our results, and introducing suggestions for future computational studies of lower mass brown dwarf atmospheres using next-generation radiative 3D GCMs.

Fig. 1 Mass-radius relation for 25 of the known highly irradiated brown dwarfs transiting either a main sequence star (Carmichael et al. 2020 and references therein) or a white dwarf (Casewell et al. 2020a; Littlefair et al. 2014; Parsons et al. 2017), with the brown dwarfs under consideration here clearly labelled. Additionally, the effective temperatures of the brown dwarf hosting stars are indicated by the colour of each brown dwarfs circle, while the size of each coloured circle is proportional to the total incident surface irradiation. Also shown are the Sonora-Bobcat evolutionary models of Marley et al. (2018) for 100 Myr and 600 Myr at only solar metalicity, and at 2 Gyr, 6 Gyr, and 10 Gyr for both solar metallicities (grey) and low metallicities ([Fe/H] =−0.5 – black).

2 Observations of brown dwarfs: Kepler-13Ab, KELT-1b, and SDSS1411B

To date, observations have revealed 25 highly irradiated brown dwarfs that transit a main sequence star (Carmichael et al. 2021 and references therein) and three that transit white dwarfs (Casewell et al. 2020a; Littlefair et al. 2014; Parsons et al. 2017). Figure 1, in which we plot the mass-radius-irradiation relation for 25 of these transiting brown dwarfs, reveals a somewhat similar trend to that seen for highly irradiated Jupiter-like planets: a general increase in the observed radius with increasing irradiation, albeit with a significant reduction in inflation strength towards higher masses (Thorngren et al. 2021). However, this trend is not universal with two of the most highly irradiated brown dwarfs included in Fig. 1 showing little to no sign of radius inflation; that is to say that the observations closely match the radii predicted by evolutionary models. So, we consider what separates these brown dwarfs from their compatriots and we find that one major factor appears to be that they are situated in short orbits around white dwarfs.

To date, no mechanism has been suggested that can explain either the general trend of increasing inflation with irradiation (albeit tempered by mass) or the exception to this trend for very highly irradiated brown dwarfs orbing hot white dwarfs.

Here, we intend to investigate whether potential temperature advection can (via a simple physical mechanism) finally provide this explanation. To that end, we intend to explore the steady state deep atmospheres of three representative brown dwarfs: Kepler-13Ab, KELT-1b, and SDSS1411B. Studies of Kepler-13Ab and KELT-1b reveals a significant discrepancy between the observed and theoretical radii for both brown dwarfs, albeit to a slightly smaller degree than might be expected for a similarly located hot Jupiter thanks to the attenuating effect of their comparatively higher masses (Thorngren et al. 2021). Whereas observations of SDSS1411b show little to no signs of radius inflation despite its very short orbit around a very hot white dwarf, resulting in a strong surface irradiation that would otherwise imply the presence radius inflation even when accounting for the attenuating effects of mass.

All three objects studied are assumed to be tidally locked due to their proximity to their host star. Host star and brown dwarf parameters can be found in Table 1.

2.1 Kepler-13Ab

Kepler-13Ab is a low-mass brown dwarf detected by the Kepler space telescope (Borucki et al. 2011; Szabó et al. 2011). It has an inflated radius of 1.512 ± 0.035 R_J, a mass of 9.28 ± 0.16 M_J, and an orbital period of ~1.76 days (Esteves et al. 2015). The observations and models suggest that Kepler-13Ab has an average night-side (brightness) temperature of ~ 2500 K (Shporer et al. 2014; Esteves et al. 2015) and a day-side (brightness) temperature of between 2750 K (Shporer et al. 2014)and 3490 K (Esteves et al. 2015). For our models, we take an average of these two results and set the day-night temperature difference to 600 K at 0.01 bar. Finally, calculations suggest that the zero albedo equilibrium temperature of this brown dwarf is 2570 K.

The spectra of Kepler-13Ab reveal an absorption feature due to water (Beatty et al. 2017), whereas 1D models would suggest that any water feature should either be in emission (because of a temperature inversion) or muted (from an isothermal structure: see Lothringer et al. 2018). A number of explanations for this feature have been put forth (e.g. Lothringer & Casewell 2020), including the suggestion that the water absorption feature is seen because observations are probing the internal adiabat. However, in 1D models, this requires that the internal temperature be increased from 1000 K, corresponding to a 1D radius of ~ 1.195 R_J, to at least 1750 K– an increase which is incompatible with the age and mass of Kepler-13Ab. Here, we investigate whether, instead, the vertical advection of potential temperature can explain such a hot interior (Sect. 4.1.1).

The host star, Kepler-13A, is a chromospherically moderately rapidly rotating (v sin i = 78 ± 15 km s⁻¹) A0V star in a close (a = 0.410 AU) triple system (consisting of Kepler-13A and the binary subsystem Kepler-13B and Kepler-13C), with an effective temperature of 7650 ± 250 K (Table 1: Shporer et al. 2014).

2.2 KELT-1b

KELT-1b was is the first object detected by the Kilodegree Extremely Little Telescope-North (KELT-North) transit survey (Siverd et al. 2012). It has an inflated radius of ${1.116}_{-0.029}^{+0.038}{R}_{\text{J}}$ and a mass of 27.38 ± 0.93 M_J, which is three times that of Kepler-13Ab. Its orbital period is 1.21 days (Siverd et al. 2012), slightly shorter than that of Kepler-13Ab. The day-side of KELT-1b has a blackbody (brightness) temperature of 3340 ± 110 K, while the night-side has a more poorly constrained (brightness) temperature of ${1820}_{-1150}^{+640}\text{K}$ (Beatty et al. 2020). Due to the uncertainty in these values, our models use a day-night temperature difference of 1200 K at 0.01 bar (with additional tests run with a 1600 K temperature difference – tests which reveal similar dynamics with just a slight increase in jet speed and deep heating rate). Finally, calculations suggest that the zero albedo equilibrium temperature of this brown dwarf is 2426 K.

Comparisons with the 1D atmospheric and stellar evolution models of Phillips et al. (2020) suggest that in order for the models to match the observed radius, the internal temperature of the brown dwarf must be increased from ~ 700 K (corresponding to a 1D radius of ~0.97 R_J) to ~ 1600 K, an increase which would require that the brown dwarf is either significantly younger or more massive than observations of the system suggest. Again, we aim to investigate whether, instead, potential temperature advection can explain such a hot, deep adiabat in Sect. 4.1.1.

The host star, KELT-1A, is an F5V, with a slightly slower surface rotation (v sin i = 56 ± 2 km s⁻¹) than Kepler-13A (Table 1: Siverd et al. 2012; Beatty et al. 2014). Due to its spectral type, KELT-1A is expected to have low levels or no activity (Siverd et al. 2012), although von Essen et al. (2021) have detected light curve modulation, which may be due to surface spots.

2.3 SDSS1411B

SDSS1411B (also know as SDSS-J1411+2009) is a highly irradiated brown dwarf in a very short orbit (~2 h) around a white dwarf (Beuermann et al. 2013; Littlefair et al. 2014; Casewell et al. 2018). It has a radius of 0.072 ± 0.004 R_⊙, consistent with model results, and a mass of 0.050 ± 0.002 M_⊙, the highest mass brown dwarf in our study. H band observations suggest that SDSS1411B has a night-side brightness temperature of ${1540}_{-70}^{+90}\text{K}$ and a day-side brightness temperature of 1730 ± 70 K (Casewell et al. 2018), giving a day-night temperature difference of 200 K at 0.01 bar. Finally, calculations suggest that the zero albedo equilibrium temperature of this brown dwarf is 1367 K.

At first glance, it might appear that the lack of radius inflation seen here is because SDSS1411B is significantly cooler than either KELT-1b or Kepler-13Ab, however, this is not the case. Firstly, SDSS1411B is not significantly cooler than numerous inflated hot Jupiters, for instance, HD 209458b, and, as Sainsbury-Martinez et al. (2019) show, this object exhibits significant radius inflation. Secondly, Thorngren & Fortney (2018) suggests that planets at ~ 1500–1750 K are often highly inflated and that the inflation efficiency (or, rather, the fraction of incident flux that goes in to inflating the planet) appears to peak at these temperatures. Thus, this is why we are so interested in finding out what separates SDSS1411B from other inflated gaseous planets as well as finding out if the advective heating mechanism can reproduce this difference.

The primary, SDSS-J1411A, is a hydrogen rich (DA) white dwarf with an effective temperature of 13 000 ± 300 K, a slightly below-average mass of 0.53 ± 0.03 M_⊙ and an age greater than 3 Gyr (Table 1: Littlefair et al. 2014).

3 Method

In order to investigate whether the advection of the potential temperature can explain the observed radius-irradiation-host trend observed for highly irradiated brown dwarfs, we must explore the steady state dynamics of the deep atmospheres of said objects. However, even if we start close to the equilibrium state (a deep adiabat based on the horizontal convergence pressure-temperature point – FS19), to reach a steady-state at all pressures is a slow process. Typically, thanks to the longer timescales of deep dynamics (see Rauscher & Menou 2010; Mayne et al. 2014b; Sainsbury-Martinez et al. 2019 etc.), these models require an order of magnitude more simulation time (~ 300 Earth years) than would be used for equivalent models which focus only upon reaching equilibrium in the outer atmosphere (e.g. ten Earth years to reach day-side radiative equilibrium at pressures of < 10 bar for the brown dwarf WD0137 – Lee et al. 2020).

This, along with the limits on available computational resources, precludes us from using any of the current generation of fully radiative GCMs (such as the MetOffice GCM, Exo-FMS, etc...) for our models. As such, we turned to the highly computationally efficient next-generation GCM DYNAMICO (Sect. 3.1), which uses a Newtonian Cooling approach to model the radiative forcing (Sect. 3.2) that we parameterise using day-side, night-side, and equilibrium 1D models calculated using the PHOENIX code (Sect. 3.3).

3.1 DYNAMICO

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 spherical grid (Dubos et al. 2015). It is being developed as a next-generation dynamical core for Earth and planetary climate studies at the Laboratoire de Météorologie Dynamique and is publicly available¹. Recently it has been used by FS19² to model the deep atmosphere of the hot Jupiter HD209458b over a very long-timescale, and by Spiga et al. (2020) to model the atmosphere of Saturn at a very high resolution.

In brief, DYNAMICO takes an energy-conserving Hamiltonian approach to solving the primitive equations of meteorology. This has been shown to be suitable for modelling the atmospheres of both hot Jupiters and brown dwarfs (see Showman et al. 2020 and references therein), although it may not be valid in other planetary atmospheres, such as small Neptunes and super Earths, (Mayne et al. 2019). 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 grid) for which the total number of horizontal cells, N, is defined by the number of subdivisions, d, of each edge of the main spherical icosahedral³: $$N=10d^2+2.$$(1)

In all the models considered here, we set the number of subdivisions to 30, which results in a total horizontal resolution of 9002 cells. This corresponds to an angular resolution of approximately 2.5°.

As for the vertical grid, DYNAMICO uses a pressure coordinate system whose levels can be defined by the user at runtime. In our models, this means 44 pressure levels that are linearly spaced in $\log \left(P\right)$ space. Finally, the boundaries of our simulations are closed and stress-free with zero-energy transfer (i.e. the only means of energy injection and removal are the Newtonian cooling relaxation scheme – described in Sect. 3.2 – and the horizontal numerical dissipation required to stabilise the system – see below). We note that unlike some other GCM models of gaseous giants (e.g. Schneider & Liu 2009; Liu & Showman 2013; Showman et al. 2019; Tan & Showman 2021), we do not include an additional frictional (i.e. Rayleigh) drag scheme at the bottom of our simulation domain, instead relying on the hyperviscosity (see below), the extended deep atmosphere (with P_max = 1000 bar), and the 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 hyperdiffusion 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). 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 model and somewhat accounts for sub-grid-scale dynamics. This approach is known as the large eddy approximation and has long been standard practice in both the stellar (e.g. Miesch 2005) and planetary (e.g. Cullen & Brown 2009) atmospheric modelling communities.

In a series of benchmark cases, Heng et al. (2011) have shown that both spectral and finite-difference-based dynamical cores which implement horizontal hyperdiffusion filters can produce differences of the order of tens of percent in the temperature and velocity fields when varying the dissipation strength. A similar sensitivity to the dissipation timescale was found and explored by FS19 for models of HD209458b. The solution enacted by FS19 was to calibrate the dissipation timescale by minimising unwanted small-scale numerical noise as well as replicating published benchmark results. As a result of these calibration tests, they eventually settled on a dissipation timescale of τ_dissip = 2500 s for their low resolution runs and τ_dissip = 1250 s for their mid resolution runs (which match the horizontal resolution of our brown dwarf models here). Due to the similarity of of our model setup, as well as the fact that we are both modelling highly irradiated gaseous planets, it would seem natural to use the same dissipation timescale for our models (since we lack benchmark models and detailed observationsto directly compare against). However, as we expect at least some of our brown dwarf models to contain an even stronger super-rotating jet than in FS19 models of HD209458b, we chose to retain the weaker dissipation (τ_dissip = 2500 s) used in their low resolution model. We note that a limited number of test cases with τ_dissip = 1250 s reveal very similar dynamics, including smaller differences in the jet speed than found when varying the temporal averaging period.

Finally, since DYNAMICO (like many other GCMs) does not include a dynamic time-step, a benchmark model was run for each brown dwarf in order to calculate the maximum stable time-step (which both balances the need to correctly reproduce the physics, while also minimising the computational cost of the long-timescale models). For both Kepler-13Ab and KELT-1b, this resulted in a 120-s time-step being chosen, whereas the more turbulent flows of SDSS1411B (Sect. 4.3) necessitated a shorter 30-s time-step, quadrupling this model’s relative computational cost.

3.2 Radiative forcing via Newtonian cooling

In our simulations of brown dwarf atmospheres using DYNAMICO, we do not directly model either the incident thermal radiation onthe day-side or the thermal emission on the night-side of the exoplanet. This would be prohibitively computationally expensive for the large array of long-timescale simulations we present here. Instead we use a simple thermal relaxation scheme to model these radiative effects, with a spatially varying equilibrium temperature profile, T_eq, and a radiative relaxation timescale, τ_rad, that increases with pressure throughout the outer atmosphere. Specifically, we model the radiation by adding a source term to the temperatureevolution equation that takes the form $$\frac{\partial T\left(P,\theta ,\varphi \right)}{\partial t}=-\frac{T\left(P,\theta ,\varphi \right)-T_{\text{eq}}\left(P,\theta ,\varphi \right)}{{\tau }_{\text{rad}}\left(P\right)}.$$(2)

This method, known as Newtonian cooling, has long been applied within the 3D GCM exoplanetary community (i.e. Guillot & Showman 2002; Showman et al. 2008; Rauscher & Menou 2010; Showman & Polvani 2011; Mayne et al. 2014a; Guerlet et al. 2014; Mayne et al. 2014b; or FS19), although it is gradually being replaced by coupling with simplified but much more computationally expensive radiative transfer schemes (e.g. Showman et al. 2009; Rauscher & Menou 2012 or Amundsen et al. 2016)due to its many limitations (particularly when studying the outer atmosphere), some of which we discuss in Sect. 4.1.

For the majority of our models, $T_{\text{eq}}\left(P,\theta ,\varphi \right)$ is calculated using the combination of a simple linear in $\log \left(P\right)$ space fit to a 1D equilibrium model (see Sect. 3.3 for details about how these 1D models where calculated) with a pressure dependent day-night temperature difference: $$\begin{array}{lll}{T}_{\text{eq}}\left(P,\theta ,\varphi \right)\hfill & =\hfill & T_{\text{eq-1D}}\left(P\right)-\frac{\Delta T(P)}{2}\hfill \\ \hfill & \hfill & +\Delta T(P)\cos \left(\theta \right)\max \left[0,\cos (\varphi -\pi )\right],\hfill \end{array}$$(3) $$\begin{array}{lll}\Delta T(P)\hfill & =\hfill & \text{}\{\begin{array}{ll}\Delta T_0\hfill & \text{if}P<P_{\text{low}}\hfill \\ \Delta T_0\log (P/P_{\text{low}})\hfill & \text{if}{P}_{\text{low}}<P<P_{\text{high}}\hfill \\ 0\hfill & \text{if}P>P_{\text{high}}\hfill \end{array}\text{}\text{},\hfill \end{array}$$(4)

where ΔT₀ is based on observations (see Sect. 2) and remains constant for all models of a particular brown dwarf (see Table 2), P_low = 0.01 bar for all models, and both P_high and $T_{\text{eq-1D}}\left(P\right)$, are model dependent. P_high is simply the pressure at which the day-night temperature difference goes to zero (which we also define as the horizontal convergence pressure – i.e. it is the pressure at which we expect mixing to horizontally homogenise the atmosphere). As for $T_{\text{eq-1D}}\left(P\right),$ not only does this depend upon the horizontal convergence pressure (i.e. P_high), it also depends upon both the brown dwarf in question as well as the presence of any deep radiative dynamics (which we explore via isothermal deep forcing – see Sect. 4.2). The full $T_{\text{eq-1D}}\left(P\right)$ profiles that form the basis for each models individual profile are plotted on the left-hand side of Fig. 2, and a full parametrisation of each models individual $T_{\text{eq-1D}}\left(P\right)$ profile can be found in Table 4.

Likewise, τ_rad is represented by a linear in $\log \left(P\right)$ fit to a 1D equilibrium model for each brown dwarf (see the right-hand side of Fig. 2), although larger differences exist between the profile used in different models here, since adjusting τ_rad allows for a simple means to control the relative strength of the radiative forcing (in the deep atmosphere). In our ‘core’ models, we used an infinite timescale to explore the formation of the adiabat deep profile by potential temperature advection alone, then we added a deep forcing (towards a deep isotherm) and verify how strong it needs to be to destroy the deep adiabatic structure that would be responsible for any observed radius inflation (in the spirit of what has been done in FS19; see Sect. 4.2).

3.3 1D models of brown dwarf atmospheres

The 1D models which are parametrised as part of the Newtonian Cooling approach to radiative dynamics in our 3D models have been calculated using PHOENIX (Hauschildt et al. 1997, 1999; Barman et al. 2001, 2011). PHOENIX is a well-tested self-consistent atmospheric model that has been used to study the atmospheres of stellar and sub-stellar objects for decades, and which has previously been adapted to model the atmosphere of brown dwarfs (e.g. Allard et al. 2011; Lothringer & Casewell 2020).

In brief, PHOENIX works by iteratively calculating pressure-temperature profiles for a given 64-layer P-T structure on an log-spaced optical depth grid that extends from τ = 10⁻⁸ to τ = 10^2.5 (which generally corresponds to pressures of between around 10⁻⁶ to 10² bar). To start, the model first calculates chemical equilibrium using solar metalicity elemental abundances (Asplund et al. 2006) and by assumingthat the atmosphere is in local thermodynamic equilibrium (LTE – i.e. collisional effects dominate the determination of the atomic and molecular compositions). Next, the opacity is calculated for each point using a wavelength grid that extends from 10 to 10⁶ Å with direct-Opacity-Sampling (dOS – Schweitzer et al. 2000). The model includes over 130 different molecular opacity sources, plus both the atomic opacity for elements up to and including uranium, as well as many continuous opacity sources, including collision-induced-absorption and H-opacity. Radiative dynamics are then calculated including any irradiation by the host star which is defined using the closest matching stellar spectrum. We note that for the main sequence stars, we used the PHOENIX models from Husser et al. (2013) and for the white dwarf, SDSS1411, we used the LTE DA white dwarf model from Koester (2010). Finally, this can be used to calculate the vertical flux and update the pressure-temperature structure, using a modified Unsöld-Lucy method (Hauschildt et al. 2003), thus bringing the model slightly closer to radiative-equilibrium. This iterative process continues until the maximum temperature change for any layer is less than 0.5 K.

For each of the brown dwarfs under consideration, we calculated three 1D models: an equilibrium profile, a sub-stellar point profile, and a non-irradiated profile. We note that the latter two of these profiles were used to represent a highly irradiated day-side and a non-irradiated night-side, respectively, in Sect. 4.1, whilst the equilibrium profiles form the basis of the rest of our models.

On top of the external irradiation by the host star (for the equilibrium and sub-stellar point profiles), each 1D model also requires that we set an internal temperature for the brown dwarf. Following Baraffe et al. (1995), we selected suitable warm internal temperatures that should be representative of the uninflated planets and checked for consistency against the models of Phillips et al. (2020). We note that we ended up using a slightly hotter internal temperature for Kepler-13Ab than the aforementioned models suggest, since the model was originally run as part of a test to see if 1D models alone could explain the atmosphere water feature of Kepler-13Ab. However, this should not be a significant problem because, in the way our models are setup, the exact value of the internal temperature has little effect on our primary results. This is because we primarily use the 1D profiles to set the forcing in the outer atmosphere, a region in which T_int has little to no impact due to the local dominance of stellar radiative heating on the atmospheric dynamics. However, we still chose to use suitable values so as to both minimise errors and to ensure the accuracy of the day–night models to calculate the horizontal convergence pressures.

For the brown dwarfs under consideration here, we set T_int = 700 K for KELT-1b, T_int = 1000 K for Kepler-13Ab, and T_int = 800 K for SDSS1411B. Finally, the two irradiated profiles for each brown dwarf include heat redistribution (with f = 0.25 and f = 1.0 for the equilibrium and sub-stellar point profiles respectively). Using these internal temperatures, the planetary parameters (given in Table 2), and suitable stellar spectra, we are then able to calculate the aforementioned T–P and radiative timescale profiles. We plot the resulting equilibrium profiles in Fig. 2 and the full set of 1D models can be found as part of the opendata archive⁴.

4 Results

In order to explore the deep atmospheres of Kepler-13Ab, KELT-1b, and SDSS1411B, we have run a large array of models of each brown dwarf designed to explore not only at what pressure the deep adiabat might be able to develop (Sect. 4.1), but also how stable this advective adiabat is against deep radiative dynamics (driven by Newtonian cooling – Sect. 4.2), as well the link between its formation and the outer and deep atmosphere flows and circulations (Sect. 4.3). The base parameters for each set of brown dwarf models are given in Table 2, and details of each models individual equilibrium temperature (T_eq) and radiative timescale (τ_rad) profiles are given in Tables 3 and 4, respectively.

Fig. 2 Fits to the 1D equilibrium temperature-pressure (Teq – left) and Newtonian cooling relaxation timescale (τrad – right) profiles, which form the basis of the profiles used in our models of Kepler-13Ab (top), KELT-1b (middle), and SDSS1411B (bottom). In all cases, the equilibrium 1D profile is shown in light grey, and the linear in $\log \left(P\right)$ space fit is shown in orange. Note that depending upon the model parameters, including the pressure at which day-night temperaturedifferences go to zero and presence of deep radiative dynamics, these fits may be truncated or modified as required. For full details of the exact fit used in each model, see both Tables 3 and 4.

4.1 Calculating the radiative and advective boundaries for Kepler-13Ab and KELT-1b

When modelling HD 209458b, FS19 had a relatively easy time setting the 3D atmospheric temperature-pressure (T–P) profile of their models thanks to the wealth of data available for this well explored exoplanet (see, for example, Iro et al. 2005; Heng et al. 2011; Mayne et al. 2014b). However this is not the case for the brown dwarfs under consideration here. These objects have undergone much less study and so there exists no consensus on what form the 3D atmospheric T–P profiles should take. Fortunately observations can help to constrain the low pressure surface temperature profile, including the day-night temperature contrast at the top of our models (see Sect. 2). However, the same cannot be said for the deeper regions of the atmosphere, particularly because it is in these regions that we expect advective heat transport to play a critical role, leading to significant changes in the local dynamics when compared to the 1D models. Yet it is the structure of this ‘mid’ atmosphere (~ 0.1 → ~ 10∕100 bar) that is essential to improving our understanding of where a deep adiabat might be able to develop. FS19 showed that the formation of a deep adiabat is strongly linked to the horizontal homogenisation of the deep atmosphere.

In Sect. 4.1.1, we show that by simply adjusting the horizontal convergence pressure (and setting the deep radiative effects to be dynamically unimportant), we are able to form a deep adiabat at almost any location, including non-convective adiabats that are hot enough to explain observed features (such as the radius, etc.; see Sect. 2). However, to do so, we needed to make a number of assumptions about both the deep radiative forcing as well as the 3D temperature structure. In Sect. 4.1.2, we show how we attempted to remove this second assumption by using 1D models of a sub-stellar-point day-side and a non-irradiated night-side to modify our 3D models radiative cooling profile. The aim of this being to constrain the pressure range at which we might expect advection to dominate over radiative effects and hence horizontally homogenise the deep atmosphere.

4.1.1 Freely varying the convergence pressure

We start by exploring models of each brown dwarf with very loose limits on the location of the horizontal convergence pressure. More specifically, we simply require that the convergence pressure should be less than or equal to that required to sustain the corresponding, uninflated 1D models’ convective adiabat, while also being high enough that it remains compatible with observations of the day-night temperature difference.

To illustrate this, in Fig. 3, we plot three exemplary models of KELT-1b with horizontal convergence pressures of 1 bar (top), 10 bar (middle), and 40 bar (bottom). In all three models, the atmospheres are initialised with an adiabat based on the convergence pressure-temperature point given in Table 4.

For models of both KELT-1b and Kepler-13Ab, we find that the deep adiabat then starts to slowly, but significantly, heat. However, the rate of heating slows as we move towards shallower (i.e. lower) horizontal convergence pressures. This results in a smaller fractional increase in the deep adiabat’s temperature, compared to the adiabat used to initialise it: a ~ 3.7% temperature increase at 1 bar, ~4.9% at 10 bar, and ~6.1% at 40 bar. This reduction in the deep heating rate can likely be linked to both the increase in thermal mass of the ‘deep’ atmosphere, as well a reduction in the downward heating rate as the model’s outer atmosphere becomes smaller and the equatorial jet becomes shallower. The latter can be attributed to the equatorial jet being primarily driven by day-night temperaturedifferences.

This sensitivity of the deep heating rate to the depth of the convergence pressure is even more pronounced for SDSS1411B. While deeper convergence pressure models do show signs of very slight deep heating, said heating appears to slow much more quickly as we move towards shallower convergence pressures. We explore how the deep heating rate changes with convergence pressure in more detail in Sect. 4.3.

Given this slightly tighter constraint on the horizontal convergence pressure (i.e. it must be deep enough for a significant, vertical advection driving, equatorial jet to develop), we next explored whether our models for Kepler-13Ab and KELT-1b can maintain advective adiabats that might explain the observed phenomena; for example, the observed water absorption feature of Kepler-13Ab or the inflated radius of KELT-1b. We note that we excluded SDSS1411B from this part of the study since observationssuggest that it should not be significantly inflated and, indeed, the observed radius is consistent with evolutionary models (Littlefair et al. 2014).

For Kepler-13Ab, explaining the observed water absorption features (without invoking an additional thermal inversion) requires at minimum that the atmosphere is adiabatic for all pressures greater than around 0.6 bar. As for KELT-1b, in order to match the observed radius, comparisons with evolutionary and 1D atmospheric models with altered (i.e. increased) internal temperatures suggests that the internal adiabat should develop at pressures less than around 10 bar. This is significantly shallower and, hence, hotter than the convective adiabat that develops at around 40 bar in our reference, radiative-convective, 1D model.

In Fig. 4, we plot (latitudinally to the left and longitudinally to the right) the averaged temperature-pressure profiles at various longitudes or latitudes, respectively, for the models of Kepler-13Ab (top) and KELT-1b (bottom). Here we find that when using a suitable (i.e. shallow) horizontal convergence pressure, both brown dwarf models are able to maintain an advective adiabat that agrees with the observational criteria. That is to say both longitudinally and, to a slightly lesser extent, latitudinally, the deep atmospheres have both converged onto an advective deep adiabat that is hot enough to explain the observed phenomena. We note that the difference in convergence rate between the longitudinal and latitudinal profiles can be linked to a combination of the relative strengths of zonal and meridional flows, and the amount of time that the simulations have been run for. More specifically, zonal dynamics tend to be much stronger than meridional dynamics thanks to the super rotating equatorial jet. As a result, after ~ 200 Earth years of simulation time, the models have converged longitude wise but are still very slowly converging latitude wise. However, since the profiles are already close to being converged at all longitudes, we do not expect this ongoing slow latitudinal convergence to drastically impact the final temperature of the deep advective adiabat.

Even though these deep advective adiabats are a good a match to the observations, two key questions remain before we can ascertain if this result is truly valid. Firstly, we examine wether these are suitable pressures to use when setting the day-night temperature difference to zero. Secondly, we consider what happens when we include radiative forcing in the deep atmosphere.

Fig. 3 Equatorially averaged (i.e. the zonal-mean at the equator) T–P profiles (orange) for three models of KELT-1b with different horizontal convergence pressures (Pconverge), and no radiative forcing, via Newtonian cooling, for P > Pconverge. From top to bottom: the convergence pressure of each model is 1 bar, 10 bar, and 40 bar, respectively. In these plots, the equilibrium Newtonian cooling profile is shown in grey, the day-side profile is shown in red, and the night-side profile in shown in blue.

Fig. 4 Snapshots of the latitudinally (left) or longitudinally (right) averaged T–P profile for two evolved models of Kepler-13Ab (top) and KELT-1b (bottom), both with a very shallow horizontal convergence pressure, at various longitudes or latitudes, respectively. These models have been selected in order to match either the observed radius of the brown dwarf (KELT-1b), or the minimum deep atmosphere temperatures required to explain a water absorption feature in the brown dwarf spectra (Kepler-13Ab); see Sect. 2 for more details.

4.1.2 Using 1D day-night models to constrain the convergence pressure

In order to better constrain the horizontal convergence pressures in our brown dwarf models, we next modify our models to remove any assumptions about where the day-night temperature difference goes to zero.

To do this, we must go beyond using just the equilibrium 1D models and instead look at day-side and night-side profiles explicitly. To that end, we next use a sub-stellar point profile to represent the day-side of our model and a non-irradiated planet to represent the cold night-side. We plot these 1D model profiles for all three of the brown dwarfs under consideration here in Fig. 5. Here, the sub-stellar day-side profiles are shown in shades of red and the non-irradiated night-side profiles in shades of blue. Further each brown dwarfs profiles uses a different line-style: Solid for Kepler-13Ab, dashed for KELT-1b, and dotted for SDSS1411B.

Immediately these 1D profiles reveal a stark difference between Kepler-13Ab and KELT-1b, and SDSS1411B. For the former two brown dwarfs we find that the 1D models do not converge at any pressure. Whereas for SDSS1411B, we find that the sub-stellar day and non-irradiated night 1D models converge at ~ 100 bar. Since we do not expect SDSS1411B to significantly inflate, and hence we do not expect to find significant differences between 1D models and observations, we set 100 bar as the preliminary horizontal convergence pressure for models of SDSS1411B (see Sect. 4.2 for results).

As for Kepler-13Ab and KELT-1b, at first glance it might appear that we can simply use these 1D profiles to set the day-side and night-side Newtonian Cooling forcing profiles in our models and then evolve said models to a steady state. However, doing this would be not only a failure to understand the distinction between the local nature of the (i.e. sub-stellar or non-irradiated) 1D models and the equilibrium nature of Newtonian cooling, but also the fundamental meaning of the term radius inflation: an increase in the observed radius that cannot be explained by 1D models alone. The 1D models lack horizontal advective mixing, such as the super-rotating jet which advectively cools the day-side and heats the night-side. As a result, when we look at the temperature difference between the 1D sub-stellar day and non-irradiated night profiles in the outer atmosphere, we find a day-night contrast that is significantly higher than what has been observed. On the other hand, equilibrium Newtonian cooling requires that (as the name might imply) the relaxation temperature profiles represent the equilibrium temperature state. Thus, if we use the non-advective and non-redistributive T–P profiles from the 1D models with our Newtonian cooling scheme, while we might find that our models develop the required super-rotating jet, any advective cooling or heating on the day or night sides would be eliminated, as the relaxation forces the atmosphere back towards the imposed, non-equilibrium, profiles.

In summary, the use of Newtonian cooling means that our outer atmosphere temperature profile must be prescriptive rather than descriptive. We note that, of course, using a GCM that includes a full description of radiative effects would eliminate this problem. However, the vast increase in computational cost means that such a model was not suitable for this initial exploration of the deep atmospheres of brown dwarfs (see Sect. 5 for a discussion of how next-generation radiative GCMs might be put to work to explore deep brown dwarf atmospheres in future studies).

However, the above does not mean that these 1D sub-stellar day and non-irradiated night profiles are not still useful for constrainingthe pressure at which the atmosphere horizontally converges. To that end, we ran a pair of models of Kepler-13Ab and KELT-1b with a Newtonian cooling profile based on the sub-stellar day and non-irradiated night 1D models in the deep atmosphere, and the equilibrium 1D model modified by the observed day-night temperature difference inthe outer atmosphere. This combined forcing profile allows us to investigate how the deep atmosphere (and, in particular, the advective dynamics) responds to the presence of a (somewhat) physically motivated deep day-night temperature difference. It also ensures that the outer atmosphere dynamics and, hence, the resulting advective downflows, are appropriate for both comparisons with observations and with our ‘core’ equilibrium forcing models.

The resulting combined day (red dashed), night (blue dashed) and ‘equilibrium’ (grey dashed) temperatureforcing profiles can be seen for both Kepler-13Ab (top) and KELT-1b (bottom) as part of Fig. 6. In this figure, we plot the latitudinally (left) and longitudinally (right) averaged T–P profiles at various longitudes or latitude, respectively, for both brown dwarf models. Here, we see that both brown dwarf models show significant signs of horizontal convergence, particularly longitudinally. However latitudinally, the models are still slowly converging, an effect which can once again be linked to the differences in strength between zonal and meridional flows.

An analysis of the longitudinal variation profiles suggests that Kepler-13Ab starts to converge between about 1 and ~ 4–5 bar, and KELT-1b starts to converge between about 10 and ~40–50 bar. The latitudinal variation profiles, at the current simulation time of approximately 300 Earth years, tend to support this, particularly for KELT-1b, where we start to see some convergence around 10 bar, which gets stronger as we shift deeper towards 40 bar.

As discussed and shown above, we have explored a number of models with convergence pressures close to or within this range and shown that a hot, deep adiabat is able to develop. However, one question remains as to whether these advective adiabats are stable when we introduce deep radiative forcing back into the system. We explore this question in the following section.

Fig. 5 1D models of sub-stellar ‘day-side’ (reds) and non-irradiated ‘night-side’ (blues) for the three brown dwarfs under consideration here: Profiles for Kepler-13Ab are plotted with solid lines, with dashed lines used for KELT-1b, and dotted lines for SDSS1411B.

4.2 Sensitivity of deep advective adiabats to the inclusion of deep radiative forcing

At this point, we have a tighter constraint on the pressure at which we might expect the deep atmospheres of Kepler-13Ab, KELT-1b, and SDSS1411B to horizontally converge. In addition, we have confirmed that for models with no deep radiative forcing, all three brown dwarfs can maintain a deep advective adiabat at these convergence pressures. We next explore how stable these advective adiabats are when radiative effects are no longer limited to the outer atmosphere.

To achieve this, we took already evolved models of all three brown dwarfs and explored their evolution when we introduced various strengths of deep radiative forcing (controlled via the radiative forcing timescale) towards an isotherm that has been set to the same temperature as that found at the horizontal convergence pressure. We note that the analysis below is split into two sections representing the two observed dynamical regimes: Kepler-13AB and KELT-1b, which both behave similarly to HD 209458b and show significant signs of deep heating, and SDSS1411B, which appears to struggle to heat the deep atmosphere significantly above the initial convective adiabat from the 1D models, even when radiative effects are limited to the outer atmosphere.

Fig. 6 Snapshots of the latitudinally (left) or longitudinally (right) averaged T–P profile for two evolved models of Kepler-13Ab (top) and KELT-1b (bottom) at various longitudes or latitudes respectively. These models include radiative forcing, via Newtonian cooling, at all pressures. This cooling is based upon a combination of the equilibrium 1D radiative timescale profile (grey) including the observed day-night temperature difference in the outer atmosphere, and 1D sub-stellar and non-irradiated 1D models in the deep atmosphere. The overall result is a combined Newtonian cooling model with a day-side (at the sub-stellar point) forcing profile that corresponds to the red dashed line and a night-side forcing profile that corresponds to the blue dashed line.

4.2.1 Kepler-13Ab and KELT-1b

We start by exploring models of Kepler-13Ab and KELT-1b with horizontal convergence pressures that fall towards the shallower end of the calculated convergence pressures (Sect. 4.1): 1 bar for Kepler-13Ab and 10 bar for KELT-1b.

Figure 7 shows the T–P profiles (left) after ~250 Earth years of simulation time for three models of both Kepler-13Ab (top) and KELT-1b (bottom), each with increasingly long radiative cooling timescales (right). For all deep radiative forcing models shown here, the base model used for initialisation is shown in dark grey and the temperature profile they are being cooled or forced towards is shown in light grey.

Here, we find that, for both brown dwarfs, the models with the strongest deep forcing (labelled A in Fig. 7) that most closely matches the radiative timescales from the 1D radiative-equilibrium models show significant signs of ongoing deep coolingtowards the forced isotherm.

For Kepler-13Ab, as we increase the radiative timescale of the deep forcing (B and C), we find a corresponding decrease in the rate of deep cooling. This effect is clearest in the weakest deep forcing model (C), which is cooling at an increasingly slow rate that suggests that it is asymptotically approaching a steady state in which at least some form of an advective adiabat is maintained. However, this stable deep adiabat is significantly cooler than either the initialisation profile or the convective adiabat that develops in the 1D equilibrium models (or which would develop in these models if we included the physics required for convective motions). This raises some doubts as to whether 1 bar is the correct horizontal convergence pressure for Kepler-13Ab.

As for KELT-1b, due to differences in the radiative timescale parametrisation between the fastest (A) deep forcing model and the slower deep forcing models (B and C – which do not include the 1D models region of near constant τ_rad between 10 and 100 bar), we find that the slower deep forcing models maintain steady-state advective adiabats that are much hotter than found in the strong deep forcing case (A). In particular, the slowest deep forcing model (C) reveals a steady advective adiabat that is only slightly cooler than that of the base model, and thus is both close to the temperature that would be required to explain the radius inflation of KELT-1b (see Fig. 4), as well as being notably hotter than the 1D convective adiabat.

This apparent sensitivity of models with shallower horizontal convergence pressures to deep radiative forcing raises some questions about the suitability of such convergence pressures. As such, we next explore models of Kepler-13Ab and KELT-1b with horizontal convergence pressures that fall towards the deeper end of the calculated convergence pressure range (Sect. 4.1): 4.3 bar for Kepler-13Ab and 40 bar for KELT-1b.

Figure 8 shows the T–P profiles (left) after ~300 Earth years of simulation time for three models of both Kepler-13Ab (top) and KELT-1b (bottom) with increasingly slow deep radiative forcing towards an isotherm (right). However, unlike in the models shown in Fig. 7, none of the deep forcing models shown here exhibit signs of ongoing cooling.

Starting with Kepler-13Ab, we find that all models show at least some signs of cooling from the initial adiabat. However this cooling is limited for models with longer radiative timescales (B and C), with both maintaining steady-state advective adiabats that are at least slightly hotter than the 1D models’ convective adiabat. Only the very strongest deep forcing model (A) shows signs of a significant cooling of the deep adiabat, and even then it remains hotter than the weaker deep forcing model (B) with a 1 bar convergence pressure. This apparent slowdown in the deep cooling rate is likely linked to differences in vertical heat transport. We explore these differences in more detail in Sect. 4.3.

Moving onto KELT-1b, we find that the differences between models with and without deep radiative forcing have become even smaller. Themodels with slower deep forcing (B and C) maintain steady-state advective adiabats that are at least slightly hotter than the 1D models’ convective adiabat. And as for the model with stronger forcing (A – which includes the near constant radiative forcing timescale between approximately 10 and 100 bar), here we also find a steady advective adiabat that is only slightly cooler (by ~200 K at 300 bar – around a 6% difference) than the 1D model.

These results suggest that models of Kepler-13Ab or KELT-1b with a 4.3 bar or 40bar convergence pressure, respectively, should be able to maintain a stable advective adiabat that is hotter than the convective adiabat that a 1D equilibrium model would suggest. That is to say they exhibit at least some form of radius inflation. We note that while all the models above exhibit at least some level of cooling from the initial state,low levels of cooling are not unexpected. Not only was a similar effect noted by FS19 for models of HD 209458b with low levels of deep radiative forcing, but this forcing was designed to be somewhat exaggerated in order to emphasise its effect. After all we do not expect the true equilibrium state of the deep atmosphere to be an isotherm.

However, this does not mean that the trend revealed by the deep forcing is not real, as even when we switch to forcing the deep atmosphere towards the 1D models’ convective adiabat, the models with a deeper convergence pressure are much more stable to the effects of deep radiative dynamics, and the models that are significantly impacted by deep radiative forcing remain so. We explore why this might be the case in Sect. 4.3.

Fig. 7 Snapshots of the T–P profile (left)for a series of test models investigating the effect of deep radiative forcing, via deep isothermal Newtonian cooling, on already evolved (dark grey) models of both Kepler-13Ab (top) and KELT-1b (bottom) with a shallower horizontal convergence pressure. For each brown dwarf, we explore models with increasingly weak deep forcing, as shown by the matching profile on the right (where the full Newtonian cooling timescale profile is found by combining the dark grey profile at lower pressures and the matching coloured profile at higher pressures).

Fig. 8 Snapshots of the T–P profile (left)for a series of test models investigating the effect of deep radiative forcing, via deep isothermal Newtonian cooling, on already evolved (dark grey) models of both Kepler-13Ab (top) and KELT-1b (bottom) with a deeper horizontal convergence pressure. For each brown dwarf, we explored models with increasingly weak deep forcing, as shown by the matching profile on the right, where the full Newtonian cooling timescale profile is found by combining the dark grey profile at lower pressures and the matching coloured profile at higher pressures.

4.2.2 SDSS1411B

Moving on to SDSS1411B, Fig. 9 plots T–P profiles for two 100 bar convergence pressure model atmospheres. We show the base model in dark grey, which, after ~70 Earth years of simulation time has only very slightly warmed from its initial adiabat, and a deep forcing model (orange) which is still steadily cooling (at a near constant rate) after an additional ~ 30 Earth years of simulation time. We note that we only included a single model of SDSS1411B with deep radiative forcing due to similarities in the results. Unless we set the deep radiative timescale to be drastically slower than the 1D models predict, for instance, τ_rad = 10^18→20 at 1000 bar, all deep forcing models show signs of ongoing deep cooling, albeit at steadily slower rates as the deep radiative timescale is lengthened.

We discuss why this might be the case in Sect. 4.3. However, briefly, it appears that downward advective transport is much weaker in SDSS1411B atmosphere than in the atmospheres of Kepler-13Ab and, KELT-1b. This has the effect of making any advective adiabat that develops in SDSS1411B much more sensitive to deep radiative effects.

These results suggest that (as observations require) the formation of a stable hot and deep adiabat is unlikely to lead to any significant inflation of SDSS1411B atmosphere. This can be understood as occurring because not only does advection struggle to heat the deep atmosphere significantly, but even weak radiative effects cool the advective adiabat below the convective adiabat 1D models suggest would develop if we included convective physics in our models.

Fig. 9 Snapshot of the still-cooling T–P profile for a variant of the 100 bar convergence pressure model of SDSS1411B (dark grey) that includes radiative forcing, via Newtonian cooling, at all pressures (orange – deep forcing, DF). This forcing is towards an equilibrium profile in the outer atmosphere (P < Pconverge) and towards an isotherm in the deep atmosphere (P > Pconverge). We note that we only include a single model with deep radiative forcing here since models with even weaker deep cooling behave similarly, just with a different evolution and cooling timescale.

4.3 Differences in flows and circulations between Kepler-13Ab/KELT-1b and SDSS1411B

In Sect. 4.2 we introduced two key results: (1) models of Kepler-13Ab and KELT-1b with a deeper convergence pressure (4.3 and 40 bar respectively) can maintain a hot advective adiabat even in the presence of low levels of deep radiative forcing. At the same time, models of SDSS1411B with a 100 bar convergence pressure struggle to maintain a hot advective adiabat in the presence of deep radiative forcing unless said deep radiative effects are very weak (and even then the deep advective adiabat is not significantly warmed compared to the 1D models’ convective adiabat); (2) models of Kepler-13Ab and KELT-1b with a shallower convergence pressure (i.e. 1 and 10 bar, respectively) appear to have a have harder time maintaining a hot advective adiabat in the presence of deep radiative effects than models with a deeper convergence pressure. Explanations for both of these results can be found by looking at both the flows and circulations that drive the vertical advective heating responsible for the formation of a hot and deep advective adiabat, as well as the associated energy fluxes.

4.3.1 Zonal flows and meridional circulations

We start, in Fig. 10, by exploring the zonally and temporally averaged zonal wind profiles (left) and meridional circulation profiles (right) for models of Kepler-13Ab (4.3 bar – top), KELT-1b (40 bar – middle), and SDSS1411B (100 bar – bottom) with deep convergence pressures. We note that the meridional circulation profile is also known as the streamfunction of the mass-flux on the meridional (i.e. latitudinal-pressure) plane.

Here, we once again find that our models fall into two distinct regimes: Kepler-13Ab and KELT-1b both exhibit zonal wind profiles that are dominated by a broad equatorial super-rotating jet that extends down to around the convergence pressure, and which tends to be flanked by much weaker, and even broader, higher latitude counterflows. Whereas SDSS1411B exhibits a very narrow jet that extends to the bottom of the atmosphere and which is flanked by similarly narrow, weaker, counterflows (which themselves give way to even weaker easterly flows etc.). We note that very similar flow structures for very short period brown dwarfs orbiting white dwarfs were found by Lee et al. (2020) and Tan & Showman (2020).

This difference in the latitudinal extent of both the equatorial jets, as well as the flanking counterflows, can more easily be seen in Fig. 11, which plots slices of the zonally and temporally averaged zonal wind at 0.1 bar (left) and 1 bar (right). Here, we can clearly see that at lower pressures, both the super-rotating easterly jet and the flanking westerly counter-flows of SDSS1411B are significantly narrower than the equivalent flows for either Kepler-13Ab or KELT-1b. Furthermore, we find that the equatorial jet that develops in the atmosphere of Kepler-13Ab is clearly weakening and narrowing as we get closer to the horizontal convergence pressure (and hence the ‘removal’ of the day-night temperature difference that drives the jet). A similar effect was found for the jet of KELT-1b as we approach 40 bar, and to a lesser extent for the narrow jet of SDSS1411B as we approach 100 bar.

Evidence for the dynamical differences between our brown dwarf models can also be seen in the meridional circulation profiles (Fig. 10, right side). As might be expected, due to the super-rotating jet being the dominating driving force behind the flows and circulations in highly irradiated, tidally locked, planets, the meridional circulations profiles once again fall into two distinct regimes: the circulation profiles of both Kepler-13Ab and KELT-1b are very reminiscent of those found for HD209458b (e.g. FS19). We find a clear equatorial downflow, slightly higher-latitude upflows, and weaker polar (or as the case is for Kepler-13Ab, near polar) downflows. These flows can all be linked back to the super-rotating equatorial zonal jet.

For example, in the meridional circulation profile of Kepler-13Ab we can see a clear difference between the near equator circulations in the outer atmosphere and in the deep atmosphere. In the outer atmosphere, the circulations cells extend horizontally to around the same latitude as the equatorial jet, whereas in the deep atmosphere the circulations become much narrower as rather than being driven flows associated with the zonal winds, they are being driven by the adiabatic redistribution of thermal energy that the outer atmosphere circulations have advected downwardly.

However, the best example of how significantly the zonal winds can affect the meridional circulations is SDSS1411B. Here, we see an almost completely different circulation profile than that found in the models of Kepler-13Ab, KELT-1b, or even HD 209458b. Instead, it is much more similar to the profile found by Tan & Showman (2020) for a rapidly rotating hot Jupiter atmosphere. The narrow equatorial jet, along with the dynamics that drive it, lead to a meridional circulation profile that is dominated by relatively weak narrow circulations cells, which do drive a downflow at the equator, while also driving significant (cooling) low latitude upflows. This circulation structure may help to explain why SDSS1411B struggles to maintain an advective adiabat even when very weak radiative effects are included in the deep atmosphere (see Sect. 4.2).

Fig. 10 Zonally and temporally averaged (over a period of ~3 yr) zonal wind profiles (left) and meridional circulation streamfunctions (right) for Kepler-13Ab (top), KELT-1b (middle), and SDSS1411B (bottom). In the zonal wind profiles easterly winds are positive and westerly winds are negative, and in the meridional circulation profile, we plot the stream function using a logged scale in order to clearly illustrate the full circulation profile. Here, clockwise circulations are shown in red and anti-clockwise in blue – these circulations combine in all models to reveal an equatorial downflow at all pressures.

4.3.2 Influence of rotation on the flows and circulations

Thus, we consider what is driving the above differences in flows and circulations. One possible explanation, other than the day-night temperature contrast, is the relative influence of rotational effects when compared to advective dynamics.

To explore this balance, we make use the of the Rossby number. This is a measure of the relative strength of advective forces to Coriolis accelerations: $$R_0=\frac{U_z}{2\Omega \sin \left(\theta \right)L},$$(5)

where U_z is the zonal velocity, Ω is the angular rotation rate, and L is the characteristic length scale of the flows. We note that we used the zonal velocity since we are interested in how rotational effects affect the zonal flows, and we approximated the planetary radius as the characteristic length scale (L = R_BD) since we are interested in the global scale effect. When the Rossby number is very small (Ro ≪ 1), rotational effects dominate over the advective dynamics; whereas when the Rossby number is large (Ro ≫ 1), rotational effects play a minor role in the dynamical balance.

Calculating the Rossby number for our three main deep-convergence brown dwarf models reveals that, as expected, the Rossby number near the equator (i.e. in the jet) tends to be significantly greater than 1: at 0.01 bar, R₀ = 16.6 for Kepler-13Ab, R₀ = 28.8 for KELT-1b, and R₀ = 22.9 for SDSS1411B. However, the Rossby number then drops rapidly as we move to higher latitudes, especially for SDSS1411B: at mid latitudes (i.e. around 45°), typical Rossby numbers are R₀ = 0.127 for Kepler-13Ab, R₀ = 0.130 for KELT-1b, and R₀ = 0.012 for SDSS15411b. The latter is similar to off-equator Rossby numbers found by Tan & Showman (2020) for a rapidly rotating HD 209458b-like model, and by Lee et al. (2020) in models of WD0137-349B, another brown Dwarf that is in a close orbit around a white dwarf.

This suggests that off-equator flows are much more susceptible to rotational effects in SDSS1411B than in either Kepler-13Ab or KELT-1b, which may help to explain the very different zonal wind and meridional circulations profiles seen here. We intend to explore the effects of changing rotation rate on the atmospheric dynamics, and the resulting deep steady-state, for highly irradiated exoplanets in much more detail as part of a future study.

Fig. 11 Temporally and zonally averaged, over a period of ~3 yr, profiles of the zonal wind at 0.1 and 1 bar for deep horizontal convergence pressure models of Kepler-13Ab (green), KELT-1b (blue), and SDSS1411B (red). Here, easterly zonal winds are positive and westerly zonal winds are negative.

4.3.3 Exploring the vertical enthalpy and radiative fluxes

To understand how the differences in flows and circulations affect the formation of a deep advective adiabat, we next explore the vertical advective heating rate and compare it to the deep radiative flux. To achieve this, we turn to the vertical enthalpy flux $$F_H=\rho c_{p}T{U}_r,$$(6)

where U_r is the vertical velocity component. For KELT-1b, we find a mean vertical enthalpy flux of − 1.11 × 10⁹ erg s⁻¹ cm⁻², which is significantly higher than the average deep radiative flux of 5.4 × 10⁶ erg s⁻¹ cm⁻². Similarly, for Kepler-13Ab, we find a mean vertical enthalpy flux of − 2.68 × 10⁹ erg s⁻¹ cm⁻², which is again significantly higher than the average deep radiative flux of 2.86 × 10⁷ erg s⁻¹ cm⁻². However, the same cannot be said for SDSS1411B, where we find a mean vertical enthalpy flux of − 4.92 × 10⁶ erg s⁻¹ cm⁻², which is almost equivalent to average deep radiative flux of 2.49 × 10⁶ erg s⁻¹ cm⁻². We note that for all comparisons above, the deep radiative flux is taken from a 1D radiative-convective model that closely matches the 3D models’ hot advective adiabat.

This explains both why the deep adiabat in our 100 bar convergence pressure model of SDSS1411B is so slow to heat and why it is so sensitive to deep radiative effects, namely: the highly rotationally influenced flows that develop in SDSS1411B lead to a significantly reduced vertical advective heating rate. The result is a heating rate that is too slow to significantly heat the deep atmosphere and, hence, to cause radius inflation.

The vertical enthalpy flux also helps to explain why the shallower convergence pressure models of Kepler-13Ab and KELT-1b have been shown to be more sensitive to deep radiative effects than models run at a deeper convergence pressures. For both Kepler-13Ab and KELT-1b, not only does the slightly shallower jet drive slightly weaker equatorial downflows, but the shallower convergence pressure means that the average deep radiative flux, particularly near the convergence pressure, is higher. As such, it is harder for the model atmospheres to maintain an advective adiabat when deep radiative effects are included. For KELT-1b, our 10 bar convergence pressure model has a mean vertical enthalpy flux of − 1.07 × 10⁸ erg s⁻¹ cm⁻² and an average deep radiative flux of 1.86 × 10⁷ erg s⁻¹ cm⁻². Similarly, for Kepler-13Ab, our 1 bar convergence pressure model has a mean vertical enthalpy flux of − 2.21 × 10⁹ erg s⁻¹ cm⁻² and an average deep radiative flux of ~1.1 × 10⁸ erg s⁻¹ cm⁻².

A similar effect has been found for test models of SDSS1411B with shallower horizontal convergence pressures, further reinforcing our analysis and conclusion that SDSS1411B should not be significantly inflated relative to 1D models.

5 Discussion and conclusions

In this work, we present a series of models of Kepler-13Ab, KELT-1b, and SDSS1411B, designed with the explicit purpose of exploring whether the observed radius inflation relation for brown dwarfs can be understood using the same mechanism that was shown to be applicable to hot Jupiters (FS19; PT17). That is the vertical advection of potential temperature from a highly irradiated outer atmosphere to the deep atmosphere, leading to the formation of a hot and deep advective adiabat that is at a higher entropy than 1D ‘radiative-convective’ models suggest (and, hence, is inflated relative to the 1D models).

However, compared to the hot Jupiter explored by FS19 (HD 209458b), the thermal structures of the outer atmospheres of the brown dwarfs under consideration here are not as well-documented. This presented the first major hurdle to this work as an accurate parametrisation of the outer atmosphere temperature-pressure (T–P) profile, including the day-night temperature difference at all pressures, is key to the Newtonian cooling approach (see Sect. 3) to radiative dynamics used here. An approach that remains necessary as the next generation of high performance, exascale, radiative GCMs, such as DYNAMICO-LMDZ⁵ or LFRic⁶, are still in active development.

Although the observations are able to explore the day-night temperature contrast at very low pressures, they do not provide details about the full vertical structure of the T–P profile. In particular, observations do not help to constrain the horizontal convergence pressure at which horizontal homogenisation by advective flows and circulations has reduced the day-night temperature difference to zero. In Sect. 4.1.1, we demonstrate that freely varying the horizontal convergence pressure in models with no deep (i.e. P >P_converge) radiative forcing allows us to form an advective adiabat at almost any pressure (e.g. Fig. 3). This includes advective adiabats that were hot enough to explain either the observed radius or observed atmospheric features (Sect. 2) of KELT-1b or Kepler-13Ab respectively (see Fig. 4).

Nonetheless, while the above models easily allow us to reproduce observations, they do not reveal any information about the true properties of the deep atmospheres of our brown dwarfs. After all, we can generate almost any deep advective adiabat in these no deep radiative forcing models just by adjusting the horizontal convergence pressure, and we cannot just run a single shallow convergence pressure model with deep radiative forcing because (as we discuss in Sect. 4.3) the location of the day-night convergence affects the deep heating rate.

In order to try and resolve this issue, in Sect. 4.1.2 we next explored a series of models designed to remove any assumptionsabout where the day-side and night-side temperature profiles converge. More specifically, we ran a series of models with Newtonian cooling profiles based on the combination of an equilibrium 1D model paired with the observed day/night temperature difference inthe outer atmosphere, and a deep atmosphere profile which used a sub-stellar-point 1D model to represent its day-side and a non-irradiated 1D model to represent its night-side. For SDSS1411B, both the 1D models, as well as our 3D models, confirm that the deep atmosphere becomes both latitudinally and longitudinally converged at around 100 bar. that is, at essentially the same pressure that a convective adiabat is expected to develop in all of its fixed T_int 1D models. On the other hand, while the 1D models of both Kepler-13Ab and KELT-1b (Fig. 6) do not converge at any pressure, the same cannot be said for the 3D models. These models reveal that the atmosphere starts to converge longitudinally (and more slowly longitudinally) at 1 bar and 10 bar for Kepler-13Ab and KELT-1b, respectively, with the profile almost completely (longitudinally) converging at around 4.3 bar and 40 bar, respectively.We note, however, that although these models can help us to identify where the deep atmosphere convergences, the use of a somewhat unphysical deep Newtonian cooling profile with an imposed deep horizontal temperature gradients prevents us from using them to explore radius inflation in general.

As such, given that we already knew that advective adiabats could form within these convergence pressure ranges when we excluded radiative dynamics from the deep atmosphere (Sect. 4.1.1), we next explored (in Sect. 4.2) what happens when this truncation is no longer in effect. To that end, we ran a series of models of all three brown dwarfs with horizontal convergence pressures based those found above, outer atmosphere forcing based on 1D equilibrium models and the day-night temperature difference, along with varying levels of deep radiative forcing towards an isotherm.

Starting with Kepler-13Ab, we found that even as we vary the strength of the forcing, all test models with a shallower convergence pressure (1 bar) are highly sensitive to the presence of deep radiative effects. However, this is not the case for models with a slightly deeper convergence pressure (4.3 bar). Here, while some cooling compared to the no-deep-forcing models does occur, the majority of our models maintain advective adiabats that are hotter than the convective adiabats found in our 1D models. Furthermore, it is important to note that this deep isothermal forcing is likely to over-exaggerate the impact of radiative cooling on the deep atmosphere – test models with forcing towards the 1D models’ convective adiabat revealed the same trend seen above, albeit to a much weaker extent, which is harder to quantify. As such, it is clear that potential temperature advection from the highly irradiated outer atmosphere to the deep atmosphere can provide an explanation of at least some of the observed radius inflation of Kepler-13Ab. Furthermore it may, depending upon the exact nature of the deep radiative dynamics, be able to explain the observed atmospheric water features (Beatty et al. 2017).

A similar story holds true for KELT-1b. Here, models with a shallower convergence pressure (10 bar) fall into one of two regimes depending upon the exact form that the radiative timescale profile takes. Models with forcing profiles that include the region of near-constant τ_rad between ~10 and ~100 bar (Fig. 2) generally exhibit significant deep cooling resulting in an advective adiabat that is notably cooler than the 1D models’ convective adiabat. On the other hand, model instead extrapolate the 1 bar timescale profile down to 1000 bar tend to exhibit much weaker cooling. This results in advective adiabats that remain hot enough to almost fully explain the observed radius inflation. Furthermore, moving onto models with a deeper convergence pressure (40 bar), we found that the importance of the form that the radiative timescale profile takes has been diminished. Instead, all the 40 bar deep forcing models explored here retain advective adiabats that are, at most, only slightly cooled compared to their initial state. Thus, as is the case for Kepler-13Ab above, it is clear that potential temperature advection can explain at least some of the observed radius inflation of KELT-1b. And in the right circumstances can explain all of the observed inflation.

Finally we looked at the stability of an advective adiabat in a model of SDSS1411B with a 100 bar convergence pressure. Here we found that, for all but the weakest of deep radiative forcing, any initial adiabat rapidly cools such that, after only 10 Earth years of simulation time, the deep, isothermally forced, T–P profile is already significantly cooler than a 1D model’s convective adiabat. This suggests that the vertical advection of potential temperaturedoes not lead to any significant radius inflation for SDSS1411B. This result with agrees with both the observations of SDSS1411B and brown dwarfs orbiting white dwarfs in general (Fig. 1 and Sect. 2).

To answer the question of why this difference in deep advective heating occurs, we next looked at the flows and circulations that drive this heating (Sect. 4.3). Here, we once again found that our brown dwarf models fall into one of two distinct regimes: the zonal wind and meridional circulation profiles of both Kepler-13Ab and KELT-1b behave very similarly to that of HD209458; that is to say, we find a super-rotating zonal jet that drives a strong equatorial downflow which, in turn, drives a strong downward enthalpy flux (− 1.11 × 10⁹ erg s⁻¹ cm⁻² for KELT-1b and − 2.68 × 10⁹ erg s⁻¹ cm⁻² for Kepler-13Ab).

On the other hand, SDSS1411B reveals a very different jet structure that consists of narrow alternating flows that weaken as they move to high latitudes. This zonal wind profile drives a very different circulation profile than seen in either Kepler-13Ab or KELT-1b, which consists of very horizontally-narrow and vertically-extended circulation cells that drive a weak net downward enthalpy flux of a strength on the order of the deep radiative flux (− 4.92 × 10⁶ erg s⁻¹ cm⁻² versus 2.49 × 10⁶ erg s⁻¹ cm⁻²), thus reducing the ability of advection to overcome radiative effects and heat the deep atmosphere of SDSS1411B.

These differences in the flows and circulations are likely linked to the relative strength of mid to high latitude rotational effects in the two regimes. Analyses of the zonal flows reveal that the mid-latitude Rossby number of SDSS1411B (R₀ = 0.012) is an order of magnitude smaller than that of either KELT-1b or Kepler-13Ab (R₀ ~ 0.13), suggesting that, for SDSS1411B, off-equator flows are significantly affected by the Coriolis effect, reducing the relative efficiency of advective heat transport.

A similar affect, albeit with a slightly different underlying cause, is also observed in the convective zones of stars: as the rotation rate is changed, the influence of rotation on the angular momentum flux balance changes. As a result, the meridional circulation profile shifts from a single cell in each hemisphere (slower, anti-solar, rotation) to multiple-cells aligned with the rotation axis (faster, solar-like, rotation – here, cells become aligned with the rotation axis thanks to the effects of convection and the Taylor-Proudman constraint – Taylor & Lamb 1917; Proudman & Lamb 1916). This means that there is a similar transition in circulation dynamics at play as what is seen here; see, for example, Miesch (2005) or Miesch et al. (2011) for a review of meridional circulation in solar convective zones.

The above findings suggest that the differences in flows and circulations that likely occur due to the differences in rotational influence between the slower rotating brown dwarfs with a main-sequence host star and the much more rapidly rotating brown dwarf orbiting a white dwarf can help to explain why the former show significant radius inflation and the latter do not. This simple, physical, explanation for the differences in brown dwarf inflation rates reinforces the idea that irradiation-induced and rotationally influenced advection of potential temperature appears to be a robust mechanism for resolving the radius-inflation problem in both hot Jupiters (FS19) and brown dwarfs.

To further confirm that this is the case, we suggest that this work should be followed by two critical future studies. The first should explore how rotational effects influence the vertical advection of potential temperature and verify that this mechanism shuts down (or at least significantly slows down) in highly rotationally influenced atmospheres. The second should, once the next generation of high performance (i.e. exascale) radiative GCMs are available, explore the steady state deep atmospheres of brown dwarfs models that include self-consistent radiative dynamics (and, hence, make fewer assumptions about the outer atmospheres’ temperature structure, including the convergence pressure and the day-night temperature difference) and, thus, are able to obtain a more reliable measurement of how much of the observed radius inflation can be adequately explained by advective heating. Additionally we suggest that WD1032b (Casewell et al. 2020a) should be included in future studies of the radius inflation of brown dwarfs in order to investigate whether the slightly lower irradiation (cooler host star) and rotational influence (slower rotation) in the atmosphere can explain why this brown dwarfs orbiting a white dwarf might exhibit radius inflation when SDSS1411B does not.

Acknowledgements

F. Sainsbury-Martinez and P. Tremblin would like to acknowledge and thank the ERC for funding this work under the Horizon 2020 program project ATMO (ID: 757858). The authors also wish to thank Idris, CNRS, and MDLS for access to the supercomputer Poincare, without which the long-timescale calculations featured in this work would not have been possible. Additionally this work was granted access to the HPC resources of IDRIS (Jean-Zay) and CEA-TGCC (Irene/Joliot-Curie) under the 2020/2021 allocation – A0080410870 made as part of the GENCI Dari A8 call. Finally the authors with to thank the referee (and editor) for useful comments, questions, and suggestions.

References

All Tables

All Figures

Fig. 1 Mass-radius relation for 25 of the known highly irradiated brown dwarfs transiting either a main sequence star (Carmichael et al. 2020 and references therein) or a white dwarf (Casewell et al. 2020a; Littlefair et al. 2014; Parsons et al. 2017), with the brown dwarfs under consideration here clearly labelled. Additionally, the effective temperatures of the brown dwarf hosting stars are indicated by the colour of each brown dwarfs circle, while the size of each coloured circle is proportional to the total incident surface irradiation. Also shown are the Sonora-Bobcat evolutionary models of Marley et al. (2018) for 100 Myr and 600 Myr at only solar metalicity, and at 2 Gyr, 6 Gyr, and 10 Gyr for both solar metallicities (grey) and low metallicities ([Fe/H] =−0.5 – black).

In the text

Fig. 2 Fits to the 1D equilibrium temperature-pressure (Teq – left) and Newtonian cooling relaxation timescale (τrad – right) profiles, which form the basis of the profiles used in our models of Kepler-13Ab (top), KELT-1b (middle), and SDSS1411B (bottom). In all cases, the equilibrium 1D profile is shown in light grey, and the linear in $\log \left(P\right)$ space fit is shown in orange. Note that depending upon the model parameters, including the pressure at which day-night temperaturedifferences go to zero and presence of deep radiative dynamics, these fits may be truncated or modified as required. For full details of the exact fit used in each model, see both Tables 3 and 4.

In the text

Fig. 3 Equatorially averaged (i.e. the zonal-mean at the equator) T–P profiles (orange) for three models of KELT-1b with different horizontal convergence pressures (Pconverge), and no radiative forcing, via Newtonian cooling, for P > Pconverge. From top to bottom: the convergence pressure of each model is 1 bar, 10 bar, and 40 bar, respectively. In these plots, the equilibrium Newtonian cooling profile is shown in grey, the day-side profile is shown in red, and the night-side profile in shown in blue.

In the text

Fig. 4 Snapshots of the latitudinally (left) or longitudinally (right) averaged T–P profile for two evolved models of Kepler-13Ab (top) and KELT-1b (bottom), both with a very shallow horizontal convergence pressure, at various longitudes or latitudes, respectively. These models have been selected in order to match either the observed radius of the brown dwarf (KELT-1b), or the minimum deep atmosphere temperatures required to explain a water absorption feature in the brown dwarf spectra (Kepler-13Ab); see Sect. 2 for more details.

In the text

	Fig. 5 1D models of sub-stellar ‘day-side’ (reds) and non-irradiated ‘night-side’ (blues) for the three brown dwarfs under consideration here: Profiles for Kepler-13Ab are plotted with solid lines, with dashed lines used for KELT-1b, and dotted lines for SDSS1411B.
In the text

Fig. 6 Snapshots of the latitudinally (left) or longitudinally (right) averaged T–P profile for two evolved models of Kepler-13Ab (top) and KELT-1b (bottom) at various longitudes or latitudes respectively. These models include radiative forcing, via Newtonian cooling, at all pressures. This cooling is based upon a combination of the equilibrium 1D radiative timescale profile (grey) including the observed day-night temperature difference in the outer atmosphere, and 1D sub-stellar and non-irradiated 1D models in the deep atmosphere. The overall result is a combined Newtonian cooling model with a day-side (at the sub-stellar point) forcing profile that corresponds to the red dashed line and a night-side forcing profile that corresponds to the blue dashed line.

In the text

Fig. 7 Snapshots of the T–P profile (left)for a series of test models investigating the effect of deep radiative forcing, via deep isothermal Newtonian cooling, on already evolved (dark grey) models of both Kepler-13Ab (top) and KELT-1b (bottom) with a shallower horizontal convergence pressure. For each brown dwarf, we explore models with increasingly weak deep forcing, as shown by the matching profile on the right (where the full Newtonian cooling timescale profile is found by combining the dark grey profile at lower pressures and the matching coloured profile at higher pressures).

In the text

Fig. 8 Snapshots of the T–P profile (left)for a series of test models investigating the effect of deep radiative forcing, via deep isothermal Newtonian cooling, on already evolved (dark grey) models of both Kepler-13Ab (top) and KELT-1b (bottom) with a deeper horizontal convergence pressure. For each brown dwarf, we explored models with increasingly weak deep forcing, as shown by the matching profile on the right, where the full Newtonian cooling timescale profile is found by combining the dark grey profile at lower pressures and the matching coloured profile at higher pressures.

In the text

Fig. 9 Snapshot of the still-cooling T–P profile for a variant of the 100 bar convergence pressure model of SDSS1411B (dark grey) that includes radiative forcing, via Newtonian cooling, at all pressures (orange – deep forcing, DF). This forcing is towards an equilibrium profile in the outer atmosphere (P < Pconverge) and towards an isotherm in the deep atmosphere (P > Pconverge). We note that we only include a single model with deep radiative forcing here since models with even weaker deep cooling behave similarly, just with a different evolution and cooling timescale.

In the text

Fig. 10 Zonally and temporally averaged (over a period of ~3 yr) zonal wind profiles (left) and meridional circulation streamfunctions (right) for Kepler-13Ab (top), KELT-1b (middle), and SDSS1411B (bottom). In the zonal wind profiles easterly winds are positive and westerly winds are negative, and in the meridional circulation profile, we plot the stream function using a logged scale in order to clearly illustrate the full circulation profile. Here, clockwise circulations are shown in red and anti-clockwise in blue – these circulations combine in all models to reveal an equatorial downflow at all pressures.

In the text

	Fig. 11 Temporally and zonally averaged, over a period of ~3 yr, profiles of the zonal wind at 0.1 and 1 bar for deep horizontal convergence pressure models of Kepler-13Ab (green), KELT-1b (blue), and SDSS1411B (red). Here, easterly zonal winds are positive and westerly zonal winds are negative.
In the text