Early Time Small-scale Structures in Hot Exoplanet Atmosphere Simulations

A hallmark of nonlinear dynamical systems is their sensitivity to initial condition (H. Poincaré 1914; E. N. Lorenz 1963). In such systems, infinitesimal perturbations at early times are quickly amplified during their evolution, leading to a loss of predictability in certain variables and measures (E. N. Lorenz 1964). This phenomenon, often referred to as the "butterfly effect," lies at the heart of chaos and underscores the inherent unpredictability of many complex systems. As a highly nonlinear system, the hot exoplanet atmosphere also exhibits this paradigmatic feature. Indeed, numerical simulations of hot exoplanet atmospheres are sensitive to their initial states and their ability to represent flows across an adequate range of dynamically significant scales (e.g., J. Y-K. Cho et al. 2008, 2015; H. T. Thrastarson & J. Y-K. Cho 2010; J. W. Skinner & J. Y-K. Cho 2021).

Early efforts to simulate hot, 1:1 spin–orbit synchronized exoplanet atmospheres have utilized a very simple setup for the initial and forcing conditions (e.g., A. P. Showman & T. Guillot 2002; J. Y-K. Cho et al. 2003, 2008; C. S. Cooper & A. P. Showman 2006; I. Dobbs-Dixon & D. N. C. Lin 2008; A. P. Showman et al. 2008; K. Menou & E. Rauscher 2009; E. Rauscher & K. Menou 2010). In this setup, an atmosphere initially at rest is set in motion by "relaxing" the temperature field to a prescribed spatial distribution, on a specified timescale. Despite its simplicity, the setup is useful and provides valuable insights—especially in the absence of detailed information about the atmosphere. For this reason, simple setups continue to be utilized in modeling studies today (e.g., F. Debras et al. 2020; V. G. A. Böning et al. 2024).

A salient feature that can be studied with the simple setup for hot exoplanets is the wide range of thermal relaxation timescales in the currently observable region of their atmospheres. In particular, on hot exoplanets the timescale can be very short (i.e., much shorter than the advective timescale) above the ∼10⁵ Pa altitude level. However, it has long been known that such quick adjustments from rest lead to a flow which is very sensitive to a wide range of modeling parameters (e.g., J. Y-K. Cho et al. 2003, 2015; H. T. Thrastarson & J. Y-K. Cho 2010, 2011; I. Polichtchouk & J. Y-K. Cho 2012). For example, the numerical resolution, initial flow state, thermal relaxation timescale, strength and form of numerical dissipation, and altitude of peak heat absorption all affect the predictions (e.g., J. Y-K. Cho et al. 2003, 2008; K. Heng et al. 2011; I. Polichtchouk et al. 2014; J. W. Skinner & J. Y-K. Cho 2021; M. Hammond & D. S. Abbot 2022; J. W. Skinner et al. 2023).

In this paper, we highlight the profound effect of small-scale structures that arise at early times on the late-time flow. This "early time sensitivity" has not been explicitly called to attention before. Past studies have generally focused on the state of the flow a long time after the start of the simulation (often referred to as the "equilibrated" state),³ generally employing only low to moderate resolution. Little focus has been given to the transient, small-scale flow structures that occur during the first ∼10 days of the simulation. The tacit assumption in the past has been that the flow evolution would "forget" the initial condition and head inexorably to the same statistically steady state. Here simulations are performed at high resolution with low dissipation (to be elucidated below) to more accurately capture the small-scale dynamics than have been done in the past.

The governing equations, planetary parameters, numerical model, and physical setup in this work are same as those in J. Y-K. Cho et al. (2021) and J. W. Skinner & J. Y-K. Cho (2022). Therefore, only a brief summary is presented here; we refer the reader to the above works for more details—as well as to J. W. Skinner & J. Y-K. Cho (2021), I. Polichtchouk et al. (2014), and J. Y-K. Cho et al. (2015) for extensive convergence tests and intermodel comparisons. As in the aforementioned works, the hydrostatic primitive equations (PEs) are solved here to study the three-dimensional (3D) atmospheric dynamics. The dissipative PEs, with pressure p serving as the vertical coordinate, read:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \zeta }{\partial t}={\boldsymbol{k}}\cdot {\boldsymbol{\nabla }}\times {\boldsymbol{n}}+{{ \mathcal D }}_{\zeta },\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \delta }{\partial t}={\boldsymbol{\nabla }}\cdot {\boldsymbol{n}}-{\nabla }^{2}\left(\displaystyle \frac{1}{2}{{\boldsymbol{v}}}^{2}\,+\,{\rm{\Phi }}\right)+{{ \mathcal D }}_{\delta },\end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}{\rm{\Theta }}}{{\rm{\partial }}t}=-{\boldsymbol{\nabla }}\cdot ({\rm{\Theta }}\,{\boldsymbol{v}})-\displaystyle \frac{{\rm{\partial }}\,}{{\rm{\partial }}p}(\omega \,{\rm{\Theta }})+\displaystyle \frac{\dot{{ \mathcal Q }}}{{\rm{\Pi }}}\,+{{ \mathcal D }}_{{\rm{\Theta }}},\end{eqnarray} \tag{ 1c }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}{\rm{\Phi }}}{{\rm{\partial }}{\rm{\Pi }}}=-{\rm{\Theta }},\end{eqnarray} \tag{ 1d }$$

$$\begin{eqnarray}&&\frac{\partial \omega }{\partial p}=-\delta .\end{eqnarray} \tag{ 1e }$$

In Equations (1), ζ(x, t) ≡ k · ∇ × v is the vorticity and δ(x, t) ≡ ∇ · v is the divergence, where ${\boldsymbol{x}}\in {{\rm{{\mathbb{R}}}}}^{3}$, k is the local vertical direction, v(x, t) is the horizontal velocity, and ∇ is the gradient, which operates along an isobaric (constant p) surface; Θ(x, t) ≡ (c_p/Π)T is the potential temperature, where c_p is the constant specific heat at constant p, ${\rm{\Pi }}\equiv {c}_{p}{(p/{p}_{\,\rm{ref}\,})}^{\kappa }$ is the Exner function with $\kappa \equiv { \mathcal R }/{c}_{p}$, ${ \mathcal R }$ is the specific gas constant, p_ref is a reference p, and T(x, t) is the temperature; Φ(x, t) ≡ gz is the specific geopotential, where g is a constant and z(x, t) is the height; ω ≡ Dp/Dt is the vertical velocity, where D/Dt ≡ ∂/∂t + v · ∇ + ω∂/∂p is the material derivative. n ≡ −(ζ + f)k × v − δv − ∂(ωv)/∂p, where $f\equiv 2{\rm{\Omega }}\sin \phi$ is the Coriolis parameter with Ω the planetary rotation rate and ϕ the latitude; $\rho =p/({ \mathcal R }\,T)$ is the density; ${{ \mathcal D }}_{\zeta }(\nu ,{\mathfrak{p}})$, ${{ \mathcal D }}_{\delta }(\nu ,{\mathfrak{p}})$, and ${{ \mathcal D }}_{{\rm{\Theta }}}(\nu ,{\mathfrak{p}})$ are the (hyper)dissipations, which are dependent on the dissipation coefficient $\nu =\nu ({\mathfrak{p}})$ and order ${\mathfrak{p}}\in {{\mathbb{Z}}}^{+}$; and $\dot{{ \mathcal Q }}({\boldsymbol{x}},t;{\tau }_{r})$ is the net heating rate, where τ_r is the relaxation time parameter. The boundary conditions in this work are free-slip (i.e., ω = 0) at the top and bottom isobaric surfaces and periodic in the zonal (longitudinal) direction. Throughout this paper, the lateral coordinates are (longitude, latitude) = (λ, ϕ); time, length, pressure, and temperature are expressed in units of planetary day (τ = 3 × 10⁵ s), planetary radius (R_p = 10⁸ m), reference pressure (p_ref = 10⁵ Pa), and reference temperature (T_ref = 1500 K), respectively.

The ζ–δ–Θ formulation of the PEs in Equation (1) facilitates the use of the pseudospectral method to solve the equations accurately. Unlike other numerical methods which offer algebraic convergence (e.g., finite difference or finite element methods), the spectral method offers exponential convergence—i.e., the error decays exponentially fast with increased resolution; this leads to dramatically improved accuracy for the same or similar computational cost (e.g., J. P. Boyd 2000). In applying the spectral method to solve Equations (1), the mapping ${\boldsymbol{v}}\mapsto {\boldsymbol{v}}\cos \phi \,$ is employed because the components of v are not well suited for representation in scalar spectral expansions (A. J. Robert 1966). The formulation of PEs in p vertical coordinates also offers practical simplifications of the equations as well as clarity of presentation; a second-order finite difference scheme is used for the p-direction.

The resolution of the numerical simulations performed here is up to T341L50—corresponding to 341 total wavenumbers and 341 zonal wavenumbers in the Legendre expansion of the variables, and 50 linearly spaced vertical levels (layers) spanning the range p ∈ [0.1, 5.0] bar. Note that the simulations here may not be numerically converged to those that span a much larger p-range, especially if the density of layers is much greater than the simulations here (J. W. Skinner & J. Y-K. Cho 2022). However, the focus here is on a fundamental feature stemming from nonlinearity—acute sensitivity to small-scale structures under strong large-scale forcing. We have verified that the current resolution is adequate to lucidly and robustly demonstrate the highlighted feature.⁴ Additionally, while a variety of p-ranges have been used in past simulation studies (see, e.g., J. Y-K. Cho et al. 2003, 2015; A. P. Showman et al. 2008; E. Rauscher & K. Menou 2010; X. Tan & T. D. Komacek 2019; J. M. Mendonça 2020; J. W. Skinner & J. Y-K. Cho 2022), the range here is chosen to cover the majority of the thermally irradiated levels while permitting the highlighted dynamical effect to be demonstrated clearly.

For the time integration, the second-order leapfrog scheme is used with a time step size of Δt = 4.0 × 10⁻⁵. The Robert–Asselin filter (A. J. Robert 1966; R. Asselin 1972) with a small filter value of  = 0.02 (H. T. Thrastarson & J. Y-K. Cho 2011) is applied to suppress the growth of the computational mode arising from using the leapfrog scheme to integrate first-order (in time) equations. All simulations are initialized from rest (i.e., v = 0) and evolve under the prescribed thermal forcing $\dot{{ \mathcal Q }}({\boldsymbol{x}},t;\,{\tau }_{r})$ (see, e.g., J. W. Skinner & J. Y-K. Cho 2022, Figure 1). Equations (1) are integrated to t = 1000, much longer than the significant dynamical timescales (e.g., advective, forcing, and small-scale dissipation timescales).

The only parameters varied in the simulations presented here are ${\mathfrak{p}}$ and ${\nu }_{2{\mathfrak{p}}}$ in the hyperdissipations:

$$\begin{eqnarray}&&{{ \mathcal D }}_{\chi }(\nu ,{\mathfrak{p}})={\nu }_{2{\mathfrak{p}}}\left[{(-1)}^{{\mathfrak{p}}+1}{\nabla }_{p}^{2{\mathfrak{p}}}+{ \mathcal C }\right]\,\chi ,\end{eqnarray} \tag{ 2 }$$

where χ ∈ {ζ, δ, Θ} and ${ \mathcal C }$ is a correction term that compensates for the damping of uniform rotation (I. Polichtchouk et al. 2014). Here a ${\nu }_{2{\mathfrak{p}}}$ of 5.9 × 10⁻⁶ and 1.5 × 10⁻⁴³ (in units of ${R}_{p}^{2{\mathfrak{p}}}\,{\tau }^{-1}$) are carefully chosen for ${\mathfrak{p}}$ of 1 and 8, respectively, to ensure that the energy dissipation rate at the truncation wavenumber n_t (=341) is the same for both ${\mathfrak{p}}$-values. At T341 resolution, decreasing ${\mathfrak{p}}$ and/or increasing ${\nu }_{2{\mathfrak{p}}}$ serve to modulate the energy dissipation behavior in small-scale flow structures. No other parameterizations (e.g., radiative transfer and chemistry) and dissipations or drags (e.g., gravity wave and ion) are used; currently, these are poorly known for all hot exoplanets and their inclusion does not obviate the issue addressed here. The focus of the present study is to investigate the dynamics of well-resolved flow structures that arise under a large and constant day–night temperature contrast.

Figure 1 presents the main result of this paper. When forced by a large day–night temperature contrast ramped up on a short timescale, energetic small-scale structures quickly emerge in hot exoplanet atmospheric flows, and the preclusion or mitigation of these structures causes significant differences in the long-term flow and temperature distributions.⁵ Here by "short" we mean a period smaller than ∼1, and by "small" we mean a lateral size smaller than ∼1/10. The significant role of small-scale structures on the flow has been noted and addressed from the inception of hot exoplanet atmosphere studies by J. Y-K. Cho et al. (2003), and explicitly demonstrated to depend on viscosity and resolution in numerical simulations in subsequent studies (e.g., H. T. Thrastarson & J. Y-K. Cho 2011; J. Y-K. Cho et al. 2015; J. W. Skinner & J. Y-K. Cho 2021; J. W. Skinner et al. 2023). In this paper, we highlight the importance of elongated, sharp fronts (that subsequently roll up into long-lived vortices) and internal gravity waves (C. Watkins & J. Y-K. Cho 2010)—both of which unavoidably arise at the beginning of the simulation (as well as throughout the simulation): these structures are generated in response to the atmosphere's attempt to adjust to the applied thermal forcing.

Zoom In Zoom Out Reset image size

The figure shows the ζ(λ, ϕ) fields from two simulations (A and B) at illustrative p-levels (0.05 and 0.95) and times (0.25, 2.50, and 500). The fields are shown in Mollweide projection, centered at the planet's substellar point (λ = 0, ϕ = 0). The two simulations are identical in every way—except energy is removed more rapidly in a slightly wider range of small scales in simulation B than in simulation A for a very brief interval of time, t ∈ [0, 3). The two p-levels are chosen because the features highlighted are generic to the upper region and to the broader, lower region in the simulations (see J. W. Skinner & J. Y-K. Cho 2022).

At t = 0.25, the flows of the two simulations are essentially identical (see A and B in the left column at both p-levels). At this time, the added dissipation in B has not had a chance to act on the flow (as quantified below). However, at t = 2.5, the difference in dissipation is clearly felt by the flow. For example, sharp vorticity fronts (shear layers) in the eastern hemisphere of the dayside and near the equator are markedly different (see A and B in the center column at both p-levels). In general, sharp fronts demarcate the outer boundaries of planetary-scale hetons⁶ and the flanks of an azonal "equatorial jet;" however, unlike in B, the fronts have spawned a large number of small-scale vortices (storms) in A. At t = 500, long after the simulations have been brought back to the common dissipation condition, the flows are still different—and even more so, compared with the flows at t = 2.5. Now hetons are no longer present. Instead, a cyclonic modon⁷ is present in A (center of the frame at p = 0.95) while it is not present in B. Note that the difference in the flow is not due to a temporary "phase offset:" the difference persists over the entire duration of the simulations after t ∼ 3. We emphasize here that this difference cannot be accurately captured below the T341 resolution because the flow structures and their motions are not accurately captured in hot, synchronized-planet simulations starting from rest (J. W. Skinner & J. Y-K. Cho 2021).

Figure 2 shows the (specific) kinetic energy spectrum, E = E(n), of the flows presented in Figure 1; here n is the total wavenumber and E(n) is an average over the zonal wavenumbers at n. In Figure 2, uniform ranges of E and n are shown for ease of comparison. The top row contains the spectra of the flow from the p = 0.05 level, and the bottom row contains the spectra of the flow from the p = 0.95 level. In summary, the figure shows that the effect of a difference in viscosity, which is limited to the small scales and for only a brief period at the beginning of the simulation, spreads to large scales and persists in spectral space—long after the difference has ceased. The spreading is a fundamental property of nonlinearity of Equations (1). It also occurs in the full Navier–Stokes equations (e.g., I. Dobbs-Dixon & D. N. C. Lin 2008; N. J. Mayne et al. 2014; J. M. Mendonça et al. 2016), from which Equations (1) derive.

Zoom In Zoom Out Reset image size

At t = 0.25, the spectra for A and B are identical at both p-levels, as expected from the corresponding physical space fields in Figure 1. Clearly, the dynamics is not affected by the difference in viscosity at this time—at all scales. In contrast, at t = 2.5, a large difference can be seen between the spectra for A and B—particularly at the small scales, n ≳ 10, as expected. The difference is huge in the 20 ≲ n ≲ 300 subrange. At this time all four spectra are still evolving, but the overall shape of each one is nearly stationary after t ∼ 20. Long after the dissipation rate has been rendered identical across the entire spectrum (at t = 3), the spectra at t = 500 are still noticeably different—this time much more at the large scales (n ≲ 10), especially at p = 0.95; at p = 0.05, the difference at the large scales is significant for only select wavenumbers (e.g., n = 2 and n = 3), but the difference is significant for the entire n ≲ 20 subrange at p = 0.95. In fact, at p = 0.95, the difference is significant across essentially the entire range of well-resolved scales above the dissipation range (i.e., n ≲ 200); this is again consistent with the corresponding physical space fields in Figure 1.

Broadly, energy is accumulated in both the large-scale and small-scale subranges (heuristically defined here as n ≲ 10 and n ≳ 100, respectively).⁸ However, the accumulations are different in the two simulations. As the simulations evolve, their spectra become increasingly dissimilar for n ≲ 10, until each simulation reaches a different quasi-equilibrium state. Note that these are the scales which are directly comparable with observations as well as explicitly represented in most current numerical models. However, because of the nonlinear interaction, inclusion of n ≫ 10 in the simulation is necessary to accurately represent n ≲ 10 (e.g., J. P. Boyd 2000; J. W. Skinner & J. Y-K. Cho 2021).

The general behavior seen here explains in part why hot exoplanet simulations in which small scales have been poorly resolved, or altogether missing, would not be able to produce the result of high-resolution simulations at the large scales—as pointed out in many studies in the past (e.g., H. T. Thrastarson & J. Y-K. Cho 2011; I. Polichtchouk & J. Y-K. Cho 2012; J. Y-K. Cho et al. 2015; J. W. Skinner & J. Y-K. Cho 2021). Because underresolved and/or overdissipated models would not be able to capture the intrinsic sensitivity and complexity of the flow, they would give a specious appearance of stability or consistency in their predictions for large scales.

Figure 3 presents a clearer picture of the behavior over time, particularly for spatially averaged quantities. As already alluded to, some averaged quantities are directly important for observations. In the figure, the blackbody total emission flux, ${\mathscr{F}}(t)\equiv \sigma {\langle {T}^{\,4}\rangle }_{\,\,\rm{SS}\,}$, is shown at two p-levels. Here 〈 (·) 〉_SS represents a line-of-sight, projection-weighted (a cosine factor) average over the dayside disk centered on the substellar point ("SS"), and σ is the Stefan–Boltzmann constant; emissivity is assumed to be unity, for simplicity. The flux for each simulation is normalized by the initial flux, ${\tilde{{\mathscr{F}}}}_{0(\cdot )}\equiv {\sigma \langle \,{T}_{0(\cdot )}^{\,4}\rangle }_{\,\,\rm{SS}\,}$, so that ${\tilde{{\mathscr{F}}}}_{(\cdot )}={\tilde{{\mathscr{F}}}}_{(\cdot )}\,/\,{\tilde{{\mathscr{F}}}}_{0(\cdot )}$; for example, ${\tilde{{\mathscr{F}}}}_{\,\rm{A}\,}={\tilde{{\mathscr{F}}}}_{\,\rm{A}\,}(t)$ is the normalized flux for simulation A. The value of the normalization is same for both simulations presented and is also independent of the location of the disk center, due to the spatially uniform temperature distribution used at t = 0. The t ∈ [0, 500] duration is shown in the figure, but the general behavior is unchanged up to t = 1000. We note that an overly long duration (e.g., t ≳ 1000) is unlikely to be physically realistic, given the highly idealized initial and forcing conditions employed; however, it does illustrate the robustness of the particular behavior discussed.

Zoom In Zoom Out Reset image size

As can be seen, large differences in the averaged flux from the simulations persist over a long time and over the entire p-range. It is important to stress here again that—because of the intrinsic, ever-present imbalance and nonlinearity of the atmosphere—the effect of small scales highlighted is not removed by a simple averaging of $\tilde{{\mathscr{F}}}(t)$ vertically, or in time.⁹ The difference ranges approximately ±0.4 at p = 0.05 and ±0.1 at p = 0.95 for t ≳ 8, and continues for the entire duration of the simulations. The variations correspond to disk-averaged temperature differences of up to ±100 K at p = 0.05 and ±75 K at p = 0.95. Such differences, which stem from the acute sensitivities inherent in the PEs,¹⁰ directly impact the ability to correctly predict and/or interpret observations from current and next-generation telescopes such as the James Webb Space Telescope and Ariel (G. Tinetti et al. 2022; J. Rigby et al. 2023). As noted above, this also underscores the critical importance of accurately and consistently representing small-scale flows throughout the entirety of the simulation; the absence of such representation, even for a brief period, permanently vitiates the reliability of model predictions.

Figure 4 shows a more complete picture of the differences at high resolution. The figure presents $\tilde{{\mathscr{F}}}(t)$ at p = 0.95 for six T341L20 simulations (A)–(E)—all identical except for the $(\nu ,\,{\mathfrak{p}})$ pairs, (1.5 × 10⁻⁴³, 8) and (5.9 × 10⁻⁶, 1); see Equation (2). In each simulation, the applied dissipation is switched to another one from {∇¹⁶, ∇²} at the time indicated by the dashed line in the figure. Panels (A) and (D) present reference simulations, with no change in dissipation for all t. Panels (B) and (C) present simulations with parameters identical to those in the simulation of (A) but with a stronger dissipation rate, $(\nu ,\,{\mathfrak{p}})=(5.9\times 1{0}^{-6},\,1)$, applied at the beginning for different durations (t < 3 and t < 20, respectively). Panels (E) and (F) present simulations with parameters identical to those in the simulation of (D) but with a weaker dissipation rate, $(\nu ,\,{\mathfrak{p}})=(1.5\times 1{0}^{-43},\,8)$, applied at the beginning for different durations (t < 3 and t < 20, respectively, here as well).

Zoom In Zoom Out Reset image size

Broadly, two distinct types of $\tilde{{\mathscr{F}}}(t)$ behaviors emerge (see J. W. Skinner & J. Y-K. Cho 2021): (1) a dynamic and chaotic large-scale flow characterized by multiple, persistent large-amplitude oscillations; and, (2) persistent, regular oscillations that "kick in" after a period of small-amplitude oscillations. The two types can be seen in the left and right columns of Figure 4, respectively. The first type is caused by dynamical instability and turbulent motion of energetic, planetary-scale vortices, which ultimately migrate around the planet; these giant vortices, which may be singlets or doublets (J. W. Skinner et al. 2023), interact with a large number of small vortices during the migration in a way reminiscent of Brownian motion (panels (A)–(C)). In this case, the planet's vorticity and temperature fields are highly inhomogeneous, with strong meridional (north–south) asymmetries and vigorous mixing at large scales. The second type is caused by a long-lived, planetary-scale, meridionally (latitudinally) symmetric modon, which is weaker (lower ∣ζ∣) than the giant vortices in the first type of behavior. After a transient period of large excursions from near the substellar point at the beginning of the simulation, there is generally a period of "quiescence," when the modon's position is nearly fixed at the substellar point (e.g., 120 ≲ t ≲ 260 in panel (D)). After this quiescent period, the modon transitions to a one of westward "migration" around the planet—subsequently either transitioning back to the quiescent state (t ≈ 275 in panel (D)) or remaining in the migrating state (t ≈ 230 in panel (D)).

Notice that the quiescent state is not always present in the simulations (panel (F)). However, when it is present (panels (D) and (E)), there is nearly a fourfold reduction in $\tilde{{\mathscr{F}}}$ variations as well as a sustained high amplitude in $\tilde{{\mathscr{F}}}$ compared with $\tilde{{\mathscr{F}}}$ in the migrating state (present in all of panels (D)–(F)). Both the reduction of variation and sustenance of high amplitude occur because there is little or no heat transport away from the dayside. In contrast, when the modon migrates westward, it mixes and advects heat to the western terminator or beyond. The timescale of the mixing/advection is relatively short, evinced by the fast decay time of the regular peaks seen in panels (D)–(F); it is roughly equal to the thermal relaxation timescale for the p-level shown in the figure.

It is clear from Figure 4 that model predictions for the large scales are highly sensitive to the dissipation rate of the small scales. It follows that the sensitivity would be entirely missed if the small scales are not represented in simulations. In the figure, the right column shows that the sensitivity is active even when the small scales are heavily suppressed throughout the majority of the flow's evolution, as long as the small scales are represented; see left column. In addition, the evolutions in (E) and (F) are noticeably different, despite the simulations being identical except for the small difference in the duration of the reduced dissipation at the beginning. Less dramatic, but still noticeable, behaviors are seen in the opposite situation, in which the small scales are more heavily suppressed for the first 3 and 20 days (only) of the evolution ((B) and (C), respectively). In (C), a long-period variation not seen in (A) appears in the evolution; in (B), a suggestive transition to a "quiescent-like" state is observed (see (A)). Finally, it is important to also note that the temporal mean values of $\tilde{{\mathscr{F}}}$ (hence 〈T〉_SS) vary among the simulations at quasi equilibration—even though the thermal forcing is identical in all of them.¹¹

In this paper, we have presented results from high-resolution atmospheric flow simulations with a setup (initial, boundary, and forcing conditions) commonly employed in hot exoplanet studies. Our simulations explicitly demonstrate that the behavior of the atmosphere at large scales is highly sensitive to the rate of energy loss at small scales—the loss both intentional and not. Surprisingly, the sensitivity is present even if the increase in the rate is operating only for a very brief period. Hence, deviation from high-resolution results is fully expected when the small scales are poorly resolved or altogether missing, as has been demonstrated by J. W. Skinner & J. Y-K. Cho (2021). As in that study, the sensitivity is comprehensively illustrated in the physical, spectral, and temporal spaces in this study. Here we clearly demonstrate that high resolution is necessary for generating accurate predictions, especially for hot synchronized planets, since the small scales nontrivially affect the evolution of flow and temperature at the large scales.

More broadly, high resolution (as well as an accurate algorithm) is critical for understanding ageostrophic (unbalanced) turbulence,¹² in general. It is found that the presence—or preclusion—of the small-scale structures, which appear almost immediately in the flow (t ≲ 1), leads the hot exoplanet atmosphere simulations to settle into qualitatively different quasi-equilibrium states. The small-scale flow structures generated at early times of the evolution are (i) sharp, elongated fronts that roll up into energetic vortices and (ii) radiated internal gravity waves. These form in response to the atmosphere's attempt to adjust to the applied forcing—unrelated to the degree to which the radiative process is idealized; the aggressive response is due to the rapidity of the thermal relaxation to a high day–night temperature gradient, leading to large Rossby and Froude numbers (see, e.g., J. Y-K. Cho et al. 2008) for the flow. The need for high resolution and balancing to address such atmospheres has been noted since the beginning of exoplanet atmospheric dynamics studies by J. Y-K. Cho et al. (2003); in that study, nonlinear balancing and slow lead-up to the full thermal forcing have been employed at T341 resolution.

This work has significant implications for general circulation and climate modeling of hot exoplanets. Models that use any combination of low-order dissipation, low spatial or temporal resolution, high explicit viscosity coefficients, or strong basal drags to force numerical stabilization are at risk of generating inaccurate and/or unphysical solutions. This is because all of the above expediencies prevent small-scale flows from being adequately captured throughout the simulation's evolution. In our view, it is unlikely that the state of the hot exoplanet atmosphere can be usefully captured in such models—as the dynamics, which forms the backbone on which physical parameterizations hang, would itself be questionable.

Accurately simulating exoplanet atmospheres is a very complex and difficult problem, requiring meticulous assessment and reduction of uncertainty at every level of the model hierarchy—from the equations solved, to the dynamical core that generate the solutions, to the parameterizations that enhance the dynamics as well as rely on it. This is the case even for the Earth, for which detailed observation-derived (referred to as "analyzed" in numerical weather and climate predictions) inputs to the numerical models are available and the dynamics are geostrophically balanced (small Rossby and Froude numbers). Importantly, effects of small scales very similar to those described here are well-known in Earth climate simulation studies (e.g., J. A. Rial et al. 2004; R. L. Sriver et al. 2015; C. Deser et al. 2020). However, as expected, their effects are weaker compared to those for hot exoplanets, which are far from geostrophic balance. Despite this, the effects seem not to have garnered much attention in exoplanet studies thus far. The preclusion of small-scale structures poses a particularly critical problem in hot exoplanet "radiative transfer and/or chemistry coupled with dynamics" simulations. The high sensitivity of the overall flow to small-scale structures, which arise early in the simulation, means that considerable care must be taken to (i) accurately represent structures of wide-ranging scales throughout the entire duration of the simulation, and (ii) sensibly initialize the simulation.

The authors thank Joonas Nättilä and Quentin Changeat for helpful discussions. J.W.S. is supported by NASA grant 80NSSC23K0345 and a Simons Foundation Pivot Fellowship awarded to Albion Lawrence. This work used high-performance computing at the San Diego Super Computing Centre through allocation PHY230189 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. This work also used high-performance computing awarded by the Google Cloud Research Credits program GCP19980904.