Predicting frequency changes of global-scale solar Rossby modes due to solar cycle changes in internal rotation

Context. Large-scale equatorial Rossby modes have been observed on the Sun over the last two solar cycles. Aims. We investigate the impact of the time-varying zonal flows on the frequencies of Rossby modes. Methods. A first-order perturbation theory approach is used to obtain an expression for the expected shift in the mode frequencies due to perturbations in the internal rotation rate. Results. Using the time-varying rotation from helioseismic inversions we predict the changes in Rossby mode frequencies with azimuthal orders from m = 1 to m = 15 over the last two solar cycles. The peak-to-peak frequency change is less than 1 nHz for the m = 1 mode, grows with m, and reaches 25 nHz for m = 15. Conclusions. Given the observational uncertainties on mode frequencies due to the finite mode lifetimes, we find that the predicted frequency shifts are near the limit of detectability.


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 nearsurface 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 . 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 varia-tions 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 Section 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 Section 4. Finally, the discussion and conclusions are given in Section 5.

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, Ω(r, θ), is taken to be the unweighted mean rotation profile over the 23-year time interval. Variations about this mean profile over time are thus δΩ(r, θ, t) = Ω(r, θ, t) − Ω(r, θ). (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).

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).

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 with 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. d p = p + ξ · ∇p 0 ). 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: Using this relationship we can rewrite the operators that depend on u as Substituting these into Eq.

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 ξ * 0 on the left, integrating over V, and utilising the self-adjointness of the operators (cf. Lynden-Bell & Ostriker 1967) gives whereẑ is the unit vector along the rotation axis,ẑ =r cos θ − θ sin θ. 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 S (ξ) to contain the centrifugal and δL terms, we have with 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 S (ξ) 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 δΩ.

Rossby wave eigenfunction and resulting kernel
From here forwards we consider reference eigenfunctions ξ 0 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 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 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 ω 0 of mode m due to a change in rotation is given by where the sensitivity kernels K are given by For the displacement eigenfunction given in Eq. (17) and a real-valued L(θ), we have and ξ * 0 · D(ξ 0 ) =R 2 (r) 2mΩ 0 (r, θ) 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 we take a sectoral spherical harmonic ( = m) solution with no radial node: 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 This kernel does not depend on the uniform rotation rate Ω 0 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.

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, Ω 0 ) and the time-averaged differential rotation of the Sun (Ω(r, θ)) 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 .
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.

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.

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.

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 The kernel given should reduce to this result in the case of uniform rotation (Ω 0 = 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 Utilising the identity I p = π 0 (sin θ) p dθ = [(p − 1)/p]I p−2 , we obtain which is the expected result.

Appendix B: Latitudinal form of the kernel
The figures in Section 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 + 2)(cos θ) 2 (sin θ) 2m−1 ; therefore, they are equal and opposite at θ s = acos m m 2 + 2 .
(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.

Appendix C: Explaining the behaviour of δω for high and low m values
The top three panels of Fig. B.2 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). 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 SC 24. 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. B.2 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 .

β-plane
Here we perform a comparison with recent direct numerical calculations of frequencies of viscous Rossby modes in the equatorial β-plane  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 U(y, t) = R 0 r sin θ [Ω(r, θ, t) − Ω eq (r)] r 2(m+1) ρdr 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.