Nonlinear Rossby Wave–Wave and Wave–Mean Flow Theory for Long-term Solar Cycle Modulations

The solar cycle observed in sunspot number time series is approximately periodic in nature, with the main period being around 11 yr. On closer inspection, however, the solar cycle exhibits both long-term modulations and fluctuations in the 11 yr period on the same timescales. It has been suggested that several periods exist for these modulations, including 100, 220, and 1000 yr periods (Usoskin 2017).

It has long been noted that an apparent connection between the duration and magnitude of the cycle does exist (Waldmeier 1936), suggesting an inverse correlation between the amplitude of the solar activity (maximum number of sunspots at the peak of the cycle) and the duration of the cycle. Possibly, one of the most remarkable manifestations of this relationship is the period of the cycle during the historical minimum of solar activity (Maunder minimum), which was found to be increased by as much as twice the usual 11 yr period. Such a relation has been used as an empirical prediction method of the maximal activity at the peak of the cycle as a function of the increased rate of the activity during the ascending phase of the cycle (Pipin et al. 2012).

All of the characteristics mentioned above impose a remarkable challenge to dynamo models of practical importance, since the improvement of such modeling issues may lead to better predictions of the solar activity. Although the theory behind such long-term variations of the solar cycle remains elusive, dynamo models are able to reproduce some long-term features of the solar cycle by the inclusion of a stochastically fluctuating alpha effect with slow variations. The physical basis behind such stochastic fluctuations, however, remains unclear, where one of the possible suggestions is that they are related to turbulent fluctuations in vortex sizes and their turnover timescales (Pipin et al. 2012). One particular class of models that has been able to reproduce long-term variation patterns of solar activity refers to alpha-omega dynamos with stochastic parameters (Hoyng 1993). This class of dynamo models yields wavelike solutions of the induction equation that propagate toward the equator, therefore reproducing the butterfly pattern of solar activity, with the activity beginning at midlatitudes, around 35°–40° of latitude, and gradually migrating toward the equator. Long-term activity variations can also be achieved in flux-transport dynamo models with a prescribed meridional flow (Dikpati et al. 2005).

Another feature of solar activity that is believed to exhibit fluctuations on slow timescales is the differential rotation, or the zonal flow profile. Although it might be difficult to detect long-term variations of the differential rotation at the tachocline or inside the convection zone, variations in differential rotation are observed at the solar surface and might be associated with hemispheric asymmetries in sunspot activity (see Zhang et al. 2015). Such variations on the mean zonal flow are apparently in anticorrelation with amplitudes of solar activity, suggesting a possible coupling between smaller-scale processes, such as the ones that may lead to sunspot activity, and the global-scale differential rotation profile.

Raphaldini & Raupp (2015) analyzed the weakly nonlinear interaction theory of magnetohydrodynamic (MHD) Rossby waves embedded in a constant toroidal magnetic field background state and showed that the periodic fluctuations of the wave amplitudes associated with the resonant triad coupling occur on the same timescale of the solar cycle (typically around 1 order of magnitude greater than the linear period of the waves). Consequently, the authors suggested that the temporal energy modulations of such MHD Rossby waves due to nonlinear interaction might be related to the periodic nature of the solar magnetic activity.

Recently, Dikpati et al. (2018) analyzed a shallow-water MHD model for the solar tachocline and highlighted that the nonlinear interaction involving magnetic Rossby waves, the differential rotation profile, and toroidal magnetic fields might be responsible for the so-called quasiperiodic tachocline nonlinear oscillations. Here we will argue that a similar mechanism might also give rise to long-term modulations in the solar magnetic activity.

In a recent article, Bustamante et al. (2014) proposed a novel mechanism of nonlinear wave systems that might produce long-term fluctuations in the wave amplitudes in an intermediate-amplitude regime. The basis of this mechanism relies on a resonance between the nonlinear oscillations of one wave triad and the fluctuations in the wave phases of an adjacent triad. In this scenario, the phases' oscillations will be strongly influenced by the frequency mismatch among the waves. The precession resonance mechanism is shown to be responsible for strong energy transfers throughout the whole nonlinear wave system in several contexts (Bustamante et al. 2014). Also, another interesting feature of this mechanism is that it allows the mean zonal flow (represented by eigenmodes having both zero time frequency and zero zonal wavenumber) to exchange energy with the waves, which is not possible in the weak turbulence limit (small amplitude) in that the zonal flow acts as a catalyst for the energy exchanges between the waves.

Here we augment the nonlinear interaction theory of MHD Rossby waves at the solar tachocline developed by Raphaldini & Raupp (2015) by accounting for effects similar to the precession resonance that allows significant energy transfer throughout different triads, as well as the interaction between Rossby waves and the differential rotation. For this purpose, we first search for sets of three interacting waves such that the mismatch among the waves' eigenfrequencies is close to one of the harmonics of the solar cycle, that is,

$$\begin{eqnarray*}&&{\omega }_{1}+{\omega }_{2}-{\omega }_{3}\sim j\pi /22\,{\mathrm{yr}}^{-1},j=1,2,3....\end{eqnarray*}$$

The triads mentioned above contain a mode with zero zonal wavenumber and eigenfrequency that mimics the solar differential rotation effects. Then, we have analyzed a representative example of such triads in which the triad is connected via one wave mode to a second triplet. If one of these triads dominates the initial energy of the system, the initial excitation of the secondary triplet can be explained by a linearized theory through a mechanism reminiscent of modulational instability (Connaughton et al. 2010). However, we have demonstrated that when the nonlinearity associated with the secondary wave triad is restored, the maximum efficiency of intertriad energy exchange is attained in the precession resonance regime with the secondary triplet having a nonlinear frequency of amplitude modulation near jπ/22 yr⁻¹, which refers to the frequency mismatch among the waves of the primary triad. The resulting energy exchanges between the two wave triads yield modulations on a timescale longer than the 22 yr cycle, which corresponds to the period of the amplitude oscillations of the second triad.

In addition, as a consequence of the Manley–Rowe invariant (Bustamante & Kartashova 2011), the nonlinear oscillation period of the mode amplitudes of a wave triad is inversely proportional to the square root of the energy of two modes of the triad, that is,

$$\begin{eqnarray*}&&T(I)\propto 1/\sqrt{I},\end{eqnarray*}$$

where I is a weighted sum of the squares of the mode amplitudes.

The above equation provides an inverse relationship between the amplitude and nonlinear period of the triad, which is remarkably similar to the aforementioned Waldmeier law. Also, the zonal flow mode amplitude modulations are found to be approximately in opposite phase with the amplitude oscillations of the second triad, which is supposedly related to the Schwabe cycle, according to our theoretical model. Therefore, we argue here that the precession resonance involving MHD Rossby wave triads might be a possible mechanism behind the long-term modulations of the solar cycle observed in sunspot number time series.

In Section 2 we introduce the model equations, which refer to a simplified version of the quasi-geostrophic MHD equations derived by Zeitlin (2013), but augmented to take into account the effects of spherical geometry. Section 2 also revisits the linear theory of the model equations for a resting and constant toroidal magnetic field background state, as well as the reduced dynamics of a single triad of nonlinearly interacting waves. In Section 2 we also show that there is a considerable amount of triads of large-scale Rossby–Hauwirtz modes whose frequency mismatch among the waves is comparable to the typical main frequency of the solar cycle. We then consider in the following sections a representative example with a set of four modes coupled to a zonal flow mode (zero eigenfrequency and zonal wavenumber mode). The solutions of the system are analyzed in Sections 3 and 4 for the conservative and forced-dissipative cases, respectively. Section 5 discusses how the precession resonance mechanism explored here might operate in the full model equations. The main conclusions are presented in Section 6. Further details of the calculations and mechanisms explored here are presented in the appendixes. Appendix A presents details of the calculation of the nonlinear coupling coefficients. Appendix B reviews the general ideas regarding the integration of the three-wave interaction equations. Appendix C provides a further explanation of the two processes of excitation that are relevant to this study: the modulational instability and the aforementioned precession resonance mechanism. Appendix D details the evaluation of the damping coefficients.

2.1. Basic Equations and Linear Theory

As the simplest context supporting the existence of MHD Rossby wave disturbances, we consider the barotropic vorticity equation in the rotating MHD case (Gilman 2000; Zaqarashvili et al. 2007; Raphaldini & Raupp 2015). This equation can be regarded as the asymptotic limit of the quasi-geostrophic MHD equations for high equivalent depth. The quasi-geostrophic MHD equations, in turn, have been derived by Zeitlin (2013) as a distinguished limit of the shallow-water MHD equations in the regime of strong rotation. The advantage of adopting this simplified model instead of the original MHD shallow-water system is that the quasi-geostrophic equation filters out the inertia-gravity waves and, therefore, allows only the rotational modes that are relevant to our analysis. We assume that there is a global structure background magnetic field in the toroidal direction given by $\overline{B}(\theta )={B}_{0}\cos \theta$. As will be shown later, this choice simplifies our mathematical analysis by yielding linear eigenmodes with a purely oscillatory nature associated with a separable linear operator. Indeed, Equations (1) and (2) are no longer separable when one considers more realistic banded structures for the background toroidal magnetic field (Zaqarashvili et al. 2007)³ . With the considerations described above, the evolution equations of the system in spherical coordinates can therefore be written in terms of the absolute vorticity q, stream function ψ, and magnetic potential A as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial q}{\partial t}+{ \mathcal J }(\psi ,q)={ \mathcal J }(A,{{\rm{\nabla }}}^{2}A),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial A}{\partial t}+{ \mathcal J }(\psi ,A)=0,\end{eqnarray} \tag{ 2 }$$

where $q={{\rm{\nabla }}}^{2}\psi +2{\rm{\Omega }}\sin \theta$ and Ω refers to the rigid body rotation rate of the Sun. In the spherical coordinate system adopted here, the Jacobian and Laplacian operators take the form

$$\begin{eqnarray}&&{ \mathcal J }(f,g)=\displaystyle \frac{1}{{a}^{2}\cos \theta }\left(\displaystyle \frac{\partial g}{\partial \theta }\displaystyle \frac{\partial f}{\partial \phi }-\displaystyle \frac{\partial f}{\partial \theta }\displaystyle \frac{\partial g}{\partial \phi }\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}f=\displaystyle \frac{1}{{a}^{2}\cos \theta }\left[\displaystyle \frac{1}{\cos \theta }\displaystyle \frac{{\partial }^{2}f}{\partial {\phi }^{2}}+\displaystyle \frac{\partial }{\partial \theta }\left(\cos \theta \displaystyle \frac{\partial f}{\partial \theta }\right)\right]\end{eqnarray} \tag{ 4 }$$

for any differentiable functions f and g, where ϕ is the longitude, θ is the latitude, and a is the tachocline solar radius. In order to linearize the equations around a background state at rest and with a zonally symmetric toroidal magnetic field in the zonal direction, $(\overline{B}(\theta ),0)$, we set $A=\overline{A}+A^{\prime}$, ψ = ψ' and discard the terms arising from products of perturbations, where $\overline{A}$ is chosen such that the mean toroidal magnetic field assumes the special form $\overline{B}={B}_{0}\cos \theta$. With the assumptions described above, the linearized perturbation equations can be written in a vector form as

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left[\begin{array}{c}{{\rm{\nabla }}}^{2}\psi \\ A\end{array}\right]=\left[\begin{array}{cc}-\tfrac{2}{{a}^{2}}{\rm{\Omega }}\tfrac{\partial }{\partial \phi } & \tfrac{{B}_{0}}{a{\mu }_{0}\rho }\tfrac{\partial }{\partial \phi }\left(\tfrac{1}{{a}^{2}}-{{\rm{\nabla }}}^{2}\right)\\ \tfrac{{B}_{0}}{a}\tfrac{\partial }{\partial \phi } & 0\end{array}\right]\left[\begin{array}{c}\psi \\ A\end{array}\right].\end{eqnarray} \tag{ 5 }$$

In the equation above, we have omitted the primes when referring to the perturbations for simplicity. Equation (5) can be solved by spherical harmonics (see Zaqarashvili et al. 2007 for similar treatment). Briefly, we assume a plane wave ansatz,

$$\begin{eqnarray}&&\begin{array}{rcl}\left[\begin{array}{l}\psi \\ A\end{array}\right] & = & {\rm{\Lambda }}{Y}_{n}^{m}(\phi ,\theta ){e}^{-i\omega t}{\boldsymbol{R}}\\ & = & {\rm{\Lambda }}{N}_{n}^{m}{P}_{n}^{m}(\sin \theta ){e}^{{im}\phi -i\omega t}{\boldsymbol{R}},\end{array}\end{eqnarray} \tag{ 6 }$$

with Λ being an arbitrary constant and the associated Legendre functions satisfying the following orthogonality relation:

$$\begin{eqnarray}&&{\int }_{-1}^{1}{P}_{{n}_{1}}^{m}{P}_{{n}_{2}}^{m}{dz}=2\displaystyle \frac{(n+m)!}{(n-m)!(2n+1)}{\delta }_{{n}_{1}{n}_{2}},\end{eqnarray} \tag{ 7 }$$

with ${\delta }_{{n}_{1}{n}_{2}}=1$ if n₁ = n₂ and zero otherwise. The normalization constant N_n^m is given by

$$\begin{eqnarray*}&&{N}_{n}^{m}={\left(\displaystyle \frac{(n-| m| )!(2n+1)}{(n+| m| )!}\right)}^{\tfrac{1}{2}}.\end{eqnarray*}$$

Inserting the ansatz (Equation (6)) into the linearized Equation (5) yields the following eigenvalue problem:

$$\begin{eqnarray}&&{{ \mathcal L }}^{* }{\boldsymbol{R}}=i\omega {\boldsymbol{R}},\end{eqnarray} \tag{ 8 }$$

where the matrix ${{ \mathcal L }}^{* }$ refers to the symbol of the linear operator of Equation (5), that is,

$$\begin{eqnarray}{{ \mathcal L }}^{* }=\left[\begin{array}{cc}-2{\rm{\Omega }}\tfrac{m}{n(n+1)} & \tfrac{{B}_{0}m}{{\mu }_{0}\rho }\left(\tfrac{1}{{a}^{2}}-\tfrac{2}{n(n+1)}\right)\\ \tfrac{{B}_{0}m}{a} & 0\end{array}\right].\end{eqnarray} \tag{ 9 }$$

The two branches of the characteristic equation of Equation (8) are defined according to

$$\begin{eqnarray}&&\begin{array}{l}{\omega }_{-}(m,n)=\displaystyle \frac{1}{2}\left(-\displaystyle \frac{2{\rm{\Omega }}m}{n(n+1)}\right.\\ \quad -\left.\,\sqrt{{\left(\displaystyle \frac{2{\rm{\Omega }}m}{n(n+1)}\right)}^{2}+\displaystyle \frac{4{B}_{0}^{2}{m}^{2}}{{\mu }_{0}\rho {a}^{2}}\left(1-\displaystyle \frac{2}{n(n+1)}\right)}\right),\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\omega }_{+}(m,n)=\displaystyle \frac{1}{2}\left(-\displaystyle \frac{2{\rm{\Omega }}m}{n(n+1)}\right.\\ \quad +\left.\,\sqrt{{\left(\displaystyle \frac{2{\rm{\Omega }}m}{n(n+1)}\right)}^{2}+\displaystyle \frac{4{B}_{0}^{2}{m}^{2}}{{\mu }_{0}\rho {a}^{2}}\left(1-\displaystyle \frac{2}{n(n+1)}\right)}\right).\end{array}\end{eqnarray} \tag{ 11 }$$

The branch ω₋ refers to the fast hydrodynamic mode, while ω₊ represents the slow magnetic mode (Zaqarashvili et al. 2007). In fact, ω₋ reduces to the classical hydrodynamic Rossby–Hauwirtz wave dispersion relation for B₀ = 0. Note also that for n = 1, the magnetic effects cancel out, and there is only one branch corresponding to a Rossby–Hauwirtz wave mode. This agrees with the dispersion relation obtained in Zaqarashvili et al. (2007). The corresponding right eigenvectors ${\boldsymbol{R}}$ are given by

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{\pm }(m,n)=\left[\begin{array}{c}\displaystyle \frac{{\omega }_{\pm }(m,n)}{m}\\ \displaystyle \frac{{B}_{0}}{a}\end{array}\right].\end{eqnarray} \tag{ 12 }$$

2.2. Nonlinear Theory of Wave Interactions

Restoring the nonlinear terms in the perturbation approach described in the previous subsection, Equation (5) now reads

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left[\begin{array}{c}{{\rm{\nabla }}}^{2}\psi \\ A\end{array}\right]={ \mathcal L }\left[\begin{array}{c}\psi \\ A\end{array}\right]={ \mathcal B }\left(\left[\begin{array}{c}\psi \\ A\end{array}\right],\left[\begin{array}{c}\psi \\ A\end{array}\right]\right).\end{eqnarray} \tag{ 13 }$$

In Equation (13), the linear operator ${ \mathcal L }$ is the same as in Equation (5), and the nonlinear (bilinear) operator ${ \mathcal B }$ is given by

$$\begin{eqnarray}&&{ \mathcal 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}-{ \mathcal J }(\psi ,{{\rm{\nabla }}}^{2}\psi )+\tfrac{1}{{\mu }_{0}\rho }{ \mathcal J }(A,{{\rm{\nabla }}}^{2}A)\\ { \mathcal J }(\psi ,A)\end{array}\right].\end{eqnarray} \tag{ 14 }$$

We now consider the following ansatz in the form of a linear combination of a few linear eigenmodes:

$$\begin{eqnarray}&&\begin{array}{rcl}\left[\begin{array}{c}\psi \\ A\end{array}\right] & = & \displaystyle \sum _{k=1}^{K}{{\rm{\Lambda }}}_{k}(t){Y}_{{n}_{k}}^{{m}_{k}}(\phi ,\theta ){{\boldsymbol{R}}}^{(k)}\\ & = & \displaystyle \sum _{k=1}^{K}{{\rm{\Lambda }}}_{k}(t){N}_{{n}_{k}}^{{m}_{k}}{P}_{{n}_{k}}^{{m}_{k}}(\sin \theta ){e}^{{{im}}_{k}\phi }{{\boldsymbol{R}}}^{(k)},\end{array}\end{eqnarray} \tag{ 15 }$$

where Λ_k(t) now represents the complex valued mode amplitudes. Let us first consider the case of a single-wave triplet (K = 3). In this case, we insert the ansatz given by Equation (15) into Equation (13) and proceed to obtain the time evolution equations for the mode amplitudes Λ_k, k = 1, 2, 3. In order to do so, we make use of the orthogonality relation (Equation (7)), the orthogonality of the ${e}^{{{im}}_{k}\phi }$ components in the [0, 2π] interval for different k, and the orthogonality of the eigenvectors ${{\boldsymbol{R}}}^{(k)}$ regarding the inner product,

$$\begin{eqnarray}&&\langle {\boldsymbol{X}},{\boldsymbol{Y}}\rangle ={X}_{1}^{* }{Y}_{1}+\displaystyle \frac{1}{{\mu }_{0}\rho }{X}_{2}^{* }{Y}_{2},\end{eqnarray} \tag{ 16 }$$

where ${\boldsymbol{X}}=({X}_{1},{X}_{2})$ and ${\boldsymbol{Y}}=({Y}_{1},{Y}_{2})$ represent two arbitrary elements of the ${{\mathbb{C}}}^{2}$ vector space, and the superscript * denotes complex conjugation. Therefore, the evolution equations for the complex valued mode amplitudes are the so-called triad equations:

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{1}}{{dt}}=i{\omega }_{1}{{\rm{\Lambda }}}_{1}+{C}_{\mathrm{1,2,3}}{{\rm{\Lambda }}}_{2}^{* }{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{2}}{{dt}}=i{\omega }_{2}{{\rm{\Lambda }}}_{2}+{C}_{\mathrm{2,3,1}}{{\rm{\Lambda }}}_{1}^{* }{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{3}}{{dt}}=i{\omega }_{3}{{\rm{\Lambda }}}_{3}+{C}_{\mathrm{3,1,2}}{{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{2}.\end{eqnarray} \tag{ 19 }$$

In the equations above, C_1,2,3, C_2,1,3, and C_3,1,2 are the interaction coefficients among the mode components of the triad, given by

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\mathrm{1,2,3}} & = & \displaystyle \frac{-i}{2a[{n}_{1}({n}_{1}+1)]| {{\boldsymbol{R}}}_{1}{| }^{2}}{K}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}}\\ & & \times \left[{I}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}}+{L}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}}\right],\end{array}\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&{I}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}}={n}_{2}({n}_{2}+1)-{n}_{3}({n}_{3}+1),\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{L}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}}={\left(\displaystyle \frac{{V}_{{\rm{A}}}}{a}\right)}^{2}\left(\displaystyle \frac{{\omega }_{2}}{{m}_{2}}-\displaystyle \frac{{\omega }_{3}}{{m}_{3}}\right),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}\begin{array}{rcl}{K}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}} & = & {N}_{{n}_{1}}^{{m}_{1}}{N}_{{n}_{2}}^{{m}_{2}}{N}_{{n}_{3}}^{{m}_{3}}{\displaystyle \int }_{-1}^{1}{P}_{{n}_{1}}^{m1}(z)\\ & & \times \left({m}_{2}{P}_{{n}_{2}}^{{m}_{2}}(z)\displaystyle \frac{{{dP}}_{{n}_{3}}^{{m}_{3}}(z)}{{dz}}-{m}_{3}{P}_{{n}_{3}}^{{m}_{3}}(z)\displaystyle \frac{{{dP}}_{{n}_{2}}^{{m}_{2}}(z)}{{dz}}\right){dz},\end{array}\end{eqnarray} \tag{ 23 }$$

with ${V}_{{\rm{A}}}=\tfrac{{B}_{0}}{\sqrt{{\mu }_{0}\rho }}$ indicating the Alfvén wave speed and z = sin θ. A detailed derivation of the interaction coefficients C_1,2,3, C_2,1,3, and C_3,1,2 is presented in Appendix A. These coefficients are nonzero provided the mode indexes (m₁, n₁), (m₂, n₂), (m₃, n₃) satisfy the following selection rules:

$$\begin{eqnarray}&&\begin{array}{l}{m}_{2}+{m}_{1}={m}_{3}\\ {({m}_{1})}^{2}+{({m}_{3})}^{2}\ne 0\\ {n}_{2}{n}_{1}{n}_{3}\ne 0\\ {n}_{2}+{n}_{1}+{n}_{3}{\rm{}}\,{\rm{is}}\,{\rm{odd}}\\ ({n}_{2}^{2}-| {m}_{2}{| }^{2})+({n}_{1}^{2}-| {m}_{1}{| }^{2})\gt 0\\ | {n}_{2}-{n}_{1}| \lt {n}_{3}\lt {n}_{2}+{n}_{1}\\ ({m}_{2},{n}_{2})\ne ({m}_{3},{n}_{3}),({m}_{1},{n}_{1})\ne ({m}_{3},{n}_{3}).\end{array}\end{eqnarray} \tag{ 24 }$$

The coupling coefficients given by Equation (20) can also be explicitly calculated in terms of Wigner 3j symbols (see Jones 1985 for details). To focus on the nonlinear terms only, the triad Equations (17–19) can be rewritten using the change of variables ${B}_{k}={{\rm{\Lambda }}}_{k}(t)\exp (-i{\omega }_{k}t)$, k = 1, 2, 3, resulting in

$$\begin{eqnarray}&&\displaystyle \frac{{{dB}}_{1}}{{dt}}={C}_{\mathrm{1,2,3}}{B}_{2}^{* }{B}_{3}{e}^{i{\rm{\Delta }}\omega t},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dB}}_{2}}{{dt}}={C}_{\mathrm{2,1,3}}{B}_{1}^{* }{B}_{3}{e}^{i{\rm{\Delta }}\omega t},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dB}}_{3}}{{dt}}={C}_{\mathrm{3,1,2}}{B}_{1}{B}_{2}{e}^{-i{\rm{\Delta }}\omega t},\end{eqnarray} \tag{ 27 }$$

where Δω = ω₃ − ω₂ − ω₁ is the mismatch among the triad eigenfrequencies. When Δω = 0, the triad is said to be resonant. In the weakly nonlinear regime, in which the wave amplitudes are assumed to be small, the contribution of nonresonant triads for the nonlinear evolution of the system is usually neglected. The justification for this approach is essentially based on the highly truncated three-wave problem dynamics described by Equations (25)–(27). In fact, in the weakly nonlinear regime, the mode amplitudes evolve in a timescale longer than that associated with the linear wave phases. Consequently, the factor ${e}^{\pm i{\rm{\Delta }}\omega t}$, in general, makes the right-hand side of Equations (25)–(27) highly oscillatory in time, so that the average contribution of the nonlinearity for the time evolution of the wave amplitudes is rather small. The exception occurs when the triad is resonant, or nearly so; in this case, the nonlinearity yields significant energy exchanges among the triad components. The predominance of resonant triads in the weakly nonlinear regime can also be demonstrated by near-identity transformation (Zakharov et al. 2012). Raphaldini & Raupp (2015) analyzed the solutions of Equations (25)–(27) for resonant triads of MHD Rossby modes. They showed that in the amplitude regime in which all three modes undergo significant energy modulations, the nonlinear amplitude oscillations have a period compatible with the Schwabe cycle. Consequently, they argued that these temporal energy modulations of MHD Rossby waves due to resonant triad interaction might be related to the periodic nature of the solar magnetic activity.

Nevertheless, for a system of several connected wave triplets, Bustamante et al. (2014) demonstrated that even if the linear frequency mismatch Δω is large, strong energy transfer can occur between different wave triads, provided the linear frequency mismatch of a wave triad resonates with the characteristic nonlinear frequency of the system (for example, the frequency of energy exchange associated with an adjacent wave triad). This novel mechanism of nonlinear wave systems is called precession resonance and occurs in an energy level that is not sufficiently small to neglect nonresonant wave triplets but still small enough to make the linear frequency mismatch Δω to dominate the time evolution of the combination of the phases of the complex valued mode amplitudes Λ_j(t), j = 1, 2, 3. Consequently, it was shown that this mechanism is able to promote events of strong energy exchanges between waves even if they are nonresonant. A more detailed description of the precession resonance mechanism for a system of two connected wave triplets is presented in Appendix C.

Therefore, in order to investigate the potential role of the precession resonance mechanism in promoting significant energy transfer between different wave triads in our MHD Rossby wave context, as well as to analyze its potential role in yielding long-term modulations in the Schwabe cycle of solar magnetic activity, we have sought triads whose linear frequency mismatch among the modes is close to one of the harmonics of the main solar cycle frequency, that is, Δω ≈ jπ/22 yr⁻¹, j = 1, 2, 3 ... . Another condition imposed for the intended triads is that they contain a mode with zero zonal wavenumber (and, consequently, zero eigenfrequency) that mimics the zonal flow effects of solar differential rotation. This condition is based on the observational results of Zhang et al. (2015), who showed that modulations of the Schwabe cycle exhibit a significant negative correlation with observed variations of the solar differential rotation strength. The number of such wave triplets as a function of the Alfvén wave speed is displayed in Figure 1 for resonances of 2:1 (Δω ≈ 4π/22 yr⁻¹) and 4:1 (Δω ≈ 8π/22 yr⁻¹).

Zoom In Zoom Out Reset image size

In the next section, we shall analyze the reduced dynamics in which a representative example of the triads mentioned above is connected with a nearly resonant triad whose nonlinear interaction period is close to the Schwabe cycle period, which in turn is one order of magnitude longer than the linear period of the waves. As we will demonstrate, strong energy transfer between the two wave triads occurs due to a modulational type of instability and the precession resonance mechanism mentioned above, which enhances the efficiency of energy transfer between two adjacent wave triplets in the unstable regime. This intertriad energy exchange yields significantly longer period modulations on the intratriad energy exchanges that exhibit a period compatible with long time modulations of the Schwabe cycle.

Let us consider now the ansatz (Equation (15)) for K = 5, in which modes 1, 2, and 3 satisfy the conditions of Equation (24), and modes 3, 4, and 5 satisfy the same selection rules, apart from the resonance condition for their eigenfrequencies (ω₃ ≈ ω₄ + ω₅). In this way, substituting Equation (15) into the nonlinear perturbation Equation (13) yields

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{1}}{{dt}}=i{\omega }_{1}{{\rm{\Lambda }}}_{1}+{C}_{\mathrm{1,2,3}}{{\rm{\Lambda }}}_{2}^{* }{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{2}}{{dt}}=i{\omega }_{2}{{\rm{\Lambda }}}_{2}+{C}_{\mathrm{2,3,1}}{{\rm{\Lambda }}}_{1}^{* }{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{3}}{{dt}}=i{\omega }_{3}{{\rm{\Lambda }}}_{3}+{C}_{\mathrm{3,1,2}}{{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{2}+{C}_{\mathrm{3,4,5}}{{\rm{\Lambda }}}_{4}^{* }{{\rm{\Lambda }}}_{5},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{4}}{{dt}}=i{\omega }_{4}{{\rm{\Lambda }}}_{4}+{C}_{\mathrm{4,3,5}}{{\rm{\Lambda }}}_{5}{{\rm{\Lambda }}}_{3}^{* },\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{5}}{{dt}}=i{\omega }_{5}{{\rm{\Lambda }}}_{5}+{C}_{\mathrm{5,3,4}}{{\rm{\Lambda }}}_{4}{{\rm{\Lambda }}}_{3}.\end{eqnarray} \tag{ 32 }$$

Equations (28)–(32) have been integrated numerically by using an explicit eighth-order Runge–Kutta discretization. A representative example of two wave triplets constituting the five-wave model (Equations (28)–(32)) is displayed in Table 1. The triplets are labeled by "a" (modes 1, 2, and 3) and "b" (modes 3, 4, and 5) and composed of magnetic-branch (slow) modes having spherical harmonics (0, 2), (1, 10), (1, 9) and (1, 9), (1, 12), (2, 10), respectively. Recall that triad a has a frequency mismatch of around 5.5 yr, whereas triad b is a nearly resonant wave triplet with a mismatch among the mode eigenfrequencies of the order of 15 yr.

Figure 2 displays the result of a numerical integration of the five-wave model (Equations (28)–(32)) for the representative example illustrated in Table 1. In this integration, the initial mode amplitudes of triad b are set to match the precession resonance regime, in which the characteristic frequency of amplitude modulation of this triad exhibits a 2:1 resonance with the eigenfrequency mismatch of triad a. From the time evolution of the mode energies presented in Figure 2, one can notice the shortest period of energy oscillation of the order of 10 yr, which is associated with the energy exchanges within each triad. This refers to half of the period of the corresponding mode amplitude oscillation. Apart from the main ≈10 yr period associated with intratriad coupling, the time evolution of the mode energies displayed in Figure 2 exhibits the alternation of small and large peaks, indicating that the magnitude of the ≈10 yr cycle is modulated on a longer timescale. Comparing the time evolution of the energies of modes (0, 2) and (1, 12), which pertain to different wave triads, shows that their large and small energy peaks are approximately in opposite phase with each other. This points out that the longer timescale modulation of the main ≈10 yr cycle is a result of the intertriad energy transfers. Therefore, Figure 2 shows that, in the precession resonance regime, a strong energy transfer between the adjacent wave triads takes place, yielding long-term modulations of the main ≈10 yr cycle associated with intratriad energy exchanges.

Zoom In Zoom Out Reset image size

There are two processes of excitation that are relevant to this study. The first process is based on instability and considers small initial amplitudes for modes 1 and 2, pertaining to triad a. The process is reminiscent of modulational instability (Connaughton et al. 2010). By linearizing the system around a quasiperiodic solution of the isolated triad b (modes 3, 4, and 5), we calculate the maximal Lyapunov exponent of the full set of equations, yielding a growth rate of 0.29 yr⁻¹ (see Appendix B). The second process is fully nonlinear and based on precession resonance (Bustamante et al. 2014). There, as the system explores several amplitude levels due to forcing and dissipation, at certain mode amplitudes, a balance can be struck between a linear combination of the frequency mismatches of the two triads and the nonlinear frequency broadening stemming from the finiteness of the amplitudes. At such amplitudes, a low-frequency oscillation is generated that can lead to strong energy transfers across modes. When some of the triads are quasi-resonant, the amplitudes at which precession resonance occurs can be quite small and therefore attainable in real situations. In Appendix C we briefly study this mechanism for Equations (28)–(32), yielding energy transfer efficiencies of up to 34%.

In order to better quantify the periods involved in the time evolution of the mode energies presented in Figure 2, Figure 3 shows the power spectrum referred to the time series corresponding to the energy of mode (1, 9) (mode 3). The spectrum has been calculated by using Welch's method (Welch 1936). From Figure 3, one observes a main peak at around the 10 yr period and a secondary peak at around 120 yr, which is reasonably compatible with the Gleissberg cycle (Usoskin 2017).

Zoom In Zoom Out Reset image size

Recall that in the case of a single-wave triplet, whose dynamics is described by the three-wave Equations (25)–(27), the time evolution of the wave amplitudes (energies) is exactly periodic in time, with the solutions being described in terms of Jacobi elliptic functions, as described in Appendix B (see also Raphaldini & Raupp 2015 and references therein). In addition, as a consequence of the Manley–Rowe invariants, the nonlinear oscillation period of the mode amplitudes of a wave triad (triad b, for instance) is inversely proportional to the square root of the energy of two modes of the triad, that is,

$$\begin{eqnarray*}&&T(I)\propto 1/\sqrt{I},\end{eqnarray*}$$

where $I=| {{\rm{\Lambda }}}_{3}{| }^{2}+| {{\rm{\Lambda }}}_{4}{| }^{2}$ is the Manley–Rowe invariant (see Bustamante & Kartashova 2011 for details). However, in the five-wave problem described by Equations (28)–(32), the quantity I of a wave triad becomes variable in time due to its coupling to the adjacent wave triplet, and so does the period of triad amplitude oscillation T(I). Consequently, in the five-wave model, the dynamics of the characteristic interaction period T of a wave triad is decreased (increased) during the periods of large (small) energy peaks of the correspondent wave triad. As the characteristic interaction period T of triad b is compatible with the timescale of the Schwabe cycle of solar magnetic activity, this relationship between the interaction period of a wave triad and the correspondent triad energy level is consistent with the well-known Waldmeier law for the solar cycle. To verify this relation, we have computed the instantaneous frequency of the spectral amplitudes Λ_k(t) by applying the Hilbert transform. Figure 4 shows the time evolution of the instantaneous frequency of the amplitude of mode (1, 9), which shows that the frequency increases (decreases) during periods of high (low) amplitudes of the 10 yr cycle.

Zoom In Zoom Out Reset image size

As waves in the solar tachocline are subjected to the action of forcing and dissipation, here we investigate how these effects can modify the five-wave dynamics in the precession resonance regime. As argued in Raphaldini & Raupp (2015), the forcing acting on barotropic Rossby waves comes from the horizontal divergence of both two-dimensional velocity and magnetic fields. This approach of considering a prescribed zero-mean horizontal divergence field as a Rossby wave source is usual in studies of Rossby waves in Earth's atmosphere (Hoskins & Karoly 1981 and references therein). In the solar tachocline, the horizontal divergence of the velocity field stems from different physical processes, such as baroclinic instability (Gilman & Dikpati 2014), gravity waves, and nonhomogeneous thermal forcings at the top of the radiative zone. In this context, the generalization of the nonlinear perturbation Equation (13) for the forced-dissipative case is given by

$$\begin{eqnarray}&&\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}\left[\begin{array}{c}{{\rm{\nabla }}}^{2}\psi \\ A\end{array}\right] & = & { \mathcal L }\left[\begin{array}{c}\psi \\ A\end{array}\right]\\ & = & { \mathcal B }\left(\left[\begin{array}{c}\psi \\ A\end{array}\right],\left[\begin{array}{c}\psi \\ A\end{array}\right]\right)+{\boldsymbol{F}}+{ \mathcal D }\left(\left[\begin{array}{c}\psi \\ A\end{array}\right]\right),\end{array}\end{eqnarray} \tag{ 33 }$$

where the prescribed forcing vector ${\boldsymbol{F}}$ and damping operator ${ \mathcal D }$ are

$$\begin{eqnarray}&&{\boldsymbol{F}}=\left[\begin{array}{c}-2{\rm{\Omega }}\sin \theta {D}_{u}(\phi ,\theta ,t)\\ 0\end{array}\right],\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}{ \mathcal D }=\left[\begin{array}{cc}\nu {{\rm{\nabla }}}^{2}({{\rm{\nabla }}}^{2}) & 0\\ 0 & \eta {{\rm{\nabla }}}^{2}({{\rm{\nabla }}}^{2})\end{array}\right],\end{eqnarray} \tag{ 35 }$$

with ${D}_{u}(\phi ,\theta ,t)$ indicating the prescribed horizontal divergence of the velocity field. We have omitted the nonlinear terms involving D_u, as well as the nonlinear terms involving the divergence of the two-dimensional magnetic field on the right-hand side of the magnetic potential equation, assuming them to be small. In Equation (35), the parameters ν and η are the coefficients of viscous and magnetic diffusivity, respectively. Here we have utilized the values of ν = 2.7 × 10 and η = 4.1 × 10² cm² s⁻¹, as suggested in Gough (2007). We have further assumed that the horizontal divergence field (actually, 2Ω sin θD_u) has the same spatial structure as mode 3, which is the mode that couples the two triads, and the time dependence of the forcing resonates with this mode. The effect of these assumptions is to yield a constant forcing coefficient only on the right-hand side of the time evolution equation of mode 3 amplitude.

With the above considerations, the resulting generalization of the five-wave Equations (28)–(32) with the inclusion of the forcing and damping is

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{1}}{{dt}}=i{\omega }_{1}{{\rm{\Lambda }}}_{1}+{C}_{\mathrm{1,2,3}}{{\rm{\Lambda }}}_{2}^{* }{{\rm{\Lambda }}}_{3}-{d}_{1}{{\rm{\Lambda }}}_{1},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{2}}{{dt}}=i{\omega }_{2}{{\rm{\Lambda }}}_{2}+{C}_{\mathrm{2,3,1}}{{\rm{\Lambda }}}_{1}^{* }{{\rm{\Lambda }}}_{3}-{d}_{2}{{\rm{\Lambda }}}_{2},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{3}}{{dt}}=i{\omega }_{3}{{\rm{\Lambda }}}_{3}+{f}_{3}+{C}_{\mathrm{3,1,2}}{{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{2}+{C}_{\mathrm{3,4,5}}{{\rm{\Lambda }}}_{4}^{* }{{\rm{\Lambda }}}_{5}-{d}_{3}{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{4}}{{dt}}=i{\omega }_{4}{{\rm{\Lambda }}}_{4}+{C}_{\mathrm{4,3,5}}{{\rm{\Lambda }}}_{5}{{\rm{\Lambda }}}_{3}^{* }-{d}_{4}{{\rm{\Lambda }}}_{4},\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{5}}{{dt}}=i{\omega }_{5}{{\rm{\Lambda }}}_{5}+{C}_{\mathrm{5,3,4}}{{\rm{\Lambda }}}_{4}{{\rm{\Lambda }}}_{3}-{d}_{5}{{\rm{\Lambda }}}_{5},\end{eqnarray} \tag{ 40 }$$

where

$$\begin{eqnarray}{d}_{i}=\left\langle {{\boldsymbol{R}}}^{(i)},\left[\begin{array}{cc}\nu {\left({n}_{i}({n}_{i}+1)\right)}^{2} & 0\\ 0 & \eta {\left({n}_{i}({n}_{i}+1)\right)}^{2}\end{array}\right]{{\boldsymbol{R}}}^{(i)}\right\rangle \end{eqnarray} \tag{ 41 }$$

for i = 1, 2, 3, 4, and 5, and the coefficient f₃ is a constant that depends on the magnitude of the divergence forcing. A more thorough derivation of the damping coefficients is presented in Appendix D.

The results of the numerical integration of Equations (36)–(40) are presented in Figures 5 and 6. Figure 5 displays the time evolution of the mode energies during a 1000 yr period, whereas Figure 6 shows the time evolution of mode 3 energy correspondent to the same numerical solution as Figure 5 but for a longer period of integration (5000 yr). As in the conservative case, from Figure 5, one notices oscillations in the mode energies on a decadal timescale superposed by modulations of the energy peaks on a timescale of centuries. Again, to better quantify the main oscillation periods involved in the time evolution of the mode energies, Figure 7 shows the power spectrum computed from the mode 3 energy time series (i.e., the mode with spherical harmonic (1, 9)). In comparison with the spectrum obtained in the conservative system, one can notice in the forced-damped case a broader main spectral peak, which is also slightly shifted to a period of 7–9 yr. Also, apart from this broad spectral peak band around 8 yr, Figure 7 shows a spectral peak corresponding to long-term modulations with period around 230 yr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One remarkable feature of the five-wave model with forcing and dissipation is that the evolution of the mode energies presents some periods of suppressed activity that resemble the Maunder minimum. This fact can be more clearly illustrated in the longer time integration presented in Figure 6. One can observe in Figure 6 the appearance of several periods with very low activity lasting several decades. Other integrations with different values of the forcing parameter f₃ show that the duration of such periods is highly dependent on the magnitude of the divergence forcing (figures not shown).

So far we have provided a simplified theoretical description with a reduced model of only five wave modes to illustrate how the precession resonance mechanism associated with MHD Rossby–Haurwitz triads might operate to generate long-term modulations of the solar cycle. The reduced models presented in the previous sections reproduce several aspects of the solar magnetic activity, and the qualitative features of the results do not depend on the particular modes chosen but rather on the type of triads that they form (the size of the mismatch and the relative size of the interaction coefficients that will determine the instability properties of the system). It is important, however, to explicate the relevance of the precession resonance mechanism in the full set of partial differential equations (PDEs) that govern the dynamics of MHD Rossby waves at the solar tachocline. The key point elucidated with the reduced five-wave system is that the precession resonance enhances the efficiency of the spectral broadening of the energy injected into the system. In this context, in the reduced model with only two interacting wave triplets, the concept of energy spectral broadening is restricted to the energy transfer between the two interacting triads, and, as the intertriad energy exchanges occur on a longer timescale than the interaction involving the modes of the primary triad, this spectral broadening leads to long-term modulations of the primary cycle. However, despite the rather oversimplified setting of only five-wave modes, the key point elucidated by the reduced model still applies to the whole spectral representation of the governing equations; that is, the reduced five-wave model (Equations (28)–(32)) should be seen as a building block for a complete understanding of the full system (Equation (13)). Indeed, as can be inferred from the results of Bustamante et al. (2014) in the full PDE governing the dynamics of Rossby/drift waves, the amplitude (energy) regime at which the interaction between different wave triads becomes the most efficient one is a regime of intermediate amplitudes, in which the concept of "intermediate" acquires a precise meaning: a regime in which there is a balance between the characteristic nonlinear timescale of the whole system and the timescale associated with the mismatch among the linear mode frequencies. Once the system is found in this regime, a synchronization of triad phases across several scales will take place, and, depending on how the energy is injected into the system by forcing mechanisms (like the one considered in Section 4), this forcing will select not only one triad but a group of interacting wave triplets, thus allowing efficient energy transfer paths along the whole spectral network. In this context, the cyclic nature of the full model equations comes from the fact that a turbulent system in the precession resonance regime operates in the vicinity of an unstable periodic orbit (UPO; Bustamante et al. 2014). The corresponding period is then set when the system chooses a particular UPO, which will depend on the energy amount injected into the system: the higher the energy level associated with the forcing injection, the faster the cycle referred to the UPO, and vice versa. Thus, one may define a precession resonance scaling between a characteristic nonlinear frequency, Ω_T, and a typical value associated with the mismatch among the linear eigenfrequencies, δω_T, to divide a turbulent system into three categories:

$$\begin{eqnarray}&&\delta {\omega }_{T}\gt {{\rm{\Omega }}}_{T}({\rm{w}}{\rm{e}}{\rm{a}}{\rm{k}}\,{\rm{t}}{\rm{u}}{\rm{r}}{\rm{b}}{\rm{u}}{\rm{l}}{\rm{e}}{\rm{n}}{\rm{c}}{\rm{e}}\,{\rm{r}}{\rm{e}}{\rm{g}}{\rm{i}}{\rm{m}}{\rm{e}}),\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&\delta {\omega }_{T}\approx {{\rm{\Omega }}}_{T}({\rm{i}}{\rm{n}}{\rm{t}}{\rm{e}}{\rm{r}}{\rm{m}}{\rm{e}}{\rm{d}}{\rm{i}}{\rm{a}}{\rm{t}}{\rm{e}}\,{\rm{t}}{\rm{u}}{\rm{r}}{\rm{b}}{\rm{u}}{\rm{l}}{\rm{e}}{\rm{n}}{\rm{c}}{\rm{e}}\,{\rm{r}}{\rm{e}}{\rm{g}}{\rm{i}}{\rm{m}}{\rm{e}}),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\delta {\omega }_{T}\lt {{\rm{\Omega }}}_{T}({\rm{s}}{\rm{t}}{\rm{r}}{\rm{o}}{\rm{n}}{\rm{g}}\,{\rm{t}}{\rm{u}}{\rm{r}}{\rm{b}}{\rm{u}}{\rm{l}}{\rm{e}}{\rm{n}}{\rm{c}}{\rm{e}}\,{\rm{r}}{\rm{e}}{\rm{g}}{\rm{i}}{\rm{m}}{\rm{e}}),\end{eqnarray} \tag{ 44 }$$

where the balance associated with the intermediate turbulence regime (Equation (43)) is closely related to the concept of critical balance in turbulence (Goldreich & Sridhar 1995; Nazarenko & Schekochihin 2011).

Another important issue that arises from our analysis is about how the mechanism presented here would apply to the magnetic cycle of other stars. First, it is important to mention that observations of magnetic activity in other stars are much more difficult, with a very limited time span of observations (e.g., see Baliunas et al. 1995 and Saar & Brandenburg 1999 for further discussions). Conversely, there seems to be a consensus that, at least for Sun-like stars, younger and more rapidly rotating stars tend to exhibit shorter cycles of magnetic activity, with some of them sometimes even exhibiting irregular cycles and strong activities. On the other hand, older and more slowly rotating stars, such as the Sun, tend to have longer cycles with weaker activity (Lorenzo-Oliveira et al. 2018; Radick et al. 2018).

Nevertheless, there are several factors that need to be taken into account when comparing our arguments with these basic trends in the state of the art of stellar magnetic activity. For instance, the intensity of both the toroidal magnetic field and the rotation of the star may alter the linear wave frequencies, as well as the typical MHD Rossby wave amplitudes that will set the corresponding characteristic nonlinear frequency of the system. Despite these difficulties, the overall tendency of stronger activity being associated with a shorter cycle of star magnetic activity seems to be compatible with both the behavior of the nonlinear frequencies in the precession resonance regime and the Waldmeier law. Finally, the other possibility is that the magnetic activity of a star operates in a regime of even higher amplitudes, so that the nonlinear timescale is faster than the typical linear wave frequency mismatch of the system. In this regime, the mechanism presented here should not be relevant, and the behavior of the corresponding star magnetic activity would probably appear to be more random, without a clear cycle, which may be the case for some young stars. Alternatively, a star operating in the weakly turbulent regime would probably have much lower levels of activity characterized by longer cycles.

Here we have augmented the nonlinear interaction theory of MHD Rossby waves in the solar tachocline developed by Raphaldini & Raupp (2015) to take into account the effect of the precession resonance mechanism that allows significant energy transfer throughout different wave triads, as well as the interaction between Rossby waves and modes having zero zonal wavenumber and eigenfrequency, which are believed to contribute to the zonal flow profile associated with the solar differential rotation. For this purpose, we have sought interacting triads containing a zero zonal wavenumber mode yielding unstable solutions (in the modulational sense according to Connaughton et al. 2010).

Consequently, we have analyzed a representative example of such triads in which the triad is connected via one wave mode to a second triplet that is nearly resonant. Numerical integrations of the five-wave system show that the energy transfers between the two wave triplets allowed by modulational type instability yield long-term modulations on the main approximately 11/22 yr cycle associated with intratriad energy exchanges. In addition, the zonal flow mode amplitude modulations are found to be approximately in opposite phase with the amplitude oscillations of the second triad, which is supposedly related to the Schwabe cycle, according to our theoretical model. This result is consistent with the observational work of Zhang et al. (2015), who showed that modulations of the Schwabe cycle exhibit significant negative correlation with observed variations of the solar differential rotation strength.

When analyzing the dynamics of the five-wave system in the presence of a divergence forcing and dissipation, a remarkable resemblance is found between the time evolution of the wave amplitudes and the observed long-term variations of the solar cycle, with an 11/22 yr cycle being modulated at timescales 1 order of magnitude longer (100 yr), as well as the emergence of periods of suppressed wave activity lasting several decades that resemble the grand minimum states. In addition, we have demonstrated that the amplitude of the Rossby wave "activity" is inversely proportional to the instantaneous period of nonlinear energy exchange. Similarly, observations of the solar cycle point out that the amplitude of the cycle, which is commonly measured by the number of sunspots at the peak phase of the cycle, is inversely proportional to the duration of the cycle. This relation between the strength and duration of the solar cycle is described by the so-called Waldmeier law. Therefore, we argue here that the modulation-like instabilities involving MHD Rossby wave triads might be a possible mechanism behind the long-term modulations of the solar cycle observed in sunspot number time series.

It was shown recently by Raphaldini & Raupp (2015) that the propagation of the magnetic branch of MHD Rossby modes is confined to an equatorial belt extending from −35° to +35° in latitude and refracted toward the equator, similar to the sunspot evolution during the solar cycle depicted by the butterfly diagram. We believe that combining the arguments of the present paper with Raphaldini & Raupp (2015) makes MHD Rossby waves strong candidates to play a major role in the dynamics of solar activity. Therefore, we have to speculate on the possible link between Rossby waves and solar magnetic activity.

A dynamo model provided by Rossby wave motions was first suggested by Gilman (1969a, 1969b) by using a two-layer quasi-geostrophic model. As further argued recently by the same author (Gilman & Dikpati 2014), baroclinic Rossby waves and instability should be able to account for a dynamo mechanism, since they provide both vorticity and small vertical motions, which constitute the necessary physical ingredients for the alpha effect, which is known to be associated with helicity. Such a combination of vorticity and vertical motions could also be able to amplify the poloidal component of the magnetic field at the expense of the toroidal one. Smaller-scale instabilities could also create ascending filaments of magnetic field associated with sunspots. Possible extensions of the present study include the analysis of larger clusters of nonlinearly interacting triads and the possibility of the emergence of self-organized synchronized states, such as in Chian et al. (2010), that could explain the approximately periodic nature of the solar dynamo.

B.R. would like to thank the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP; grant No. 2017/23417-5), the FAPESP-PACMEDY project (grant No. 2015/50686-1), and CNPq (grant No. 235203/2014-9/PDE). The work of C.F.M.R. and A.S.T. was supported by CAPES IAG/USP PROEX (grant No. 0531/2017). C.F.M.R. also acknowledges support from the FAPESP-PACMEDY project (grant No. 2015/50686-1). We would like to thank Diego Lorenzo-Oliveira for the references on star activity.

In Section 2 we introduced the coupling coefficients that arise from the nonlinear terms in the perturbation equations, which have the form of a Jacobian operator ${ \mathcal J }(.,.)$. Here we provide a more detailed description of the derivation of such coefficients. The definition of the interaction coefficient C_j,k,l involving three arbitrary modes, j, k, and l, is the projection, in terms of the pseudoenergy norm, of the nonlinear term applied to modes l, k onto the first mode j:

$$\begin{eqnarray}\begin{array}{rcl}{C}_{j,k,l} & = & \displaystyle \frac{1}{{E}_{j}}{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{-\tfrac{\pi }{2}}^{\tfrac{\pi }{2}}\langle { \mathcal B }({{\boldsymbol{u}}}_{{\boldsymbol{k}}},{{\boldsymbol{u}}}_{{\boldsymbol{l}}})\\ & & +\ { \mathcal B }({{\boldsymbol{u}}}_{{\boldsymbol{l}}},{{\boldsymbol{u}}}_{{\boldsymbol{k}}}),{{\boldsymbol{u}}}_{{\boldsymbol{j}}}\rangle {a}^{2}\cos \theta d\theta d\phi ,\end{array}\end{eqnarray} \tag{ 45 }$$

where the vector ${{\boldsymbol{u}}}_{k}$ describes the spatial structure of a particular eigensolution of the linear perturbation equations (similarly for the jth and lth modes)

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{k}(\phi ,\theta )=\left[\begin{array}{c}{\psi }_{k}(\phi ,\theta )\\ {A}_{k}(\phi ,\theta )\end{array}\right],\end{eqnarray} \tag{ 46 }$$

with

$$\begin{eqnarray}&&{\psi }_{k}(\phi ,\theta )={N}_{{n}_{k}}^{{m}_{k}}{P}_{{n}_{k}}^{{m}_{k}}(\sin \theta ){e}^{{{im}}_{k}\phi }{R}_{1}^{(k)},\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}&&{A}_{k}(\phi ,\theta )={N}_{{n}_{k}}^{{m}_{k}}{P}_{{n}_{k}}^{{m}_{k}}(\sin \theta ){e}^{{{im}}_{k}\phi }{R}_{2}^{(k)}.\end{eqnarray} \tag{ 48 }$$

In the equations above, ${R}_{1}^{(k)}=\tfrac{{\omega }_{k}}{{m}_{k}}$ and ${R}_{2}^{(k)}=\tfrac{{B}_{0}}{a}$ refer to the components of the corresponding eigenvector ${{\boldsymbol{R}}}^{(k)}$ defined by Equation (12), and E_j is the pseudoenergy norm of the jth eigenmode, given by

$$\begin{eqnarray}&&{E}_{j}={\int }_{0}^{2\pi }{\int }_{-\tfrac{\pi }{2}}^{\tfrac{\pi }{2}}\left[| {\rm{\nabla }}{\psi }_{j}{| }^{2}+\displaystyle \frac{1}{{\mu }_{0}\rho }| {\rm{\nabla }}{A}_{j}{| }^{2}\right]{a}^{2}\cos \theta d\theta d\phi .\end{eqnarray} \tag{ 49 }$$

Therefore, to evaluate the nonlinear coupling coefficient ${C}_{j,k,l}$, one needs to obtain the pseudoenergy norm from the nonlinear perturbation equation (Equation (13)). As the first equation is written for ${{\rm{\nabla }}}^{2}\psi$ and the second for A, one must multiply the first equation by ${\psi }^{* }$ and the second by ${{\rm{\nabla }}}^{2}{A}^{* }$ and integrate by parts to yield the pseudoenergy norm, where the superscript * denotes the complex conjugate. Consequently, to be consistent with the pseudoenergy norm, one projects the nonlinear term onto the corresponding adjoint eigensolution ${{\boldsymbol{u}}}_{j}^{\dagger }$ referred to the jth mode, given by

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{j}^{\dagger }=\left[\begin{array}{c}{\psi }_{j}^{* }\\ {{\rm{\nabla }}}^{2}{A}_{j}^{* }\end{array}\right].\end{eqnarray} \tag{ 50 }$$

In this way, substituting the ansatz of Equation (15) into the nonlinear perturbation given by Equation (13), multiplying the resulting equations by the adjoint eigensolution of mode j given above, integrating by parts the resulting equations, and using the boundary conditions (periodic solutions in ϕ and regularity at the poles), as well as the orthogonality relations, we get

$$\begin{eqnarray}&&\begin{array}{l}{E}_{j}\displaystyle \frac{d{{\rm{\Lambda }}}_{j}}{{dt}}-{E}_{j}i{\omega }_{j}{{\rm{\Lambda }}}_{j}={{\rm{\Lambda }}}_{k}(t){{\rm{\Lambda }}}_{l}(t)\times {\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{-\tfrac{\pi }{2}}^{\tfrac{\pi }{2}}\langle { \mathcal B }({{\boldsymbol{u}}}_{{\boldsymbol{k}}},{{\boldsymbol{u}}}_{{\boldsymbol{l}}})\\ +{ \mathcal B }({{\boldsymbol{u}}}_{{\boldsymbol{l}}},{{\boldsymbol{u}}}_{{\boldsymbol{k}}}),{{\boldsymbol{u}}}_{{\boldsymbol{j}}}\rangle {a}^{2}\cos \theta d\theta d\phi ,\end{array}\end{eqnarray} \tag{ 51 }$$

where the nonlinear operator is given by

$$\begin{eqnarray}&&{ \mathcal B }({\boldsymbol{u}},{\boldsymbol{u}})=\left[\begin{array}{c}-{ \mathcal J }(\psi ,{{\rm{\nabla }}}^{2}\psi )+\tfrac{1}{{\mu }_{0}\rho }{ \mathcal J }(A,{{\rm{\nabla }}}^{2}A)\\ { \mathcal J }(\psi ,A)\end{array}\right].\end{eqnarray} \tag{ 52 }$$

Consequently, the nonlinear terms in the equations above can be explicitly written in terms of the eigenmodes as

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal J }({\psi }_{k},{{\rm{\nabla }}}^{2}{\psi }_{l})+{ \mathcal J }({\psi }_{l},{{\rm{\nabla }}}^{2}{\psi }_{k})\\ \quad =\,\left({{im}}_{k}{P}_{{n}_{k}}^{{m}_{k}}\displaystyle \frac{{{dP}}_{{n}_{l}}^{{m}_{l}}}{d\theta }-{{im}}_{l}{P}_{{n}_{l}}^{{m}_{l}}\displaystyle \frac{{{dP}}_{{n}_{k}}^{{m}_{k}}}{d\theta }\right)\displaystyle \frac{{\omega }_{k}}{{m}_{k}}\displaystyle \frac{{\omega }_{l}}{{m}_{l}}\displaystyle \frac{1}{{a}^{2}\cos \theta }\\ \quad \quad \times \,\left({n}_{l}({n}_{l}+1)-{n}_{k}({n}_{k}+1))\right){e}^{i({m}_{k}+{m}_{l})\phi },\end{array}\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{{\mu }_{0}\rho }[{ \mathcal J }({A}_{k},{{\rm{\nabla }}}^{2}{A}_{l})+{ \mathcal J }({A}_{l},{{\rm{\nabla }}}^{2}{A}_{k})]\\ \quad =\,\left({{im}}_{k}{P}_{{n}_{k}}^{{m}_{k}}\displaystyle \frac{{{dP}}_{{n}_{l}}^{{m}_{l}}}{d\theta }-{{im}}_{l}{P}_{{n}_{l}}^{{m}_{l}}\displaystyle \frac{{{dP}}_{{n}_{k}}^{{m}_{k}}}{d\theta }\right)\displaystyle \frac{{B}_{0}^{2}}{{\mu }_{0}\rho {a}^{2}}\displaystyle \frac{1}{{a}^{2}\cos \theta }\\ \quad \quad \times \,\left({n}_{l}({n}_{l}+1)-{n}_{k}({n}_{k}+1))\right){e}^{i({m}_{k}+{m}_{l})\phi },\end{array}\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal J }({\psi }_{l},{A}_{k})+{ \mathcal J }({\psi }_{k},{A}_{l})=\left({{im}}_{k}{P}_{{n}_{k}}^{{m}_{k}}\displaystyle \frac{{{dP}}_{{n}_{l}}^{{m}_{l}}}{d\theta }-{{im}}_{l}{P}_{{n}_{l}}^{{m}_{l}}\displaystyle \frac{{{dP}}_{{n}_{k}}^{{m}_{k}}}{d\theta }\right)\\ \quad \quad \times \displaystyle \frac{{B}_{0}^{2}}{{a}^{2}}\left(\displaystyle \frac{{\omega }_{k}}{{m}_{k}}-\displaystyle \frac{{\omega }_{l}}{{m}_{l}}\right)\displaystyle \frac{1}{{a}^{2}\cos \theta }{e}^{i({m}_{k}+{m}_{l})\phi }.\end{array}\end{eqnarray} \tag{ 55 }$$

Evaluating the integrals in Equation (51), as well as the inner product according to Equation (16), it follows that Equation (51) becomes

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{j}}{{dt}}-i{\omega }_{j}{{\rm{\Lambda }}}_{j}={{\rm{\Lambda }}}_{k}(t){{\rm{\Lambda }}}_{l}(t){C}_{j,k,l},\end{eqnarray} \tag{ 56 }$$

with the coupling coefficient ${C}_{j,k,l}$ being expressed according to

$$\begin{eqnarray}&&{C}_{j,k,l}=\displaystyle \frac{-i}{2a}{K}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}}\left[{I}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}}+{L}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}}\right],\end{eqnarray} \tag{ 57 }$$

where the constants ${I}_{{n}_{j}{n}_{k}{n}_{l}}^{{m}_{j}{m}_{k}{m}_{l}}$, ${L}_{{n}_{j}{n}_{k}{n}_{l}}^{{m}_{j}{m}_{k}{m}_{l}}$ and the coupling integral ${K}_{{n}_{1}{n}_{2}{n}_{3}}^{{m}_{1}{m}_{2}{m}_{3}}$ are given by

$$\begin{eqnarray}&&{I}_{{n}_{j}{n}_{k}{n}_{l}}^{{m}_{j}{m}_{k}{m}_{l}}=\left[({n}_{k}({n}_{k}+1))-({n}_{l}({n}_{l}+1))\right],\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&{L}_{{n}_{j}{n}_{k}{n}_{l}}^{{m}_{j}{m}_{k}{m}_{l}}={\left(\displaystyle \frac{{V}_{{\rm{A}}}}{a}\right)}^{2}\left(\displaystyle \frac{{\omega }_{l}}{{m}_{l}}-\displaystyle \frac{{\omega }_{k}}{{m}_{k}}\right),\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&\begin{array}{l}{K}_{{n}_{j}{n}_{k}{n}_{l}}^{{m}_{j}{m}_{k}{m}_{l}}={N}_{{n}_{j}}^{{m}_{j}}{N}_{{n}_{k}}^{{m}_{k}}{N}_{{n}_{l}}^{{m}_{l}}\\ \quad \qquad \qquad \times {\displaystyle \int }_{-1}^{1}{P}_{{n}_{j}}^{{mj}}\left({m}_{k}{P}_{{n}_{k}}^{{m}_{k}}\displaystyle \frac{{{dP}}_{{n}_{l}}^{{m}_{l}}}{{dz}}-{m}_{l}{P}_{{n}_{l}}^{{m}_{l}}\displaystyle \frac{{{dP}}_{{n}_{k}}^{{m}_{k}}}{{dz}}\right){dz}.\end{array}\end{eqnarray} \tag{ 60 }$$

Here we review some basic features of the dynamics of the three-wave equations. Typically, in nonlinear wave problems with quadratic nonlinearities, the nonlinear interactions involving the normal modes of the linear system are described by a complex chain of three-wave systems of the form

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{1}}{{dt}}=i{\omega }_{1}{{\rm{\Lambda }}}_{1}+{C}_{1}{{\rm{\Lambda }}}_{2}^{* }{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{2}}{{dt}}=i{\omega }_{2}{{\rm{\Lambda }}}_{2}+{C}_{2}{{\rm{\Lambda }}}_{1}^{* }{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{3}}{{dt}}=i{\omega }_{3}{{\rm{\Lambda }}}_{3}+{C}_{3}{{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{2},\end{eqnarray} \tag{ 63 }$$

where Λ_i denotes the complex valued amplitude of the ith wave, ω_i is the corresponding eigenfrequency, and C_i is the corresponding coupling coefficient. The dynamical system described above has three independent conserved quantities, namely, the Hamiltonian

$$\begin{eqnarray}&&H=\mathrm{Im}({{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{2}{{\rm{\Lambda }}}_{3}^{* }),\end{eqnarray} \tag{ 64 }$$

as well as two quantities called Manley–Rowe relations:

$$\begin{eqnarray*}\begin{array}{rcl}{I}_{12} & = & | {{\rm{\Lambda }}}_{1}{| }^{2}+| {{\rm{\Lambda }}}_{2}{| }^{2}\\ {I}_{13} & = & | {{\rm{\Lambda }}}_{1}{| }^{2}+| {{\rm{\Lambda }}}_{3}{| }^{2}.\end{array}\end{eqnarray*}$$

These three conserved quantities make the system integrable, as will be discussed below. In fact, using the polar representation ${{\rm{\Lambda }}}_{j}={A}_{j}{e}^{i{{\rm{\Phi }}}_{j}}$, the complex Equations (61)–(63) can be rewritten as four equations describing the time evolution of the real amplitudes and the combination of the phases represented by Φ:

$$\begin{eqnarray}&&\displaystyle \frac{{{dA}}_{1}}{{dt}}={C}_{1}{A}_{2}{A}_{3}\cos {\rm{\Phi }},\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dA}}_{2}}{{dt}}={C}_{2}{A}_{1}{A}_{3}\cos {\rm{\Phi }},\end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dA}}_{3}}{{dt}}=-{C}_{3}{A}_{1}{A}_{2}\cos {\rm{\Phi }},\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Phi }}}{{dt}}={\rm{\Delta }}\omega +{A}_{1}{A}_{2}{A}_{3}\left(\displaystyle \frac{{C}_{1}}{{A}_{1}^{2}}+\displaystyle \frac{{C}_{2}}{{A}_{2}^{2}}+\displaystyle \frac{{C}_{3}}{{A}_{3}^{2}}\right),\end{eqnarray} \tag{ 68 }$$

where Δω = ω₃ − ω₂ − ω₁ is the mismatch among the mode eigenfrequencies. With the conserved quantities described above, these equations are integrable by quadrature, with the solutions being expressed in terms of Jacobi elliptic functions (see Bustamante & Kartashova 2011 for details):

$$\begin{eqnarray}\begin{array}{rcl}{A}_{1}^{2} & = & -\mu {\left(\displaystyle \frac{2K(\mu )}{{ZT}}\right)}^{2}{{sn}}^{2}\left(\displaystyle \frac{2K(\mu )(t-{t}_{0})}{T},\mu \right)\\ & & +\displaystyle \frac{{I}_{13}}{3}\left(2-\rho +2\sqrt{(}1-\rho -{\rho }^{2})\cos (\alpha /2)\right),\end{array}\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}\begin{array}{rcl}{A}_{2}^{2} & = & -\mu {\left(\displaystyle \frac{2K(\mu )}{{ZT}}\right)}^{2}{{sn}}^{2}\left(\displaystyle \frac{2K(\mu )(t-{t}_{0})}{T},\mu \right)\\ & & +\displaystyle \frac{{I}_{13}}{3}\left(2\rho -1+2\sqrt{(}1-\rho -{\rho }^{2})\cos (\alpha /2)\right),\end{array}\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}\begin{array}{rcl}{A}_{3}^{2} & = & \mu {\left(\displaystyle \frac{2K(\mu )}{{ZT}}\right)}^{2}{{sn}}^{2}\left(\displaystyle \frac{2K(\mu )(t-{t}_{0})}{T},\mu \right)\\ & & +\displaystyle \frac{{I}_{13}}{3}\left(\rho +1-2\sqrt{(}1-\rho -{\rho }^{2})\cos (\alpha /2)\right).\end{array}\end{eqnarray} \tag{ 71 }$$

In the equations above, sn stands for the elliptic sine, and K(μ) is the elliptic integral of the first kind, given by

$$\begin{eqnarray}&&K(\mu )={\int }_{0}^{\pi /2}\displaystyle \frac{d\theta }{\sqrt{1-\mu \sin \theta }},\end{eqnarray} \tag{ 72 }$$

where the argument μ of the elliptic integral above is

$$\begin{eqnarray}&&\mu =\displaystyle \frac{\cos (\alpha /3+\pi /6)}{\cos (\alpha /3-\pi /6)}.\end{eqnarray} \tag{ 73 }$$

The solutions described above are periodic in time, with period T given by

$$\begin{eqnarray}&&T=\displaystyle \frac{\sqrt{2\sqrt{3}}K(\mu )}{Z(1-\rho -{\rho }^{2})\sqrt{{I}_{13}\cos \left(\alpha /3-\pi /6\right)}},\end{eqnarray} \tag{ 74 }$$

where ρ is the ratio between two Manley–Rowe constants, ρ = I₁₃/I₂₃, and the angle $\alpha \in [0,\pi ]$ is defined by

$$\begin{eqnarray}&&\cos \alpha =\displaystyle \frac{\sqrt{2\sqrt{3}}K(\mu )}{Z(1-\rho -{\rho }^{2})\sqrt{{I}_{13}\cos \left(\alpha /3-\pi /6\right)}}.\end{eqnarray} \tag{ 75 }$$

Likewise, the solution for the phases is given by

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}(t)=\mathrm{sign}({{\rm{\Phi }}}_{0})\mathrm{arccot}\left(\displaystyle \frac{\mu }{| H| }{\left(\displaystyle \frac{2K(\mu )}{{ZT}}\right)}^{3}{sn}\left(\displaystyle \frac{2K(\mu )(t-{t}_{0})}{T},\mu \right)\right.\\ \left.\times \ {cn}\left(\displaystyle \frac{2K(\mu )(t-{t}_{0})}{T},\mu \right){dn}\left(\displaystyle \frac{2K(\mu )(t-{t}_{0})}{T},\mu \right)\right),\end{array}\end{eqnarray} \tag{ 76 }$$

with dn and cn indicating the other Jacobi elliptic functions.

Let us now couple two triads of nonlinearly interacting waves through one mode, resulting in the five-wave model

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{1}}{{dt}}=i{\omega }_{1}{{\rm{\Lambda }}}_{1}+{C}_{123}{{\rm{\Lambda }}}_{2}^{* }{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{2}}{{dt}}=i{\omega }_{2}{{\rm{\Lambda }}}_{2}+{C}_{231}{{\rm{\Lambda }}}_{1}^{* }{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{3}}{{dt}}=i{\omega }_{3}{{\rm{\Lambda }}}_{3}+{C}_{312}{{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{2}+{C}_{345}{{\rm{\Lambda }}}_{4}^{* }{{\rm{\Lambda }}}_{5},\end{eqnarray} \tag{ 79 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{4}}{{dt}}=i{\omega }_{4}{{\rm{\Lambda }}}_{4}+{C}_{435}{{\rm{\Lambda }}}_{5}^{* }{{\rm{\Lambda }}}_{3},\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{5}}{{dt}}=i{\omega }_{5}{{\rm{\Lambda }}}_{5}+{C}_{534}{{\rm{\Lambda }}}_{3}{{\rm{\Lambda }}}_{4},\end{eqnarray} \tag{ 81 }$$

where ω_j and C_ijk can be read off from Table 1, namely,

$$\begin{eqnarray*}\begin{array}{rcl}{\omega }_{1} & = & 0,{\omega }_{2}=1.78236\times {10}^{-7},{\omega }_{3}=1.72695\times {10}^{-7},\\ {\omega }_{4} & = & 1.85859\times {10}^{-7},{\omega }_{5}=3.56473\times {10}^{-7},\end{array}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{C}_{123}=-0.200293i,{C}_{231}=-1.75195i,{C}_{312}=-2.15463i,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{C}_{345}=0.620163i,\quad {C}_{453}=0.27978i,\quad {C}_{534}=0.904184i.\end{eqnarray*}$$

In general, Equations (77)–(81) are not integrable, but it is easy to show that the following quadratic functions of the dependent variables ${\{{{\rm{\Lambda }}}_{j}\}}_{j=1}^{5}$ are constants of the motion (Manley–Rowe invariants):

$$\begin{eqnarray}&&I=| {{\rm{\Lambda }}}_{3}{| }^{2}+\displaystyle \frac{| {C}_{345}| }{| {C}_{534}| }| {{\rm{\Lambda }}}_{5}{| }^{2}+\displaystyle \frac{| {C}_{312}| }{| {C}_{123}| }| {{\rm{\Lambda }}}_{1}{| }^{2},\end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray}&&J=\displaystyle \frac{| {{\rm{\Lambda }}}_{1}{| }^{2}}{| {C}_{123}| }-\displaystyle \frac{| {{\rm{\Lambda }}}_{2}{| }^{2}}{| {C}_{231}| },\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&K=\displaystyle \frac{| {{\rm{\Lambda }}}_{4}{| }^{2}}{| {C}_{453}| }+\displaystyle \frac{| {{\rm{\Lambda }}}_{5}{| }^{2}}{| {C}_{534}| }.\end{eqnarray} \tag{ 84 }$$

The constancy of J has a direct interpretation: $| {{\rm{\Lambda }}}_{1}{| }^{2}$ and $| {{\rm{\Lambda }}}_{2}{| }^{2}$ are directly coupled, as evidenced by the plots of the energies of modes with spherical wavenumbers (0, 2) and (1, 10) in Figures 2 and 5. Similarly, the constancy of K means the direct coupling of $| {{\rm{\Lambda }}}_{4}{| }^{2}$ and $| {{\rm{\Lambda }}}_{5}{| }^{2}$, corresponding to spherical wavenumbers (1, 12) and (2, 10) in the same figures. Finally, the constancy of I represents the coupling between triads a and b and means that the energies $| {A}_{1}{| }^{2},| {A}_{3}{| }^{2},| {A}_{5}{| }^{2}$ must lie on a certain spheroid. Notably, from the fact that the numerical factors in these formulae take finite values, it follows that all modes' energies are bounded from above.

Below, we briefly discuss two robust mechanisms of energy transfer between the triads.

Modulational instability. Note that, in case the amplitude and mismatch frequencies are commensurable, a particular solution of triad b alone defines a periodic orbit of the system of five waves. This is done by setting initial conditions Λ₁ = Λ₂ = 0 at t = 0, and the first and second modes' amplitudes will remain zero for all times. We can, in principle, linearize the system around a periodic solution of triad b and analyze the stability of the system to small perturbations on the first and second amplitudes. Such instability is reminiscent of the modulational instability explored in Connaughton et al. (2010) in the case of Rossby waves. In the periodic case, this instability can be studied by using Floquet analysis (Hale 1969). However, in general, the triad equations are quasiperiodic, so the instability analysis can be done by calculating the largest Lyapunov exponent of the system. In order to do this, we use the procedure of Benettin et al. (1976) and the implementation available in Datseris (2018). In the case of Equations (77)–(81), with initial amplitudes for triad b ${{\rm{\Lambda }}}_{3}(0)=9.56(1+i)\times {10}^{-9}$, ${{\rm{\Lambda }}}_{4}(0)={10}^{-2}{{\rm{\Lambda }}}_{3}(0)$, and ${{\rm{\Lambda }}}_{5}(0)=\tfrac{3}{2}{{\rm{\Lambda }}}_{3}(0)$ (and infinitesimally small initial amplitudes for modes Λ₁, Λ₂), the largest Lyapunov exponent is 9 × 10⁻⁹, corresponding to a growth rate of 0.28 yr⁻¹, which seems compatible with the timescales associated with the solar cycle. We illustrate the instability by plotting the real and imaginary parts of mode Λ₄ in a 100 yr integration in Figure 8.

Zoom In Zoom Out Reset image size

Precession resonance. To explain this mechanism, let us consider a simple instance whereby triad (1, 2, 3) initially has low energy in comparison with triad (3, 4, 5). Let us also assume, for simplicity of exposition, that triad (3, 4, 5) is resonant or quasi-resonant (Δω₂ ≡ ω₃ + ω₄ − ω₅ ≈ 0). In this case, as presented before, the time evolution of the mode amplitudes of triad (3, 4, 5) is periodic in time, with period T₂ inversely proportional to the wave amplitudes according to Equation (74). Consequently, in this appropriate amplitude regime, in which the Δω₁ term dominates the corresponding time evolution equation of the relative phase of triad (1, 2, 3), if $\tfrac{2\pi }{{T}_{2}}\approx {\rm{\Delta }}{\omega }_{1}$, the last term on the right-hand side of Equation (79) will act as a resonant forcing for triad (1, 2, 3), making the energy of this triad grow in time. This resonance between the linear frequency mismatch of one triad and the nonlinear frequency of the energy oscillation of the other one was called precession resonance by Bustamante et al. (2014). In this resonant case, Bustamante et al. (2014) showed that there is a strong energy transfer between different wave triads, even in the case of nonresonant interactions (${\rm{\Delta }}{\omega }_{1}\ne 0$). To illustrate this effect, we consider a more appropriate set of initial conditions: ${{\rm{\Lambda }}}_{1}=1.48\alpha \exp (1.53i)\times {10}^{-13}$, ${{\rm{\Lambda }}}_{2}=1.71\alpha \exp (0.304i)\times {10}^{-13}$, ${{\rm{\Lambda }}}_{3}=1.03\alpha \exp (1.11i)\times {10}^{-8}$, ${{\rm{\Lambda }}}_{4}=1.22\alpha \exp (5.06i)\times {10}^{-8}$, and ${{\rm{\Lambda }}}_{5}=1.17\alpha \exp (3.95i)\times {10}^{-8}$, where α is a real scale parameter (of order 1) used to search for the resonance. The phases of these initial conditions were randomly generated uniformly over [0, 2π]; apart from the scale parameter α, the amplitudes $| {{\rm{\Lambda }}}_{1}| ,| {{\rm{\Lambda }}}_{2}|$ were randomly generated uniformly over the domain [1, 2] × 10⁻¹³, and the amplitudes $| {{\rm{\Lambda }}}_{3}| ,| {{\rm{\Lambda }}}_{4}| ,| {{\rm{\Lambda }}}_{5}|$ were randomly generated uniformly over the domain [1, 2] × 10⁻⁸. In order to find the resonance, a simulation of Equations (77)–(81) is done for selected choices of scale parameter α. For each simulation, the quotient

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal Q }(t)\equiv \displaystyle \frac{| {C}_{312}| \,| {{\rm{\Lambda }}}_{1}(t){| }^{2}}{| {C}_{123}| \,I},\quad 0\leqslant { \mathcal Q }(t)\leqslant 1\quad \mathrm{for}\ \mathrm{all}\quad t\geqslant 0,\end{array}\end{eqnarray} \tag{ 85 }$$

where I is defined in Equation (82), is plotted over a long time range (about 1524 yr in this work), and its maximum value over that time range is recorded, giving the so-called efficiency ${ \mathcal E }(\alpha )\equiv {\max }_{t}{ \mathcal Q }(t)$. Plots of ${ \mathcal Q }(t)$ for selected values of the scale parameter (α = 0.45, 0.55, and 0.70) are shown in Figure 9, and a plot of efficiency ${ \mathcal E }(\alpha )$ over an extended range 0.1 ≤ α ≤ 14.0 is shown in Figure 10. Remarkably, the peak of efficiency at α = 0.55 gives 34% efficiency, significantly larger than the calculated efficiency at unfeasibly higher amplitudes α ≫ 1. The time series plotted in Figure 9 show three timescales: a fast one corresponding to the typical nonlinear timescale of the order of 10 yr, an intermediate one corresponding to the envelopes' widths of the order of 100 yr, and a slow one corresponding to the distance between the envelopes, which can vary between 100 and 1000 yr. as we could measure in our extended sweep over values of α between 0.1 and 14. The energy share of this long time range, as a function of α, has a marked peak in the vicinity of the peak at α = 0.55 and also in the vicinity of the transition point at α = 3.40 (figure not shown). These low-frequency peaks provide evidence that precession resonance is the mechanism behind the observed strong energy transfers toward zonal modes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Due to the fact that the dissipation coefficients are different for the velocity and magnetic fields, we have to project the resulting effects onto the eigenvectors corresponding to each wave mode in order to obtain the dissipation coefficients. Consider an equation for a freely decaying vector field,

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\boldsymbol{V}}={ \mathcal D }{\boldsymbol{V}},\end{eqnarray} \tag{ 86 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{V}}=\left[\begin{array}{c}{{\rm{\nabla }}}^{2}\psi \\ A\end{array}\right]\end{eqnarray} \tag{ 87 }$$

and the linear dissipation operator is given by

$$\begin{eqnarray}{ \mathcal D }=\left[\begin{array}{cc}\nu {{\rm{\nabla }}}^{2}({{\rm{\nabla }}}^{2}) & 0\\ 0 & \eta {{\rm{\nabla }}}^{2}({{\rm{\nabla }}}^{2})\end{array}\right].\end{eqnarray} \tag{ 88 }$$

As described before, in the equations above, ν and η are the coefficients of viscous and magnetic diffusivity, respectively. In order to obtain the spectral amplitude equations associated with the freely decaying vector field, one multiplies the first component equation by ψ and the second by ∇²A and integrates by parts, obtaining in the spectral space the following equation for the amplitudes Λ_j:

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Lambda }}}_{j}}{{dt}}=\displaystyle \sum _{i}{{ \mathcal D }}_{{ij}}{{\rm{\Lambda }}}_{i},\end{eqnarray} \tag{ 89 }$$

where i, j denote the particular eigenmode, characterized by a spherical harmonic (m, n) and one of the mode types (slow magnetic or fast hydrodynamic branch). Because of the orthogonality relations of the spherical harmonics, a given mode can be influenced only by itself and the corresponding mode with the same wavenumber but in the opposite branch (other type of wave), since the nonlinear interaction of one mode with another with the same wavenumber is forbidden by the Ellsaesser rules (Ellsaesser 1966). Therefore, we consider only diagonal interactions, whose coefficients are given by Equation (41).