Simulating Solar Near-surface Rossby Waves by Inverse Cascade from Supergranule Energy

Over the last several decades nonlinear hydrodynamic shallow-water models have been applied for studying global oceanic and atmospheric dynamics (see, e.g., p. 387 of Pedlosky 1987). Such models, as well as the generalized versions, including magnetohydrodynamics, have been extensively applied for studying global solar dynamics (Gilman 2000; Zaqarashvili et al. 2010; Klimachkov & Petrosian 2017; Dikpati et al. 2018b). A full set of nonlinear hydrodynamic and magnetohydrodynamic shallow-water equations, applied to study tachocline instabilities, have been presented in an inertial frame of reference in spherical coordinates (see, e.g., Dikpati & Gilman 2001; Gilman & Dikpati 2002). Here we focus primarily on the hydrodynamic shallow-water system, and present the equations in the rotating reference frame, corotating with the core-rotation rate.

Basic formulation of a global nonlinear shallow-water model, including physical assumptions, approximations, and boundary conditions, applied to the Sun have already been presented in various papers (see, e.g., Dikpati et al. 2018b), which also provide detailed formulations of energy integrals in closed form. We do not repeat them here. In brief, shallow-water models are quasi-3D models with a much larger horizontal extent than thickness. The variations in the horizontal direction are much larger than in the vertical direction. Due to the thinness of the fluid layer compared to radius of the astrophysical body divergence of the radii and the density variation are ignored in the momentum and mass continuity equations. Traditionally, shallow-water models are hydrostatic in the vertical because the horizontal spatial scale of the motions is much larger than the thickness, but it has recently been shown (Jalali & Dritschel 2021) that such models can be extended to include nonhydrostatic motions whose horizontal scale is comparable to the thickness of the layer. We have not employed that extension here, but it may be useful in the future.

Before we write the full set of nonlinear hydrodynamic shallow-water equations in the frame rotating with core-rotation rate (ω_c ), we discuss the application of the shallow-water model in our case here. The supergranulation layer of the Sun is thin, with thickness no more than a few percent of the solar radius. Therefore, for the simulations that follow, we use a global thin spherical shell. Its undisturbed thickness varies with latitude according to whether the shell is rotating, differentially rotating, or not rotating. In our case, the variations in thickness are all small, but the dynamical perturbations to the thickness are very important because they allow for gravity waves that can absorb energy put in at the smallest scales of the model, without giving rise to energy transfer to excite global-scale Rossby waves. Thus, in this model the nonlinear mode-mode interactions can choose to inverse cascade and/or to get absorbed in gravity waves.

If we denote λ as longitude, ϕ as latitude, t as time, u as the longitudinal velocity in the rotating frame, v the latitudinal velocity, and h the height deformation of the top surface, then the dimensionless governing shallow-water equations are given by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial u}{\partial t} & = & \displaystyle \frac{v}{\cos \phi }\left[\displaystyle \frac{\partial v}{\partial \lambda }-\displaystyle \frac{\partial }{\partial \phi }(u\cos \phi )\right]\\ & & -\,\displaystyle \frac{1}{\cos \phi }\displaystyle \frac{\partial }{\partial \lambda }\left(\displaystyle \frac{{u}^{2}+{v}^{2}}{2}\right)-G\displaystyle \frac{1}{\cos \phi }\displaystyle \frac{\partial h}{\partial \lambda }+2{\omega }_{c}\,v\,\sin \phi ,\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial v}{\partial t} & = & -\displaystyle \frac{u}{\cos \phi }\left[\displaystyle \frac{\partial v}{\partial \lambda }-\displaystyle \frac{\partial }{\partial \phi }(u\cos \phi )\right]\\ & & -\,\displaystyle \frac{\partial }{\partial \phi }\left(\displaystyle \frac{{u}^{2}+{v}^{2}}{2}\right)-G\displaystyle \frac{\partial h}{\partial \phi }-2{\omega }_{c}\,u\,\sin \phi ,\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}(1+h) & = & -\displaystyle \frac{1}{\cos \phi }\displaystyle \frac{\partial }{\partial \lambda }\left((1+h)u\right)\\ & & -\,\displaystyle \frac{1}{\cos \phi }\displaystyle \frac{\partial }{\partial \phi }\left((1+h)v\cos \phi \right),\end{array}\end{eqnarray} \tag{ 3 }$$

In Equations (1) and (2) above, the parameter G is no longer a measure of the actual effective gravity of the supergranule layer, but rather a parameter to set the characteristic frequencies of gravity waves present in the model. The dimensional parameters within G, such as the radius, layer thickness, and gravity, have somewhat different values than were originally chosen for the application of the model to the tachocline (radius larger, gravity smaller, thickness in the same range), but the appropriate G values to include gravity waves that are relatively high frequency compared to Rossby waves are in the same range as they were for the tachocline case. Gravity waves in this model have frequencies proportional to G^1/2, so the smaller G is, the lower the frequency. We take a G value large enough that the model gravity waves propagate much faster than Rossby waves, so they are quite distinct. We take G = 10 for all cases studied. To evaluate G, we use an estimate of the thickness of the supergranule layer. In order to separate the timescales of Rossby and gravity waves, G just has to be large enough. In this study G has been specified to be of value 10. By following the definition of equivalent height as given in Pedlosky (1987), the thickness of the supergranule layer would be somewhat larger than the equivalent height, since the latter is close to ∼1 scale height. So if we consider an equivalent height of that layer, the value of G would be greater than 10, for which the dynamics of the model will be similar because the behavior of modes was found to change the regime from high G to low G below about G = 1 (Dikpati 2012).

Nonlinear shallow-water equations have been solved by Dikpati and colleagues using pseudo-spectral schemes. Since the details can be found in Dikpati (2012) (see also Dikpati et al. 2018b), here we briefly describe the major steps, and describe certain small differences in the present simulations compared to that of TNOs (tachocline nonlinear oscillations).

To solve Equations (1)–(3) we perform the following steps: (i) compute u, v, and h in a Gaussian grid of ϕ − λ; certain nonlinear terms are computed in real space and then are converted to spectral space; (ii) for time evolution, first kick-start using the fourth-order Runge–Kutta method, then apply the Adams–Bashforth method; (iii) momentum is checked and balanced every few steps; (iv) in order to avoid ringing/Gibb's phenomena, a small spectral viscosity is added; (v) a semi-implicit scheme is used to integrate the effect of very high-frequency gravity waves. For the length of integrations we have done, the small spectral viscosity added makes no measurable difference in the constancy of the total energy of the system with time. Velocities do not decay in total, which we will see later in the results section while discussing our simulation results.

In the TNO study, Dikpati et al. (2018b) focused on low longitudinal wavenumbers m up to m = 10 or so, and the evolution of the model-system containing spot-producing magnetic fields was advanced for several months. However the present study requires all possible l and m in spherical harmonics representations of u, v, and h. With m ≤ l, we consider all possible m for l up to 21 and 42, respectively, in the T21 and T42 models. The computation is done in a 4.5 petaflop supercomputer "Cheyenne" at the NCAR-Wyoming supercomputer center, using OpenMP parallelization.

To perform a simulation run for 1 day, 500 cpu hours are required in a T21 model. But in a T42 model, this time increases by about a factor of 10. Therefore, we demonstrate only one case of our simulations in a T42 model, and the rest we do using a T21 model. Performing simulations of one case in both resolutions also confirms the validity of the numerical scheme. The T42 model truly captures the supergranule scale, but due to being too expensive, we do almost all our simulations in a T21 model. T21 is not fine enough to truly include supergranule scale horizontal motions; however, if similar results are obtained using both the T42 and T21 models, we can safely use the T21 model for exploring various numerical experiments.

For the initial conditions in all cases below we put kinetic energy in almost all l, m modes to ensure nonlinear interaction, but with the spectrum in m heavily weighted toward the highest wavenumbers m allowed in the model; this spectrum peaks at m = 35 for the T42 model and m = 15 for the T21 version, above which it falls off sharply to the highest m allowed, namely, m = 42 and m = 21, respectively. We also include initial potential energy in the form of undulations of the top surface to further enhance the ability of the model to stimulate gravity waves on all spatial scales.

We compute the radial component of vorticity and horizontal divergence present in the initial motions. We know that Rossby wave velocities are almost all vortex motions with very little divergence. By contrast, supergranule horizontal motions are mostly divergent as estimated from results in Langfellner et al. (2015) and Böning et al. (2015). We select the initial velocities in such a way as to reflect these characteristics of supergranulation. The vorticity ζ and horizontal divergence δ in our initialization are given, respectively, by

$$\begin{eqnarray}&&\zeta =\displaystyle \frac{\partial v}{\partial \lambda }-\displaystyle \frac{\partial }{\partial \phi }(u\cos \phi )\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&\delta =\displaystyle \frac{1}{\cos \phi }\left(\displaystyle \frac{\partial u}{\partial \lambda }+\displaystyle \frac{\partial }{\partial \phi }(v\cos \phi )\right).\end{eqnarray} \tag{ 5 }$$

Since energy is a scalar quantity it is straightforward to compute the kinetic energy in the initial state as a function of m. But the vorticity and divergence, being or containing vectors, can be computed at each grid point of the domain. To compute the mean-square vorticity and divergence as function m at t = 0, we integrate the square of these quantities over the entire domain for each m. A similar methodology is used for the calculation of mean-square vorticity, which is often called enstrophy in fluid mechanics (see, e.g., Weiss 1991; Umurhan & Regev 2004).

How the system responds to the initial kinetic energy at small scales critically depends on the magnitude of the energy. Small initial energy means the nonlinear advection terms in Equations (1) and (2) will always be small compared to horizontal pressure gradients and Coriolis forces, and even compared to local accelerations (∂u/∂t, ∂v/∂t). There will be no preference for Rossby waves over gravity waves modestly influenced by rotation. That is clearly not the situation in the Sun's supergranule layer, where the kinetic energy of horizontal supergranulation flow is three or four orders of magnitude larger than the Rossby waves observed. If the initial kinetic energy in the system is large enough to be responsible for producing significant nonlinear advection terms, mode-mode nonlinear interactions in small-scale motions can lead to a different quasi-stable state in the system. While the nonlinear advection terms in the equations of motion couple modes with different l and different m together to excite other modes with other l, m values, in rotating spherical geometry there is a second mechanism for mode coupling, namely, the Coriolis forces, which couple modes of the same m but differing l values. In the general case, both mode coupling mechanisms operate.

As previewed above, we first carry out a simulation with the high resolution T42 model, which truly captures supergranule scale. Figure 1(a) shows the initial kinetic energy spectrum for this simulation, which peaks around m = 35, similar to where supergranulation energy peaks. Initially, there is essentially no energy at low m. Figure 1(b) shows the initial distribution of squared vorticity and horizontal divergence by wavenumber m, which shows that at the highest m, which represent supergranule scale horizontal motions, the motions have little vorticity compared to their divergence, consistent with observations. In fact, all but the very lowest m modes have more initial divergence than vorticity. And both are small for m < 15 or so compared to higher m. This ensures that initial conditions are not favored for producing large vortical motions.

Zoom In Zoom Out Reset image size

Can the Coriolis forces acting on the small-scale divergent motion generate large vortical motions as the dynamics evolve? Can mode-mode nonlinear interactions carry kinetic energy in wavenumber space from large to small m, where Coriolis forces have much greater influence on the flow?

Three snapshots of the evolving kinetic energy spectra in l, m space are shown in Figure 2. The movie shows their continuous evolution (attached in the caption of Figure 2). The total calculation shown is for 1.7 days. We see that the energy very quickly migrates from high m and l to low m, l, to modes of less than 6,6 or so. Each of these modes contains far more kinetic energy than do high individual m's at the start of the simulation. Thus, we see that an inverse cascade of energy via mode-mode interactions can quickly take energy from supergranule scales and deposit it in global low m modes that are likely to be Rossby waves.

Zoom In Zoom Out Reset image size

These low m modes now contain mostly vorticity, which is characteristic of Rossby waves. Having obtained the basic inverse-cascade phenomena in a T42 model, several numerical experiments as well as analysis need to be performed to understand the features of these low m modes. However, the T42 model is extremely expensive: producing a 1.7 day run uses 36 cores × 12 hr per run × 45 runs, which is approximately 20,000 core hours. Therefore, for the feasibility of completing the necessary simulations, we do all subsequent simulations reported below with the T21 version of the model, which is roughly 10 times faster. The T21 model can retain all the characteristics of the T42 model, just with less scale separation between the smallest resolved scales and global scales.

Figure 3 shows the initial kinetic energy spectrum for a simulation at T21 summed over l for each m. For this case the spectrum peaks around m = 15, in contrast to the the T42 case shown in Figure 1, where it peaks near m = 35. We also see in Figure 3 that the vorticity and divergence mean-square amplitudes are similar to those for the T42 case, but of course peaking at a lower wavenumber. For the T21 simulations, the horizontal divergence is initially much larger than the vorticity for all longitudinal wavenumbers.

Zoom In Zoom Out Reset image size

Figure 4 shows a sequence of energy spectra during an evolving solution. Panel (a) shows the spectra shortly after the initialization (about 0.18 day). Here kinetic energy is spread over a wide range of wavenumbers, but concentrated at the high wavenumber end. As time progresses (panels (b), (c)) one can clearly see energy migrating to both low l and low m, and by panel (d) the highest energy concentrations are in l, m at 6,6 and below. By about 2.4 days and onward (see the bottom two panels, (e) and (f), the spectrum is stabilizing with dominant energies present in m = 1−4 with corresponding $l^{\prime} s$. These are most likely global-scale Rossby waves. In the accompanying video, it is possible to see waves of kinetic energy migrating from high to low m with time. This appearance is partly due to potential energy at high m in the initial state being converted to kinetic energy of the same m as the calculation proceeds. This initial potential energy is providing a source of energy more extended over time, but it is also limited by the constraint that the total energy, kinetic plus potential, is pretty accurately conserved in these simulations, very much like what we saw in the TNO simulations (see, e.g., Figure 3(a) of Dikpati et al. 2017).

Zoom In Zoom Out Reset image size

It is to be noted that the energy spectra cannot reach a fully stable state because of the presence of large nonlinearity. The spectra reach a quasi-stable state, in which the dominant m = 1−4 modes continue nonlinear exchange of energies among themselves in a slow rate (compare panels (e) and (f) of Figure 4), but no evidence of the energy cascading back to higher wavenumbers is seen even in the longest simulation of 25 day duration. A small difference between the solid rotation simulations performed using a T21 model and T42 model is that in the T21 model the Rossby waves are even more concentrated at the lowest m's, but with very similar overall evolution and final spectral profiles.

Figure 5 shows how kinetic (panel (a)) and potential (panel (b)) energies migrate in m-space as a function of time. In Figure 5(a) we present the kinetic energy for each longitudinal wavenumber m for all possible spherical harmonic indices l and display the evolution of the kinetic energy in a contour plot, in which m is the vertical axis and t the horizontal axis. The constant energy contours in panel (a) show how the kinetic energy is being transferred from high m to low m. The peak in kinetic energy, denoted by yellow in this PARULA color map (dark blue–sky blue–green–yellow), migrates smoothly from wavenumber m = 15 down to a band of energy in wavenumbers 1–4, with a peak at m = 2,3. Note also that only a small or negligible amount of kinetic energy ends up in the m = 0 band, which implies latitudinal differential rotation. This result suggests that in the solar supergranule layer, the inverse cascade of energy from supergranule scales to larger scale does not play a significant role in the generation and sustenance of observed latitudinal rotation gradient. But it has been known for many years that vertical momentum transfers by supergranulation may be responsible in significant measure for the near-surface shear layer observed there (Gilman & Foukal 1979; Hotta et al. 2015; Matilsky et al. 2019). Finally, we note that there appear to be weak bands of kinetic energy at higher wavenumbers that show migration toward both lower and higher m with time. These represent some combination of gravity waves and nonlinear mode-mode interactions that at times are reflecting back toward higher m with time.

Zoom In Zoom Out Reset image size

We can see in the bottom panel (b) of Figure 5 that the potential energy also moves from high to low m in spectral space in about the same time interval, but with a small shift in time by about a few to several hours with respect to the movement of the kinetic energy spectra. Note that the shell thickness variations that determine it are closely spatially correlated with the horizontal velocities. This is because Rossby waves are nearly geostrophic, with vortical velocities swirling round high and low pressure centers, consistent with a near balance between horizontal pressure gradient and Coriolis forces.

Figure 6 shows that the total energy, kinetic plus potential, is conserved over the length of the T21 simulation, despite the presence of a small amount of spectral viscosity at the highest wavenumbers. The simulations presented here are essentially for an ideal HD setup. We see several interesting features revealed from Figure 6. First, some out of phase oscillations with time between kinetic and potential energies are evident as the simulation proceeds. This is expected because the kinetic and potential energy reservoirs can nonlinearly exchange energies between them, here in an approximate timescale of gravity waves. We also see that the kinetic energy slowly grows with time at the expense of the potential energy of the system. This is also a clear evidence that the Rossby waves are getting generated in the system. Much longer simulations reported in Dikpati (2012) and Dikpati et al. (2017) show the same precision of total energy conservation of the system.

Zoom In Zoom Out Reset image size

We bring out further properties of the emerging Rossby waves by constructing synoptic maps of the flow and height fields, and by measuring the phase speed of the Rossby wave patterns on these maps. Figure 7 shows a sequence of three synoptic maps of the horizontal flow vectors and the variations in the thickness of the single layer of the shallow-water model. We can see clearly that as time advances, the initial pattern, dominated by small scales, evolves over time to one with global organization of flow and thickness into N-S bands at a much lower wavenumber. These bands are the manifestation of low longitude wavenumber Rossby waves. Over time, there are slow variations in the total pattern, as well as migration of the pattern in the retrograde direction (right to left in Figure 7), which is characteristic of classical Rossby waves. A more detailed analysis reveals that there are nonlinear interactions among these Rossby waves, contributing to time variations in the amplitude of each wave. From such data it may be possible to define lifetimes of individual Rossby waves, in a way analogous to what is done in helioseismological studies of solar acoustic waves.

Zoom In Zoom Out Reset image size

Our initial flow pattern, concentrated at high m, contained much more horizontal divergence than vertical vorticity. What has changed as a result of the inverse cascade? Figure 8 shows the mean-square vorticity and divergence after about 4 days of simulation. We see that now there is much more vorticity than divergence, especially at low m. This has occurred as a result of Coriolis forces acting on the flow at all m, coupled with the inverse cascade of kinetic energy, which inevitably results in vorticity being concentrated where in m space the velocity amplitudes are the largest. High vorticity compared to divergence is a fundamental characteristic of Rossby waves. Since Figure 8 is a snapshot of vorticity and divergence rather than an average, we should expect the envelope of vorticity amplitudes will not be smooth at any given time. This accounts for the somewhat smaller amplitude in the m = 5 mode compared to its immediate neighbors. At another time, m = 5 might be larger than its neighbors. The nonlinear interactions among Rossby waves with neighboring wavenumbers will continually fluctuate the amplitudes relative to each other.

Zoom In Zoom Out Reset image size

We can also see in Figure 8 that there is a small amount of vorticity and divergence in the m = 0 mode. The divergence is associated with the time dependent meridional circulation, obviously much weaker than latitudinal flow for m > 0 modes. Coriolis forces acting on this circulation generates a small amount of vorticity in the small induced differential rotation. This small effect is concentrated in middle and high latitudes where the Coriolis forces are stronger. As can be seen in Figure 8, this feature is quite weak compared to the Rossby waves induced by the inverse cascade.

In this section, we simulate the mode-mode interactions when the T21 model is initialized with a differential rotation of the form as derived from helioseismic observations (see, e.g., Charbonneau et al. 1999):

$$\begin{eqnarray}&&\omega ={s}_{0}-{s}_{2}{\sin }^{2}\phi -{s}_{4}{\sin }^{4}\phi .\end{eqnarray} \tag{ 6 }$$

So for these cases we have a prescribed initial differential rotation velocity for the longitudinal motion u for the m = 0 mode, but then we allow the system to evolve as it will according to the nonlinear dynamics, still conserving total energy (kinetic plus potential energies) of the system.

For G = 10, we know that the differential rotation given by Equation (7) is hydrodynamically stable to perturbations unless s₄ is large compared to s₂. When we simulated the nonlinear interactions among large m modes, including a pole-to-equator differential rotation of 28% equally contributed by s₂ and s₄ terms, we get very similar results to those with no differential rotation because this differential rotation is stable. A very similar inverse cascade of kinetic energy occurs in this case, but the Rossby wave longitudinal phase speeds are altered somewhat by the presence of the differential rotation.

The case with unstable differential rotation is more interesting. For this case we distributed the same 28% pole-to-equator differential rotation contributed by s₂ and s₄ terms as s₂ = 0.04 and s₄ = 0.24. This profile is unstable with perturbation e-folding growth time of about 5 months. Figure 9 shows the kinetic energy spectrum plotted in l, m space, analogous to that shown in Figure 4 for the case of solid rotation. We again see evidence of an inverse energy cascade, but one that appears less vigorous than without differential rotation present, in that kinetic energy is somewhat less concentrated at the very lowest wavenumbers, and with somewhat more energy remaining in the high wavenumbers. What is happening here is that the kinetic energy is coming into low l, m modes both from higher l, m and from the differential rotation itself, resulting in some oppositely directed energy cascades, leading to peak energy for modes l, m of 3 and 4, rather than 1 and 2 as in Figure 4. In some sense, the low m modes 1 and 2, which are unstable, push back against the inverse energy cascade.

Zoom In Zoom Out Reset image size

We can estimate the phase speed in longitude for each longitudinal wavenumber m from synoptic maps made at a sequence of times, and tracking the distance each peak in vorticity or thickness moves. We do this by creating images of velocities filtered for each m mode, in which all l's are present. Because the system is nonlinear, the movement is not always at the same speed, and there are amplitude changes as well, arising from various mode-mode interactions due to the nonlinear fluid advection terms in the equations of motion. Gravity waves also cause occasional rapid movements in latitude, during which the tracked vortex pattern may fade out and reappear. Therefore, a phase-speed estimate is a formidably difficult task. However, we manage this by tracking simultaneously a few vortices, so that all do not fade out at the same time. Thus, the estimates are reasonably stable, and consistent with each other. We compare all phase speeds measured in this way with the classical Rossby wave formula for phase speed, originally derived by Haurwitz (1940), given in dimensional form as

$$\begin{eqnarray}&&\displaystyle \frac{\omega }{m}=-\displaystyle \frac{2{{\rm{\Omega }}}_{c}}{l(l+1)},\end{eqnarray} \tag{ 7 }$$

in which Ω_c is the rotation rate of our coordinate system, and as above, l is the spherical harmonic index.

Figure 10 shows measures of the westward (retrograde) propagation speeds of Rossby waves with longitudinal wavenumbers between 1 and 10, all measured dimensionally in units of degrees of longitude per day, a familiar measure for solar rotation observations. The classical Rossby wave speed is denoted with black filled circles, together with two cases starting from solid rotation, one with large initial kinetic energy perturbations at high m (deep blue circles), the other with medium kinetic energy perturbations (sky-blue circles), and finally the phase speeds for a case with large initial kinetic energy perturbations at high m along with an initial latitudinal differential rotation, the profile of which is unstable to perturbations for m = 1 and 2 (orange circles). In order to avoid clumsiness in the plot, the error bar has only been put in one curve, namely, for the high-perturbation case with solid rotation. It is to be noted that the error is slightly more for low m modes due to their faster speeds.

Zoom In Zoom Out Reset image size

We can immediately see in Figure 10 that all cases show very similar profiles of phase velocity as functions of longitudinal wavenumber m, even the case for which the differential rotation profile used is itself unstable for m = 1 and 2. Our interpretation of the differences in phase speed is that, unlike the original Rossby–Haurwitz waves, our simulated Rossby waves are large enough in amplitude to themselves be nonlinear; just as finite amplitude acoustic waves commonly travel faster than the sound of speed, so too these Rossby waves travel faster than their linear counterparts. The results that the medium initial perturbation in the kinetic energy of high m modes are faster than the Haurwitz waves by about half the amount for the high initial kinetic energy modes supports this interpretation.

By contrast to solid rotation, for the case with unstable differential rotation present, the Rossby waves are a little slower retrograde than the Haurwitz formula (see the orange-filled circles in Figure 10). This is because the positive differential rotation in low latitudes reduces the retrograde propagation speed by a small amount. But these waves are still propagating retrograde. Herein lies a possible paradox. It is well known that in rotating shear flow in channels, differential rotation is unstable to nonaxisymmetric (m > 0) perturbations only if the phase speed of the growing wave falls within or close to the range of differential rotation; see, e.g., the so-called semicircle theorem of Drazin & Howard (1962), (Pedlosky 1987, see Chapter 7, p. 514). For a thin spherical shell like the Sun's supergranule layer, this means unstable Rossby waves must have phase speeds less than the maximum rotation rate at the equator, but not much less than the rotation rate at the poles (the limit being set by a factor proportional to the latitude derivative of the Coriolis parameter, the β-effect in meteorological parlance). Certainly the very retrograde speed for m, l = 1 modes in the Haurwitz formula are ruled out as unstable modes. So what Rossby waves are appearing in Figure 10 for this case? They are in fact a set of neutral Rossby waves, much like those for the solid rotation cases. Examination of a time sequence of synoptic maps for this case indeed shows that both neutral and unstable waves are present, with very large variations with time in phase speed and mode structure with latitude, as unstable and neutral waves compete with and interact to produce the total map. The unstable waves have their origin in the initial differential rotation, while the neutral, retrograde propagating waves are driven by the inverse energy cascade, just as in the solid rotation cases. The competition between these waves of different origin is most intense for the m = 3 and 4 modes, each supplying energy respectively from lower and higher m, leading to the peak energies seen there.

It is worth emphasizing that the difference in longitudinal phase speeds between energetically neutral Rossby waves uninfluenced by differential rotation, and Rossby waves excited by instability of differential rotation, is very large: the m = 1 neutral wave has a retrograde phase velocity equal to the rotation rate of the system itself, which means to an observer looking at the Sun from faraway the wave would appear to not be rotating at all. By contrast, the unstable wave would be rotating with the Sun with a rate less than the maximum rate at the equator, likely a rate between that and the minimum rate at the poles, or perhaps slightly less than the polar rate.

Given the phase velocity results described just above, some comments about implications of these results for interpreting the various observations of solar Rossby waves are in order. Since unstable Rossby waves must have phase velocities within the range of differential rotation, meaning they cannot have fast retrograde speeds, the fast retrograde waves observed in the supergranule layer must be energetically neutral waves, which from observations appears to be so. While this does not rule out the possibility that unstable Rossby waves are present but not yet helioseismically observed, it makes such a discovery unlikely. Rossby waves in the tachocline excited by instability of the combination of differential rotation and toroidal field there acquire a phase speed equal to the rotation rate of the latitude where the narrow-banded toroidal field, from which sunspots likely emerge, peaks. This property may explain the Rossby wave speeds detected from coronal bright points (McIntosh et al. 2017; Krista et al. 2018). Since a strong toroidal field (say 10–15 kGauss) is needed to produce this result, this is strong evidence of tachocline origin of Rossby wave properties seen in the corona, and of the strength of the toroidal field in the tachocline.

For comparison with cases with rotation, we ran a brief T21 simulation with rotation switched off. This is done by simply setting ω_c = 0 in Equations (1) and (2) so Coriolis forces are absent. With Coriolis forces present, the factor $\cos \phi$ couples modes of different spherical harmonic index l. Without Coriolis forces, the nonlinear coupling of modes comes only from the velocity advection terms. Without rotation we see no evidence of an inverse cascade that piles up energy at the lowest m. Obviously without rotation, there are no nonlinear mode-mode interactions involving different m's, and hence an inverse cascade is not possible. The energy spectra in this case fail to show an organized pattern over the duration of the simulation. This should not be surprising because without rotation the fluid does not know where the poles are. Any orientation is equally likely; the system is rather degenerate, with many patterns of modes equally likely. The organizing effect of rotation is missing.