Spontaneous Generated Convective Anticyclones in Low Latitude -- A Model for the Great Red Spot

The Great Red Spot at about latitude $22^{\circ}S$ of Jupiter has been observed for hundreds of years, yet the driving mechanism on the formation of this giant anticyclone still remains unclear. Two scenarios were proposed to explain its formation. One is a shallow model suggesting that it might be a weather feature formed through a merging process of small shallow storms generated by moist convection, while the other is a deep model suggesting that it might be a deeply rooted anticyclone powered by the internal heat of Jupiter. In this work, we present numerical simulations showing that the Great Red Spot could be naturally generated in a deep rotating turbulent flow and survive for a long time, when the convective Rossby number is smaller than a certain critical value. From this critical value, we predict that the Great Red Spot extends at least about 500 kilometers deep into the Jovian atmosphere. Our results demonstrate that the Great Red Spot is likely to be a feature deep-seated in the Jovian atmosphere.


INTRODUCTION
Vortices are ubiquitous in gas giant planets. In Jupiter, more than 500 vortices have been identified from the sequence of 70-day Cassini images with latitude spanning from 80 • S to 80 • N (Li et al. 2004). Recent observation from the Juno spacecraft has detected multiple closely-packed cyclones in the northern and southern polar regions (Adriani et al. 2018;Tabataba-Vakili et al. 2020). Anticyclones of sizes larger than 1000km were also observed in the polar regions of Jupiter, among which most were found to form a dipole configuration with a cyclone (Adriani et al. 2020). The Great Red Spot (GRS), is probably the most prominent large-scale vortex in Jupiter, which is an anticyclone centered at about 22 • S and has been observed for hundreds of years. The origin of GRS, however, remains a puzzle.
To understand the origin and existence of large-scale vortices (LSVs) in Jupiter, two scenarios of flow dynamics have been proposed based on different assumptions on the depth of the Jovian atmosphere. The first is a shallow model in which jets (Liu & Schneider 2010;Lian & Showman 2010) and vortices (Zhang & Showman 2014;O'Neill et al. 2015) are driven by the energy generated from the latent heat of moist convection. The second is a deep model in which columnar jets (Busse 1976;Christensen 2001;Heimpel & Aurnou 2007;Cai & Chan 2012;Aurnou et al. 2015) and vortices (Chan & Mayr 2013;Heimpel et al. 2016;Cai et al. 2021) are maintained by the energy deep from the interior of Jupiter.
While the depth of the Jovian atmosphere is uncertain, recent measurement on gravity field by Juno indicates that the zonal jets could possibly extend thousands of kilometers deep into the atmosphere (Kaspi et al. 2018;Guillot et al. 2018). The heating of Jupiter's upper atmosphere also reveals that the GRS should be heated from below, suggesting that it might have a deep origin (O'Donoghue et al. 2016). Although there are some evidences favoring the deep model, it is fair to say that the possibility of the shallow scenario is not entirely excluded (Kong et al. 2018).
Even if one assumes that the atmosphere is deep, it is natural to ask whether LSVs, such as GRS could be formed. Previous simulations on rotating compressible convection (Chan 2007;Käpylä et al. 2011) revealed that LSVs could be generated in a rapidly rotating turbulent flow. Later, the dynamics of LSVs has been intensively studied in a rapidly rotating Rayleigh-Bénard convection (Guervilly et al. 2014;Rubio et al. 2014;Favier et al. 2014). Despite the considerable progress achieved in these simulations, there are some limitations in these studies. First, the heat fluxes used in these simulations (Chan 2007;Käpylä et al. 2011) are larger than Jupiter's internal heat flux by several magnitudes. Second, in previous simulations on rapidly rotating convection only shear (Novi et al. 2019;Currie & Tobias 2020) or large-scale cyclones (Chan & Mayr 2013) can be found in low latitudes. The present study is motivated by the following questions. Could LSVs be generated in deep convective flow with a flux at or close to Jupiter's internal heat? Could large-scale anticyclones, such as GRS, be generated in a rapidly rotating turbulent flow at low latitudes? Results of some numerical simulations will be discussed in the following sections to address these questions. We first consider the simulations of rotating turbulent convection in a f-plane at high latitudes. The computational domain is a Cartesian box with a constant flux fed at the bottom. Large eddy simulations were performed with a semi-implicit mixed finite-difference spectral code (Cai 2016), which accounts for the effects of density stratification, compressibility, and subgrid scale turbulence (see Appendix A). We use dimensionless variables in the code, with all the variables normalized by the values of height, pressure, density, and temperature at the top of the box. For example, the velocity is normalized by (p top /ρ top ) 1/2 , and the flux is normalized by ρ top (p top /ρ top ) 3/2 , where p top and ρ top are the pressure and density at the top of box, respectively. The pressure and density at Jupiter's surface (defined at the location with 1 bar level pressure) are 10 5 Pa and 0.167kg/m 3 respectively. The internal heat flux of Jupiter is estimated to be 7.48Wm −2 (Li et al. 2018). For Jupiter, the dimensionless total flux F tot J is about 1 × 10 −7 , which is mostly carried by convection in the convection zone. Due to the low Mach number and long period of integration time, simulation at small flux is still a challenging problem in stellar or planetary convection (Kupka & Muthsam 2017). So far, only a few attempts were made on the simulations of solar-like low flux convection (Hotta 2017; Käpylä 2019). Simulations on rotating convection at low fluxes have not yet been reported.
The simulation box has an aspect ratio of Γ = 4 (lateral dimension to height), and a resolution of 512 2 × 101 grid points. By using this resolution, the simulation marginally shows the Kolmogrov inertial range and the small scale effects are approximated by a simple small scale turbulence model (see Appendix D and Fig. 6 for discussion). We simulate a deep flow with a total depth of about 3.2 pressure scale height at latitude θ = 90 • , with the dimensionless rotation rate Ω spanning from 0 to 5.12, and the dimensionless total flux F tot spanning from 10 −5 to 10 −2 . Additional parameters of the simulation cases are given in Table 1 (see Appendix D). The initial thermal structure of each simulation is a polytrope, with about 80 percent of the total flux F tot carried by radiation or conduction, and the rest 20 percent is carried by convection when convection is well developed (convective flux F c ∼ 0.2F tot ). In Jupiter, almost all the flux is transported by convection. Here we raise the proportion of flux carried by radiation or conduction to keep the Prandtl number small (the thermal conductivity κ decreases more rapidly than the turbulent kinematic viscosity ν when flux is reduced, which leads to a raise of Prandtl number; note that scalings κ ∝ F tot and ν ∝ F 1/3 tot are approximately followed in large eddy simulations if grid resolution is fixed) and to speed up the thermal relaxation. This technique has been widely used in simulations of stellar and planetary convection (Brummell et al. 2002;Rogers & Glatzmaier 2005;Cai 2020;Cai et al. 2021).In our numerical setting, the Ekman number has a scaling E ∝ ν ∝ F 1/3 tot and the flux-based Rayleigh number (see Aurnou et al. (2020) for the definition) has a scaling Ra F ∝ F tot /(νκ 2 ) ∝ F −4/3 tot . Therefore, E decreases and Ra F increases with decreasing F tot , respectively. Basically, the simulation reaches a parameter regime of higher Ra F and lower E when F tot decreases. Previous simulations (Chan 2007;Käpylä et al. 2011;Chan & Mayr 2013;Favier et al. 2014;Guervilly et al. 2014;Cai 2021) have revealed that the convective Rossby number Ro c plays an important role on the appearance of LSVs. In our cases, we compute where v ref is the averaged root-mean-square convective velocity calculated in an independent non-rotating reference case, H is the height of the box, and Ω is the rotation rate. v ref has been widely used as the characteristic velocity to evaluate Ro c in the astrophysical community (Kim & Demarque 1996;Chan 2001;Augustson et al. 2016). In the fluid dynamics community, the free-fall convective velocity is usually used as characteristic velocity to evaluate Ro c for Rayleigh-Bénard convection. Simulations on Rayleigh-Bénard convection (Cai 2021) have shown that the free-fall convective velocity is comparable to v ref in magnitude. Therefore, there is no significant difference between using either of these two velocities for Ro c .
After a small initial perturbation and a long period of evolution (at least several thermal relaxation time τ th ∼ 1/F tot ), the convective flow reaches a statistically thermal equilibrium state. Two transitions in flow patterns take place as the rotation rate increases (Fig. 1). The convective Rossby number Ro c , which measures the importance of buoyancy to rotational effects (see Appendix B), is an essential parameter to control the transitions of flow patterns (Chan & Mayr 2013;Guervilly et al. 2014;Cai 2021).
The first transition from turbulence (regime I) to large-scale cyclones (regime II, the flow is still turbulent but large-scale cyclones appear) occurs when Ro c is smaller than a critical number Ro c1 ∼ 0.11 ( Fig. 2A). It is followed by a second transition to the appearance of large-scale anticyclones (regime III, the flow is still turbulent but large-scale anticyclones appear) when Ro c further decreases to another critical number Ro c2 ∼ 0.023 ( Fig. 2A). This result is insensitive to various parameters, such as the magnitude of flux, the Prandtl number, and the Ekman number, at least in the parameter space that has been explored (10 −5 ≤ F tot ≤ 10 −3 , 0.07 ≤ P r ≤ 8.1, 2.5 × 10 −6 ≤ E ≤ 1.8 × 10 −4 ). The large-scale cyclones and anticyclones can survive for a long time. For example, the vortices in Case D7 are still alive after 1600 system rotation periods. In regime II, a pair of cyclones in Case A4 have survived for at least 100 system rotation periods till we terminated the simulation. It indicates that multiple cyclones could be generated in a deep rotating convection zone. When curvature effect is taken into account, these multiple cyclones tend to cluster around the pole (Cai et al. 2021), which may explain the driving mechanics of the circumpolar cyclones observed in Jupiter's poles (Adriani et al. 2018(Adriani et al. , 2020. In regime III, counter rotating vortex dipoles are formed in most cases, similar to the dipole configurations observed in the Jupiter's southern polar region (10 • away from the pole) (Adriani et al. 2020). Mach number also has effect on the survivorship of cyclones. Case A6 shows that a cyclone ( Fig. 1) is twisted and almost destructed when the Mach number reaches one. Similar phenomenon is observed in Case A7. In other cases, cyclones are stable since their speeds are subsonic. Mach number decreases with the convective flux ( Table 1). The vertical velocity is much smaller than the horizontal velocity in rapidly rotating flow (Fig.2C), thus the Mach number reflects the maximum wind speed of the LSVs. As the convective flux of Jupiter is lower, we expect that the LSVs in Jupiter are subsonic. The horizontal velocity reaches a peak value near the second transition state just before the large-scale anticyclones appear (Fig. 2C). Passing through this transition state, the horizontal velocity and vortex size decrease with decreasing Ro c . The sizes of LSVs are quantitatively measured and listed in the last column of Table 1. The observed decreasing vortex size with decreasing Ro c is in agreement with the theoretical analysis of Coriolis-Inertial-Archimedean (CIA) balance on rapidly rotating convection (Christensen 2010), from which the relative size of LSVs at small Ro c obeys a scaling /h ∼ Ro 1/2 c , where and h are horizontal and vertical length scales of a vortex respectively (see Appendix B). Now we can apply the simulation result to Jupiter. The convective flux F c calculated in the simulation group D (F c ∼ 2 × 10 −6 ) is close to the realistic values in gas giant planets (for Jupiter, F c J ∼ 1 × 10 −7 ). Our simulation has shown that the essential factor for determining the appearance of LSVs is the convective Rossby number Ro c . It leads to the argument that this result could be extrapolated to the Jovian atmosphere. To obtain Ro c in Jupiter, we first make an estimation on the convective velocity v of a non-rotating Jupiter. From the mixing length theory, v has a scaling relation with the convective velocity F c of v ∼ (F c /ρ) 1/3 (Böhm-Vitense 1958; Chan & Sofia 1989), where ρ is density. The 1/3 scaling law is tested against the data of the non-rotating cases (Fig.2B). The data shows a scaling relation v = 1.11(F c /ρ top ) 1/3 , which is remarkably consistent with the theoretical mixing-length prediction. With this scaling relation, the dimensionless convective velocity of a non-rotating Jupiter is predicted to be 0.0052, which corresponds to 4.02m/s. As the angular velocity of Jupiter is a constant 1.76 × 10 −4 rad/s, the appearances of large-scale cyclones or anticyclones in Jupiter solely depend on the depth of atmosphere H. For a given Ro c1 or Ro c2 , the lower limit of H can be inferred. Recently, large-scale cyclones and anticyclones have been observed on Jupiter's poles by the Juno spacecraft (Adriani et al. 2018(Adriani et al. , 2020. The condition for the emergence of LSVs can be applied to the poles. To satisfy Ro c ≤ Ro c1 , we expect that the depth of Jovian atmosphere is deeper than in the polar region. To generate anticyclone in the polar region, it is expected that the Jovian atmosphere is deeper than H c2 = v ref /(2ΩRo c2 ) ≈ 501km. As both cyclones and anticyclones exist, we predict that Jupiter's atmosphere is at least 500km deep in the polar regions. Juno's gravitational measurement predicted that the depth of atmosphere at polar regions is in the range of 0 ∼ 500km (Kaspi et al. 2018). The lower limit of our prediction coincides with the upper limit of the prediction from gravitational measurement.

How to generate GRS-like LSVs in low latitudes?
We have shown that LSVs can be generated in a rapidly rotating deep convective flow at high latitudes. Now we demonstrate by simulations that LSVs can also survive in a rotating flow at low latitudes. For simulation parameters at low latitudes, θ is fixed at 22.5 • , F tot is from 10 −5 to 10 −3 , and Ω is from 0.16 to 2.56 (see Table 2 in Appendix D). The basic settings are almost the same as the simulations at high latitudes, except that now the rotational axis is inclined with gravitational axis by an angle 67.5 • .
The appearances of LSVs in rapidly rotating flow at low latitudes are shown in Fig. 3. Also apparent is that large scale anticyclone, similar to GRS observed at Jupiter's low latitude, are formed when Ro c is small. In simulations at high latitudes, columnar circular vortices are aligned with the rotational axis. Unlike previous simulations, here we see that the axially aligned columnar vortices are elliptic on the plane perpendicular to Ω rather than circular (Fig. 4A). This can be shown more clearly by showing the flow structure at a cut plane perpendicular to the rotational axis (Fig. 4E). It demonstrates that elliptical vortex could be generated in deep convective rotating flow. The flow patterns on the meridional-zonal (x-y), zonal-vertical (x-z), and meridional-vertical (y-z) planes can be viewed as the projections of the elliptical vortex. Interestingly, the large-scale anticyclone on the x-y plane shows a circular structure (Fig. 4B), while on the x-z plane it remains elliptic (Fig. 4C). Recently, observation has shown that the GRS is changing from elliptical shape to circular shape (Simon et al. 2018), which shares some similarity with our simulation result. Tilted patterns, which are parallel to Ω as exhibited on the y-z plane (Fig.4D), demonstrating that the columnar structure is indeed aligned with Ω. As indicated in the Taylor-Proudman theory, it is not surprising to observe that the columnar structure is aligned with the rotation axis. In a rapidly rotating system, the vortices are expected to form parallel to the rotation axis regardless of the simulation domain geometry and latitude. The LSVs are stable for a very long time. For example, the vortices in Case D7b have survived for at least 400 rotating periods till we terminate the simulation. Movies S1, S2, S3 and S4 (see Appendix E) show the flow motions at different planes z = 0.1, z = 0.9, y = 3.6, and x = 1.0 for about 27 system rotating periods. In movies S1 (horizontal plane cut at z=0.1) and S2 (at z=0.9), small scale turbulent motions are observed inside the anticyclone, similar to the fluid motions observed in Juno images (Sánchez-Lavega et al. 2018). The large-scale anticyclone can be better understood by showing vortical structures. Surprisingly, small vortical tubes (SVTs) (aligned with Ω) are observed in the instantaneous vortical structure (Fig.4F). Tube-like vortices were found associated with strong down flows in deep convection (Brummell et al. 2002). Movie S5 on vortical structure shows that LSVs are associated with the organized movements of these SVTs. The distributions of positive and negative SVTs are asymmetric. The positive SVTs are almost evenly distributed in the fluid domain (Fig. 5B). However, the distribution of negative SVTs are uneven (Fig. 5A), with more distributed inside the anticyclone and less distributed inside the cyclone. This is further demonstrated in Fig. 4G, which clearly shows the large-scale vortical columns in the time averaged vortical structure.
There are a number of processes that are responsible for causing the asymmetric distributions of positive/negative SVTs in the positive/negative LSVs. First, both laboratory experiments (Vorobieff & Ecke 2002) and numerical simulations (Guervilly et al. 2014; Cai 2021) of rotating convection have shown that ejection of thermal plumes is more frequent than injection at thermal boundary layers. The more frequent ejection of thermal plumes, the more positive SVTs are created. Second, the rotational effects are different in the cyclones and anticyclones. Cyclone spins in the same direction as the system rotation, while anticyclone spins in the opposite way. Thus, the effective rotation rate (the summation of the spinning rate of the LSV and the system rotation rate) is larger inside a cyclone than an anticyclone. As a result, the suppressing effect of rotation on convection is stronger inside a cyclone than an anticyclone (Chan & Mayr 2013). The suppressing effect reduces the numbers of both injected and ejected plumes, but more are reduced inside the cyclone. As the number of injected plumes (negative SVTs) is fewer than that of ejected plumes (positive SVTs), a significant reduction of injected plumes can lead to the sparse distribution of negative SVTs inside the cyclone. Third, two like-signed vortices tend to rotate around each other when they are strong and close enough (Boubnov & Golitsyn 1986;Hopfinger & Van Heijst 1993). The patch of the two like-signed vortices has stronger intensity, and can attract more like-signed vortices to rotate together (Yasuda & Flierl 1995). Finally, LSVs can be formed by the clustering of these like-signed vortices. Inside the large-scale cyclone, the intensities of positive SVTs are stronger and their distributions are denser than those of negative SVTs. The clustering of these strong positive SVTs expels negative SVTs to form an anticyclonic region (Guervilly et al. 2014). When the negative SVTs are strong enough and their distances are short enough, a large-scale anticyclone can be formed in the anticyclonic region.
Compared with the cases at high latitudes, the sizes of LSVs at low latitudes are larger. Since the dominant size of LSVs is in the x-y plane, we expect that the vortex size follows the relation /h ∼ [ v ref /(2ΩH sin θ)] 1/2 = (Ro c / sin θ) 1/2 at the low Rossby regime, and thus, the vortex size at θ = 22.5 • is about 1.6 times larger than that at θ = 90 • . To fit LSVs in the box at a small θ, it requires that either the box size is wide enough or Ro c is small enough. Otherwise, only vague vortices or shear flows exist (the first column of Fig.3). The shear flows might be closely related to the strong prograde and retrograde jets observed in Jupiter's equatorial region. Near the equatorial region, θ is very small. Since the vortex size is proportional to (1/ sin θ) 1/2 , LSVs are unlikely to be formed in the meridional-zonal plane. Alternatively, shear flows will be favored in this region. We expect that, if the box is large enough, the critical Rossby numbers for the appearances of LSVs would remain the same as at high latitudes. Thus we predict that the Jovian atmosphere at low latitudes is deeper than 500km.

DISCUSSION
Our results suggest that long-lived LSVs could be generated in a rotating deep convective flow at both high and low latitudes, when the convective Rossby number Ro c is smaller than certain critical values. In particular, anticyclonic LSVs similar to the GRS can be found in low latitudes for certain critical values of this number. We find that there are positive and negative SVTs inside the LSVs. The distributions of negative SVTs are nonuniform inside the LSVs, with more distributed inside the large-scale anticyclone and less distributed inside the large-scale cyclone. We have proposed three mechanisms to explain it. First, the preference of ejection of thermal plumes in the thermal boundary layers tend to create more positive SVTs than negative ones. Second, the effective rotation rate is larger in the large-scale cyclone than anticyclone, which reduces more negative SVTs in the large-scale cyclone. Third, the clustering of positive SVTs expels negative SVTs in the cyclone, but the clustering of negative SVTs is not strong enough to significantly expel positive SVTs in the anticyclone. We also found that the critical values of Ro c for the appearances of LSVs are insensitive to other parameters, such as the magnitude of normalized flux and the Ekman number. This finding may have certain implication in understanding astrophysical or geophysical flows. In stars or planets, the normalized flux and Ekman number are usually too small to simulate directly. The convective Rossby number, a viscous-free parameter, however, is not as small as the normalized flux and Ekman number, and may be reachable in numerical simulations. It can be anticipated that using similar values of convective Rossby number in the simulations of stars or planets may provide some insights on understanding the dynamics of these astrophysical or geophysical flows. For the Sun, one can estimate that Ro c is in the range of 0.02 to 0.2. The estimation is based on the physical values that the depth of the solar convection zone H ≈ 2.1 × 10 8 m, the rotation rate Ω ≈ 2.6 × 10 −6 s −1 , and the convective velocity v ≈ 20 − 200m/s (Miesch et al. 2008). On the basis of our calculation, we speculate that LSVs might exist at the bottom of solar convection zone, where the convective Rossby number is at around 0.02. For fast rotating solar-like stars, such as HII 296, BD-16351, and TYC 5164-567-1 (Folsom et al. 2016), the convective Rossby numbers are about 10 times smaller than that of Sun. Therefore, LSVs could probably be formed in the surface convection zones of these solar type stars. However, these estimations do not consider the effect of magnetic field. It has been found that a strong magnetic field can inhibit the formation of LSVs in rapidly rotating convection (Guervilly et al. 2015). Application of the result to stars requires further investigation on the effect of magnetic field.
The microwave radiometer in the Juno spacecraft can see up to 350 km deep into the GRS in Jupiter. The Juno gravity measurements can also provide information on the depth of GRS (Parisi et al. 2020;Galanti et al. 2019). Recently, the Juno microwave radiometer data (Bolton et al. 2021) and the gravity data (Parisi et al. 2021) reveal that the GRS has a root which could possibly extend down to 500 km deep. It has to be mentioned that our model is based on the theory of rotating convection on a f -plane. We did not exclude the possibility that LSVs could be formed by other mechanisms, such as the formation of LSVs by opposite zonal jets (Lemasquerier et al. 2020) where LSVs can be formed in a shallower atmosphere. Also, we did not include the curvature β-effect in our model. When β-effect is taken into account, the conditions for the appearance of LSVs at low latitudes could be more restricted. It is challenging to set an appropriate boundary condition in the meridional direction for a box simulation if β-effect and zonal jets are taken into account. Global simulations (Cai & Chan 2012; will be more suitable for such kind of problems. We defer the discussion of global simulations of LSVs to future works. We thank the reviewer for the helpful comments on the manuscript. This work was partially supported by NSFC (Nos. 12173105,11503097) where ρ is density, M is momentum, Σ is stress tensor, p is pressure, g is gravity, Ω is angular velocity, E is total energy, F d is diffusive or conductive flux. We solve the hydrodynamic equations by a mixed finite-difference spectral model (Cai 2016). This method has been used in the study of compressible turbulent convection (Cai 2018) and convective overshooting (Cai 2020). In this paper, additional Coriolis term has been added in the momentum equation. In the model, all the linear terms, including the Coriolis term, are integrated semi-implicitly. An advantage of this model is that the numerical time step is not restricted by the CFL conditions associated with speeds of acoustic wave, gravity wave, or inertial wave. Thus the model is suitable for simulating rapidly rotating compressible turbulent convection at small flux. In this paper, the initial thermal structure is assumed to be in a polytropic state: where T is the temperature, η is a ratio measuring the depth of fluid, m is the polytropic index, and the subscript top denotes the corresponding value at the top of the computational domain. All the physical variables are normalized by the height, density, pressure, and temperature at the top. The adiabatic polytropic index is set to be m ad = 1.5. The fluid is convectively unstable when m < m ad , and stable when m > m ad . In our settings, we consider the rapidly rotating convection in a pure convection zone with η = 4 and m = 1. The numerical experiments are performed with a large eddy simulation method, in which the small scale effects are modelled by a sub-grid-scale (SGS) turbulent model. In the SGS model (Smagorinsky 1963), the turbulent kinematic viscosity is evaluated by where c d = 0.28 is the Deardorff coefficient, σ is the strain rate tensor, and ∆x, ∆y, and ∆z are the mesh grid sizes in the x-, y-, and z-directions, respectively. The turbulent viscosity is decomposed into two components: a time-dependent horizontal mean component and a fluctuating component. In the numerical calculations, the higher order viscous terms associated with the fluctuating component are disregarded. Impenetrable and stress-free boundary conditions are used for velocities at the top and bottom of the box. The temperature is set to be a constant at the top, and a constant flux F tot = κη is supplied at the bottom. Here the thermal conductivity κ has a constant value throughout the domain. In the horizontal directions, periodic boundary conditions are used for all the thermodynamic variables. The computational domain is a f-plane with vertical z-direction inclined to the rotational axis by an angle 90 • − θ, where θ is the latitude. The horizontal x-and y-directions of the f-plane are corresponding to the west-to-east zonal and south-to-north meridional directions in the spherical coordinates, respectively. Fourier spectral expansions are used in the horizontal directions, and finite difference method is used in the vertical direction. We use density ρ, pressure p, total energy E, divergence of horizontal momentum δ = ∇ where µ is the dynamic viscosity, Ω z is the vertical component of Ω, and the overline denotes horizontal average.

B. SCALINGS OF CIA BALANCE
For a compressible flow, the vorticity equation can be obtained by taking curl on both sides of the moment equation: where u is the velocity, ω = ∇ × u is the vorticity, Σ is the stress tensor, and p and ρ are perturbations from their equilibrium states. In rapidly rotating convection, CIA balance among the Coriolis (C) term, the inertial (I) term, and the Archimedean (A) buoyancy term is achieved in the vorticity equation (Christensen 2010): A scale analysis of the terms in the CIA balance gives (Aurnou et al. 2020) where U is the characteristic velocity of the LSVs, U c is the characteristic velocity of non-rotating convection, and h and are the vertical and horizontal length scales of LSVs, respectively. It quickly obtains As a result, the relative size of LSVs at small Ro c obeys a scaling /h ∼ Ro 1/2 c . It is consistent with the prediction of Aurnou et al. (2020), where they used the terminology system-scale Rossby number.

C. EXTRACTION OF VORTICES
We extract vortices by the λ 2 -criterion (Jeong & Hussain 1995) that the quantity λ 2 is defined as the second largest eigenvalue of the 3 × 3 martrix A = S 2 + ω 2 , where S = (J + J T )/2 and ω = (J − J T )/2, and the velocity gradient tensor J = ∇u. The vortex shape is then identified by choosing a specific negative value of λ 2 and plotting the corresponding isosurface.

D. SIMULATION PARAMETERS
To explore how flux affects the flow pattern in rapidly rotating flow, we performed simulations on four groups of fluxes with dimensionless values F tot ∈ {10 −2 , 10 −3 , 10 −4 , 10 −5 } at latitude θ = 90 • . We then choose eight different angular velocities Ω in each group to investigate the effect of rotation, giving a total of 32 simulation cases. The nonrotating cases are computed for reference. The detailed parameters of each simulation case are given in Table 1. F tot is the total flux (about 80% of F tot is transported by adiabatic temperature gradient); Ω is the angular velocity; v is the root-mean-square (rms) velocity; v z is the rms vertical velocity; v h is the rms horizontal velocity; P r = ρc p ν/κ is the averaged Prandtl number, where c p is the specific heat capacity at constant pressure, ν is the kinematic viscosity, and κ is the thermal conductivity; Re = v H/ν is the averaged Reynolds number; Re z = v z H/ν is the averaged vertical Reynolds number; M a = max(v/c s ) is the time averaged maximum Mach number at the top layer; Ro c = v /(2ΩH) is the averaged convective Rossby number; E = ν/2ΩH 2 is the averaged Ekman number. The symbol denotes that the average is taken both temporarily and spatially throughout the computational domain. The symbol overline denotes that the average is taken temporarily. The symbols 'n', 'c', and 'ac' denote no LSV, large-scale cyclone, and large-scale anticylone, respectively. The last column gives the corresponding radii of LSVs. The radius of a LSV is measured by the distance from the center of the LSV to the maximum value of averaged tangential velocity around the center at the plane z = 0.5.

Case
Ftot    To investigate whether LSVs could be formed at low latitudes, we performed a total of 9 simulations with different fluxes and angular velocities in the tilted f-plane at the latitude θ = 22.5 • . The basic settings of B5b-B7b, C5b-C7b, and D5b-D7b are all the same to the simulation Cases B5-B7, C5-C7, and D5-D7 at the latitude θ = 90 • . The only difference is that the angular velocity has two components in the tilted f-plane: the vertical component Ω sin θ along the z-direction, and the horizontal component Ω cos θ along the y-direction. The symbol 's' in the last column denotes shear flow. The detailed parameters of each simulation case are given in Table 2. Other parameters and symbols have the same meanings as mentioned in Table 1. At a given z, the velocity v can be decomposed into the summation of a Fourier series of v mn . Then the two-dimensional kinetic energy spectrum (Chan & Sofia 1996;Cai 2018) can be evaluated as where the symbol star denotes the conjugate of a variable. Fig. 6 shows the time averaged compensated power spectral density of the kinetic energy kP 2 (k) as a function of k for case D7b. The simulation marginally shows the Kolmogrov inertial range.  Figure 6. Time averaged compensated power spectral density of the kinetic energy kP2(k) as a function of wavenumber k at z = 0.5 for case D7b. The averaged period is about 14 system rotation periods. The Kolmogrov k −5/3 scaling is shown with the dashed line for reference.
positive (red color) and negative (blue color) vertical vorticity of Case D7b. The time period for each movie is about 27 system rotation periods.