NONLINEAR DYNAMICS OF MAGNETOHYDRODYNAMIC ROSSBY WAVES AND THE CYCLIC NATURE OF SOLAR MAGNETIC ACTIVITY

The magnetic activity of the Sun exhibits several approximately periodic cycles, the 11/22 yr cycle (Schwabe/Hale cycle) being the most pronounced one. This cycle is associated with both reversals in the magnetic polarity of the dipolar mode of the Sun's magnetic field and the sunspots. Other cycles with a wide range of periodicities have been investigated in the literature, including 2 yr, 80 yr, and 200 yr cycles (Jager & Dahau 2011). A fundamental question is, how do these well-defined periods emerge from the turbulent flows in the solar interior?

In this context, modeling of solar activity has largely been based on kinematic dynamo models, in which the feedback of the Lorentz force on the velocity field is neglected. This approach has successfully simulated the Sun's periodic cycles, being capable of reproducing the spatiotemporal behavior of the solar activity depicted in the butterfly diagram (Charbonneau 2010). On the other hand, global simulations of the full magnetohydrodynamic equations have failed to reproduce the main features of solar magnetic activity (Guerrero et al. 2013). From this point of view, a fundamental theoretical description of these reversals is still elusive, and a fully magnetohydrodynamic basis for describing the solar activity is still needed.

In the kinematic interpretation of the solar dynamo, the trapping of sunspot activity and its respective migration toward the equator, represented by the butterfly diagram, can be understood as a result of the propagation of dynamo waves (Schmitt & Schüssler 2004). However, in the system of coupled magnetohydrodynamic equations, dynamo waves no longer exist as linear eigenmodes. In contrast, the purely magnetic oscillation of dynamo waves in this case should be coupled to the hydrodynamic component of the model, yielding a magnetohydrodynamic wave mode. Hence, a relevant question in the transition from a kinematic to a magnetohydrodynamic interpretation of the solar activity is, what kind of magnetohydrodynamic wave does inherit the equatorward migration properties of dynamo waves? This is a problem addressed in Section 3 of the present paper.

Another main observational constraint, which is relevant for the modeling of the cyclic behavior of the solar activity, is that cyclic oscillations are reflected in several spherical harmonic coefficients of the solar magnetic field, both axisymmetric (m = 0) and nonaxisymmetric (Knaack & Stenflo 2005). According to Knaack & Stenflo (2005), fluctuations in the spherical harmonic coefficients have shown a dominant periodicity of 22 yr as well as second (11 yr) and third (7 yr) harmonics, but a pronounced periodicity of 2–2.5 yr is also observed in some coefficients. In addition, Knaack's harmonic analysis also shows that oscillations of different harmonics of the magnetic field are out of phase and have the same period, indicating a possible nonlinear energy exchange among different harmonics of the magnetic field. Here, we present a possible theoretical explanation for these observed periodic modulations of nonaxisymmetric harmonic coefficients of the solar magnetic field based on resonant nonlinear interactions among magnetohydrodynamic Rossby waves. This theoretical description of the cyclic activity of the Sun's magnetic field may open up a new perspective to construct a more comprehensive theory of the solar cycle.

The magnetic activity of the Sun is believed to have its origin in the solar tachocline, a thin layer in the base of the convection zone and above the rigid interior of the Sun (Spiegel & Zahn 1992) at about 0.7 solar radii. Because the thickness of the tachocline is small compared to the radius of the Sun, a suitable approximation for the dynamics of this layer is the so-called shallow-water magnetohydrodynamic (MHD) approximation (Gilman 2000). This approximation is valid for motions whose vertical scales are small compared with their horizontal ones, and, consequently, magnetic and velocity fields are assumed to be two- dimensional vector fields. This is the MHD analog to the well-known shallow-water approximation applied to the Earth's ocean and atmosphere (Pedlosky 1986). The only difference is the inclusion of the magnetic field effects.

The shallow-water MHD approximation has been widely used in the study of both waves and instabilities in the solar tachocline (Gilman & Dikpati 2002; Dikpati & Gilman 2005). In particular, Rossby-like waves arise as one of the possible eigensolutions of the governing equations either in spherical geometry or in the beta-plane approximation (Zaqarashvili et al. 2007, 2010). The effect of Rossby wave motion in the context of the Sun was studied by Gilman (1969), who analyzed a kinematic MHD model where the hydrodynamic component is considered as a prescribed background flow for the magnetic field evolution. He obtained oscillatory behavior with periods on the order of a few years. However, Zaqarashvili et al. (2007), by including the feedback effect of the magnetic field on the hydrodynamic component of the model represented by the Lorentz force, showed the existence of a lower-frequency branch that was not present in Gilman's analysis. Zaqarashvili et al. (2007) studied MHD Rossby waves through the shallow-water MHD equations linearized around a resting and constant toroidal magnetic field background state, and Zaqarashvili et al. (2010) augmented this linear theory to include a latitudinal variation of the background magnetic field and a zonal mean flow represented by a differential rotation in order to analyze the possible role of wave instability in Rieger-type periodicities. The possible role of Rossby waves as a mechanism for the dynamo alpha effect in thick shells because of its helical structure, was highlighted by Avalos-Zuniga (2009). Here, we extend the MHD Rossby wave theory of Zaqarashvili et al. (2007) in order to take into account the effect of nonlinearity. For this purpose, the weakly nonlinear Wentzel–Kramers–Brillouin (WKB) theory for dispersive waves has been applied to MHD Rossby waves embedded in a resting and constant toroidal magnetic field. In order to further demonstrate the potentiality of the present nonlinear theory in describing the observed spatial–temporal variability of solar activity, we have also analyzed a particular case of the linear theory of MHD Rossby waves embedded in a latitudinally varying background state presented by Zaqarashvili et al. (2010) and investigated possible waveguides of MHD Rossby waves.

The weakly nonlinear theory of dispersive waves has been extensively studied in the context of plasma physics (Davidson 1972) and geophysical fluid dynamics (Majda 2003). According to this theory, the effect of nonlinearity is to break down the independence of wave eigenmodes, allowing their amplitudes to evolve slowly in time through energy exchanges among resonant modes. For quadratic nonlinearities, as is usually the case in most geophysical fluid dynamics and MHD problems, wave coupling takes place in a triad configuration, where the energy exchanges occur only if the modes satisfy the well-known resonant triad interaction conditions in which the wavenumbers and time frequencies of the triad components both add to zero. In this context, the nonlinear energy exchanges are described by a set of coupled ordinary differential equations in the case of a discrete spectrum of wave eigenmodes and by a set of coupled nonlinear Schrödinger equations in the case of wave packets (see, for instance, Majda 2003; Craik 1985 and references therein).

Here we apply the weakly nonlinear WKB method to the shallow-water MHD equations and analyze the reduced dynamics of a single resonant triad of MHD Rossby wave modes. In addition, we have assumed a discrete spectrum associated with periodic boundary conditions so that the wave amplitudes are described by a Hamiltonian set of three coupled ordinary differential equations, the so-called triad equations. In the absence of forcing, these equations are integrable and can be solved analytically in terms of Jacobian elliptic functions (Bretherton 1964; Craik 1985). In this case, the period of energy modulations depends on the coupling constants and the initial mode amplitudes. The coupling constants, in turn, represent the intrinsic nonlinear coupling among the triad components and depend on the background magnetic field, the wavenumbers, and the eigenfrequencies of the modes composing the resonant triad and on the density and magnetic permeability. For a representative example of a resonant triad explored here, which is composed of two magnetic and one hydrodynamic Rossby modes, more expressive energy exchanges among the triad components were found to occur for a parameter regime of the solution in which the period of the energy modulations is close to the 11 yr period or even longer periods of the solar cycle. Thus, we argue that the wave interaction mechanism explored here might be relevant for the nearly 11 yr cycle as well as for longer timescale cycles of solar magnetic activity. In the presence of convective forcing, the time evolution of the wave amplitudes of the resonant triad exhibits a growth/decay asymmetry similar to that observed in the magnetic field during the solar dynamo cycle. Also, the linear divergence terms of both the vorticity and magnetic density current equations that work as parametric forcings for MHD disturbances may potentially balance the dissipative terms and, consequently, appear to be crucial to sustain the quasiperiodic fluctuations associated with the resonant triad interactions against dissipative actions. Thus, we believe that convective forcings, possibly related to thermal inhomogeneities throughout the tachocline's boundaries, might be important in triggering and sustaining the resonant triad interactions explored here.

The analysis of the linear MHD Rossby wave theory, including a latitudinally varying mean flow performed here, points out the existence of a waveguide between 35°S and 35°N for the magnetic branch, with a maximum of the refraction index at the equator. Therefore, according to our theoretical results, the observed evolution of solar activity emerging at midlatitudes and migrating toward the equator, which is often represented by the so-called butterfly diagram, is associated with a waveguide at the equator that is due to the differential rotation of the Sun. Beyond the WKB limit, the solar differential rotation, as well as the background toroidal magnetic field of the solar tachocline, may be unstable and consequently excite planetary-scale MHD Rossby mode disturbances. Thus, apart from convective forcings related to thermal inhomogeneities along the tachocline's boundaries, global instabilities of both solar differential rotation and toroidal magnetic field can also be accounted for by another energy source that prevents resonant triad interactions of MHD Rossby modes from declining because of viscous and ohmic dissipative actions. Therefore, we argue here that the magnetic branch of Rossby waves might be a good candidate for replacing the dynamo waves in a magnetohydrodynamic interpretation. In addition, our results suggest that the observed time modulation of sunspots with a nearly 11 yr period might be associated with energy exchanges among MHD Rossby modes in resonant triad interactions.

The shallow-water equation model governs the dynamics of a homogeneous and incompressible fluid under hydrostatic balance. It constitutes an appropriate approximation of the full Navier–Stokes equations for quasi-horizontal motions in a thin layer and has long been used in studies of Earth's ocean and atmosphere (Pedlosky 1986). Its magnetohydrodynamic extension was introduced by Gilman (2000). For the sake of simplicity, we make use of Cartesian coordinates with the so-called beta-plane approximation, in which the dependence of the Coriolis parameter on latitude is given by a linear function (Pedlosky 1986). In this context, the shallow-water MHD equations are written as

$$\begin{equation} \frac{{D {\boldsymbol u}}}{{Dt}} + f {\boldsymbol u}^{\bot} = - g\nabla { h} + \frac{1}{{\mu _0 \rho }} {\boldsymbol b} \cdot \nabla {\boldsymbol b} + \nu \nabla ^2 {\boldsymbol u}\end{equation} \tag{ 1a }$$

$$\begin{equation} \frac{{D {\boldsymbol b}}}{{Dt}} = {\boldsymbol b} \cdot \nabla {\boldsymbol u} + \eta \nabla ^2 {\boldsymbol b}\end{equation} \tag{ 1b }$$

$$\begin{equation} \frac{{Dh}}{{Dt}} = h\, {\rm div}\, {\boldsymbol u}\end{equation} \tag{ 1c }$$

$$\begin{equation} {\rm div} ({{ h}{\boldsymbol b}} ) = 0, \end{equation} \tag{ 1d }$$

where ${\boldsymbol u}\left({x,y,t} \right) = \left({u_1,u_2 } \right)$ is the two-dimensional velocity field, ${\boldsymbol u}^ \bot \left({x,y,t} \right) = \left({ - u_2,u_1 } \right)$, ${\boldsymbol b}\left({x,y,t} \right) = \left({b_1,b_2 } \right)$ represents the two-dimensional magnetic field, and the scalar field h(x, y, t) refers to the thickness of the fluid layer; g is the gravity acceleration, ρ is the density, μ₀ refers to the magnetic permeability, and f(y) = f₀ + βy corresponds to the Coriolis parameter, with both f₀ and β being assumed to be constant. The parameters η and ν are magnetic diffusivity and viscosity, respectively. Equation (1d) is the constraint associated with the requirement that the magnetic field be parallel to the upper free surface (Gilman 2000).

Because we aim to adopt model Equations (1) to study MHD Rossby waves, it is convenient to take the curl of Equations (1a) and (1b), where the two-dimensional vorticity and the magnetic density current are defined as ζ = (∂u₂/∂x) − (∂u₁/∂y) and j = (∂b₂/∂x) − (∂b₁/∂y), respectively. Therefore, we obtain the following scalar equations for the vorticity and the magnetic density current:

$$\begin{equation} \frac{{D\left({f + \zeta} \right)}}{{Dt}} + \left({f + \zeta} \right){\rm div}\, {\boldsymbol u} = \frac{1}{{\mu _0 \rho }}({{\boldsymbol b} \cdot \nabla j + j{\rm div}\, {\boldsymbol b}} ) + \nu \nabla ^2 \zeta \end{equation} \tag{ 2a }$$

$$\begin{equation} \frac{{Dj}}{{Dt}} = {\boldsymbol b}\cdot \nabla \zeta + \zeta {\rm div}\, {\boldsymbol b} + j {\rm div}\, {\boldsymbol u} + \eta \nabla ^2 { j}{\rm.}\end{equation} \tag{ 2b }$$

Because our analysis focuses on planetary-scale perturbations (Rossby waves), to a first approximation, both velocity and magnetic fields could be assumed as divergenceless. However, we will retain the divergence of both velocity and magnetic fields in the equations above as small forcing terms associated with thermal anomalies in the boundaries of the solar tachocline. In fact, because in an incompressible fluid the continuity equation may be approximated as

$$\begin{equation*} {\rm div}\, {\boldsymbol u} = - \frac{{\partial u_3 }}{{\partial z}},\end{equation*}$$

with u₃ representing the vertical velocity field, an upwelling fluid parcel in a confined layer will result in horizontal divergence; hence horizontal divergence can be used as a measure of thermal convection. In order to keep the effect of thermal convection without dealing with the thermodynamic equations, we define $D_u \left({x,y,t} \right) = {\rm div}\, {\boldsymbol u}$ as a prescribed term representing a Rossby wave convective source. This assumption of keeping the horizontal divergence as a prescribed Rossby wave source has been extensively used in geophysical fluid dynamics in the study of teleconnection patterns in the atmosphere (see, for instance, Hoskins & Karoly 1981; Hoskins & Ambrizzi 1993 and many others). In addition, these numerical studies of atmospheric teleconnection patterns using the vorticity equation have disregarded the term $\zeta {\rm div}\, {\boldsymbol u}$ based on arguments of scale analysis. However, as shall be discussed in Section 5, this parametric forcing term of the vorticity balance is important in our nonlinear analysis because it can potentially balance the dissipative terms and, consequently, sustain the nonlinear oscillations in the presence of dissipation. Thus, we will consider here the effect of both prescribed ($f{\rm div}\, {\boldsymbol u}$) and parametric ($\zeta{\rm div}\, {\boldsymbol u}$) forcing terms of the vorticity equation.

Similarly, the divergence of the horizontal magnetic field in Equations (2a)–(2b) will also be assumed as a prescribed function representing magnetic anomalies associated with inhomogeneities in the boundaries of the solar tachocline. In fact, from Equation (1d) the divergence of the horizontal magnetic field can be written as

$$\begin{equation*} D_B \left({x,y,t} \right) = {\rm div}\, {\boldsymbol b}\left({x,y,t} \right) = \frac{{{\boldsymbol b}\cdot \nabla { h}}}{{h}},\end{equation*}$$

which shows that this prescribed forcing can be particularly important in the case of a variable tachocline thickness (see Guerrero et al. 2007 for a discussion on the effects of a prolate tachocline).

Therefore, with the assumptions mentioned above, Equations (2) become

$$\begin{equation} \frac{{D\left({f + \zeta} \right)}}{{Dt}} - \frac{1}{{\mu _0 \rho }}{\boldsymbol b}\cdot \nabla j = - {fD}_{u} - \zeta{D}_{u} + \frac{{{jD}_{B} }}{{\mu _0 \rho }} + {\rm }\nu \nabla ^2 \zeta\end{equation} \tag{ 3a }$$

$$\begin{equation} \frac{{Dj}}{{Dt}} - {\boldsymbol b}.\nabla \zeta = \zeta {D}_{B} + {jD}_{u} + \eta \nabla ^2 {j}, \end{equation} \tag{ 3b }$$

where D_u(x, y, t) = ${\rm div}\, {\boldsymbol u}$ and D_B(x, y, t) = ${\rm div}\, {\boldsymbol b}\left({x, y, t} \right)$.

3.1. Resting and Constant Magnetic Field Background State

MHD Rossby waves are eigensolutions of the linearized and unforced version of Equations (3). Let us consider, for instance, the linearized and unforced version of (3) for a background state with no motion and a constant toroidal magnetic field. Note that this background state represents an exact stationary solution of Equation (3). Thus, considering infinitesimal amplitude perturbations around a background state at rest and with a constant toroidal magnetic field (i.e., ${\boldsymbol B} = \left({B_0,0} \right)$), the linearized and unforced (D_u = D_B ≡ 0) problem governing the time evolution of these perturbations can be written in terms of stream function ψ and vector potential A as

$$\begin{equation} \frac{\partial }{{\partial t}}\left({\begin{array}{*{20}c} {\nabla ^2 \psi } \\ {\nabla ^2 A} \\ \end{array}} \right) = {\cal L}\left({\begin{array}{*{20}c} \psi \\ A \\ \end{array}} \right)\end{equation} \tag{ 4a }$$

where the linear differential operator ${\cal L}$ is given by

$$\begin{equation} {\cal L} = \left[ {\begin{array}{@{}cc@{}} {\beta \frac{\partial }{{\partial x}}} & {\frac{1}{{\mu _0 \rho }}B_0 \frac{{\partial \nabla ^2 }}{{\partial x}}} \\ {B_0 \frac{{\partial \nabla ^2 }}{{\partial x}}} & 0 \\ \end{array}} \right]\end{equation} \tag{ 4b }$$

and the streamfunction ψ and vector potential A satisfy the relations ∇²ψ = ζ and ∇²A = j, respectively. Taking the Fourier transform of ${\cal L}$, we obtain the algebraic operator

$$\begin{equation} {\cal L}^* = i\left[ {\begin{array}{@{}cc@{}} {\beta k} & {\frac{1}{{\mu _0 \rho }}B_0 k\left| k \right|^2 } \\ {B_0 k\left| k \right|^2 } & 0 \\ \end{array}} \right],\end{equation} \tag{ 5 }$$

whose eigenvalues correspond to the two possible branches of the MHD Rossby wave dispersion relation and are determined by the characteristic equation

$$\begin{equation} det\left({ - iI\omega - \frac{1}{{\left| {{\boldsymbol k}} \right|^2 }}{\cal L}^* \left({i{\boldsymbol k}} \right)} \right) = 0.\end{equation} \tag{ 6 }$$

In the equations above, ${\boldsymbol k} = \left({k,l} \right)$ is the vector wavenumber. The two solutions of the characteristic Equation (6) are given by

$$\begin{equation} \omega _ + = \frac{{ - \left({\beta + \sqrt {\beta^2 + 4v_A^2 \left| {{\boldsymbol k}} \right|^4 } } \right)k}}{{2\left| {{\boldsymbol k}} \right|^2 }}\end{equation} \tag{ 7a }$$

and

$$\begin{equation} \omega _ - = \frac{{ - \left({\beta - \sqrt {\beta^2 + 4v_A^2 \left| {{\boldsymbol k}} \right|^4 } } \right)k}}{{2\left| {{\boldsymbol k}} \right|^2 }},\end{equation} \tag{ 7b }$$

where $v_A = ({{B_0 }}/{{\sqrt {\mu _0 \rho } }})$ refers to the Alfvén wave speed. The branches ω₊ and ω₋ refer to the fast hydrodynamic mode and the slow magnetic mode, respectively (Zaqarashvili et al. 2007). In fact, for small wavenumbers such that ${{4v_A^2 \left| {{\boldsymbol k}} \right|^4 }}/{{\beta }} \ll 1$, the dispersion relations of the two branches can be approximated by $\omega _ + \approx - ({{\beta k}}/{{\left| {{\boldsymbol k}} \right|^2 }})$ and $\omega _ - \approx v_A^2 ({{\left| {{\boldsymbol k}} \right|^2 k}}/{\beta })$. Thus, the asymptotic form of ω₊in the long-wave limit corresponds to the dispersion relation of purely hydrodynamic Rossby waves, whereas the asymptotic form of ω₋ in the same limit corresponds to the dispersion relation of Alfvén waves. For this reason we shall label mode "−" as the magnetic mode and mode "+'' as the hydrodynamic mode. In the calculations displayed in this paper we have set the parameters β = 0, 95.10⁻¹³ m⁻¹ s⁻¹, and v_A = 200 m s⁻¹.

The corresponding right eigenvectors of matrix ${1}/{{\left| {{\boldsymbol k}} \right|^2 }}{\cal L}^*$ are given by

$$\begin{equation} {\boldsymbol R}^ - \left({{\boldsymbol k}_i } \right) = \left({\begin{array}{@{}c@{}} {\omega _ - \left({{\boldsymbol k}_i } \right)k_i } \\ {B_0 } \\ \end{array}} \right)\end{equation} \tag{ 8a }$$

and

$$\begin{equation} {\boldsymbol R}^ + \left({{\boldsymbol k}_i } \right) = \left({\begin{array}{@{}c@{}} {\omega _ + \left({{\boldsymbol k}_i } \right)k_i } \\ {B_0 } \\ \end{array}} \right).\end{equation} \tag{ 8b }$$

For each wavenumber ${\boldsymbol k}$, the eigenvectors defined above form an orthogonal basis in ${{\mathbb C}} ^2$ space under a suitable inner product defined by

$$\begin{equation} \langle {\boldsymbol u}_1,{\boldsymbol u}_2 \rangle _{{\mathbb C}{^2 }} = \left({\bar u_1 u_2 + \frac{1}{{\mu _0 \rho }}\bar v_1 v_2 } \right),\end{equation} \tag{ 9 }$$

where ${\boldsymbol u}_1 = \left({u_1,v_1 } \right)$ and ${\boldsymbol u}_2 = \left({u_2,v_2 } \right)$ represent arbitrary elements of $\mathbb{C}^2$ vector space and the over bar denotes the complex conjugate. The norm defined by this inner product represents the energy of a given eigenmode.

3.2. Latitudinally Varying Mean Zonal Flow and Waveguides for MHD Rossby Waves

Schmitt & Schüssler (2004) have shown that dynamo waves, propagating oscillations that arise as eigensolutions of the induction equation, can reproduce the equatorward migration of the solar activity represented in the butterfly diagram. Here we argue that the magnetic branch of MHD Rossby waves might be a good magnetohydrodynamic counterpart of the dynamo waves because they are confined in a belt around the equator and are refracted toward the equator, as will be shown in this section.

Waveguides associated with westerly jet streams play a central role in the propagation of Rossby wave disturbances and, consequently, in establishing the main teleconnection patterns in Earth's atmosphere (Hoskins & Ambrizzi 1993; Hoskins & Karoly 1981). In addition, latitudinally varying shear flows constitute an essential ingredient of the Sun's dynamics, and, as a consequence, solar MHD Rossby waves are seldom embedded in a resting background state. Thus, in this subsection we examine the linear theory of MHD Rossby waves embedded in a latitudinally varying background state presented by (Zaqarashvili et al. 2010) to investigate possible MHD Rossby waveguides in the solar dynamo. For simplicity, we consider the special case of a latitudinally varying zonal flow and a constant toroidal magnetic field.

In this context, let us now consider infinitesimal amplitude disturbances embedded in a background state characterized by a constant magnetic field in the x direction and by a slowly varying zonal flow $\bar U\left(y \right)$. In this case, according to the WKB approximation, the MHD Rossby wave dispersion relation is given by

$$\begin{equation} \omega ({ {\boldsymbol k}, y} ) = \bar U\left(y \right)k - \frac{{\left({\beta ^* \pm \sqrt {\beta ^{*2} + 4v_A^2 \left| {{\boldsymbol k}} \right|^4 } } \right)k}}{{2\left| {{\boldsymbol k}} \right|^2 }}, \end{equation} \tag{ 10 }$$

where $\beta ^* = \beta - \bar U \left(y \right)^{\prime\prime}$ is the meridional gradient of absolute vorticity.

For simplicity in exposition, we consider the asymptotic forms of the two branches in the long-wave limit introduced in the last section. For the hydrodynamic branch $\omega _ + \approx \bar Uk - ({{\beta ^* k}}/{{\left| {{\boldsymbol k}} \right|^2 }})$, and the total wavenumber is given by

$$\begin{equation} K_ + = \sqrt {\frac{{\beta ^* }}{{\bar U - c_ + }}}.\end{equation} \tag{ 11a }$$

On the other hand, for the magnetic branch $\omega _ - \approx \bar Uk - v_A^2 ({{\left| {{\boldsymbol k}} \right|^2 k}}/{{\beta ^* }})$, and the total wavenumber is

$$\begin{equation} K_ - = \sqrt {\frac{{({\bar U - c_ - } )\beta ^* }}{{v_A^2 }}}.\end{equation} \tag{ 11b }$$

In the equations above, c_± = ω_±/k refers to the zonal phase speed of the modes. Wave propagation is allowed only where the total wavenumber K is real because in the latitudes where the total wavenumber K is imaginary, the solution is evanescent.

In general, the phase speed of Rossby waves is smaller in magnitude than the velocity of the background flow. Consequently, a usual approximation of the equations above is the assumption of a stationary wave (c = 0). In this case, the total wavenumber defined above is replaced by the stationary wavenumber K_s. The concept of stationary wavenumber K_s has been extensively used in geophysical fluid dynamics problems to determine possible Rossby waveguides in the presence of realistic background zonal flows and has been successfully used to interpret some model results on the atmospheric response to heat or orographic forcing in a qualitative fashion (e.g., Hoskins & Karoly 1981; Hoskins et al. 1977; Hoskins & Ambrizzi 1993).

Figure 1 shows the latitudinal profile of K_s² referred to the magnetic branch for a zonal flow $\bar U$(y) obtained from a theoretical profile of the solar differential rotation taken from MacGregor & Charbonneau (1997). One notices in Figure 1 the existence of critical latitudes at around 80° and 35° in both hemispheres. These critical latitudes refer to regions where the solution changes from a wave propagation (K_s² > 0) regime to an evanescent (K_s² < 0) regime. This is in agreement with observations of solar magnetic activity, which show that sunspots are present in the tropical belt. Furthermore, for the magnetic branch, the stationary wavenumber K_s is maximal at the equator. According to the WKB theory of Rossby wave propagation in a slowly varying medium (Bühler 2009; Held 1983; Hoskins & Ambrizzi 1993, Section 2), Rossby waves are refracted toward regions of maximum stationary wavenumber such that these regions work as waveguides. Thus, according to the result displayed in Figure 1, in the solar dynamo context the equator acts as a waveguide for the magnetic branch of MHD Rossby waves such that waves generated at the critical latitudes around ±35° are refracted toward this latitude of maximum stationary wavenumber. This feature is also consistent with observations of the solar activity whose latitude versus time structure is depicted by the butterfly diagram.

Zoom In Zoom Out Reset image size

The behavior described above has a striking agreement with the structure of the activity zones (drifting groups of sunspots) as described by Joy's law (Hale et al. 1919). Joy's law predicts that the tilt of a bipolar activity zone will vary from approximately 15° as they appear around 35° of latitude and will decrease as the activity zone approaches the equator, becoming parallel to the equator at 0° of latitude. This suggests that activity zones might be traveling MHD Rossby waves and Joy's law might represent the ray of the Rossby wave as it travels in the slowly varying medium due to differential rotation.

In this section we present the nonlinear MHD Rossby wave interaction theory as a fundamental theoretical description of observed periodic cycles of nonaxisymmetric components of solar magnetic activity. For simplicity, we consider wave disturbances embedded in the background state at rest and with a constant toroidal magnetic field ${\boldsymbol B} = \left({B_0,0} \right)$. Let us consider, for instance, finite amplitude perturbations around the background state described above in the unforced and inviscid version of model Equations (3). In this case, the governing equations for these perturbations are obtained by restoring the quadratic nonlinear terms into Equations (4):

$$\begin{equation} \frac{\partial }{{\partial t}}\left({\begin{array}{@{}c@{}} {\nabla ^2 \psi } \\ {\nabla ^2 A} \\ \end{array}} \right) + {\cal L}\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right) = {\cal B}\left({\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right),\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right)} \right).\end{equation} \tag{ 12a }$$

In Equation (12), ${\cal L}$ is the linear operator given by Equation (4b) and ${\cal B}$ corresponds to the quadratic nonlinear operator, given by

$$\begin{equation} {\cal B}\left({\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right),\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right)} \right) \!\,{=}\,\! \left[ {\begin{array}{@{}c@{}} { - ({{\boldsymbol \nabla} ^ \bot \psi } ).{\boldsymbol \nabla} \nabla ^2 \psi + \frac{1}{{\mu _0 \rho }} ({{\boldsymbol \nabla} ^ \bot A} ).{\boldsymbol \nabla} \nabla ^2 A} \\ { - ({{\boldsymbol \nabla} ^ \bot \psi } ).{\boldsymbol \nabla} \nabla ^2 A + ({{\boldsymbol \nabla} ^ \bot A} ).{\boldsymbol \nabla} \nabla ^2 \psi } \\ \end{array}} \right].\end{equation} \tag{ 12b }$$

The standard asymptotic treatment of Equation (12) in the weakly nonlinear regime requires elimination of secular terms that are due to resonances arising in the right-hand side of Equation (12a). To eliminate secular solutions, we enforce conditions for the slow evolution of the wave amplitudes given by the Fredholm alternative (Majda 2003, Chapter 5). For a set of three modes whose wavenumbers ${\boldsymbol k}_1$, ${\boldsymbol k}_2$, and ${\boldsymbol k}_3$and corresponding eigenfrequencies $\omega _{\alpha _1 } \left({{\boldsymbol k}_1 } \right)$, $\omega _{\alpha _2 } \left({{\boldsymbol k}_2 } \right)$, and $\omega _{\alpha _3 } \left({{\boldsymbol k}_3 } \right)$ satisfy the resonance conditions

$$\begin{equation*} {\boldsymbol k}_1 + {\boldsymbol k}_2 + {\boldsymbol k}_3 = 0\end{equation*}$$

and

$$\begin{equation*} \omega _{\alpha _1 } \left({{\boldsymbol k}_1 } \right) + \omega _{\alpha _2 } \left({{\boldsymbol k}_2 } \right) + \omega _{\alpha _3 } \left({{\boldsymbol k}_3 } \right) = 0,\end{equation*}$$

the respective wave amplitudes Ω₁, Ω₂, and Ω₃ evolve in time according to

$$\begin{equation} \frac{{d\Omega _1 }}{{dt}} = {\rm\mathbb C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)} \overline \Omega _2 \overline \Omega _3 \end{equation} \tag{ 13a }$$

$$\begin{equation} \frac{{d\Omega _2 }}{{dt}} = {\rm\mathbb C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)} \overline \Omega _1 \overline \Omega _3 \end{equation} \tag{ 13b }$$

$$\begin{equation} \frac{{d\Omega _3 }}{{dt}} = {\rm\mathbb C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({\alpha _3,\alpha _1 \alpha _2 } \right)} \overline \Omega _1 \overline \Omega _2.\end{equation} \tag{ 13c }$$

In the equations above, α_i, i = 1, 2, 3, refers to the branch (α_i = + for hydrodynamic branch and α_i = − for magnetic branch), the over bar indicates the complex conjugate, and ${\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)}$, ${\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)}$, and ${\rm C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({,\alpha _3,\alpha _1 \alpha _2 } \right)}$ are the nonlinear interaction coefficients, given by

$$\begin{eqnarray} {\rm C}_{({{\boldsymbol k},{\boldsymbol k}_1,{\boldsymbol k}_2 } )}^{\left({\alpha,\alpha _1,\alpha _2 } \right)} &=& \langle {{\boldsymbol L}^\alpha ({{\boldsymbol k}_1 + {\boldsymbol k}_2 } ),{\cal B}^* ({{\boldsymbol R}^{\alpha _1 } ({{\boldsymbol k}_1 } ),{\boldsymbol R}^{\alpha _2 } ({{\boldsymbol k}_2 } )} )} \rangle \nonumber\\ & =& \left\lbrace {\left[ {\left({{\boldsymbol k}_2^ \bot .{\boldsymbol k}_1 } \right)\left({\left| {{\boldsymbol k}_2 } \right|^2 - \left| {{\boldsymbol k}_1 } \right|^2 } \right)} \right]}\right. \nonumber\\ &&\left. \times {\left({\omega _{\alpha _1 } \left({{\boldsymbol k}_1 } \right)\omega _{\alpha _2 } \left({{\boldsymbol k}_2 } \right) + v_A^2 } \right)} \right\rbrace\omega _{\alpha _3 } \left({{\boldsymbol k}_1 + {\boldsymbol k}_2 } \right), \quad\quad \end{eqnarray} \tag{ 14 }$$

where ${\boldsymbol L}^\alpha$ stands for the left eigenvector of the linear operator${\cal L}$ and 〈〉 represents the inner product (9).

Equations (13) are integrable and can be solved analytically in terms of Jacobian elliptic functions (Craik 1985). In general, the solution of Equation (13) is periodic, except for the special case of m = 1, as will be shown below.

To obtain the analytical solution of system (13), let us consider mode 1 as the one with the highest absolute value interaction coefficient and mode 3 as having the smallest absolute value interaction coefficient, that is,

$$\begin{equation*} \left| {{\rm C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({,\alpha _3,\alpha _1 \alpha _2 } \right)} } \right|{\rm < }\left| {{\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)} } \right|{\rm < }\left| {{\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)} } \right|.\end{equation*}$$

From Equation (14) one notices that the interaction coefficients are proportional to their respective eigenfrequencies, so the condition above implies that

$$\begin{equation*} \left| {\omega _{\alpha _3 } \left({{\boldsymbol k}_3 } \right)} \right| < \left| {\omega _{\alpha _2 } \left({{\boldsymbol k}_2 } \right)} \right| < \left| {\omega _{\alpha _1 } \left({{\boldsymbol k}_1 } \right)} \right|.\end{equation*}$$

The conditions above, or equivalently $({{{\rm C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({,\alpha _3,\alpha _1 \alpha _2 } \right)} }}/{{\omega _{\alpha _3 } \left({{\boldsymbol k}_3 } \right)}}) = ({{{\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)} }}/{{\omega _{\alpha _2 } \left({{\boldsymbol k}_2 } \right)}}) = ({{{\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)} }}/{{\omega _{\alpha _1 } \left({{\boldsymbol k}_1 } \right)}})$, represent a characteristic feature of conservative resonant triads (Craik 1985). In this context, for the initial condition Ω₁(t = 0) = 0, Ω₂(t = 0) = Ω₂(0), and Ω₃(t = 0) = Ω₃(0),¹ the solution can be written as (Domaracki & Loesch 1977)

$$\begin{equation} {\rm \Omega }_1 \left(t \right) = {\rm \Omega }_2 \left(0 \right)\left| {\frac{{{\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)} }}{{{\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)} }}} \right|^{\frac{1}{2}} sn\left({\frac{{{ \tilde T}}}{m}} \right)\end{equation} \tag{ 15a }$$

$$\begin{equation} {\rm \Omega }_2 \left(t \right) = {\rm \Omega }_2 \left(0 \right)cn\left({\frac{{{ \tilde T}}}{m}} \right)\end{equation} \tag{ 15b }$$

$$\begin{equation} {\rm \Omega }_3 \left(t \right) = {\rm \Omega }_3 \left(0 \right)dn\left({\frac{{{\tilde T}}}{m}} \right),\end{equation} \tag{ 15c }$$

where sn, cn, and dn are the Jacobian elliptic functions, with the argument corresponding to the rescaled time ${ \tilde T}$, given by

$$\begin{equation} {\tilde T} = {\rm \Omega }_3 \left(0 \right)\left| {{\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)} {\rm C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({\alpha _3,\alpha _1 \alpha _2 } \right)} } \right|t.\end{equation} \tag{ 16 }$$

Also, the parameter m of the Jacobian elliptic functions above is written as

$$\begin{equation} m = \left| {\frac{{{\rm C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({\alpha _3,\alpha _1 \alpha _2 } \right)} }}{{{\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)} }}} \right|\left({\frac{{{\rm \Omega }_2 \left(0 \right)}}{{{\rm \Omega }_3 \left(0 \right)}}} \right)^2.\end{equation} \tag{ 17 }$$

In the limiting case of m = 0, the solution (15) reduces to

$$\begin{equation} {\rm \Omega }_1 \left(t \right) = {\rm \Omega }_2 \left(0 \right)\left| {\frac{{{\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)} }}{{{\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)} }}} \right|^{\frac{1}{2}} {\rm sin}({{ \tilde T}} )\end{equation} \tag{ 18a }$$

$$\begin{equation} {\rm \Omega }_2 \left(t \right) = {\rm \Omega }_2 \left(0 \right){\rm cos}({{ \tilde T}} )\end{equation} \tag{ 18b }$$

$$\begin{equation} {\rm \Omega }_3 \left(t \right) = {\rm \Omega }_3 \left(0 \right) = {\rm cte}{\rm.}\end{equation} \tag{ 18c }$$

The m = 0 case occurs for resonant triads whose interaction coefficient of the third component (13c) is zero, which in turn refers to the case where $\omega _{\alpha _3 } \left({{\boldsymbol k}_3 } \right) = 0$ (zero-frequency wave). For such triads, the amplitude of the third wave does not change in time. Consequently, this mode acts as a catalyst for the interaction of the other two waves, which exchange energy periodically in time, with the period of the energy exchange varying linearly with the amplitude of the third wave. In this case, the higher the amplitude of mode 3, the shorter the interaction period.

According to the dispersion relation of MHD Rossby waves (Equation (7)), the modes with zerofrequency are those having an axisymmetric x structure (k = 0). For instance, dipole modes (i.e., the two branches of the MHD Rossby mode with zonal wavenumber k = 0 and meridional wavenumber l = 1) constitute an example of such stationary modes. Thus, the present theory points out that the dipole meridional structure MHD Rossby modes, which might contribute to the dipole harmonic of the Sun's magnetic field, can act as a catalyst for nonaxisymmetric propagating MHD Rossby modes to interact and exchange energy with each other. An example of such a resonant triad is composed of hydrodynamic Rossby waves with corresponding wave numbers (0, 1) (zero-frequency mode) and (1, 8) and (1, 9), both with frequency of approximately 8.45 × 10⁻⁶ Hz.

According to Equation (17), Equations (18) can also be obtained for resonant triads containing only propagating modes as the asymptotic limit of Equations (15) when |Ω₂(0)| ≪ |Ω₃(0)|. However, in this case, Equation (18) shows that modes 1 and 2 exchange energy, but their amplitudes remain small compared with the amplitude of mode 3, that is, |Ω_{1, 2}(t)| ≪ |Ω₃(0)| for all t > 0. Consequently, the triad interaction in this case does not yield significant changes in the dynamical fields ${\boldsymbol u}$ and ${\boldsymbol b}$. Thus, we believe that the MHD Rossby wave interactions in this regime, whose interaction period is on the order of a few years or shorter, might not be relevant for the observed cycles of solar magnetic activity.

In the other limit (m = 1), the amplitude of mode 1 grows rapidly in time at the expense of the amplitude decrease of modes 2 and 3. After mode 1 achieves its maximum amplitude level, the solution attains a stationary regime where the amplitudes of the three waves do not change in time. Hence, for resonant triads involving three propagating MHD Rossby modes, the period of the wave amplitude modulations can be arbitrarily long for m → 1 or short for m → 0, with the parameter m depending crucially on the interaction coefficients and on the initial condition.

In general, the initial amplitudes and the interaction coefficients are such that 0 < m < 1, with the three modes undergoing significant energy modulations in this regime (Domaracki & Loesch 1977). In fact, Raupp et al. (2008) have shown that the resonant triad interactions among equatorial atmospheric waves in this parameter regime are capable of yielding significant amplitude modulations in the wind and pressure fields. Thus, we illustrate our analysis with a representative example of a resonant triad in the present model for the initial amplitudes of modes 2 and 3 chosen to fall into this regime. This representative example of a resonant triad in the present model is displayed in Table 1. This triad is composed of two hydrodynamic modes and one magnetic mode. For this triad illustrated in Table 1, one notices that the coupling coefficients do satisfy the relation ${\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)} + {\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)} + {\rm C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({,\alpha _3,\alpha _1 \alpha _2 } \right)} = 0$, which ensures that the modes exchange energy but the total energy of the triad remains constant. Figure 2 displays the time evolution of mode amplitudes for the triad illustrated in Table 1 composed of three propagating MHD Rossby modes. In this resonant triad interaction, mode 1 is a hydrodynamic mode with zonal wavenumber k = 1 and meridional wavenumber l = 1, essentially representing a planetary-scale wave mode. Although mode 1 does not have an axisymmetric x structure, it exhibits a dipole meridional structure. Hence, this mode might contribute to the observed dipole mode activity of the Sun's magnetic field. The remaining modes of the triad have more complex meridional structures, contributing to the energetics of higher harmonics of the Sun's magnetic field; namely, mode 2 is a hydrodynamic mode with zonal wavenumber 2 and meridional wavenumber l = 18; mode 3 is a magnetic mode with meridional wavenumber 19 and zonal wavenumber 3. For this triad interaction involving two hydrodynamic and one magnetic modes, the planetary-scale hydrodynamic mode is the highest frequency mode and, consequently, the unstable mode of the triad. Consequently, its energy (amplitude) always changes in time in the opposite direction of the remaining triad components; that is, this mode always supplies energy to or receives energy from the remaining two waves. When the initial state is projected only onto modes 2 and 3 (Figure 2), the planetary-scale hydrodynamic wave is excited, and its amplitude (energy) grows in time at the expense of the remaining modes of the triad. When the planetary-scale hydrodynamic Rossby mode attains its maximum energy level, this mode begins to excite the shorter-scale waves again; that is, its amplitude begins to decrease in time at the expense of the amplitude growth of modes 2 and 3. After mode 1 attains its minimum energy level again (Ω₁(t) = 0), its amplitude grows again but in the opposite polarity, in detriment of the energy loss of modes 2 and 3. The complete cycle, including the two polarities (positive and negative amplitudes of the dipolar wave mode), has a period on the order of 20 yr, whereas the period of the energy (which is proportional to the amplitude squared) modulation is one-half of the period of the amplitude modulation.

Zoom In Zoom Out Reset image size

The scenario described above may explain in a simplified fashion certain global magnetic patterns on the Sun that are part of the dynamo, such as the so-called active longitudes (Dikpati & Gilman 2005 and references therein). The active longitudes refer to some specific longitude intervals where the cyclic magnetic activity of the Sun is more pronounced. A possible conjecture for the existence of such specific longitude intervals of pronounced periodic signals is that nonaxisymmetric modes of magnetic activity may be responsible for or amplify (some pre-existing) magnetic oscillations. This theoretical conjecture is supported by the observational evidence of Knaack & Stenflo (2005), who found significant periodic oscillations in nonaxisymmetric spherical harmonic coefficients of the solar magnetic field. According to the theory described here, these observed periodic modulations of nonaxisymmetric spherical harmonic coefficients of the solar magnetic field presented by Knaack & Stenflo (2005) are due to resonant triad interactions containing nonaxisymmetric MHD Rossby waves. In fact, according to the results displayed in Figure 2, when the dipole harmonic (l = 1) of the meridional structure of the solar magnetic field is at its maximum amplitude level, the higher meridional harmonics reach their minimum amplitude level and vice versa.

5.1. Constant Convective Forcing

Let us consider now the same finite amplitude analysis carried out in the previous section but with the inclusion of the prescribed divergence term ($f{\rm div}\, {\boldsymbol u}$) in the vorticity balance equation. This approach of considering only ($f{\rm div}\, {\boldsymbol u}$) as a prescribed Rossby wave source is in the spirit of theoretical and modeling studies on the impact of thermal or topographic forcings on the large-scale atmospheric dynamics, as previously discussed. In this context, Equation (12a) is replaced by

$$\begin{equation} \frac{\partial }{{\partial t}}\left({\begin{array}{@{}c@{}} {\nabla ^2 \psi } \\ {\nabla ^2 A} \\ \end{array}} \right) + {\cal L}\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right) = {\cal B}\left({\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right),\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right)} \right) + F,\end{equation} \tag{ 19 }$$

where the operators ${\cal L}$ and ${\cal B}$ are given by Equations (4b) and (12b), respectively; F corresponds to the prescribed forcing vector defined according to

$$\begin{equation} F = \left[ {\begin{array}{@{}c@{}} {{fD}_u \left({x,y,t} \right)} \\ 0 \\ \end{array}} \right].\end{equation} \tag{ 20 }$$

In the equations above, the forcing ${fD}_u = f{\rm div}\, {\boldsymbol u}$ is required to have a zero average in order to ensure conservation of angular momentum, meaning that there will always be a region of divergence and a region of convergence. By expanding F in a Fourier basis, it is readily seen that this prescribed forcing will be composed of a series of wavepatterns, each one with a positive and a negative region of divergence.

A possible source of divergence in the tachocline is the baroclinic instability, an important process for the generation of vorticity through conversion of potential energy (Ghil & Childress 1987). Baroclinic instability was recently shown to occur on modes with a global scale (wavenumbers 1 and 2; Dikpati & Gilman 2014). Different from the barotropic Rossby waves described here, baroclinic Rossby waves are associated with horizontally divergent motions, so the prescribed divergence term in Equation (20) could be seen as a representation of the net effect of baroclinic motions on barotropic modes, especially motions associated with unstable modes that will mostly contribute to horizontal divergence.

The asymptotic treatment of Equation (19) in the weakly nonlinear and weak forcing regime yields the following reduced equations for the amplitudes of three resonantly interacting waves:

$$\begin{equation} \frac{{d{\rm \Omega }_1 }}{{dt}} = {\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)} \overline {\Omega} _2 \overline {\Omega} _3 + {\rm \Gamma }_1 \end{equation} \tag{ 21a }$$

$$\begin{equation} \frac{{d{\rm \Omega }_2 }}{{dt}} = {\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)} \overline {\Omega} _1 \overline {\Omega} _3 + {\rm \Gamma }_2 \end{equation} \tag{ 21b }$$

$$\begin{equation} \frac{{d{\rm \Omega }_3 }}{{dt}} = {\rm C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({,\alpha _3,\alpha _1 \alpha _2 } \right)} \overline {\Omega} _1 \overline {\Omega} _2 + {\rm \Gamma }_3, \end{equation} \tag{ 21c }$$

where the coupling coefficients ${\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)}$, ${\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)}$, and ${\rm C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({,\alpha _3,\alpha _1 \alpha _2 } \right)}$are given by Equation (14) and Γ₁, Γ₂, and Γ₃ correspond to the projection of forcing vector F onto modes 1, 2, and 3, respectively, meaning that only the special scales of the divergence corresponding to the spatial scale of waves 1, 2, and 3 will be relevant to our analysis.

Figure 3 shows the time evolution of the mode amplitudes of the same resonant triad and initial condition as in Figure 2, but with the inclusion of a Rossby wave source F = [fD(x, y, t), 0]^T that projects only onto the unstable mode of the triad (mode 1, the large-scale dipole mode). According to the dynamics of periodically forced resonant triads, this situation where the forcing resonates only with the unstable mode is the condition for total energy to remain bounded in time (Raupp & Dias 2009; Harris et al. 2012). Apart from producing a period doubling in the time evolution of the amplitude of mode 3, from Figure 3 one notices that the effect of the forcing on this resonant interaction is to yield an asymmetry between the growth and decay phases of the mode amplitudes. This growth and decay asymmetry is also an observed characteristic of the time evolution of solar magnetic activity. Thus, the results suggest that a vorticity source on the most unstable mode may lead sunspots to alternate small and large energy peaks during the solar cycle.

Zoom In Zoom Out Reset image size

5.2. Combined Effect of Linear Forcing and Dissipation

In a realistic scenario, the flows in the solar interior are subjected to dissipation, which may contribute to the amplitude modulations associated with the nonlinear wave interactions discussed thus far to die out with time. However, as mentioned in Section 2, the fluid at the tachocline is also heated from below, which may act as a forcing mechanism for MHD Rossby waves, preventing them from vanishing through dissipative losses. Thus, in what follows, we describe the effects of the interaction of the complete divergence forcing and dissipation on the nonlinear MHD Rossby wave interactions.

For this purpose, we consider the same finite-amplitude perturbation analysis performed thus far but add the linear divergence and dissipative terms of Equation (3). In this way, the governing equations for the MHD disturbances are

$$\begin{equation} \frac{\partial }{{\partial t}}\left({\begin{array}{@{}c@{}} {\nabla ^2 \psi } \\ {\nabla ^2 A} \\ \end{array}} \right) + {\cal L}\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right) = {\cal B}\left({\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right),\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right)} \right) + F + {\rm \Sigma }\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right) + {\rm \Lambda }\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right).\end{equation} \tag{ 22 }$$

In Equation (22), the operators ${\cal L}$, ${\cal B}$, and F are defined according to Equations (4b), (12b), and (20), respectively, and Σ and Λare the linear operators representing the parametric divergence forcing and dissipation, respectively, given by

$$\begin{equation} {\rm \Sigma } = \left[ {\begin{array}{@{}cc@{}} {D_u \left({x,y,t} \right)\nabla ^2 } & {D_B \left({x,y,t} \right)\nabla ^2 } \\ {D_B \left({x,y,t} \right)\nabla ^2 } & {D_u \left({x,y,t} \right)\nabla ^2 } \\ \end{array}} \right]\end{equation} \tag{ 23 }$$

$$\begin{equation} {\rm \Lambda } = \left[ {\begin{array}{@{}cc@{}} {\nu \nabla ^4 } & 0 \\ 0 & {\eta \nabla ^4 } \\ \end{array}} \right].\end{equation} \tag{ 24 }$$

Consequently, the general form of the wave amplitude equations in a forced-dissipative resonant triad of model Equations (22) is

$$\begin{equation} \frac{{d{\rm \Omega }_1 }}{{dt}} = {\rm C}_{\left({{\boldsymbol k}_1,{\boldsymbol k}_2,{\boldsymbol k}_3 } \right)}^{\left({\alpha _1,\alpha _2,\alpha _3 } \right)} {\rm \Omega }_2 \overline {{\rm \Omega }_3 } + {\rm \Gamma }_1 + \mathop \sum \limits_{i = 1}^3 \left({F_{i1} - {\rm \kappa }_{i1} } \right){\rm \Omega }_i \end{equation} \tag{ 25a }$$

$$\begin{equation} \frac{{d{\rm \Omega }_2 }}{{dt}} = {\rm C}_{\left({{\boldsymbol k}_2,{\boldsymbol k}_3,{\boldsymbol k}_1 } \right)}^{\left({\alpha _2,\alpha _3 \alpha _1 } \right)} {\rm \Omega }_1 \overline {{\rm \Omega }_3 } + {\rm \Gamma }_2 + \mathop \sum \limits_{i = 1}^3 \left({F_{i2} - {\rm \kappa }_{i2} } \right){\rm \Omega }_i \end{equation} \tag{ 25b }$$

$$\begin{equation} \frac{{d{\rm \Omega }_3 }}{{dt}} = {\rm C}_{\left({{\boldsymbol k}_3,{\boldsymbol k}_1,{\boldsymbol k}_2 } \right)}^{\left({,\alpha _3,\alpha _1 \alpha _2 } \right)} {\rm \Omega }_1 \overline {{\rm \Omega }_2 } + {\rm \Gamma }_3 + \mathop \sum \limits_{i = 1}^3 \left({F_{i3} - {\rm \kappa }_{i3} } \right){\rm \Omega }_i ,\end{equation} \tag{ 25c }$$

where

$$\begin{equation} F_{ij} = \left\langle\left({k + {\rm l}} \right)^2 {\boldsymbol L}_i,\left[ {\begin{array}{@{}cc@{}} {D\left({x,y,t} \right)} & {B\left({x,y,t} \right)} \\ {B\left({x,y,t} \right)} & {D\left({x,y,t} \right)} \\ \end{array}} \right]{\boldsymbol R}_j\right\rangle \end{equation} \tag{ 26 }$$

$$\begin{equation} {\rm \kappa }_{{\rm ij}} = \left\langle\left({k + {\rm l}} \right)^4 {\rm {\boldsymbol L}}_{\rm i} \left[ {\begin{array}{@{}cc@{}} \nu & 0 \\ 0 & \eta \\ \end{array}} \right]{\boldsymbol R}_j\right\rangle.\end{equation} \tag{ 27 }$$

One notices that the linear divergence terms in both the vorticity balance and induction equations compete with the dissipative terms so that the last term in the right-hand side of Equation (25) can work in different ways for the three-wave interacting dynamics. If F_ij − κ_ij has positive eigenvalues, it will lead to instability. In this context, the pattern of the convective cells will determine the most excited modes, which in turn will be the most energetic waves. On the other hand, if F_ij − κ_ij has only negative eigenvalues, it will work as an effective damping. Finally, if the eigenvalues of F_ij − κ_ij are purely imaginary, the last term in the right-hand side of Equation (25) yields a resonant coupling between two of the three modes of the triad, contributing to the oscillatory behavior of the wave amplitudes. Consequently, the different ways by which the competition between the linear forcing and dissipation can interact with nonlinearity makes the dynamics of system (25) rather complex, exhibiting a rich variety of behaviors depending on the parameters, including transitions between regular (i.e., periodic) and chaotic regimes. For a comprehensive account of the bifurcation structure of these equations, see (Wersinger et al. 1980a, 1980b). Chian et al. (2000) analyzed a special case of system (25) in the context of plasma physics, where the combination of the last two terms in the right-hand side of Equation (25) acts as a linear parametric forcing for the unstable mode (mode 1) and as a dissipation mechanism for the remaining modes of the triad. They showed that the three-wave interacting dynamics in this scenario undergoes transitions from regular to chaotic behaviors via either the period-doubling route or the type I Pomeau–Manneville intermittency route. Here we stick to the case where the dissipative terms are slightly greater than linear divergence forcing in order to demonstrate the potentiality of these thermal convection parametric forcing terms to sustain, against dissipative losses, the amplitude oscillations related to the nonlinear coupling among the MHD Rossby modes.

In fact, the dynamical system governed by Equation (25) can be conservative provided it preserves volume in phase space, which is accomplished if Equation (25) satisfies the Liouville property:

$$\begin{equation*} \mathop \sum \limits_{i = 1}^3 \frac{{\partial {\rm \dot \Omega }_i }}{{\partial {\rm \Omega }_i }} = \mathop \sum \limits_{i = 1}^3 \left({F_i - {\rm \kappa }_i } \right) = 0.\end{equation*}$$

An unbalance between forcing and dissipation will lead the phase space volume to shrink (dissipation is greater than forcing), or to expand in the opposite situation.

Figure 4 shows an example of a numerical solution of system (25), for the same triad presented previously, representing the situation in which the forcing is greater than dissipation on the dipole scale, and dissipation is stronger on a small scale. This is reasonable because dissipation varies linearly with the wavenumber, being therefore more effective for the small-scale disturbances. We use parameters (F₁ − κ₁) = 0.1, (F₂ − κ₂) = −0.4, (F₃ − κ₃) = 0.45, discarding diagonal terms. In order to obtain this numerical solution, we have adopted the fourth-order Runge–Kuta numerical scheme. In Figure 4, we see that this combination of forcing and dissipation will alter the frequency of oscillation, where in this particular case the solution obtained seems to represent an attractor corresponding to a periodic orbit.

Zoom In Zoom Out Reset image size

5.3. Interaction with Differential Rotation and Nonuniform Zonal Magnetic Field

Apart from the convective forcing that excites MHD Rossby waves through the divergent circulation associated with the vertical motion, instabilities of either the differential rotation or the toroidal magnetic field might constitute another energy source for MHD Rossby waves in the solar tachocline. Recall that the meridionally varying background flow associated with solar differential rotation can trap the magnetic branch of MHD Rossby waves within the latitudinal belt between around ±35° and refract them toward the equator, as demonstrated earlier. However, here we shall perform an alternative analysis of the effect of a latitudinally varying background state on MHD Rossby waves in order to elucidate how instabilities associated with either solar differential rotation or toroidal magnetic fields of a solar tachocline might interact with resonant nonlinear interactions among the MHD Rossby waves discussed in this paper. For this purpose, let us analyze here how a meridionally varying background flow combined with a meridionally varying toroidal magnetic field can modify the eigenmodes of the homogeneous operator ${\cal L}$ defined by Equation (4b). Thus, for clarity in exposition, we consider the linearized and unforced (D_u = D_B = 0) version of Equations (3) for perturbations around a background state defined according to

$$\begin{equation} {\boldsymbol V} = \left[ {U\left(y \right),0} \right]\end{equation} \tag{ 28a }$$

$$\begin{equation} {\boldsymbol B} = \left[ {B_0 + B\left(y \right),0} \right],\end{equation} \tag{ 28b }$$

with B₀ being a constant. In this context, the time evolution of the infinitesimal amplitude perturbations around the background state defined above can be written in terms of stream function ψ and vector potential A as

$$\begin{equation} \frac{\partial }{{\partial t}}\left({\begin{array}{@{}c@{}} {\nabla ^2 \psi } \\ {\nabla ^2 A} \\ \end{array}} \right) = \left[ {{\cal L} + \bar {\cal L}\left(y \right)} \right]\left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right),\end{equation} \tag{ 29 }$$

where ${\cal L}$ is the homogenous linear operator given by Equation (4b) and $\bar {\cal L}\left(y \right)$ is the linear differential operator with variable coefficients, given by

$$\begin{equation} {\bar {\cal L}\left(y \right) = \left[\! {\begin{array}{@{}cc@{}} { - U\left(y \right)\displaystyle\frac{{\partial \nabla ^2 }}{{\partial x}}} &\!\!\! {B^{\prime\prime}\left(y \right) \displaystyle\frac{\partial }{{\partial x}}} \\ {\left({B^{\prime\prime}\left(y \right) \displaystyle\frac{\partial }{{\partial x}} + 2B^{\prime}\left(y \right) \displaystyle\frac{{\partial ^2 }}{{\partial x\partial y}}} \right)} &\!\!\! {\left({U\left(y \right) \displaystyle\frac{{\partial \nabla ^2 }}{{\partial x}} - U^{\prime\prime}\left(y \right) \displaystyle\frac{\partial }{{\partial x}} - 2U^{\prime}\left(y \right) \displaystyle\frac{{\partial ^2 }}{{\partial x\partial y}}} \right)} \\ \end{array}} \!\right].} \end{equation} \tag{ 30 }$$

To analyze how the nonhomogeneous operator $\bar {\cal L}\left(y \right)$ modifies the eigenmodes of the homogeneous operator ${\cal L}$, let us assume that the perturbed fields in Equation (29) are composed of a single mode:

$$\begin{equation} {\boldsymbol W} = \left({\begin{array}{@{}c@{}} \psi \\ A \\ \end{array}} \right) = {\rm \Omega }_{j} e^{i\left({{\boldsymbol k}_j.{\boldsymbol x} - \omega _j t} \right)} {\boldsymbol R}^{\alpha _j } ({{\boldsymbol k}_j } ),\end{equation} \tag{ 31 }$$

where Ω_j refers to an arbitrary constant amplitude and ${\boldsymbol R}^{\alpha _j } ({{\boldsymbol k}_j } )$ corresponds to the normalized right eigenvector of ${\cal L}$. Inserting the ansatz (31) into Equation (29) yields

$$\begin{equation} \omega _j ({{\boldsymbol k}_j } ) = {\tilde \omega} _{j } ({{\boldsymbol k}_j } ) + iC_{jj}, \end{equation} \tag{ 32a }$$

where

$$\begin{equation} C_{jj} = \int_{T} e^{ - ily} \langle\bar {\cal L}^* \left(y \right){\boldsymbol R}^{\alpha _j } ({{\boldsymbol k}_j } ),{\boldsymbol R}^{\alpha _j } ({{\boldsymbol k}_j } )_{{\mathbb C}_{^2 } }\rangle dy.\end{equation} \tag{ 32b }$$

In Equation (32b), l is the meridional wavenumber, $\langle\rangle _{{\mathbb C}_2}$ refers to the inner product defined according to Equation (9), $\bar {\cal L}^* (y )$ is the operator obtained from $\bar {\cal L}(y )$ by replacing ∂/∂x by ik and ∂/∂y by il, and ${\tilde \omega} _{j } ({{\boldsymbol k}_j } )$ corresponds to the eigenvalue of the homogeneous operator ${\cal L}$ associated with eigenvector ${\boldsymbol R}^{\alpha _j } ({{\boldsymbol k}_j } )$ and is given by Equation (7). From Equation (32) one can notice that the nonhomogeneous operator $\bar {\cal L}^* (y )$ can either modify the time frequency of MHD Rossby modes from their eigenvalues (7) or generate instabilities or exponential growth of these modes. In this context, the imaginary part of coefficient C_jj is responsible for the Doppler-shift effect on the propagation of MHD Rossby modes, whereas the real part of C_jj yields either exponential growth (Re(C_jj) < 0) or decay (Re(C_jj) > 0) of the amplitude of these modes.

Thus, according to the linear dynamics represented by Equation (29), if a particular infinitesimal amplitude mode j is unstable (Re(C_jj) < 0),its amplitude will grow exponentially in time. However, after a certain time, the amplitude growth of the unstable mode leads the linear approximation to break down, and the quadratic nonlinearities must be restored for a satisfactory description of system (3) in this regime. In this case, the role of the quadratic nonlinearities is to yield saturation for the amplitude growth of the unstable mode j. In this scenario, we believe that the nonlinear saturation is mostly due to resonant triad interactions in which mode j participates as the pump mode (the triad member with the highest absolute value of nonlinear coupling coefficient). In fact, if we restore the quadratic nonlinear and dissipative terms in Equation (29), the resulting dynamical system governing the resonant three-wave interacting dynamics in this case will be identical to system (25), but with the coefficients F_ij as given by Equation (26) being replaced by Re(C_ij) as given by Equation (32b). In this case, the modes composing the resonant triad would match the usual resonant triad interaction conditions for their wavenumbers and the real part of their time–frequencies, whereas the imaginary part of their time–frequencies would appear as linear forcing coefficients in the right-hand side of Equation (25).

A linear instability analysis of MHD shallow-water equations describing the dynamics of a solar tachocline has been investigated by Dikpati & Gilman (2005) for the profiles U(y) and B(y), representative of solar differential rotation and toroidal magnetic field bands in the solar tachocline, respectively. They show that the most unstable modes represent planetary-scale disturbances with a wavenumber 1–2 spatial structure. As demonstrated in Section 4, the hydrodynamic branch of wavenumber 1–2 MHD Rossby modes, in turn, can work as a pumpwave in resonant triad interactions with shorter-wavelength MHD Rossby modes. On the other hand, dissipation is believed to be small for planetary-scale disturbances and to have its effect enhanced as wavenumber increases. As a consequence, in the presence of a background state represented by Equation (28), with U(y) and B(y) representing the solar differential rotation and toroidal magnetic field of the solar tachocline, respectively, for resonant triad interactions involving a hydrodynamic MHD Rossby mode with zonal and meridional wavenumbers one or two and shorter spatial-scale MHD Rossby modes, we would expect that dissipation is weaker than linear forcing for the planetary-scale pump mode and stronger than linear forcing for the shorter spatial-scale modes. Hence, for this kind of resonant triad of MHD Rossby modes, the competition between linear forcing and dissipation yields an effective linear forcing for the pump mode and an effective damping for the remaining modes of the triad. As demonstrated by Chian et al. (2000) for the nonlinear triad interaction involving Langmuir, whistler, and Alfvén waves in the context of planetary magnetospheres, the nonlinear triad interaction in this case is capable of saturating the exponential growth of the unstable mode, with the solution of the three-wave interaction equations undergoing either periodic or chaotic behaviors. Thus, we believe that the resonant triad interactions described in Section 4, when embedded in a background state characterized by a differential rotation mean zonal flow and the toroidal magnetic field of the solar tachocline, might be capable of saturating the exponential growth of the unstable planetary-scale mode and yielding nonaxisymmetric periodic modulations in the solar magnetic activity. Saturation of global instabilities by nonlinear waves at the solar tachocline is also suggested by numerical simulations of the full hydrodynamic shallow-water equations performed by Dikpati (2012).

Therefore, apart from convective sources related to thermal inhomogeneities along a tachocline's boundaries, we have pointed out in the present part of this section that global instabilities of both solar differential rotation and toroidal magnetic field can also be accounted for by potentially important energy sources that prevent resonant triad interactions of MHD Rossby modes from declining through viscous and Ohmic dissipative actions. In addition, we believe that such instabilities might be the natural triggering for the resonant interactions of the MHD Rossby waves analyzed in Section 4.

6.1. Generation of Waves Due to Meridional Circulation

6.2. Meridional Flows Due to Rossby Waves

In view of the success of transport-flux dynamos, we briefly discuss a possible connection with Rossby waves. As discussed in Schmitt & Schüssler (2004), it is difficult to distinguish between a dynamo wave mechanism and an advection explanation for the structure of the butterfly diagram, a distinction that could only be made by helioseismological observation of the properties of the meridional flow at tachocline depths, but Zhao et al. (2013) suggest that the equatorial return of meridional circulation occurs at depths shallower than the tachocline. In Section 3.2 we provide a mechanism that can be considered as an MHD generalization of the dynamo waves that keep the desired properties. Here we show that Rossby waves might also help explain longitudinally dependent structures in the currently prevalent view of equatorward flux transport. Our arguments could be adapted to near-surface poleward advection of poloidal magnetic fields.

Letting$\psi = A{\rm exp} [ {i ({{\boldsymbol k}.{\boldsymbol x} - \omega t} )} ]$ be the stream function associated with a Rossby wave, the corresponding flow will be

$$\begin{equation*} {\boldsymbol u} = \left({\begin{array}{@{}c@{}} u \\ v \\ \end{array}} \right) = \nabla ^ \bot \psi = i\left({\begin{array}{@{}c@{}} { - {\rm l}} \\ {{ k}} \\ \end{array}} \right)Ae^{i\left({{\boldsymbol k}.{\boldsymbol x} - \omega t} \right)},\end{equation*}$$

where we see that for all Rossby waves with nonzero zonal wavenumber k (nonaxisymmetric) there will be a meridional flow associated with them. Supposing that the horizontal velocity field ${\boldsymbol u}_H$ is composed of a superposition of resonantly interacting Rossby waves, in a particularly interesting case one of the waves is axisymmetric, providing a zonal flow, and one is not axisymmetric, providing a meridional flow. A cubic mechanism is needed to couple such waves; see Burzlaff et al. (2008) for details. When the amplitude of the axisymmetric wave is maximal, the zonal flow (differential rotation) is maximal, and when the amplitude of the nonaxisymmetric wave is maximal, the meridional flow has maximum velocity. Therefore, resonantly coupled Rossby waves can provide a mechanism for oscillatory behavior of both zonal and meridional flows; the oscillation can be in or out of phase depending on the modes involved and on their respective amplitudes at a given time (i.e., initial conditions). This mechanism can be associated with observed variations in the meridional flow throughout a cycle (Javaraiah et al. 2006 see also Dikpati et al. 2004 and references therein for details).

Because of the continuity of the radial shear, the meridional flow of the Rossby waves must match the meridional flow of the layers above. An interaction of Rossby waves at the tachocline and the meridional flow at the convection zone will occur; this combination will affect the preexisting meridional flow by adding a longitudinal structure to it. This effect might combine with the traditional interpretation of equatorward flux transport to explain longitudinally dependent structures, such as preferred longitudes (Dikpati & Gilman 2005).

Here we have presented a possible fundamental theoretical description of the cyclic oscillations of solar magnetic activity based on resonant nonlinear interactions among MHD Rossby waves in the solar tachocline. According to our theory, planetary-scale modes having a meridional structure of a nonaxisymmetric dipole mode of the Sun's toroidal magnetic field can act as a pump wave in resonant triad interactions, supplying energy to and receiving energy from shorter-scale modes in a periodic fashion, and the energy exchanges resulting from this parametric instability yield periodic oscillations of the polarity of the pump mode with a period close to the nearly 11 yr period of the solar cycle. These periodic oscillations are in phase with the maximum energy level of the remaining shorter-scale waves of the triad. To demonstrate a closer fidelity of our theory with observations, we have also examined the linear theory of MHD Rossby waves embedded in a latitudinally varying zonal flow mimicking the solar differential rotation. In this context, we have computed the parameters associated with wave propagation and found that in the solar dynamo the equator can act as a waveguide for the magnetic branch such that waves generated in midlatitudes are refracted toward the equator. In a more realistic model that takes into account instabilities of the background state, as these waves propagate they may cause disturbances on the surrounding magnetic field, which would possibly trigger instabilities. This could account for the appearance of sunspots as the waves propagate throughout the equatorial waveguide. The path followed by the waves could explain the format of the butterfly diagram, in which there is magnetic activity appearing at around 35° and propagating toward the equator. This feature is in agreement with Joy's law for activity sunspot zones, suggesting that activity zones could be determined by Rossby waves. The modulation of this wave activity in a nearly 11 yr cycle observed in the butterfly diagram could be explained by resonant nonlinear interactions of these waves with larger-scale modes representing the dipole meridional structure of the magnetic field.

We would like to stress the fact that nearly periodic coherent phenomena are rather common in the context of turbulent systems. This was perhaps first noticed in the classic work of Fermi and collaborators on a chain of coupled oscillators (Fermi et al. 1955), the so-called Fermi–Pasta–Ulam paradox. The counterintuitive result of Fermi and collaborators is that despite the high dimensionality and nonlinearity of their problem, there is a violation of ergodicity, and the final states of the system are characterized by the energy being concentrated in a few modes that exchange energy periodically through nonlinear resonance. In our view, the situation of the solar cycle is analogous, as may become evident by examining the evolution of the amplitudes of spherical harmonic coefficients of the solar magnetic field (Knaack & Stenflo 2005). The present article is therefore a first attempt to put the periodic nature of the solar cycle in this context.

In summary, we believe that this "quasi-geostrophic MHD" picture could help in clarifying the interpretation of several aspects of the magnetic activity of the Sun, including periodicities and spatio-temporal evolution, particularly the features with longitudinal structure.

In order to develop a complete theory explaining both the butterfly structure of sunspots and its decadal modulation, we need to put the arguments presented in Section 5.3 on resonant MHD Rossby wave triads embedded in a latitudinally varying background state into a quantitative form, as well as augment this theory for a continuous spectrum of wave packets. We intend to develop this theory in future work. Further generalization of this model to spherical geometry might be important as well because, in the spherical geometry, Rossby waves can generate zonal jets at higher orders of amplitude, as suggested by Burzlaff et al. (2008), permitting exchange of energy with zero-frequency modes. A closer fidelity of our results with the observed temporal evolution of solar activity might also be achieved by augmenting our nonlinear theory to include coupled resonant triads as well as parametric forcings that could yield either exponential growth or a linear resonant coupling. These points will be investigated in a future paper.

The work of Carlos F. M. Raupp is supported by FAPESP, and Breno Raphaldini was supported by CAPES. We also thank the anonymous reviewer for his or her fruitful comments for improving this manuscript.