Issue 
A&A
Volume 673, May 2023



Article Number  A124  
Number of page(s)  19  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202245666  
Published online  23 May 2023 
Interaction of solar inertial modes with turbulent convection
A 2D model for the excitation of linearly stable modes
^{1}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: philidet@mps.mpg.de
^{2}
Institut für Astrophysik, GeorgAugustUniversität Göttingen, 37077 Göttingen, Germany
Received:
12
December
2022
Accepted:
28
March
2023
Context. Inertial modes have been observed on the Sun at low longitudinal wavenumbers. These modes probe the dynamics and structure of the solar convective zone down to the tachocline. While linear analysis allows the complex eigenfrequencies and eigenfunctions of these modes to be computed, it gives no information about their excitation nor about their amplitudes.
Aims. We tested the hypothesis that solar inertial modes are stochastically excited by the turbulent motions entailed by convection. Unlike the acoustic modes, which are excited by vertical turbulent motions, the inertial modes are excited by the radial vorticity of the turbulent field.
Methods. We have developed a theoretical formalism where the turbulent velocity fluctuations provide the mechanical work necessary to excite the modes. The modes are described by means of a 2D linear wave equation with a source term, under the β plane approximation. This wave equation restrained to a spherical surface is relevant for the quasitoroidal inertial modes that are observed on the Sun. Latitudinal differential rotation is included in the form of a parabolic profile that approximates the solar differential rotation at low and mid latitudes. The turbulent vorticity field underlying the source term is treated as an input to the model and is constrained by observations of the solar surface. The solution to the linear inhomogeneous wave equation is written in terms of a Green function, which is computed numerically.
Results. We obtain synthetic power spectra for the wave’s latitudinal velocity, longitudinal velocity, and radial vorticity, with azimuthal orders between 1 and 20. The synthetic power spectra contain the classical equatorial Rossby modes, as well as a rich spectrum of additional modes. The mode amplitudes are found to be of the same order of magnitude as observed on the Sun (∼1 m s^{−1}). There is a qualitative transition between low and high azimuthal orders: the power spectra for m ≲ 5 show modes that are clearly resolved in frequency space, while the power spectra for m ≳ 5 display regions of excess power that consist of many overlapping modes.
Conclusions. The general agreement between the predicted and observed inertial mode amplitudes supports the assumption of stochastic excitation by turbulent convection. Our work shows that the power spectra are not easily separable into individual modes, thus complicating the interpretation of the observations.
Key words: waves / turbulence / Sun: oscillations / Sun: interior / Sun: helioseismology
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model.
Open access funding provided by Max Planck Society.
1. Introduction
Multiple types of waves can propagate in the interior of a star. In the case of a nonrotating star, these modes are spheroidal; they are the pmodes (or acoustic modes), the fmodes (or surfacegravity modes), and the gmodes (or gravity modes). The first two have been observed on the Sun for a long time (Leighton et al. 1962; Deubner 1975) and are used to probe the equilibrium structure of the solar interior (see ChristensenDalsgaard 2002, for a review). Gravity modes, by contrast, are evanescent throughout the solar convective envelope, precluding us from using them as probes of the solar radiative interior.
The Sun, however, is a rotating star, and the inclusion of rotation entails the possibility of additional modes of oscillation. In a uniformly rotating star, theory predicts the existence of quasitoroidal modes, known as rmodes (Papaloizou & Pringle 1978), with the Coriolis force as the restoring force. In the rotating frame, they propagate in the retrograde direction and have frequencies comparable to the rotation rate. They are similar to the Rossby waves that are ubiquitous in the atmosphere of the Earth (Rossby 1939) and other planets (e.g. Allison 1990; SánchezLavega et al. 2014). Equatorial Rossby modes were recently observed on the Sun (Löptien et al. 2018; Liang et al. 2019) with a dispersion relation close to that of the theoretical sectoral (l = m) modes. A much richer spectrum of modes, collectively referred to as inertial modes, was subsequently reported by Gizon et al. (2021); in addition to the equatorial Rossby modes, these include highlatitude modes and critical latitude modes, allowed by latitudinal differential rotation. Furthermore, Hanson et al. (2022) report the observation of highfrequency waves of vorticity that are antisymmetric with respect to the equator.
These inertial modes allow the convective zone to be probed in a way that complements pmode helioseismology, especially when it comes to the superadiabatic temperature gradient and the turbulent viscosity. To do so, it is necessary to develop a theoretical understanding of the physics of inertial modes. Gizon et al. (2020) carried out a linear analysis of viscous modes of a parabolic shear flow in the β plane. This analysis, which applies to toroidal modes, was subsequently extended to viscous modes on a sphere by Fournier et al. (2022), thus allowing for a more realistic differential rotation profile and a treatment of the lowest azimuthal orders. A linear analysis of the 3D solar convection zone was carried out by Bekki et al. (2022b), who identified several of the modes reported by Gizon et al. (2021). Likewise, Triana et al. (2022) proposed an identification of the modes reported by Hanson et al. (2022).
Such linear analyses give information about the frequencies, the stability, and the eigenfunctions of the inertial modes, and allow for the identification of some of them in the observations. They also reveal that some of these modes can be linearly unstable, as a result of strong latitudinal differential rotation (Fournier et al. 2022) or a baroclinic instability due to a latitudinal entropy gradient (Bekki et al. 2022b). However, for latitudinal differential rotation profiles that are not too pronounced, most of the modes are predicted to be linearly stable, meaning that they are likely excited by turbulent convection. Understanding the excitation process of the linearly stable solar inertial modes would not only allow us to put stronger constraints on the dynamics of the solar convective zone, but would also help us predict which modes are expected to be visible and identifiable. This would go a long way towards helping us interpret the observational data at our disposal. In the present paper, we do not address the case of unstable inertial modes.
In the absence of a destabilising mechanism, the turbulent motions associated with solar convection provide an excitation mechanism, as in the case of solar pmodes (Lighthill 1967; Goldreich & Keeley 1977) or gravitoinertial waves (Mathis et al. 2014; Augustson et al. 2020). Nonlinear 3D simulations of the convection zone can help us assess the importance of this mechanism (Bekki et al. 2022a; Dikpati et al. 2022). In this paper, we follow a different approach and present a theoretical model for the turbulent stochastic excitation of purely toroidal inertial modes in 2D in order to test the hypothesis that this mechanism is indeed responsible for the amplitude level at which inertial modes are observed on the Sun. Because the analysis is done in 2D, it is only relevant for the predominantly toroidal modes, such as the equatorial Rossby modes and the other inertial modes that have been observed. Furthermore, we place ourselves in the equatorial βplane approximation, similarly to Gizon et al. (2020). We assume that the inertial modes are excited by turbulent emission, meaning that the nonlinear advection term in the momentum equation plays the role of a source term in the linear wave equation. This is also in accordance with the commonly accepted picture for pmodes (e.g. Samadi & Goupil 2001).
The paper is organised as follows. We present the stochastic excitation formalism in Sect. 2. The formalism requires two main ingredients, namely the Green function associated with the homogeneous wave equation (which is computed numerically, as described in Sect. 2.2), and the convective turbulent spectrum underlying the source term (which is treated as an input to the model and is the subject of Sect. 2.3). This model provides us with synthetic power spectra containing both the normal inertial modes of the system and the turbulent noise responsible for their generation. We focus on the expectation value of the power spectrum near the solar equator in Sect. 3.1 and discuss the latitude dependence of the power spectra in Sect. 3.2. Conclusions are drawn in Sect. 4.
2. Synthetic power spectra
2.1. Stochastic excitation by turbulence
We study the excitation of the vorticity modes observed on the Sun. These modes are quasitoroidal (characterised by their horizontal motions); we assume that the horizontal part of the wave equation can be decoupled from the radial part. In the following, we focus on the horizontal part and study the excitation of vorticity waves in a 2D shear flow mimicking the solar differential rotation (Gizon et al. 2020). We place ourselves in the equatorial β plane, where the 2D spherical coordinates λ and ϕ (denoting respectively the latitude and longitude) are transformed into
where R is the radius of the Sun. The velocity components on the sphere (v_{λ} and v_{ϕ}) can be approximated by the Cartesian components in the β plane (v_{x} and v_{y}). The equations of motion in the rotating frame and in the inviscid limit become
where ρ is the density, p is the gas pressure, f = βy = (2Ω_{eq}/R)y is the Coriolis parameter, and Ω_{eq} is the rotation rate at the equator. A linear inhomogeneous wave equation can be derived from Eqs. (3) and (4); the details are provided in Appendix A. To do so, the total flow velocity v is decomposed into a background zonal flow U ≡ U(y) e_{x}, which represents differential rotation, and a residual flow u, which contains both the waves and the turbulent noise. Working under the assumption that the residual flow is incompressible, we introduce the stream function Ψ such that
where e_{z} is the unit vector normal to the surface. This stream function is then decomposed into a contribution Ψ_{osc} for the oscillations and a contribution Ψ_{turb} for the convective turbulent noise. Linearising in terms of Ψ_{osc} while keeping all orders in Ψ_{turb}, we obtain (see Eq. (A.18))
where U″ is the second derivative of U(y), is the Laplacian operator, ν_{turb} is the turbulent viscosity, and the operator δ denotes a fluctuation around the horizontal average taken on scales larger than the turbulence scale, but shorter than the wavelength of the inertial modes (δq ≡ q − ⟨q⟩_{h} for any quantity q). The definition of this average requires a separation of scale, which is discussed in Sect. 2.3. The lefthand side of Eq. (6) governs the propagation of the inertial waves. As will be checked later, these modes are all linearly stable, because of the dissipative turbulent viscosity, ν_{turb}, which continuously pumps energy away from the modes, and because the shear induced by the differential rotation included in the model is not too strong. On the other hand, the righthand side of Eq. (6) acts as a source term that continuously injects energy into the modes. The equilibrium amplitude reached by the modes is the result of a balance between the damping and driving processes. Physically, the source term corresponds to the fluctuations of the divergence of the Reynolds stress tensor around its statistical average, and is stochastic by nature, because of the highly turbulent nature of the flow in the solar convective zone. The mechanism is similar to the traditionally accepted picture of solarlike pmode excitation (Goldreich & Keeley 1977); a key difference, however, is that whereas pmodes are mainly excited by the vertical turbulent motions, the quasitoroidal inertial modes considered here are much more sensitive to the horizontal, vorticity component of turbulence.
The inhomogeneous wave equation can be written in terms of Fourier modes exp(jωt − jk_{x}x), where j denotes the imaginary unit. We used the following convention for the Fourier transform,
where T_{obs} and X_{obs} are respectively the t and x windows over which the integral defining the Fourier transform is computed. In the Fourier domain, Eq. (6) becomes
The notation refers to the (x, t) Fourier transform of the righthand side of Eq. (6), and ℒ is the linear propagation operator, given by
where . The linear operator ℒ and the source depend on y, the angular frequency ω, and the longitudinal wavenumber k_{x}.
The solution to Eq. (8) can be expressed in terms of the Green function G(y, y_{s}), which is defined as the solution of the following differential equation (and with the same boundary conditions as the full linear problem),
where δ is the Dirac distribution, so that
This solution leads to other physical quantities, such as the azimuthal velocity , the latitudinal velocity and the radial vorticity , for which we obtain
Then we obtain the expectation of the power spectral density by forming the modulus squared of Eqs. (12)–(14) and taking the ensemble average. We assume that the spatial scale of the Green function and of the source are well separated – the validity of this assumption will be checked in Sect. 2.3. We find
where y_{obs} is the latitudinal coordinate at which the power spectrum is evaluated, and ⟨.⟩ denotes an ensemble average. The function ℐ(y_{s}) denotes the source covariance, and is defined by
We recall that the source term depends on ω and k_{x}, so that ℐ does too. The computation of the source covariance is detailed in Appendix B. We assume, in particular, that the source term is homogeneous, in the sense that its statistical properties do not depend on latitude; as such, the source autocorrelation spectrum no longer depends on y_{s}. Eventually, we find (see Eq. (B.35))
We note that the integrals span all frequencies and wavenumbers, positive and negative alike. The function ℰ_{Ψ} represents the turbulent stream function spectrum, and is given by (see Eq. (B.31))
We assumed that the turbulence is stationary and homogeneous, such that the turbulent spectrum depends on neither absolute time, T, nor on absolute space, X.
Equations (15)–(17) correspond to the contribution of the inertial modes to the velocity and vorticity power spectra. These spectra also contain a contribution from the turbulent noise, which can be expressed solely as a function of the turbulent stream function spectrum, ℰ_{Ψ}. Since this spectrum is the same quantity that appears in Eq. (19), the contribution of the inertial modes and of the turbulent noise can be modelled simultaneously. We find
Forming the sum of the inertial mode contributions (i.e. Eqs. (15)–(17)) and the noise contributions (i.e. Eqs. (21)–(23)) yields our synthetic power spectrum model, in terms of azimuthal velocity, latitudinal velocity, and radial vorticity, respectively. Only two ingredients are needed to quantify these expressions, namely (i) the Green function G(y, y_{s}) associated with the linear operator ℒ, and (ii) the turbulent stream function spectrum ℰ_{Ψ}.
2.2. The Green function
The angular frequency, ω, and azimuthal wavenumber, k_{x}, being fixed, the linear operator, ℒ, depends on (i) the differential rotation profile, U(y), (ii) the Coriolis parameter, β = 2Ω_{eq}/R, and (iii) the turbulent viscosity, ν_{turb}. Concerning the differential rotation profile, as a first step, we approximated it by a parabolic profile,
where we have implicitly placed ourselves in a frame of reference rotating at the solar equatorial rotation rate Ω_{eq}/(2π) = 453.1 nHz, and we have introduced the nondimensionalised latitudinal coordinate ξ. We chose the same value m s^{−1} as Gizon et al. (2020), which approximates the solar differential rotation at low latitudes. With this value of Ω_{eq}, the Coriolis parameter becomes β = 8.18 × 10^{−15} m^{−1}s^{−1}. Finally, the turbulent viscosity, ν_{turb}, is specified through the turbulent Reynolds number,
which we leave as a free parameter.
Once all these parameters are fixed, in order to numerically compute the Green function, we expand Eq. (10) on the basis formed by the Chebyshev polynomials of the first kind. Those are defined, for every positive integer n, by
which reduces to a polynomial expression after some algebraic manipulations. These polynomials are orthogonal to each other with respect to the following inner product,
in the sense that
where δ_{ij} is the Kronecker delta and c_{i} = 1 + δ_{i0}. We denote the column vector containing the decomposition of the Green function on the Chebyshev basis as 𝒢(ξ_{s}), so that
where we also introduced the nondimensionalised source position ξ_{s} ≡ y_{s}/R. The vector 𝒢(y_{s}) is the solution of the following linear system:
where the column vector on the righthand side comprises the projections of the Dirac distribution on the Chebyshev polynomials,
and the matrix ℳ on the lefthand side is defined by
We note that the factor 2/(πc_{i}) in both Eqs. (31) and (32) stems from the fact that the set of Chebyshev polynomials is orthogonal but not orthonormal.
The Chebyshev polynomials of the first kind prove particularly well suited for solving Eq. (9), because the ℳ_{ij} take a conveniently simple form on that basis, as was shown for example by Orszag (1971). We detail the derivation of these matrix coefficients in Appendix C. Naturally, while the matrix ℳ is of infinite dimension, it is necessary to crop it to a finite size, in order for the numerical computations to be carried out. We found that truncating the Chebyshev expansion at N = 500 was a good compromise between a reasonable computation time and an accurate representation of the Green function. We note that while the righthand side of Eq. (30) depends on the source position ξ_{s}, this is not the case of the matrix ℳ, meaning that for any given angular frequency ω and azimuthal wavenumber k_{x}, only one single N × N matrix inversion is necessary to find the Green function for all possible source positions.
Solving Eq. (30) for 𝒢(ξ_{s}) yields one Green function, corresponding to arbitrary and completely uncontrolled boundary conditions. It is therefore also necessary to enforce the correct boundary conditions:
The Chebyshev polynomials of the first kind verify T_{i}(1) = 1, T_{i}(−1) = (−1)^{i}, and , so that enforcing these boundary conditions amounts to ensuring that the solution 𝒢(ξ_{s}) of Eq. (30) verifies
These conditions can be enforced in Eq. (30) by replacing the last four lines of the matrix ℳ by
and by replacing the last four components of ℬ(ξ_{s}) by zero. This is perfectly equivalent to the τmethod applied, for example, by Orszag (1971). Stated more intuitively, this means that the highfrequency behaviour of the solution is now controlled not by the dynamical behaviour of the system, but by the mechanical constraints imposed on the boundaries. Solving this modified linear system for 𝒢(ξ_{s}) now yields the correct Green function, with the appropriate boundary conditions.
2.3. The turbulent stream function spectrum
The turbulent spectrum ℰ_{Ψ}, defined by Eq. (20), is written in terms of the turbulent fluctuation of the stream function Ψ_{turb}. On the other hand, a more common definition for the turbulent spectrum relies on the turbulent velocity u_{turb} rather than the stream function
We note that, whether it be in Eqs. (20) or (39), the angular frequency ω is not restricted to be positive, but can be of any sign. These two spectra are related through
The turbulent velocity spectrum is usually expressed as (e.g. Lesieur 2008)
If the turbulence is incompressible, then the quantity E(k, ω) is rigorously isotropic (in the sense that it does not depend on the direction of the wavevector k, but only on its modulus), and the directional information is entirely contained in the projection operator that follows. This function E(k, ω) is what is commonly referred to as the turbulent spectrum. Plugging Eq. (41) into Eq. (40), we simply find
so that knowing ℰ_{Ψ} is perfectly equivalent to knowing E.
In the following we consider that the turbulent spectrum E(k, ω) can be separated into a spatial part and a temporal part. The argument was proposed by Stein (1967) and later consistently used in the context of wave excitation by turbulent emission (see for instance Musielak et al. 1994, or, in the context of solarlike pmodes, Samadi & Goupil 2001; Chaplin et al. 2005; see also Tennekes & Lumley 1978 Chapter 8 for a textbook on the subject). The idea is to introduce the spatial spectrum
and then define the temporal part of the spectrum as
We consider that, as in the case of incompressible turbulence, the spatial spectrum E(k) is isotropic, and can therefore be rewritten as a kinetic energy per unit k instead of per unit k,
Concerning the temporal spectrum, following dimensional arguments, Stein (1967) argued that the temporal evolution of a turbulent eddy of size λ_{k} ≡ 2π/k should contain frequencies around ν_{k} ∼ u_{k}/λ_{k} (or, equivalently, angular frequencies around ω_{k} ∼ ku_{k}), where u_{k} is the typical velocity associated with these eddies. Guided by the work of Kraichnan (1965; see his Eq. (9.6) and the discussion in the paragraph above it), Stein (1967) wrote
where we have introduced the reduced frequency , and χ is now a function of frequency that does not depend on the eddy size k. Stein (1967) proved that the typical velocity u_{k} associated with these eddies is determined by the spatial spectrum through
All in all, we can write
where u_{k} is given by Eq. (47). Describing the whole turbulent spectrum requires only two ingredients, namely (i) the ωindependent spatial spectrum E_{iso} and (ii) the kindependent temporal spectrum χ. Both can be extracted from solar observations. We used measurements of the vertical vorticity deduced from granulation tracking by Langfellner et al. (2015; their Fig. 3, where the units where reassessed). We use the observations near the solar equator and adopt the same vorticity spectrum at all latitudes, for the sake of simplicity.
The spatial spectrum E_{iso}(k) is shown in the left panel of Fig. 1, and can clearly be described by three distinct power laws in three separate wavenumber regimes:
Fig. 1. Spatial turbulent spectrum and turbulent frequency. Left: Spatial turbulent spectrum, E_{iso}(k), extracted from solar observations (solid blue line). The spectrum comprises roughly three sections, each with a different power law in k^{α}. The dashed orange line shows the best fit of the data to a threepiecewise power law. The exponents and the boundaries between the three regimes are indicated on the plot. Right: Turbulent frequency ω_{k} ≡ ku_{k} as a function of wavenumber, k, where u_{k} is the typical velocity of the turbulent eddies of inverse size k, given by Eq. (47). 
where the choice of k_{ref} is completely arbitrary and can be absorbed in the factors C_{i}. The spectrum is therefore parameterised by (i) the three exponents α_{1, 2, 3}, (ii) the two wavenumber cutoffs k_{1} and k_{2}, and (iii) the total turbulent kinetic energy per unit mass ∫E_{iso}(k)dk. The best fit to the observational data is also shown in the left panel of Fig. 1. We find exponents α_{1} = 0.27, α_{2} = −1.88 and α_{3} = −7.60; in particular, the slope α_{2} of the middle section is quite close to the −5/3 power law theoretically predicted for the inertial range under the Kolmogorov hypotheses. As for the wavenumber cutoffs – also indicated in the left panel of Fig. 1 – we find k_{1}R = 140 and k_{2}R = 457. These values are significantly larger than the wavenumber associated with the solar inertial modes, thus validating the assumption made earlier concerning the scale separation.
We then use this spatial spectrum to compute the turbulent frequencies ω_{k} = ku_{k} as a function of k, where u_{k} is given by Eq. (47). Those are shown in the right panel of Fig. 1. It can be seen that the typical frequencies associated with the turbulent motions are of the order of a few tens of μHz. The frequencies of the inertial modes, by contrast, are typically much lower (of the order of ∼100 nHz), which means that there is a timescale separation between the inertial modes and the turbulence, similar to the length scale separation already mentioned above.
Finally, we checked the assumption made in writing Eq. (46) (i.e. the fact that the quantity ω_{k}χ_{k}, when plotted against the reduced frequency , collapses onto a unique, slowly varying curve, independent of k). The result, shown in Fig. 2, seems to indicate that this is indeed the case. What is more, it is also possible to determine the analytical function that best describes this curve. In the context of pmode excitation, traditional models for the turbulent temporal spectrum usually assume either a Gaussian or a Lorentzian function (e.g. Goldreich & Keeley 1977; Balmforth 1992; Samadi et al. 2007; Belkacem et al. 2010). We tried the two following models:
Fig. 2. Temporal spectrum , as a function of the reduced frequency . The data points (in orange) have been binned according to the value of (which depends on both ω and k), and only the mean over each bin is shown. These data points collapse onto a unique, slowly varying curve, which was adjusted alternatively with a Gaussian function (Eq. (50), solid black line) and a Lorentzian function (Eq. (51), solid green line). The latter is clearly the best fit, obtained for an amplitude A = 1.85 and a standard deviation σ = 0.62. 
and
In each case, the dimensionless factor σ is introduced to account for the uncertainty on the relation ω_{k} = ku_{k}. It constitutes a free parameter in each fit, but is expected to be of order unity. This parameter is akin to the parameter λ in the pmode excitation formalism of Samadi & Goupil (2001), which they also left free for the same reason. Figure 2 shows the result of each fitting procedure. The Lorentzian function clearly yields the best agreement, with an amplitude A = 1.85 and standard deviation σ = 0.62; we therefore adopted this prescription. We note that this value of σ is consistent with the values of λ found by Samadi & Goupil (2001), who constrained it using the observed excitation rate of solar acoustic modes (see their Table 1).
3. Results
3.1. Equatorial power spectrum
We first considered the synthetic spectrum as it would be observed at the solar equator, and written in terms of latitudinal velocity u_{y}. This is obtained by setting y_{obs} = 0 in Eqs. (16) and (22), which yield respectively the contribution of the inertial modes and of the turbulent noise to the overall spectrum. For the moment, however, we only consider the inertial mode contribution, and we only add the turbulent noise contribution later on. The power spectrum thus obtained corresponds to a spectral density per unit longitudinal wavevector k_{x} and angular frequency ω. It is then straightforward to transform it into a power spectral density per unit ω only, for any given azimuthal order m ≡ k_{x}R, by dividing the power spectrum by the radius, R, of the spherical domain. Naturally, our model can be applied to any real value of m, because in the scope of the equatorial βplane approximation we did not make any assumption regarding the ϕperiodicity of the velocity field. However, for the sake of comparison with observations, only integer values of m are relevant.
3.1.1. Individual mode contributions
Our model allows us to decompose the total spectrum into the individual contributions of each inertial mode. To that effect, we first needed to compute the complex eigenfrequencies and eigenfunctions of the linear homogeneous system ℒΨ = 0, where we recall that ℒ is given by Eq. (9). The eigenfrequencies ω_{n} and eigenfunctions Ψ_{n} are the solution of the following generalised boundary eigenvalue problem
where
The enforced boundary conditions are the same as those of the inhomogeneous problem (see Sect. 2.2). We solve the system numerically, as described in Sect. 2.2, by projecting the generalised boundary eigenvalue problem on the basis formed by the Chebyshev polynomials of the first kind. The discrete eigenspectrum is showcased in the complex plane, for azimuthal orders m = 1, 5 and 10 in the top panels of Figs. 3–5, where each point represents one mode. As previously mentioned, all the eigenfrequencies have a negative real part, meaning that the modes propagate in the retrograde direction, and a negative imaginary part, meaning that they are all stable, and none of them is exponentially growing. In the m = 5 and m = 10 cases, one can clearly recognise the three classical mode branches associated with Poiseuille flows (Mack 1976). Those are indicated in the top panel of Fig. 5. While the two upper branches (i.e. the A branch and the P branch) comprise a finite number of modes, all of which are shown, the vertical branch (i.e. the S branch) is infinite, and is truncated in the plots. Furthermore, in those same cases, an additional mode sticks out (it is marked as R mode in the top panel of Fig. 5), characterised by a comparatively longer lifetime (i.e. an imaginary frequency closer to zero). It is entailed by the presence of global rotation in the system (through the β factor in Eq. (9)), and can be thought of as the equivalent of an equatorial Rossby mode in Cartesian geometry. Overall, the structure of the eigenspectrum is the same as the one presented in Gizon et al. (2020), as the homogeneous part of our wave equation is the same as theirs; for more details on the matter, we therefore refer the reader to their work.
Fig. 3. Latitudinal velocity spectrum. Top: discrete inertial mode eigenspectrum for m = 3, shown in the complex plane. Each point correspond to one mode: all eigenfrequencies have a negative real part (meaning the modes are retrograde) and a negative imaginary part (meaning they are stable). The coloured dots mark the modes whose contribution to the u_{y} spectrum is most prominent. The vertical dashed lines indicate the central frequencies of each resonant peak in the u_{y} power spectrum (see bottom panel). Middle: equatorial u_{y} spectrum, obtained from Eq. (16) for azimuthal order m = 3. The solid black line shows the total spectrum, while each coloured dashed line corresponds to the individual contribution of the normal eigenmodes of the system, as described in the main text. The colours of the dashed lines match the colour scheme of the top panel. The red vertical dashed lines show the local maxima of the total spectrum. Bottom: same as the middle panel, but the vertical scale is linear. 
Fig. 5. Same as Fig. 3, but for m = 10. The three classical branches of Poiseuille flows are indicated on the plot, in addition to the R mode. 
We also need to compute the complex eigenfrequencies and eigenfunctions of the associated adjoint system ℒ^{†}Ψ = 0, where the adjoint linear operator ℒ^{†} is defined in such a way that for any function f and g satisfying the boundary conditions of the problem, we have
Performing integration by parts we find that the adjoint eigenfrequencies and eigenfunctions are the solution of the following generalised boundary eigenvalue problem (Salwen & Grosch 1981)
where
It is easy to show that . Furthermore, the eigenfunctions and adjoint eigenfunctions form a biorthogonal set, in the sense that
The eigenfunctions are normalised in such a way that the proportionality coefficient in Eq. (57) is unity. This biorthonormality relation allows us to project the Green function G(ξ, ξ_{s}) onto each of the modes alternatively, and therefore to compute the u_{y} spectrum associated with each mode separately.
The results are shown in the bottom panels of Figs. 3–5, for azimuthal orders m = 3, 5, and 10, respectively. The total spectrum is represented by the solid black line, while the coloured dashed lines represent the individual contributions of the most prominent modes in the discrete spectrum. We note that for the most part, and as is expected of a resonant mode whose driving source has a broadband frequency dependence, the spectrum contribution of each individual mode takes the form of a Lorentzian profile. For low values of m (more specifically for m ≤ 5), the spectrum is clearly dominated by two main peaks, and the frequencies of these peaks makes them clearly identifiable, as they fall very close to the real part of discrete eigenmodes of the homogeneous system, as shown in the top panels.
As m increases, the peaks become wider (because the imaginary part of the eigenfrequencies increase in modulus), and for m ≥ 6 the spectrum becomes dominated by one mode only. From the eigenspectrum shown in the top panels it can be seen that the dominant mode corresponds to the equatorial Rossby mode. Interestingly, we find that several modes have an amplitude that, if taken individually, should lead to a visible peak in the spectrum, but remain invisible in the total spectrum, where all modes are accounted for at once. This seems to suggest that the reason these modes cannot be extracted from the spectrum is not the inefficiency of the excitation process but rather the fact that they form a mutually destructive interference pattern with each other. This is possible because the modes all share the same driving source. This interference phenomenon can only occur between modes that share the same azimuthal order m, because the problem is axisymmetric, and therefore different m are completely decoupled. Furthermore, the modes must share similar eigenfrequencies (more specifically, the frequency difference must be of the order or smaller than the inverse of their lifetime).
3.1.2. Frequencies, amplitudes, and linewidths
For the resonant peaks that are sufficiently separated from each other in the synthetic spectrum, it is possible to directly infer, from the model, not only their frequency, but also their amplitude and linewidths. We define the angular frequency ω_{0} of each peak as the location of their local maximum, and their full linewidth at half maximum Γ as the angular frequency range where the power spectral density is above half the height of the peak. The amplitude is defined as
where the boundaries ω_{a} and ω_{b} should be chosen to enclose most of the peak, without overlap with the other peaks. In practice, we chose ω_{a, b} = ω_{0} ∓ Γ/2. In the case of a Lorentzian profile, this range encloses exactly half the energy of the mode, so we only have to apply a factor of to Eq. (58) to find the total mode amplitude. We note that, because of the typical linewidths of the modes, the upper boundary ω_{b} can very well be positive, despite the fact that the central frequency, ω_{0}, is systematically negative.
The frequencies of the modes are showcased, as a function of m, in Fig. 6. The coloured diamonds are identical on each panel: blue and green symbols represent each of the upper mode branches in the spectrum (see Sect. 3.1.1 for a description of these categories), while red symbols represent equatorial Rossby modes. On the background of each panel is superimposed an image of the spectrum in terms of different physical quantities (latitudinal velocity on top, azimuthal velocity in the middle, and vorticity at the bottom), in the mω plane, where each vertical slice is normalised so that the maximum is unity. One can distinguish the same transition, already mentioned above, between low azimuthal orders, for which several clearly identifiable resonant peaks can be resolved in the spectrum, and high azimuthal orders, where the distinction is no longer as clear. The theoretical Rossby mode dispersion relation is also plotted in each panel: in Cartesian coordinates, and in the presence of a background azimuthal jet U, it is given by
Fig. 6. Synthetic equatorial spectrum in the mω plane, in terms of u_{y} (top), u_{x} (middle), and ζ (bottom). Each vertical slice is normalised separately such that the maximum is unity. The diamonds show the real part of the eigenfrequencies of the linear homogeneous problem, computed as described in Sect. 3.1.1. The colour code refers to the mode categories: the blue diamonds represent the A branch, the green diamonds represent the P branch, and the red diamonds represent the Rossby modes (see Sect. 3.1.1 for a description of these branches). The solid red line shows the theoretical dispersion relation for sectoral Rossby modes in Cartesian coordinates (see Eq. (59)). 
The high m frequencies match the theoretical dispersion relation quite well. The agreement, however, is not as close for lower values of m. We also note that while the lowm equatorial Rossby modes are quite distinguishable in the u_{y} spectrum or the vorticity equatorial spectrum, they do not show in the equatorial u_{x} spectrum, due to the fact that their u_{x} eigenfunction has a node at the equator.
In the left panel of Fig. 7, we plot the linewidths Γ of the synthetic Rossby modes, as a function of m, for several values of the turbulent Reynolds number Re_{turb}. We also show, on the same plot, the observed linewidths reported by Liang et al. (2019) for solar Rossby modes at the equator. While the order of magnitude is consistent with the observations, the uncertainties associated with them do not permit us to discriminate between the different models. Of particular interest is the fact that we recover, for high m, and especially for the Re_{turb} = 300 case, the same m^{2} dependence that can be inferred from the observations. This law is consistent with the theoretical Rossby mode dispersion relation (also shown on the plot), whose imaginary part yields
Fig. 7. Equatorial Rossby mode parameters. Left: full width at half maximum of the resonant peaks that could be identified as equatorial Rossby modes, as a function of m. The diamonds show the linewidths measured in the u_{y} equatorial synthetic power spectrum, and the coloured solid lines represent the theoretical Rossby mode linewidth, obtained from the classical dispersion relation (see Eq. (60)). The colour code refers to the value of the turbulent Reynolds number: Re_{turb} = 300 (red), 700 (blue), and 1000 (green). The black line shows the mode linewidths inferred from solar observations at the equator in the latitudinal velocity spectrum, as reported by Liang et al. (2019). Error bars from the fitting procedure reported by the authors are also shown. Right: Rossby mode amplitude (coloured solid lines) in the u_{y} equatorial synthetic power spectra, as a function of azimuthal order m, defined as described in the text (see Eq. (58)). The colour code is identical to the one in the left panel. 
We note that the value of the turbulent Reynolds number that seems to give the best agreement between the theoretical dispersion relation and the observations is Re_{turb} = 300, which corresponds to a turbulent viscosity ν_{turb} = 570 km^{2} s^{−1}. This value is in accordance with the surface value inferred by Gizon et al. (2020). However, it is significantly larger than the upper limit of 100 km^{2} s^{−1} inferred by Gizon et al. (2021), which was obtained under the assumption that the turbulent viscosity is constant over the entire convection zone.
Finally, we report the amplitude of the synthetic Rossby modes in the right panel of Fig. 7, in terms of latitudinal velocity u_{y}, for several models corresponding to different values of Re_{turb}. We also superimposed the observed amplitude reported by Liang et al. (2019). The agreement that we find between the observations and the amplitudes yielded by our synthetic spectrum model, especially for Re_{turb} = 700 and 1000, is consistent with our initial hypothesis that the inertial modes observed on the Sun are stochastically excited by turbulent convection. More specifically, we find that the amplitude of the equatorial Rossby modes initially increases with m, and then reaches a plateau where it remains fairly independent of m. The low amplitude of the lowm equatorial Rossby modes explains why the m = 1 or 2 equatorial Rossby modes elude observation on the surface of the Sun. These results also show that increasing the turbulent viscosity (i.e. decreasing Re_{turb}) causes both an earlier start of the plateau and a lower value thereof. For instance, a value Re_{turb} = 300 causes the equatorial Rossby modes to reach an amplitude of 1.5 m s^{−1} after m ∼ 10, while a value Re_{turb} = 1000 causes them to reach an amplitude of 2.5 m s^{−1} after m ∼ 15.
The fact that the amplitudes predicted by our simple 2D model compare well with the solar observations deserves some discussion. This good agreement suggests that neither the radial eigenfunctions nor the source function depend strongly on radius r. Regarding the eigenfunctions, linear eigenmode computations in 3D by Bekki et al. (2022a) indicate that the dependence on r is less pronounced for the smallest m than for the larger m. For example, the velocity eigenfunction of the equatorial Rossby mode scales like the slowly varying function r^{m} for m = 3 and 4. Regarding the source of excitation, a turbulent vorticity spectrum that peaks at a scale k_{0} can only drive modes with azimuthal wavenumbers k_{x} = m/R < 2k_{0} (see Fig. 8). Since the radial vorticity for these large scales does not vary fast with depth (see e.g. Miesch et al. 2008, their Fig. 13e), it is perhaps not too surprising that our 2D model predicts the amplitudes of the lowm inertial modes correctly (in order of magnitude). Future work should nevertheless account for the radial dependence of the source properties, as well as the mode eigenfunctions and viscous damping.
Fig. 8. Geometrical argument used to identify the modes, k_{mode}, that can be excited by a single scale of turbulence, k_{0}. The integrand in Eq. (19) peaks at k′ such that both k′ + k_{mode}/2 and k′ − k_{mode}/2 are equal to k_{0}. The driving can therefore only be efficient if the two black circles in the figure intersect, i.e. if k_{mode} < 2k_{0}. 
3.1.3. Visibility of the modes
So far, we have only investigated the component of the spectrum caused by the resonant inertial modes, and we have excluded the turbulent noise component from our analysis. However, the inclusion of this noise component is necessary in order to assess whether the signal of the modes is visible above the noise level. To do this, we simply add the u_{y} noise component given by Eq. (22) to the inertial mode component given by Eq. (16). Furthermore, we would like to directly compare our results to the observed equatorial u_{θ} spectrum reported by Liang et al. (2019). The authors measured southnorth helioseismic travel times along the solar equator using datasets from the Michelson Doppler Imager (MDI) on board the Solar and Heliospheric Observatory (SoHO) and the Helioseismic and Magnetic Imager (HMI) on the Solar Dynamics Observatory (SDO). Their procedure was carried out for several latitudes uniformly distributed within the range λ = ±15°. Therefore, to mimic their procedure, we average our synthetic u_{y} power spectrum over this latitude range (that is, we average between y_{obs} = ±sin(15° ) instead of simply taking y_{obs} = 0).
The comparison is shown in Fig. 9 for azimuthal orders m = 3, 5 and 10. It can be seen that, by construction, the turbulent noise level in our synthetic power spectrum matches the noise level in the corresponding solar observations. We also point out, before diving further in the question of mode visibility, that the central frequencies in our model do not correspond to the observed mode frequencies, especially for higher values of m. This is to be expected, because the homogeneous part of the wave equation governing the frequencies of the inertial modes is somewhat simplified in our model (in particular the fact that we reduced the system to a 2D model, or the fact that we did not account for the spherical geometry of the domain). The focus of this study is not on the mode frequencies, so that these discrepancies do not constitute an issue.
Fig. 9. Latitudinal velocity spectrum, estimated at the solar equator, for azimuthal orders m = 3 (top), 5 (middle), and 10 (bottom), and for a turbulent Reynolds number Re_{turb} = 300. The solid red line represents our model, including the contribution both from the inertial modes and from the turbulent noise. The thin blue line shows the solar observations reported by Liang et al. (2019), and the thicker blue line shows the best Lorentzian fit to their data, as reported in their paper. 
We find that for m = 1 and m = 2, the amplitude of the modes is lower than the turbulent noise level. On this, our model prediction agrees with observational evidence; Liang et al. (2019) had indeed reported that no significant peak could be detected for such low values of the azimuthal order.
For m = 3 to 5, we recall that several peaks could be distinguished in our synthetic spectrum. However, as illustrated for m = 3 in the top panel of Fig. 9, it turns out that only one of these resonant peaks has enough amplitude to rise above the turbulent noise level. This is also in agreement with the observations.
For m ≥ 5, we recall that only one peak is predicted to dominate in the equatorial spectrum, corresponding to the equatorial Rossby mode. At the start of this m ≥ 5 regime, this mode is clearly visible above the noise level (as illustrated for m = 10 in the bottom panel of Fig. 9). However, we find that the visibility of the mode becomes more and more questionable as m increases. This is in agreement with the solar observations reported by Liang et al. (2019), who could not detect any significant peak in the equatorial spectrum for m ≥ 16.
Our model allows some light to be shed on this highm visibility issue. We show in Sect. 3.1.2 that, in this highm regime, the amplitudes of the equatorial Rossby modes are fairly independent of m, while their linewidths increase as m^{2}. Because the peaks become wider and wider for the same total power, their spectral height decreases as 1/m^{2}, thus explaining why they end up below the turbulent noise level after some point. We point out that, because the modes are stochastically excited by turbulent convection, their signaltonoise ratio is inherent to the excitation process, and not contingent on the observational setup. It would therefore not be possible to improve their visibility through a longer observation time period.
3.2. Power spectrum in the frequency–latitude plane
Throughout Sect. 3.1, we focused our analysis on equatorial power spectra. However, our model also gives us access to the power spectra at any given latitude. This is illustrated in Fig. 10, for the u_{y} power spectrum and for azimuthal orders m = 3, 5, and 10. The solid black curve in each panel corresponds to the critical layer where the azimuthal velocity, U, stemming from differential rotation exactly matches the azimuthal phase velocity of the inertial waves, that is, the curve y_{c}(ω) given by the following implicit relation:
Fig. 10. Power spectrum of the latitudinal velocity, u_{y}, as a function of frequency (horizontal axis) and latitude (vertical axis). Each panel corresponds to a different azimuthal order: m = 3 (top), m = 5 (middle), and m = 10 (bottom). The turbulent Reynolds number is set to Re_{turb} = 300. The solid black curve shows the critical latitudes, where the differential rotation exactly matches the phase velocity of the inertial waves, and is defined by the implicit relation Eq. (61). The dashed vertical lines show the real part of the eigenfrequencies of the homogeneous problem, with the same colour code as in Fig. 6. 
where U is given by Eq. (24).
The shape of the power spectrum clearly transitions from a lowm regime, dominated by a few, clearly identifiable and distinguishable resonant modes, to a highm regime, where the power is concentrated in a characteristic crescentshaped region along the critical latitude (or, more precisely, just below the critical latitude). The transition between the two regimes occurs at m ∼ 5 and corresponds to the same transition that we already mentioned for the equatorial spectra in Sect. 3.1. This seems to indicate that the detection of inertial modes in observational data is much more delicate for high m than for low m: the overlapping of the excess power regions associated with each mode in the frequency – latitude plane is indeed likely to make the interpretation of the observed spectra considerably more complex.
Similarly, Fig. 11 shows the u_{x} power spectrum in the same frequency – latitude plane. While the same lowm to highm regime transition can be observed, there is still a number of qualitative differences with the u_{y} power spectrum. The main difference is that, for high values of m, the region where most of the power is located passes across the critical latitude as frequency becomes more and more negative, while this region was kept confined at lower latitudes in the u_{y} spectrum. This is simply a symptom of the qualitative difference between the u_{y} and u_{x} eigenfunctions. The effect grows stronger as m increases, and indicates that, depending on the observable used to detect inertial mode, the region just above the critical latitude may be of similar interest than the region below.
We show similar results for the vorticity power spectrum in Fig. 12. Again, the main difference comes from the latitudinal structure of the eigenfunctions, which contains significantly more power at the poles than the u_{y} or u_{x} eigenfunctions. However, the lowm to highm transition described above can still be observed in the vorticity power spectra.
4. Conclusion
We have designed a model for the stochastic excitation of linearly stable, quasitoroidal solar inertial modes by turbulent convection. In order to do so, we adopted a simplified 2D framework, where inertial modes are described in an equatorial β plane close to the surface of the Sun. We included latitudinal differential rotation in the form of a parabolic, Poiseuillelike profile, with values chosen to best approximate the solar differential rotation at low latitudes. Using this model, we successfully reproduce the observed amplitude of the linearly stable, low and midlatitude inertial modes, with latitudinal velocities ranging between ∼0.1 and ∼1.5 m.s^{−1}, similar to those reported, for instance, by Liang et al. (2019). The amplitude of the linearly stable inertial modes observed in the equatorial region of the Sun is therefore consistent with a stochastic excitation by turbulent convection. However, we did not treat the case of the unstable, highlatitude inertial modes.
We also show that the power spectra in the frequency–latitude plane have a very different qualitative behaviour depending on whether m ≲ 5 or m ≳ 5. In the lowm regime, the spectra are dominated by nonoverlapping, clearly identifiable and distinguishable resonant modes. On the other hand, in the highm regime, the line profile of the modes is much wider and so the excess power region associated with each mode overlap, thus forming a single crescentshaped excess power region along the critical latitude. In this regime, it is much harder to distinguish between individual modes. This has important implications for the identification of inertial modes in solar data, as this seems to indicate that the interpretation of the observed spectra becomes increasingly complex as m increases.
In the equatorial spectra, the predicted amplitudes of the modes are such that they are only visible above the turbulent noise level for m ≥ 3, in accordance with solar observations. Between m = 3 and 5, the equatorial spectra feature several dominant peaks, whose line profiles do not overlap. By contrast, for m > 5, the power spectrum is dominated by the inertial modes with the longest lifetime, and which correspond to the Cartesian equivalent of the classical equatorial Rossby modes. We predict that the amplitude of these equatorial Rossby modes will increase with m until m ∼ 10, after which the amplitudes reach a plateau and become fairly independent of azimuthal order. By contrast, their linewidth increases with m such that their spectral height decreases in such a way that the mode stops being visible above the turbulent noise level for m ∼ 16, in accordance with observations. Interestingly, we find that some modes incur mutually destructive interference, to such a degree that their overall amplitude is negligible even when their individual amplitudes should make them visible. This is possible because all the modes share the same source of excitation.
Additionally, we find that the theoretically predicted full linewidths at half maximum of the equatorial Rossby modes agree reasonably well with the observed linewidths, provided we choose the turbulent viscosity to be ∼570 km s^{−1}. This confirms constraints previously obtained in the literature for the convectioninduced turbulent viscosity. We also show that the linewidths of the equatorial Rossby modes vary as m^{2}, which can be predicted through the classical Rossby wave dispersion relation (see Eq. (60)). Because this simple square law primarily depends on the turbulent viscosity, the value of the turbulent viscosity can be constrained throughout the solar convection zone, even potentially through the use of inversion techniques.
While the equatorial βplane approximation constitutes a drastic simplification for low azimuthal orders, m, it is not so much the case for higher m. Therefore, we do not expect our results to be overly affected should this approximation be lifted, and should the derivations be carried out in spherical geometry. This would nevertheless warrant further investigation. The suppression of the radial coordinate in the problem, by contrast, undoubtedly constitutes a more important approximation. In particular, the use of a 3D model, rather than 2D, would increase the density of modes in the eigenspectrum, and therefore the complexity of the predicted power spectrum. We find that the present 2D model makes relevant predictions in the case of quasitoroidal modes; however, going from 2D to 3D remains necessary to obtain more accurate mode amplitudes and to constrain the radial dependence of the turbulent viscosity and the source properties.
The present synthetic spectrum model may also be of interest to test mode detection pipelines. To that effect, specific spectrum realisations need to be drawn from the expected power distribution. This not only requires knowledge of the expected power spectrum (i.e. the variance associated with the complex Fourier transform for each pixel in the frequency–latitude plane), but also the correlation matrix between the different frequencies and latitudes. While different frequencies are completely uncorrelated under the hypothesis that the source of excitation is a stationary stochastic process, the correlations between different latitudes must be accounted for. The present framework would allow us to do this with only minimal adaptation. This constitutes one of the uses to which the present model will be put in the near future.
Finally, as we pointed out before, the present model is not meant to treat the case of unstable inertial modes. Some highlatitude modes can be shown to be selfexcited, either because of a steep rotational shear (Fournier et al. 2022) or by a baroclinic instability (Bekki et al. 2022a). The latter is responsible for the unstable nature of the m = 1 highlatitude mode, which is easily identifiable in the solar observations. The amplitude reached by this mode is related to a nonlinear saturation process, which is well outside the scope of the present linear study.
Acknowledgments
We are very grateful to ZhiChao Liang for providing the Rossbymode observations and Christian Baumgartner for the surface vorticity data. This work is supported in part by the ERC Synergy Grant WHOLESUN 810218.
References
 Allison, M. 1990, Icarus, 83, 282 [NASA ADS] [CrossRef] [Google Scholar]
 Augustson, K. C., Mathis, S., & Astoul, A. 2020, ApJ, 903, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Balmforth, N. J. 1992, MNRAS, 255, 639 [NASA ADS] [Google Scholar]
 Bekki, Y., Cameron, R. H., & Gizon, L. 2022a, A&A, 666, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bekki, Y., Cameron, R. H., & Gizon, L. 2022b, A&A, 662, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., et al. 2010, A&A, 522, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boussinesq, J. 1877, Essai sur la théorie des eaux courantes, Mémoires présentés par divers savants à l’Académie des sciences de l’Institut national de France (Imprimerie Nationale) [Google Scholar]
 Chaplin, W. J., Houdek, G., Elsworth, Y., et al. 2005, MNRAS, 360, 859 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J. 2002, Rev. Modern Phys., 74, 1073 [NASA ADS] [CrossRef] [Google Scholar]
 Deubner, F. L. 1975, A&A, 44, 371 [Google Scholar]
 Dikpati, M., Gilman, P. A., Guerrero, G. A., et al. 2022, ApJ, 931, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Fournier, D., Gizon, L., & Hyest, L. 2022, A&A, 664, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gizon, L., Fournier, D., & Albekioni, M. 2020, A&A, 642, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gizon, L., Cameron, R. H., Bekki, Y., et al. 2021, A&A, 652, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldreich, P., & Keeley, D. A. 1977, ApJ, 212, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Hanson, C. S., Hanasoge, S., & Sreenivasan, K. R. 2022, Nat. Astron., 6, 708 [NASA ADS] [CrossRef] [Google Scholar]
 Kraichnan, R. H. 1965, Phys. Fluids, 8, 575 [NASA ADS] [CrossRef] [Google Scholar]
 Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 581, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leighton, R. B., Noyes, R. W., & Simon, G. W. 1962, ApJ, 135, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Lesieur, M. 2008, Turbulence in Fluids (Springer) [CrossRef] [Google Scholar]
 Liang, Z.C., Gizon, L., Birch, A. C., & Duvall, T. L. 2019, A&A, 626, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lighthill, M. J. 1967, in Aerodynamic Phenomena in Stellar Atmospheres, ed. R. N. Thomas, 28, 429 [NASA ADS] [Google Scholar]
 Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat. Astron., 2, 568 [Google Scholar]
 Mack, L. M. 1976, J. Fluid Mech., 73, 497 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, S., Neiner, C., & Tran Minh, N. 2014, A&A, 565, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miesch, M. S., Brun, A. S., DeRosa, M. L., & Toomre, J. 2008, ApJ, 673, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Millionshchikov, M. 1941, Dokl. Akad. Nauk SSSR, 32, 611 [Google Scholar]
 Musielak, Z. E., Rosner, R., Stein, R. F., & Ulmschneider, P. 1994, ApJ, 423, 474 [Google Scholar]
 Orszag, S. A. 1971, J. Fluid Mech., 50, 689 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J., & Pringle, J. E. 1978, MNRAS, 182, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Rossby, C.G. 1939, J. Marine Res., 2, 38 [CrossRef] [Google Scholar]
 Salwen, H., & Grosch, C. E. 1981, J. Fluid Mech., 104, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Samadi, R., & Goupil, M. J. 2001, A&A, 370, 136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Georgobiani, D., Trampedach, R., et al. 2007, A&A, 463, 297 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 SánchezLavega, A., RíoGaztelurrutia, T., Hueso, R., et al. 2014, Geophys. Rev. Lett., 41, 1425 [CrossRef] [Google Scholar]
 Stein, R. F. 1967, Sol. Phys., 2, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, R. F., & Nordlund, Å. 1998, ApJ, 499, 914 [Google Scholar]
 Tennekes, H., & Lumley, J. 1978, in A First Course in Turbulence (MIT Press) [Google Scholar]
 Triana, S. A., Guerrero, G., Barik, A., & Rekier, J. 2022, ApJ, 934, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Xiong, D.R. 1989, A&A, 209, 126 [NASA ADS] [Google Scholar]
Appendix A: Vorticity wave equation in the equatorial β plane
A.1. Transport equation for the stream function
The total flow velocity v is decomposed into a background zonal flow U ≡ U(y) e_{x}, which represents differential rotation, and a residual flow u ≡ u_{x}(x, y)e_{x} + u_{y}(x, y)e_{y}, which contains both the waves and the convective noise. The equations of motion become
where U′=dU/dy. Assuming that the residual flow is incompressible, there exists a scalar field Ψ, the stream function, such that
where e_{z} is the radial unit vector. Explicitly,
Differentiating the xcomponent of the equation of motion with respect to y and the ycomponent with respect to x, and subtracting the two, one obtains the following equation:
where is the horizontal Laplacian operator. In the case of a nonadiabatic flow, the first term on the righthand side does not vanish, because the density and pressure streamlines are distinct. However, if we disregard nonadiabatic effects, the flow can be considered barotropic, meaning that the gas pressure is a function of density only,
Then the gas pressure gradient and density gradient are aligned, and we have
We thus obtain a purely mechanical equation,
A.2. Turbulent viscosity
The last two terms of Eq. A.9 can be split into an average quantity and a fluctuation around this average. As we will see, the former gives rise to a turbulent viscous term in the wave equation, whereas the latter will act as a source term. We define
where the notation ⟨.⟩_{h} refers to a horizontal average over scales larger than the turbulent scale, but smaller than the scale of the inertial modes. The scale separation that allows for the definition of this average is discussed in the main text (see Sect. 2.3).
The average term can be rewritten as
where the second equality stems from the assumed incompressible nature of the flow, which translates to ∇ ⋅ u = 0. It can be seen that this is the divergence of a mean flux representing the transport of the quantity ΔΨ (i.e. minus the radial vorticity) by the flow itself. We assumed that this mode of transport can be described by a diffusion process, characterised by an effective turbulent diffusion coefficient – or turbulent viscosity – ν_{turb}, so that
This is analogous to the Boussinesq (1877) approximation for the Reynolds stress tensor, which is customarily extended to other moments of the form ⟨uX⟩ in Reynoldsaveraged NavierStokes (RANS) models (see for instance Xiong 1989).
Then, Eq. A.9 becomes
A.3. Linear inhomogeneous wave equation
In order to derive a linear wave equation from Eq. A.13, we decomposed the total stream function, Ψ, into the contribution from the oscillations, Ψ_{osc}, and a contribution from the convective noise (i.e. the turbulence), Ψ_{turb}:
In the region of excitation, we assume that
This is justified in the bulk of the convective region, where smallscale convection dominates the dynamics of the star. This remains true close to the surface of the star or the tachocline: for example, in the Sun, simulations show that the typical turbulent velocities near the surface of the Sun are of the order of a few km s^{−1} (Stein & Nordlund 1998), while the typical amplitudes of inertial modes are of the order of several m s^{−1} at most. However, there must be a smooth transition between the region of wave excitation and the neighbouring radiative zone, where Ψ_{turb} is negligible. Therefore, the ordering Ψ_{osc} ≪ Ψ_{turb} cannot be valid everywhere. In the following, we work in the region of wave excitation. We also note that, by definition of the horizontal average used to introduce the turbulent viscosity, we have
Keeping only the firstorder contributions in Ψ_{osc}, but all orders in Ψ_{turb}, we obtain, in the region of wave excitation,
where we have gathered all inhomogeneous terms on the righthand side: these constitute the random forcing terms. We note that we also considered ν_{turb} to be uniform so that it can be pulled out of the gradient operator. The lefthand side of Eq. A.18 is split two ways: everything outside the brackets represents the deterministic linear operator governing the propagation of the waves, and everything inside the brackets represent the random fluctuations of the medium in which the waves propagate, and constitute a stochastic perturbation to the linear propagation operator.
In this study, we are interested in the excitation of the vorticity waves by turbulence, and therefore will not concern ourselves with the stochastic perturbation to the propagation operator. For this reason, we cast aside the bracket term on the lefthand side of Eq. A.17. Furthermore, as is usually done while dealing with pmodes (e.g. Samadi & Goupil 2001), we consider from the start that the linear forcing has a negligible effect on the inertial modes, because the frequencies and wavevectors in which the turbulence has significant power are far removed from those of the oscillations. This amounts to neglecting the first two terms on the righthand side of Eq. A.17 (which are linear in Ψ_{turb}) compared to the bracket term (which is quadratic in Ψ_{turb}). The vorticity wave equation then becomes
Appendix B: Source correlations in the frequencywavenumber domain
The contribution of the inertial modes to our synthetic power spectrum model involves the following integral (see Eq. 18)
where we recall that denotes the Fourier transform in (t, x) and is defined by Eq. 7, and the source term S(t, x, y) is defined by the righthand side of Eq. 6. To shorten notations, we note that the latter can be rewritten as
where we have adopted Einstein’s convention on repeated indices. Subtracting the horizontal average from the source term has, in fact, no effect on its autocorrelation spectrum, as is easily seen if the source term given by Eq. B.2 is explicitly expanded according to Eq. A.10; therefore, in the following, we drop the notation δ altogether. Eq. B.1 then becomes
Since the source is a quadratic function of the turbulent fluctuations of the stream function, the power spectrum ends up depending on fourthorder correlation products thereof. In the following, we make the assumption that Ψ_{turb} and its derivatives follow a multivariate normal distribution (Millionshchikov 1941), in which case the fourthorder correlation product can be expanded in terms of secondorder products only according to
Then Eq. B.3 becomes
where
and
The first integral ℐ_{a} vanishes, because it only involves onepoint, onetime correlation products, and therefore none of them depends on the time increment t − t′ or the space increment x − x′. This leaves us with only the last two integrals to consider.
First, we performed the following change of variables,
so that
Then we consider that the twopoint, twotime correlation products only depend on the time and space differences (i.e. τ, δx, and Y), and not on the absolute time and space coordinates (i.e. T, X, and y_{s}), in which case these two integrals simplify to
and
If we introduce the following functions,
where δx ≡ δxe_{x} + Ye_{y}, then these two integrals can be rewritten more compactly as
where the notation denotes the Fourier transform in (t, x, y), and is defined by
Then, expanding the Fourier transform of the products as convolution integrals, and exploiting the fact that we obtain
where k ≡ k_{x}e_{x}.
Next, we need to express the quantities , , and appearing in these integrals. Taking the Fourier transforms of each of Eq. B.14, we find
and
Because we have assumed a homogeneous and stationary turbulence, the correlation products can be rewritten in terms of Dirac distributions, yielding nonzero contributions only if ω″= − ω′ and k″ = k′. We obtain
The function ℰ_{Ψ} denotes the stream function turbulent spectrum, defined by
Plugging these into Eq. B.18 and Eq. B.19, we find
and
Flipping the indices k and l in the second integral, and exploiting the fact that ϵ_{klz} = −ϵ_{lkz}, the sum of these two integrals becomes
Finally, we explicitly expanded the index contractions. A basic property of the LeviCivita symbol is
Seeing as neither k′ nor k has a component along the zaxis, we have
Thus we obtain
This expression can be rendered more symmetric by shifting the two variables ω′ and k′ by ω/2 and k/2, respectively, which leads to the final expression for the integral ℐ as a function of the stream function turbulent spectrum, ℰ_{Ψ},
Appendix C: Chebyshev representation of the linear operator
In this Appendix, we derive the components of the linear operator ℒ, defined by Eq. 9, in the dual basis formed by the Chebyshev polynomials of the first kind. The derivation follows the work by Orszag (1971). We recall that the Chebyshev polynomials of the first kind are defined by
Using the definition of the linear operator, we can write
where
and and T_{j}⁗ denote the second and fourth derivatives with respect to ξ, respectively. We derive and in Sect. C.1, and and in Sect. C.2.
C.1. Projection of the derivatives
Given a function
we can write, for any positive integer, p,
In order to compute the coefficients and , we needed to find the recurrence relation between the and the , that is, a recurrence relation for all the derivatives of f.
For any p ≥ 1 we can write
On the other hand, from Eq. C.1, it is easily seen that
where we recall that c_{i} = 1 + δ_{i0}, and we have also introduced a new factor (defined by and ). Therefore, we can also write
Equating the T_{i} coefficients in Eqs. C.12 and C.14, we find
which constitutes a recurrence relation in terms of p. This recurrence relation, combined with the condition that we must have for any value of p, is easily solved to yield
where the notation ≡a[b] means ‘equal to a modulo b’.
For example, the coefficients of the first derivative straightforwardly read
The coefficients of the second derivative can be computed in the following way
Similarly, we find
From these expressions it directly follows that
and
C.2. Projection of the products
In order to compute the coefficients and , we need a rule for the expansion of a product on the Chebyshev basis. We now show that this rule takes the form of a convolution. We considered two arbitrary functions f and g, expanded on the Chebyshev polynomials according to
It will be useful, in the following, to introduce the following functions
so that we have
It directly follows that
Forming the product of these two expressions, and using the identity , we obtain
which can be rewritten as
Finally, using Eq. C.26 to revert back to an expansion on the T_{i}, one finally finds
The differential rotation profile in the expressions of and is given by , which is straightforwardly expanded as
so that
From this expression directly stem the following
All Figures
Fig. 1. Spatial turbulent spectrum and turbulent frequency. Left: Spatial turbulent spectrum, E_{iso}(k), extracted from solar observations (solid blue line). The spectrum comprises roughly three sections, each with a different power law in k^{α}. The dashed orange line shows the best fit of the data to a threepiecewise power law. The exponents and the boundaries between the three regimes are indicated on the plot. Right: Turbulent frequency ω_{k} ≡ ku_{k} as a function of wavenumber, k, where u_{k} is the typical velocity of the turbulent eddies of inverse size k, given by Eq. (47). 

In the text 
Fig. 2. Temporal spectrum , as a function of the reduced frequency . The data points (in orange) have been binned according to the value of (which depends on both ω and k), and only the mean over each bin is shown. These data points collapse onto a unique, slowly varying curve, which was adjusted alternatively with a Gaussian function (Eq. (50), solid black line) and a Lorentzian function (Eq. (51), solid green line). The latter is clearly the best fit, obtained for an amplitude A = 1.85 and a standard deviation σ = 0.62. 

In the text 
Fig. 3. Latitudinal velocity spectrum. Top: discrete inertial mode eigenspectrum for m = 3, shown in the complex plane. Each point correspond to one mode: all eigenfrequencies have a negative real part (meaning the modes are retrograde) and a negative imaginary part (meaning they are stable). The coloured dots mark the modes whose contribution to the u_{y} spectrum is most prominent. The vertical dashed lines indicate the central frequencies of each resonant peak in the u_{y} power spectrum (see bottom panel). Middle: equatorial u_{y} spectrum, obtained from Eq. (16) for azimuthal order m = 3. The solid black line shows the total spectrum, while each coloured dashed line corresponds to the individual contribution of the normal eigenmodes of the system, as described in the main text. The colours of the dashed lines match the colour scheme of the top panel. The red vertical dashed lines show the local maxima of the total spectrum. Bottom: same as the middle panel, but the vertical scale is linear. 

In the text 
Fig. 4. Same as Fig. 3, but for m = 5. 

In the text 
Fig. 5. Same as Fig. 3, but for m = 10. The three classical branches of Poiseuille flows are indicated on the plot, in addition to the R mode. 

In the text 
Fig. 6. Synthetic equatorial spectrum in the mω plane, in terms of u_{y} (top), u_{x} (middle), and ζ (bottom). Each vertical slice is normalised separately such that the maximum is unity. The diamonds show the real part of the eigenfrequencies of the linear homogeneous problem, computed as described in Sect. 3.1.1. The colour code refers to the mode categories: the blue diamonds represent the A branch, the green diamonds represent the P branch, and the red diamonds represent the Rossby modes (see Sect. 3.1.1 for a description of these branches). The solid red line shows the theoretical dispersion relation for sectoral Rossby modes in Cartesian coordinates (see Eq. (59)). 

In the text 
Fig. 7. Equatorial Rossby mode parameters. Left: full width at half maximum of the resonant peaks that could be identified as equatorial Rossby modes, as a function of m. The diamonds show the linewidths measured in the u_{y} equatorial synthetic power spectrum, and the coloured solid lines represent the theoretical Rossby mode linewidth, obtained from the classical dispersion relation (see Eq. (60)). The colour code refers to the value of the turbulent Reynolds number: Re_{turb} = 300 (red), 700 (blue), and 1000 (green). The black line shows the mode linewidths inferred from solar observations at the equator in the latitudinal velocity spectrum, as reported by Liang et al. (2019). Error bars from the fitting procedure reported by the authors are also shown. Right: Rossby mode amplitude (coloured solid lines) in the u_{y} equatorial synthetic power spectra, as a function of azimuthal order m, defined as described in the text (see Eq. (58)). The colour code is identical to the one in the left panel. 

In the text 
Fig. 8. Geometrical argument used to identify the modes, k_{mode}, that can be excited by a single scale of turbulence, k_{0}. The integrand in Eq. (19) peaks at k′ such that both k′ + k_{mode}/2 and k′ − k_{mode}/2 are equal to k_{0}. The driving can therefore only be efficient if the two black circles in the figure intersect, i.e. if k_{mode} < 2k_{0}. 

In the text 
Fig. 9. Latitudinal velocity spectrum, estimated at the solar equator, for azimuthal orders m = 3 (top), 5 (middle), and 10 (bottom), and for a turbulent Reynolds number Re_{turb} = 300. The solid red line represents our model, including the contribution both from the inertial modes and from the turbulent noise. The thin blue line shows the solar observations reported by Liang et al. (2019), and the thicker blue line shows the best Lorentzian fit to their data, as reported in their paper. 

In the text 
Fig. 10. Power spectrum of the latitudinal velocity, u_{y}, as a function of frequency (horizontal axis) and latitude (vertical axis). Each panel corresponds to a different azimuthal order: m = 3 (top), m = 5 (middle), and m = 10 (bottom). The turbulent Reynolds number is set to Re_{turb} = 300. The solid black curve shows the critical latitudes, where the differential rotation exactly matches the phase velocity of the inertial waves, and is defined by the implicit relation Eq. (61). The dashed vertical lines show the real part of the eigenfrequencies of the homogeneous problem, with the same colour code as in Fig. 6. 

In the text 
Fig. 11. Same as Fig. 10, but for the power spectrum of the azimuthal velocity, u_{x}. 

In the text 
Fig. 12. Same as Fig. 10, but for the power spectrum of the radial vorticity, ζ. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.