Interaction of solar inertial modes with turbulent convection A 2D model for the excitation of linearly stable modes

,


Introduction
Multiple types of waves can propagate in the interior of a star.In the case of a non-rotating star, these modes are spheroidal; they are the p-modes (or acoustic modes), the f -modes (or surfacegravity modes), and the g-modes (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 Christensen-Dalsgaard 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 r-modes (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ánchez-Lavega 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 high-latitude modes and critical latitude modes, allowed by latitudinal differential rotation.Furthermore, Hanson et al. (2022) report the observation of high-frequency waves of vorticity that are anti-symmetric with respect to the equator.These inertial modes allow the convective zone to be probed in a way that complements p-mode 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 p-modes (Lighthill 1967;Goldreich & Keeley 1977) or gravito-inertial waves (Mathis et al. 2014;Augustson et al. 2020).Non-linear 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 non-linear 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 p-modes (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 Section 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.

Stochastic excitation by turbulence
We study the excitation of the vorticity modes observed on the Sun.These modes are quasi-toroidal (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), ∆ = ∂ 2 x + ∂ 2 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 Section 2.3.The left-hand 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 right-hand 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 solar-like p-mode excitation (Goldreich & Keeley 1977); a key difference, however, is that whereas p-modes are mainly excited by the vertical turbulent motions, the quasi-toroidal 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 S refers to the (x, t) Fourier transform of the righthand side of Eq. ( 6), and L is the linear propagation operator, given by where ∆ ≡ d 2 / dy 2 − k 2 x .The linear operator L and the source S depend on y, the angular frequency ω, and the longitudinal wavenumber k x .
The solution Ψ osc 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 u x,osc , the latitudinal velocity u y,osc and the radial vorticity ζ osc , for which we obtain Then we obtain the expectation of the power spectral density by forming the modulus squared of Eqs. ( 12) to ( 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 Section 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 I(y s ) denotes the source covariance, and is defined by We recall that the source term S depends on ω and k x , so that I 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 E Ψ represents the turbulent stream function spectrum, and is given by (see Eq. (B.28)) 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) to ( 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, E Ψ .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) to ( 17)) and the noise contributions (i.e.Eqs. ( 21) to ( 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 L, and (ii) the turbulent stream function spectrum E Ψ .

The Green function
The angular frequency, ω, and azimuthal wavenumber, k x , being fixed, the linear operator, L, 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 non-dimensionalised latitudinal coordinate ξ.We chose the same value U = 244 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 δ i j 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 G(ξ s ), so that where we also introduced the non-dimensionalised source position ξ s ≡ y s /R.The vector G(y s ) is the solution of the following linear system: where the column vector on the right-hand side comprises the projections of the Dirac distribution on the Chebyshev polynomials, and the matrix M on the left-hand 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 M i j 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 M 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 right-hand side of Eq. ( 30) depends on the source position ξ s , this is not the case of the matrix M, 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 G(ξ 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 , T i (1) = i 2 and T i (−1) = (−1) i+1 i 2 , so that enforcing these boundary conditions amounts to ensuring that the solution G(ξ s ) of Eq. ( 30) verifies These conditions can be enforced in Eq. ( 30) by replacing the last four lines of the matrix M by and by replacing the last four components of B(ξ 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 G(ξ s ) now yields the correct Green function, with the appropriate boundary conditions.

The turbulent stream function spectrum
The turbulent spectrum E Ψ , 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 Eq. ( 20) or in Eq. ( 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 E Ψ 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 solar-like p-modes, 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 equation 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 k-independent 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: 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 p-mode 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: 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 p-mode 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).

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 Eq. ( 16) and Eq. ( 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.

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 LΨ = 0, where we recall that L 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 Section 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) to ( 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.
We also need to compute the complex eigenfrequencies and eigenfunctions of the associated adjoint system L † Ψ = 0, where the adjoint linear operator L † 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 ω † n and eigenfunctions Ψ † n are the solution of the following generalised boundary eigenvalue problem (Salwen & Grosch 1981) where It is easy to show that ω † n = ω * n .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) to ( 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).

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 √ 2 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.1for 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 ω 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 low-m 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 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 low-m 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 x (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.m/s after m ∼ 10, while a value Re turb = 1000 causes them to reach an amplitude of 2.5 m/s 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 low-m 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.

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 south-north 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 fo-cus of this study is not on the mode frequencies, so that these discrepancies do not constitute an issue.
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 high-m visibility issue.We show in Sect.3.1.2that, in this high-m 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 signal-to-noise 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.

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: where U is given by Eq. 24.
The shape of the power spectrum clearly transitions from a low-m regime, dominated by a few, clearly identifiable and distinguishable resonant modes, to a high-m regime, where the power is concentrated in a characteristic crescent-shaped 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 low-m to high-m 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 in- dicates 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 low-m to high-m transition described above can still be observed in the vorticity power spectra.

Conclusion
We have designed a model for the stochastic excitation of linearly stable, quasi-toroidal 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, Poiseuille-like pro- file, 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 mid-latitude 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, high-latitude inertial modes.
We also show that the power spectra in the frequencylatitude plane have a very different qualitative behaviour depend-ing on whether m 5 or m 5.In the low-m regime, the spectra are dominated by non-overlapping, clearly identifiable and distinguishable resonant modes.On the other hand, in the high-m 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 crescent-shaped 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 convection-induced 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 quasi-toroidal 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 self-excited, 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 high-latitude mode, which is easily identifiable in the solar observations.The amplitude reached by this mode is related to a non-linear saturation process, which is well outside the scope of the present linear study.
Differentiating the x-component of the equation of motion with respect to y and the y-component with respect to x, and subtracting the two, one obtains the following equation: where ∆ ≡ ∂ 2 x + ∂ 2 y is the horizontal Laplacian operator.In the case of a non-adiabatic flow, the first term on the right-hand side does not vanish, because the density and pressure streamlines are distinct.However, if we disregard non-adiabatic effects, the flow can be considered barotropic, meaning that the gas pressure is a function of density only, p = p(ρ) .We thus obtain a purely mechanical equation, (A.9) Appendix 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  .11)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 .12)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 Reynolds-averaged Navier-Stokes (RANS) models (see for instance Xiong 1989).Then, Eq.A.9 becomes 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 .15)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 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 left-hand 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 left-hand side of Eq.A.17. Furthermore, as is usually done while dealing with p-modes (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 right-hand 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 frequency-wavenumber 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 right-hand 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 fourth-order 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 fourth-order correlation product can be expanded in terms of second-order products only according to abcd = ab cd + ac bd + ad bc . (B.4) Then Eq.B.3 becomes where and The first integral I a vanishes, because it only involves one-point, one-time 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, Then we consider that the two-point, two-time 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, B.14) 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 T obs X obs Y obs dt dx dy f (t, x, y)e j(ωt−k x x−k y y) . (B.17) Then, expanding the Fourier transform of the products as convolution integrals, and exploiting the fact that g(−ω, −k) = g * (ω, k) we obtain where k ≡ k x e x .Next, we need to express the quantities f b,ik , g b, jl , f c,il and g c, jk appearing in these integrals.Taking the Fourier transforms of each of Eq.B.14, we find Finally, we explicitly expanded the index contractions.A basic property of the Levi-Civita symbol is i jz klz = δ ik δ jl δ zz +δ il δ jz δ kz +δ iz δ jk δ lz −δ ik δ jz δ lz −δ iz δ jl δ kz −δ il δ jk δ zz . (B.32) Seeing as neither k nor k has a component along the z-axis, we have Thus we obtain .34)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 I as a function of the stream function turbulent spectrum, E Ψ ,

Appendix C: Chebyshev representation of the linear operator
In this Appendix, we derive the components of the linear operator L, 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 and T j and T j denote the second and fourth derivatives with respect to ξ, respectively.We derive M (2) i j and M (3) i j in Sect.C.1, and M (4) i j and M (5) i j in Sect.C.2.In order to compute the coefficients M (2) i j and M (3) i j , we needed to find the recurrence relation between the λ (p)  i and the λ (p−1) i , that is, a recurrence relation for all the derivatives of f .
For any p 1 we can write (C.17) The coefficients of the second derivative can be computed in the following way (C.18) Similarly, we find (C.20)

Fig. 1 .
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 three-piece-wise 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).

Fig. 2 .
Fig.2.Temporal spectrum χ( ω) ≡ ω k χ k (ω), as a function of the reduced frequency ω ≡ ω/ω k .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.

Fig. 3 .
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 .
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.

Fig. 6 .
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.1for a description of these branches).The solid red line shows the theoretical dispersion relation for sectoral Rossby modes in Cartesian coordinates (see Eq. 59).

Fig. 7 .
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 Γ R = ν turb k 2x (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 byLiang 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.

Fig. 8 .
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 .

Fig. 9 .
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 byLiang et al. (2019), and the thicker blue line shows the best Lorentzian fit to their data, as reported in their paper.

Fig. 10 .
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.

Fig. 11 .
Fig. 11.Same as Fig. 10, but for the power spectrum of the azimuthal velocity, u x .

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