© C. R. Goddard et al. 2020

1. Introduction

Rossby waves are global-scale waves of vorticity which arise due to the conservation of the radial component of the absolute vorticity for rotating fluids. They are most commonly discussed in the context of planetary atmospheres and oceans (e.g. Vallis 2017). Detections of their stellar analogues, also known as r modes, have been proposed for Gamma Doradus stars (Van Reeth et al. 2016; Saio et al. 2018; Saio 2018). Recently Rossby waves have been unambiguously detected in the near-surface flows of the Sun (Löptien et al. 2018); they are characterised by a robust dispersion relation which is close to that of classical Rossby waves in a uniformly rotating fluid (e.g. Saio 1982). Follow-up studies with helioseismic methods and data have confirmed the detection of modes with azimuthal orders m = 3–15 (Liang et al. 2019; Hanasoge & Mandal 2019; Hanson et al. 2020). The measured latitudinal eigenfunctions deviate from sectoral Legendre polynomials (Proxauf et al. 2020). Calculations in the equatorial β-plane suggest that the effect of latitudinal differential rotation can account for the real part of the observed latitudinal eigenfunctions (Gizon et al. 2020). The observed Rossby waves can be described as equatorial modes trapped by latitudinal differential rotation.

The theory of r modes in the solar and stellar context has developed over the last few decades for various limiting cases (e.g. Provost et al. 1981; Wolff 1998; Zaqarashvili et al. 2007; Saio et al. 2018). Damiani et al. (2020) recently computed eigenmodes in slowly rotating polytropes and confirmed the r^m radial dependence of the displacement eigenfunctions proposed by Provost et al. (1981).

Global-scale Rossby waves are probes of the convection zone. In this paper we predict the effect of solar cycle variations in the zonal flows on their frequencies, since these frequencies are directly related to the solar rotation rate. This work is required to interpret future measurements of time variations, in particular to disentangle perturbations caused by solar cycle changes in rotation from perturbations caused by the evolving large-scale magnetic field.

Helioseismology with data from the Michelson Doppler Imager on the Solar and Heliospheric Observatory (SOHO/MDI; Scherrer et al. 1995) and the Global Oscillation Network Group (GONG; Hill et al. 1996) allowed the internal rotation profile of the Sun to be measured in the convection zone. The latitudinal differential rotation seen at the surface extends throughout the convection zone (e.g. Schou et al. 1998). Additionally, data from the Helioseismic and Magnetic Imager on the Solar Dynamics Observatory (SDO/HMI; Schou et al. 2012) confirmed these findings. All data sets revealed changes in the rotation rate that are associated with the solar cycle. These time-varying zonal flows or torsional oscillations were first noted at the surface (Howard & Labonte 1980) and have since been shown to extend throughout most of the convection zone (Vorontsov et al. 2002), although they weaken with depth (e.g. Howe 2009). The time-varying zonal flows consist of bands of slower and faster rotation that move towards the equator at low latitudes and towards the poles at high latitudes. These zonal flows vary significantly between the last two solar cycles, solar cycle 23 (SC23) and solar cycle 24 (SC24; Komm et al. 2014; Howe et al. 2018). Additionally, there is a difference in the average rotation between these two cycles. In SC24, the convection zone rotates more quickly at low latitudes, but more slowly at high latitudes (e.g. Basu & Antia 2019).

Section 2 describes the internal rotation data. In Sect. 3, we use a first-order perturbation approach (Born approximation) to obtain a kernel for the sensitivity of Rossby waves to variations of an arbitrary rotation profile, and the predictions via a series of steady problems are presented in Sect. 4. Finally, the discussion and conclusions are given in Sect. 5.

2. Solar cycle variations of internal rotation

The rotation rate data used to prescribe the zonal flows are the 2D regularised least-squares (RLS) global helioseismology inversions for solar rotation from MDI and HMI (Larson & Schou 2018). The data used here cover the period from 1996 June 06 to 2019 June 30, with a cadence of 72 days. The HMI data are used during the period of overlap between MDI and HMI. We work in spherical-polar coordinates, where r is the radial coordinate, ϕ is longitude, and θ is the co-latitude. A reference rotation profile, $\overline{\mathrm{\Omega }}(r,\theta )$, is taken to be the unweighted mean rotation profile over the 23-year time interval. Variations about this mean profile over time are thus

$$\begin{array}{c}\hfill \delta \mathrm{\Omega }(r,\theta ,t)=\mathrm{\Omega }(r,\theta ,t)-\overline{\mathrm{\Omega }}(r,\theta ).\end{array}$$(1)

The few missing time steps are dealt with via linear interpolation. The rotation residual is then binned by taking the unweighted mean for each calendar year.

Figure 1 shows the binned rotation rate residual, δΩ(r, θ, t), at r/R_⊙ = 0.90, 0.95, and 0.99 for the specified time window. The zonal flows and high latitude variations associated with the 11-year solar cycle are seen alongside the asymmetry between SC23 and SC24. As mentioned previously, in the convection zone SC24 has higher rotation rates at low latitudes, but lower rotation rates at high latitudes. This asymmetry is not due to the inconsistencies between the MDI and HMI data, and is clearly seen in GONG data (Basu & Antia 2019).

Fig. 1. Rotation rate residuals, δΩ(r, θ, t)/2π, at r = 0.9, 0.95, and 0.99 R⊙ for 1996-Jun.-06 to 2019-Jun.-30. The data is binned for each calendar year, and the colour scale is saturated at ±3 nHz. The vertical dashed grey lines denote the positions of the maxima of solar cycles 23 and 24 (2002 and 2014). The vertical dashed brown line indicates the minimum between the two cycles (2009).

3. Perturbation to Rossby mode frequencies due to variations in rotation

We utilise first-order perturbation theory to obtain a kernel for the sensitivity of Rossby wave frequencies to small variations in the rotation profile. This method is commonly used to obtain sensitivity kernels for the effect of rotation on p modes (e.g. Gough & Thompson 1990; Howe 2009; Basu 2016).

3.1. Oscillation equation

We start from the linearised equation of motion for the displacement vector ξ, written in an inertial frame (Lynden-Bell & Ostriker 1967; Unno et al. 1989). Assuming a steady rotation rate Ω(r, θ) (the wave period is shorter than the timescale of the evolution of rotation) and a displacement vector proportional to exp(−iωt), we have

$$\begin{array}{c}\hfill {\omega }^2\mathit{\xi }+\omega B(\mathit{\xi })-T(\mathit{\xi })=V(\mathit{\xi })+P(\mathit{\xi }),\end{array}$$(2)

with

$$\begin{array}{cc}\hfill B(\mathit{\xi })& =2i(\mathit{u}·\mathbf{\nabla })\mathit{\xi },\hfill \\ \hfill T(\mathit{\xi })& ={(\mathit{u}·\mathbf{\nabla })}^2\mathit{\xi },\hfill \\ \hfill V(\mathit{\xi })& =-d(\mathbf{\nabla }\psi ),\hfill \\ \hfill P(\mathit{\xi })& =d\left(\frac{1}{\rho }\mathbf{\nabla }p\right),\hfill \end{array}$$

where

$$\begin{array}{c}\hfill \mathit{u}=\mathbf{\Omega }\times \mathit{r}=rsin\theta \mathrm{\Omega }(r,\theta )\widehat{\mathit{\varphi }}\end{array}$$(3)

is the rotational velocity, r is the radial position vector, p is pressure, ψ is the gravitational potential, and d is the Lagrangian change operator (e.g. dp = p′+ξ ⋅ ∇p₀). Frictional and electromagnetic forces are not considered.

The following relation holds if we assume the azimuthal dependence of the displacement to be proportional to exp(imϕ), where m is the azimuthal order of a solar Rossby mode:

$$\begin{array}{c}\hfill (\mathit{u}·\mathbf{\nabla })\mathit{\xi }=im\mathrm{\Omega }\mathit{\xi }+\mathbf{\Omega }\times \mathit{\xi }.\end{array}$$(4)

Using this relationship we can rewrite the operators that depend on u as

$$\begin{array}{cc}\hfill B(\mathit{\xi })& =-2m\mathrm{\Omega }\mathit{\xi }+2i\mathbf{\Omega }\times \mathit{\xi },\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill T(\mathit{\xi })& =-m^2{\mathrm{\Omega }}^2\mathit{\xi }+2im\mathrm{\Omega }\mathbf{\Omega }\times \mathit{\xi }+\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathit{\xi }).\hfill \end{array}$$(6)

Substituting these into Eq. (2), and defining ℒ(ξ) = V(ξ)+P(ξ), gives

$$\begin{array}{c}\hfill {(\omega -m\mathrm{\Omega })}^2\mathit{\xi }+2i(\omega -m\mathrm{\Omega })\mathbf{\Omega }\times \mathit{\xi }-\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathit{\xi })-\mathcal{L}(\mathit{\xi })=0.\end{array}$$(7)

3.2. First-order oscillation equation

We now perturb Eq. (7) around an arbitrary reference rotation rate Ω₀(r, θ) using

$$\begin{array}{cc}\hfill \mathrm{\Omega }(r,\theta )& ={\mathrm{\Omega }}_0(r,\theta )+\delta \mathrm{\Omega }(r,\theta ),\hfill \\ \hfill \omega & ={\omega }_0+\delta \omega ,\hfill \\ \hfill \mathit{\xi }& ={\mathit{\xi }}_0+\delta \mathit{\xi },\hfill \\ \hfill \mathcal{L}& ={\mathcal{L}}_0+\delta \mathcal{L}.\hfill \end{array}$$(8)

The zero-order eigenvalue problem for (ω₀, ξ₀) is

$$\begin{array}{c}\hfill {({\omega }_0-m{\mathrm{\Omega }}_0)}^2{\mathit{\xi }}_0+2i({\omega }_0-m{\mathrm{\Omega }}_0){\mathbf{\Omega }}_0\times {\mathit{\xi }}_0-{\mathbf{\Omega }}_0\times ({\mathbf{\Omega }}_0\times {\mathit{\xi }}_0)-{\mathcal{L}}_0({\mathit{\xi }}_0)=0.\end{array}$$(9)

The first-order eigenvalue problem for (δω, δξ) is

$$\begin{array}{cc}& {({\omega }_0-m{\mathrm{\Omega }}_0)}^2\delta \mathit{\xi }+2({\omega }_0-m{\mathrm{\Omega }}_0)(\delta \omega -m\delta \mathrm{\Omega }){\mathit{\xi }}_0\hfill \\ \hfill & +2i({\omega }_0-m{\mathrm{\Omega }}_0){\mathrm{\Omega }}_0\times \delta \mathit{\xi }+2i({\omega }_0-m{\mathrm{\Omega }}_0)\delta \mathbf{\Omega }\times {\mathit{\xi }}_0\hfill \\ \hfill & +2i(\delta \omega -m\delta \mathrm{\Omega }){\mathrm{\Omega }}_0\times {\mathit{\xi }}_0-\delta \mathbf{\Omega }\times ({\mathbf{\Omega }}_0\times {\mathit{\xi }}_0)-{\mathbf{\Omega }}_0\times (\delta \mathbf{\Omega }\times {\mathit{\xi }}_0)\hfill \\ \hfill & -{\mathbf{\Omega }}_0\times ({\mathbf{\Omega }}_0\times \delta \mathit{\xi })-{\mathcal{L}}_0(\delta \mathit{\xi })-\delta \mathcal{L}({\mathit{\xi }}_0)=0.\hfill \end{array}$$(10)

3.3. First-order frequency shift

We define an inner product of two functions as ⟨η, ξ⟩ = ∫_Vη^* ⋅ ξρ dV, where V is the volume of the star, the asterisk denotes the complex conjugate, and ρ(r) is the radial density profile.

Dotting Eq. (10) with ${\mathit{\xi }}_0^{*}$ on the left, integrating over V, and utilising the self-adjointness of the operators (cf. Lynden-Bell & Ostriker 1967) gives

$$\begin{array}{cc}& 2⟨{\mathit{\xi }}_0,({\omega }_0-m{\mathrm{\Omega }}_0)(\delta \omega -m\delta \mathrm{\Omega }){\mathit{\xi }}_0⟩+2i⟨{\mathit{\xi }}_0,({\omega }_0-m{\mathrm{\Omega }}_0)\delta \mathrm{\Omega }\widehat{\mathit{z}}\times {\mathit{\xi }}_0⟩\hfill \\ \hfill & +2i⟨{\mathit{\xi }}_0,(\delta \omega -m\delta \mathrm{\Omega }){\mathrm{\Omega }}_0\widehat{\mathit{z}}\times {\mathit{\xi }}_0⟩-2⟨{\mathit{\xi }}_0,\delta \mathrm{\Omega }{\mathrm{\Omega }}_0\widehat{\mathit{z}}\times (\widehat{\mathit{z}}\times {\mathit{\xi }}_0)⟩\hfill \\ \hfill & -⟨{\mathit{\xi }}_0,\delta \mathcal{L}({\mathit{\xi }}_0)⟩=0,\hfill \end{array}$$(11)

where $\widehat{\mathit{z}}$ is the unit vector along the rotation axis, $\widehat{\mathit{z}}=\widehat{\mathit{r}}cos\theta -\widehat{\mathit{\theta }}sin\theta$. Using this expression the frequency shift due to variations in an arbitrary rotation profile may be calculated from the unperturbed (zero-order) displacement eigenfunctions. Rearranging the terms that include δω, and defining 𝒮(ξ) to contain the centrifugal and δℒ terms,

$$\begin{array}{c}\hfill \mathcal{S}({\mathit{\xi }}_0)=\delta \mathrm{\Omega }{\mathrm{\Omega }}_0\widehat{\mathit{z}}\times (\widehat{\mathit{z}}\times {\mathit{\xi }}_0)+\frac{1}{2}\delta \mathcal{L}({\mathit{\xi }}_0),\end{array}$$(12)

we have

$$\begin{array}{c}\hfill \delta \omega =\frac{⟨{\mathit{\xi }}_0,\delta \mathrm{\Omega }\mathcal{N}({\mathit{\xi }}_0)+\mathcal{S}({\mathit{\xi }}_0)⟩}{⟨{\mathit{\xi }}_0,\mathcal{D}({\mathit{\xi }}_0)⟩},\end{array}$$(13)

with

$$\begin{array}{cc}\hfill \mathcal{N}({\mathit{\xi }}_0)& =-i({\omega }_0-2m{\mathrm{\Omega }}_0)(\widehat{\mathit{z}}\times {\mathit{\xi }}_0)+m({\omega }_0-m{\mathrm{\Omega }}_0){\mathit{\xi }}_0,\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill \mathcal{D}({\mathit{\xi }}_0)& =({\omega }_0-m{\mathrm{\Omega }}_0){\mathit{\xi }}_0+i{\mathrm{\Omega }}_0\widehat{\mathit{z}}\times {\mathit{\xi }}_0.\hfill \end{array}$$(15)

We note that the denominator is more complicated than in the standard p-mode case with a non-rotating reference model. For the rest of the paper, we neglect the term 𝒮(ξ) in Eq. (13), which corresponds to neglecting the perturbation of the centrifugal acceleration and the perturbation to the pressure and gravity terms due to the change in rotation δΩ.

3.4. Rossby wave eigenfunction and resulting kernel

From here forwards we consider reference eigenfunctions ξ₀ which correspond to incompressible global Rossby waves for the case of slow uniform rotation. In this case deviations of the solar structure (i.e. oblateness) from variations in rotation have a negligible impact on the eigenfrequencies (Damiani et al. 2020). The relevant zero-order equation of motion, Eq. (9) without the centrifugal (third) term, is now

$$\begin{array}{c}\hfill {({\omega }_0-m{\mathrm{\Omega }}_0)}^2{\mathit{\xi }}_0+2i({\omega }_0-m{\mathrm{\Omega }}_0){\mathbf{\Omega }}_0\times {\mathit{\xi }}_0-{\mathcal{L}}_0({\mathit{\xi }}_0)=0.\end{array}$$(16)

This is the inertial frame equivalent of the oscillation equation given in Provost et al. (1981) and Damiani et al. (2020), under the Cowling approximation.

For the case of slow uniform rotation, the eigenfunction of the Rossby mode with azimuthal wavenumber m is of the form

$$\begin{array}{cc}\hfill {\mathit{\xi }}_0& ={\mathbf{\nabla }}_h\times \left(\mathit{r}R(r)L(\theta ){\mathrm{e}}^{\mathrm{i}(m\varphi -{\omega }_{0}t)}\right)\hfill \\ \hfill & =R(r)\left(\widehat{\mathit{\theta }}\frac{imL(\theta )}{sin\theta }-\widehat{\mathit{\varphi }}{L}^{\prime }(\theta )\right){\mathrm{e}}^{\mathrm{i}(m\varphi -{\omega }_{0}t)},\hfill \end{array}$$(17)

where ∇_h is the horizontal gradient, R(r) and L(θ) are functions to be specified for each mode, and L′(θ) = dL/dθ.

The frequency shift to the unperturbed eigenfrequency ω₀ of mode m due to a change in rotation is given by

$$\begin{array}{c}\hfill \delta \omega ={\displaystyle {\int }_0^{\pi }{\int }_0^{R_{\odot }}K(r,\theta )\delta \mathrm{\Omega }(r,\theta )\mathrm{d}r\mathrm{d}\theta ,}\end{array}$$(18)

where the sensitivity kernels K are given by

$$\begin{array}{c}\hfill K(r,\theta )=\frac{{\mathit{\xi }}_0^{\ast }·\mathcal{N}({\mathit{\xi }}_0)\rho r^{2}sin\theta }{{\int }_0^{\pi }{\int }_0^{R_{\odot }}{\mathit{\xi }}_0^{\ast }·\mathcal{D}({\mathit{\xi }}_0)\rho r^{2}sin\theta \mathrm{d}r\mathrm{d}\theta }·\end{array}$$(19)

For the displacement eigenfunction given in Eq. (17) and a real-valued L(θ), we have

$$\begin{array}{cc}\hfill {\mathit{\xi }}_0^{\ast }·\mathcal{N}({\mathit{\xi }}_0)=& m{R}^2(r)[(4m{\mathrm{\Omega }}_0(r,\theta )-2{\omega }_0)\left(\frac{L(\theta )L^{\prime }(\theta )cos\theta }{sin\theta }\right)\hfill \\ \hfill & +({\omega }_0-m{\mathrm{\Omega }}_0(r,\theta ))(\frac{m^2{L}^2(\theta )}{{sin}^2\theta }+L^{\prime }{(\theta )}^2)]\hfill \end{array}$$(20)

and

$$\begin{array}{cc}\hfill {\mathit{\xi }}_0^{\ast }·\mathcal{D}({\mathit{\xi }}_0)=& R^2(r)[2m{\mathrm{\Omega }}_0(r,\theta )\left(\frac{L(\theta )L^{\prime }(\theta )cos\theta }{sin\theta }\right)\hfill \\ \hfill & +({\omega }_0-m{\mathrm{\Omega }}_0(r,\theta ))(\frac{m^2{L}^2(\theta )}{{sin}^2\theta }+L^{\prime }{(\theta )}^2)].\hfill \end{array}$$(21)

Now we choose specific radial and latitudinal dependences for the eigenfunction given by Eq. (17), along with a corresponding reference eigenfrequency. Considering the case of uniform slow rotation

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_0(r,\theta )={\mathrm{\Omega }}_0=\text{const.},\end{array}$$(22)

we take a sectoral spherical harmonic (ℓ = m) solution with no radial node:

$$\begin{array}{cc}\hfill R(r)& =r^m\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill L(\theta )& ={(sin\theta )}^m\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\omega }_0& =m{\mathrm{\Omega }}_0-\frac{2{\mathrm{\Omega }}_0}{m+1}·\hfill \end{array}$$(25)

This simple form for the radial eigenfunctions is suggested by the study of r modes in slowly (and uniformly) rotating polytropes in Damiani et al. (2020). Corrections may be needed to account for strong radial gradients in solar rotation and for the sub-adiabatic stratification below the convection zone; however, this is beyond the scope of this paper. We note that the kinetic energy density of most of the modes we consider here (m ≥ 4) is largely restricted to the convection zone.

For this case, the kernel simplifies to

$$\begin{array}{c}\hfill K(r,\theta )=\frac{\rho r^{2(m+1)}{(sin\theta )}^{2m-1}[m-(m^2+2){cos}^2\theta ]}{{\int }_0^{R_{\odot }}{r}^{2(m+1)}\rho \mathrm{d}r{\int }_0^{\pi }{(sin\theta )}^{2m-1}(1-m{cos}^2\theta )\mathrm{d}\theta }·\end{array}$$(26)

This kernel does not depend on the uniform rotation rate Ω₀ from the zero-order model. Appendix A shows that this kernel predicts the correct answer for the case of a perturbation δΩ that is constant.

Figure 2 shows the kernel calculated with the smoothed solar-like density profile from Gizon et al. (2017). Figure 3 shows normalised latitudinal and radial cuts for azimuthal orders m = 3, 9, and 15. The frequency shifts are mostly sensitive to perturbations of the rotation rate in the low latitude convection zone for m ≥ 3. The negative higher latitude component of the kernel reduces in amplitude with increasing azimuthal order. The detected solar Rossby waves currently span m = 3–15, and hence the plots are restricted to m ≤ 15. Appendix B discusses the latitudinal dependence of the kernel in more detail.

Fig. 2. Kernels for the sensitivity of sectoral Rossby waves to variations in rotation rate for azimuthal orders m = 1, 3, 6, 9, 12, and 15 (from left to right, top to bottom). The different panels show the kernels in the meridional plane and the kernels are symmetric about the equator. The vertical axis is the rotation axis. The spatial grid is in units of R⊙, and the grided region denotes the convection zone.

Fig. 3. Radial (top) and latitudinal (bottom) dependences of the kernel for sectoral Rossby modes with azimuthal orders m = 3, 9, and 15, normalised by their maxima. The vertical dashed line in the top panel denotes the position of the base of the convection zone at r/R⊙ = 0.7.

4. Results

4.1. Predicted frequency shifts over the last two solar cycles

Here we compute the Rossby wave frequency shifts δω for each mode m, which would be expected from the mean-subtracted rotation residuals, δΩ(r, θ, t), defined in Eq. (1). The rotation residuals are binned by calendar year, as shown in Fig. 1. We do not consider the time-independent frequency shift for each m caused by the difference between our zero-order model (uniform rotation, Ω₀) and the time-averaged differential rotation of the Sun ($\overline{\mathrm{\Omega }}(r,\theta )$) used to define δΩ(r, θ, t).

The frequency-shift time series for each azimuthal order were calculated from the binned rotation residuals and Eqs. (18) and (26) using trapezoidal rule integration. The perturbation to the rotation profile δΩ(r, θ, t) is set to zero in the radiative zone (r/R_⊙ >  0.7), due to the inaccuracies in the inversions there and the lack of zonal flow signatures at these depths.

Figure 4 shows the obtained frequency shifts separately for low and high azimuthal orders, 1 ≤ m ≤ 7 and 8 ≤ m ≤ 15, respectively. For m = 1 no discernible trend is seen as this mode has little sensitivity to the strongest zonal flows in the upper convection zone. From m = 2 to m = 7 an apparent cycle with a period of > 22 years is observed: the frequency shift is negative for SC23 and positive for SC24. For m >  7 a double-peaked trend is seen, similar to the 11-year Sunspot cycle, but with the frequency during SC24 systematically higher than in SC23. This is a superposition of the 11-year cycle of the zonal flows with the asymmetry of the two cycles (as seen for low m). Appendix C discusses the relationship between the frequency shift and the rotation residuals in more detail. Appendix D shows a comparison to frequency shifts computed from direct numerical calculations in the β-plane (Gizon et al. 2020).

Fig. 4. Rossby mode frequency-shift time series due to variation in the solar rotation profile, plotted for m from 1 to 8 (top), and 9 to 15 (middle), calculated for each calendar year. The vertical grey dashed lines denote the positions of the maxima of solar cycles 23 and 24, to the nearest year. The vertical brown dashed line indicates the minimum between the two cycles. The smoothed sunspot number is plotted in the bottom panel.

As one simple way to quantify these trends, we calculate the change in frequencies between two times for each m. The upper panel of Fig. 5 shows the difference in the frequencies between 2014 and 2002 (two solar maxima) and between 2014 and 2009 (maximum and minimum). The difference between maximum and minimum highlights the 11-year component, and continues to increase for high m, reaching 23 nHz for m = 15. The difference between the cycle maxima highlights the asymmetry between the two cycles, which saturates for m >  8 at about 7 nHz.

Fig. 5. Top: predicted difference in the Rossby mode frequencies between the two solar maxima (grey curve, 2014 minus 2002) and between the maximum of cycle 24 and the preceding solar minimum (brown curve, 2014 minus 2009). The error bars are propagated errors from the scatter in the solar rotation measurements. Bottom: estimated errors in the Rossby mode frequencies due to the finite mode lifetimes as a function of observation time T. These errors can be used to determine whether the predicted frequency shifts from the top panel are detectable.

4.2. Predicted measurement uncertainties due to measured finite mode lifetimes

Let us denote by σ the uncertainty on the measured frequency ω of a Rossby mode m due to the finite mode lifetimes. Under the assumption that Rossby modes are stochastically excited, the uncertainties σ scale as T^−1/2, where T is the observation time. Specifically, we have σ = (Γ/2T)^1/2, where Γ is the mode linewidth in frequency space (Libbrecht 1992).

We used the measurements of the mode linewidths from Liang et al. (2019) to calculate the measurement uncertainties σ for various values T. The estimate of the uncertainties in the r mode frequencies depends on the measured lifetimes. We note that the various methods of helioseismology and granulation tracking give essentially the same answers for the mode linewidths (Hanson et al. 2020). The bottom panel of Fig. 5 shows this for several values of m. We note the strong dependence on m. For some azimuthal orders and T = 3 years, for example, the uncertainties are smaller than the predicted solar cycle frequency shifts. Thus, the detection of the Rossby mode frequency shifts due to the time-varying zonal flows may be possible. However, as this is near the limit of detectability some averaging over several values of m may be required to reduce noise.

5. Conclusion

We used first-order perturbation theory to obtain a kernel for the frequency shift of solar Rossby waves caused by perturbations in the rotation profile. We used these kernels together with solar torsional oscillation data to predict the time varying component of the frequency shift for solar Rossby waves. The predicted frequency shifts over the last two solar cycles have a strong dependence on the azimuthal order m. Modes with m ≤ 7 are mainly sensitive to the asymmetry between the two cycles, while those with m >  7 are mainly sensitive to the 11-year cycle of the zonal flows. The predicted frequency difference between cycle maximum and minimum reaches 23 nHz for m = 15.

Several approximations and simplifications were made in the construction of the kernel, and in the chosen latitudinal and radial dependence for the sectoral Rossby waves. These are adequate for the problem at hand, but may require further development once the frequency shifts have been observed. For example, it might be useful to consider eigenfunctions from a differentially rotating reference model. In future work we intend to investigate potential solar cycle changes in the amplitudes and the lifetimes of the modes, as well as potential changes in the eigenfunctions due to the torsional oscillations or the magnetic field.

The predicted shift of the mode frequencies due to variation in the internal rotation over the last two solar cycles is near the limit of detection. A detection of this effect would help further constrain the physics of solar Rossby modes. Should there be much larger solar cycle variations in the mode frequencies, they will have to be attributed to the effect of the solar magnetic field in the convection zone.

Acknowledgments

We thank Zhi-Chao Liang and Chris S. Hanson for useful discussions. Author contributions: LG proposed the research. CRG performed the kernel calculations and the data analysis. ACB and LG contributed to the kernel calculations. DF performed the β-plane calculations in Appendix D. CRG wrote the draft paper and all authors contributed to the final manuscript. This research is supported by ERC Synergy grant WHOLE SUN #810218. The data were processed at the German Data Center for SDO, funded by the German Aerospace Center (DLR). LG acknowledges partial support from the NYU Abu Dhabi Center for Space Science under grant G1502. The HMI and MDI helioseismology data products are used courtesy of the JSOC team. NumPy (van der Walt et al. 2011) and matplotlib (Hunter 2007) were used in this research.

References

Appendix A: Case δΩ = const.

The frequency shift in any reference frame for the described sectoral Rossby modes caused by a change in the rotation rate is

$$\begin{array}{c}\hfill \delta \omega =m\delta \mathrm{\Omega }-\frac{2\delta \mathrm{\Omega }}{m+1}·\end{array}$$(A.1)

The kernel given should reduce to this result in the case of uniform rotation (Ω₀ = const) with a uniform perturbation (δΩ = const). This has been shown via numerical integration, once converged the resulting frequency shifts agree with Eq. (A.1) to machine precision.

The kernel may also be integrated directly to obtain the frequency shift for the case of uniform rotation and perturbation. Defining the frequency shift as

$$\begin{array}{c}\hfill \delta \omega =\delta \mathrm{\Omega }{\displaystyle {\int }_0^{\pi }{\int }_0^{R_{\odot }}K(r,\theta )\mathrm{d}r\mathrm{d}\theta }\end{array}$$(A.2)

and taking the kernel given in Eq. (26), the radial component cancels from top and bottom giving

$$\begin{array}{c}\hfill \delta \omega =\frac{\delta \mathrm{\Omega }{\int }_0^{\pi }{(sin\theta )}^{(2m-1)}[m-(m^2+2){cos}^2\theta \mathrm{d}\theta ]}{{\int }_0^{\pi }{(sin\theta )}^{(2m-1)}[1-m{cos}^2\theta ]\mathrm{d}\theta }·\end{array}$$(A.3)

Utilising the identity $I_p={\int }_0^{\pi }{(\sin \theta )}^p\text{d}\theta =[(p-1)/p]I_{p-2}$, we obtain

$$\begin{array}{c}\hfill \delta \omega =\frac{(m^2+m-2)}{m+1}\delta \mathrm{\Omega }=m\delta \mathrm{\Omega }-\frac{2\delta \mathrm{\Omega }}{m+1},\end{array}$$(A.4)

which is the expected result.

Appendix B: Latitudinal form of the kernel

The figures in Sect. 3.4 show that the sign of the kernel switches at a given latitude, which decreases as the azimuthal order increases. This behaviour is due to the terms from the advection and coriolis force in the final kernel having opposite signs, and different latitudinal dependences. The remaining advective term in the kernel is proportional to m(sin θ)^2m − 1 and the remaining coriolis term is proportional to −(m² + 2)(cos θ)²(sin θ)^2m − 1; therefore, they are equal and opposite at

$$\begin{array}{c}\hfill {\theta }_s=\mathrm{acos}\sqrt{\frac{m}{m^2+2}}·\end{array}$$(B.1)

Figure B.1 shows the latitudinal dependence of the coriolis and advective terms of the kernel, and their sum, for m = 3. The vertical dashed line shows the position as which the two terms are equal and opposite, summing to zero (θ_s). As m increases this position shifts towards the equator and the relative contribution of the coriolis term decreases.

Fig. B.1. Latitudinal dependences of the coriolis and advective terms that contribute to the kernel, and their sum, for m = 3.

Appendix C: Explaining the behaviour of δω for high and low m values

The top three panels of Fig. C.1 show the result of multiplication of the rotation rate residual and the kernel, δΩ(r, θ) K(r, θ), with integration over the radial coordinate. These plots aid in understanding the connection between the rotation data and the calculated frequency shifts for the Rossby waves. We note that faster rotation than the mean profile will result in a positive frequency shift (except where the sensitivity kernel is negative).

Fig. C.1. Top three panels: rotation rate residual multiplied by the sensitivity kernel and integrated over the radial coordinate. Bottom panel: difference in the solar rotation rate in the convection zone between the maxima of cycles 24 and 23.

For m = 3, we see that both high and low latitudes have a negative contribution to the frequency shift for SC23 and positive contribution for SC24. This is due to SC23 having lower rotation rates at low latitudes and higher rotation rates at high latitudes compared to SC24, in combination with the kernel being negative for higher latitudes (see Fig. 2). This is shown in the lower panel of Fig. C.1 where the difference between the rotation residual at the maxima of SC23 and SC24 is plotted. For m = 9 we see a combination of this, and a strong contribution at low latitudes from the zonal flows, which becomes even stronger for m = 15 as the kernel becomes more constrained to low latitudes near the surface (see Fig. 2). In all three plots the switch of the kernel from positive to negative is visible as the pair of horizontal lines of zero frequency at latitudes π/2 − θ_s.

Appendix D: Comparison with calculations in the β-plane

Here we perform a comparison with recent direct numerical calculations of frequencies of viscous Rossby modes in the equatorial β-plane (Gizon et al. 2020) to validate our first-order frequency shifts. These calculations cannot inform us about the effects of the radial variations in rotation, but they can be performed for any latitudinal variations in the zonal flows. In addition, the β-plane approximation is only valid for m ≫ 1.

To make a direct comparison, we define the zonal velocity in the β-plane as a function of the coordinate y = R_⊙cosθ as

$$\begin{array}{c}\hfill U(y,t)=\frac{{\int }_0^{R_{\odot }}rsin\theta [\mathrm{\Omega }(r,\theta ,t)-{\mathrm{\Omega }}_{\text{eq}}(r)]r^{2(m+1)}\rho \mathrm{d}r}{{\int }_0^{R_{\odot }}{r}^{2(m+1)}\rho \mathrm{d}r},\end{array}$$(D.1)

where Ω_eq(r) is the equatorial rotation rate as a function of radius. This flow is used as an input to the computations described in Gizon et al. (2020). In addition, we specify the viscosity through the Reynolds number Re = 300.

The β-plane calculation is performed at the full 72-day cadence of the rotation data, and the frequency-shift time series for each m is then mean subtracted to give a result equivalent to our first-order calculations. These are then binned into calendar years.

Figure D.1 shows the resulting frequency-shift time series along with our first-order results for m = 9, 12, and 15. There is good agreement for m = 12 and 15. The agreement is still reasonable for m = 9. For even smaller values of m the agreement becomes worse as the problem is less suited to the β-plane approximation.

Fig. D.1. Frequency-shift time series from calculations in the β-plane (orange) plotted along with the results from the kernel (blue) for m = 9, 12, and 15.

All Figures

	Fig. 1. Rotation rate residuals, δΩ(r, θ, t)/2π, at r = 0.9, 0.95, and 0.99 R⊙ for 1996-Jun.-06 to 2019-Jun.-30. The data is binned for each calendar year, and the colour scale is saturated at ±3 nHz. The vertical dashed grey lines denote the positions of the maxima of solar cycles 23 and 24 (2002 and 2014). The vertical dashed brown line indicates the minimum between the two cycles (2009).
In the text

Fig. 2. Kernels for the sensitivity of sectoral Rossby waves to variations in rotation rate for azimuthal orders m = 1, 3, 6, 9, 12, and 15 (from left to right, top to bottom). The different panels show the kernels in the meridional plane and the kernels are symmetric about the equator. The vertical axis is the rotation axis. The spatial grid is in units of R⊙, and the grided region denotes the convection zone.

In the text

	Fig. 3. Radial (top) and latitudinal (bottom) dependences of the kernel for sectoral Rossby modes with azimuthal orders m = 3, 9, and 15, normalised by their maxima. The vertical dashed line in the top panel denotes the position of the base of the convection zone at r/R⊙ = 0.7.
In the text

Fig. 4. Rossby mode frequency-shift time series due to variation in the solar rotation profile, plotted for m from 1 to 8 (top), and 9 to 15 (middle), calculated for each calendar year. The vertical grey dashed lines denote the positions of the maxima of solar cycles 23 and 24, to the nearest year. The vertical brown dashed line indicates the minimum between the two cycles. The smoothed sunspot number is plotted in the bottom panel.

In the text

Fig. 5. Top: predicted difference in the Rossby mode frequencies between the two solar maxima (grey curve, 2014 minus 2002) and between the maximum of cycle 24 and the preceding solar minimum (brown curve, 2014 minus 2009). The error bars are propagated errors from the scatter in the solar rotation measurements. Bottom: estimated errors in the Rossby mode frequencies due to the finite mode lifetimes as a function of observation time T. These errors can be used to determine whether the predicted frequency shifts from the top panel are detectable.

In the text

	Fig. B.1. Latitudinal dependences of the coriolis and advective terms that contribute to the kernel, and their sum, for m = 3.
In the text

	Fig. C.1. Top three panels: rotation rate residual multiplied by the sensitivity kernel and integrated over the radial coordinate. Bottom panel: difference in the solar rotation rate in the convection zone between the maxima of cycles 24 and 23.
In the text

	Fig. D.1. Frequency-shift time series from calculations in the β-plane (orange) plotted along with the results from the kernel (blue) for m = 9, 12, and 15.
In the text