© The Authors 2023

1. 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 surface-gravity 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 quasi-toroidal 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 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 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

$$\begin{array}{cc}& x\equiv R\varphi ,\hfill \end{array}$$(1)

$$\begin{array}{cc}& y\equiv Rsin\lambda ,\hfill \end{array}$$(2)

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

$$\begin{array}{cc}& \frac{\partial v_x}{\partial t}+(\mathbf{v}·\mathbf{\nabla })v_x=-\frac{1}{\rho }\frac{\partial p}{\partial x}+f{v}_y,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial v_y}{\partial t}+(\mathbf{v}·\mathbf{\nabla })v_y=-\frac{1}{\rho }\frac{\partial p}{\partial y}-f{v}_x,\hfill \end{array}$$(4)

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

$$\begin{array}{c}\hfill \mathbf{u}=\mathbf{\nabla }\wedge \left(\mathrm{\Psi }{\mathbf{e}}_z\right),\end{array}$$(5)

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

$$\begin{array}{cc}& (\frac{\partial }{\partial t}+U\frac{\partial }{\partial x})\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{osc}}+(\beta -U^{″})\frac{\partial {\mathrm{\Psi }}_{\mathrm{osc}}}{\partial x}-{\nu }_{\mathrm{turb}}{\mathrm{\Delta }}^2{\mathrm{\Psi }}_{\mathrm{osc}}\hfill \\ \hfill & =\delta (\frac{\partial {\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}}{\partial y}-\frac{\partial {\mathrm{\Psi }}_{\mathrm{turb}}}{\partial y}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}),\hfill \end{array}$$(6)

where U″ is the second derivative of U(y), $\mathrm{\Delta }={\partial }_x^2+{\partial }_y^2$ 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 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_xx), where j denotes the imaginary unit. We used the following convention for the Fourier transform,

$$\begin{array}{c}\hfill {\displaystyle \widehat{f}(\omega ,k_x,y)\equiv \frac{1}{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}}}{\displaystyle \int \mathrm{d}}t\mathrm{d}xf(t,x,y)e^{j(\omega t-k_{x}x)},}\end{array}$$(7)

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

$$\begin{array}{c}\hfill \mathcal{L}{\widehat{\mathrm{\Psi }}}_{\mathrm{osc}}=\widehat{S}.\end{array}$$(8)

The notation $\widehat{S}$ refers to the (x, t) Fourier transform of the right-hand side of Eq. (6), and ℒ is the linear propagation operator, given by

$$\begin{array}{c}\hfill \mathcal{L}\equiv (\omega -k_{x}U)\widehat{\mathrm{\Delta }}-k_x(\beta -U^{″})-j{\nu }_{\mathrm{turb}}{\widehat{\mathrm{\Delta }}}^2,\end{array}$$(9)

where $\widehat{\mathrm{\Delta }}\equiv {\text{d}}^2/\text{d}{\mathit{y}}^2-k_x^2$. The linear operator ℒ and the source $\widehat{S}$ depend on y, the angular frequency ω, and the longitudinal wavenumber k_x.

The solution ${\widehat{\mathrm{\Psi }}}_{\mathrm{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),

$$\begin{array}{c}\hfill \mathcal{L}G(y,y_s)=\delta (y-y_s),\end{array}$$(10)

where δ is the Dirac distribution, so that

$$\begin{array}{c}\hfill {\displaystyle {\widehat{\mathrm{\Psi }}}_{\mathrm{osc}}(y)={\displaystyle {\int }_{-R}^R}\mathrm{d}{y}_{s}G(y,y_s)\widehat{S}(y_s).}\end{array}$$(11)

This solution leads to other physical quantities, such as the azimuthal velocity ${\widehat{u}}_{x,\mathrm{osc}}$, the latitudinal velocity ${\widehat{u}}_{\mathit{y},\mathrm{osc}}$ and the radial vorticity ${\widehat{\zeta }}_{\mathrm{osc}}$, for which we obtain

$$\begin{array}{cc}\hfill {\widehat{u}}_{x,\mathrm{osc}}(y)& {\displaystyle ={\displaystyle {\int }_{-R}^R}\mathrm{d}{y}_s\frac{\partial G}{\partial y}\widehat{S}(y_s),}\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill {\widehat{u}}_{y,\mathrm{osc}}(y)& {\displaystyle =-j{k}_x{\displaystyle {\int }_{-R}^R}\mathrm{d}{y}_{s}G(y,y_s)\widehat{S}(y_s),}\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill {\widehat{\zeta }}_{\mathrm{osc}}(y)& {\displaystyle ={\displaystyle {\int }_{-R}^R\mathrm{d}{y}_s(k_x^{2}G-\frac{{\partial }^{2}G}{\partial y^2})\widehat{S}(y_s).}}\hfill \end{array}$$(14)

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

$$\begin{array}{cc}\hfill \langle {\left|{\widehat{u}}_{x,\mathrm{osc}}(y_{\mathrm{obs}})\right|}^2\rangle & {\displaystyle ={\displaystyle {\int }_{-R}^R\mathrm{d}{y}_s{\left|\frac{\partial G}{\partial y_{\mathrm{obs}}}\right|}^2\mathcal{I}(y_s),}}\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill \langle {\left|{\widehat{u}}_{y,\mathrm{osc}}(y_{\mathrm{obs}})\right|}^2\rangle & {\displaystyle ={\displaystyle {\int }_{-R}^R}\mathrm{d}{y}_s{\left|k_{x}G(y_{\mathrm{obs}},y_s)\right|}^2\mathcal{I}(y_s),}\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill \langle {\left|{\widehat{\zeta }}_{\mathrm{osc}}(y_{\mathrm{obs}})\right|}^2\rangle & {\displaystyle ={\displaystyle {\int }_{-R}^R}\mathrm{d}{y}_s{\left|(k_x^{2}G-\frac{{\partial }^{2}G}{\partial {y_{\mathrm{obs}}}^2})\right|}^2\mathcal{I}(y_s),}\hfill \end{array}$$(17)

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

$$\begin{array}{c}\hfill {\displaystyle \mathcal{I}(y_s)\equiv {\displaystyle \int }\mathrm{d}Y\langle \widehat{S}(y_s){\widehat{S}}^{\ast }(y_s+Y)\rangle .}\end{array}$$(18)

We recall that the source term $\widehat{S}$ 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))

$$\begin{array}{cc}& {\displaystyle \mathcal{I}=\frac{1}{4{\pi }^3}{\displaystyle \int }\mathrm{d}{\omega }^{\prime }{\mathrm{d}}^2{\mathbf{k}}^{\prime }{k}_x^3{k}_x^{\prime }{k}_y^{\prime 2}{|{\mathbf{k}}^{\prime }+\mathbf{k}/2|}^2}\hfill \\ \hfill & \times {\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime }-\frac{\omega }{2},{\mathbf{k}}^{\prime }-\frac{\mathbf{k}}{2}){{\mathcal{E}}_{\mathrm{\Psi }}}^{\ast }({\omega }^{\prime }+\frac{\omega }{2},{\mathbf{k}}^{\prime }+\frac{\mathbf{k}}{2}).\hfill \end{array}$$(19)

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

$$\begin{array}{c}\hfill {\displaystyle {\mathcal{E}}_{\mathrm{\Psi }}(\omega ,\mathbf{k})\equiv {\displaystyle \int \mathrm{d}}\tau {\mathrm{d}}^2\mathbf{x}\langle {\mathrm{\Psi }}_{\mathrm{turb}}(T,\mathbf{X}){\mathrm{\Psi }}_{\mathrm{turb}}(T+\tau ,\mathbf{X}+\mathbf{x})\rangle e^{j(\omega \tau -\mathbf{k}·\mathbf{x})}.}\end{array}$$(20)

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

$$\begin{array}{cc}\hfill \langle {\left|{\widehat{u}}_{x,\mathrm{turb}}(y_{\mathrm{obs}})\right|}^2\rangle & {\displaystyle =\frac{1}{2\pi }{\displaystyle \int \mathrm{d}}{k}_y{\mathcal{E}}_{\mathrm{\Psi }}(\omega ,\mathbf{k})k_y^2,}\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill \langle {\left|{\widehat{u}}_{y,\mathrm{turb}}(y_{\mathrm{obs}})\right|}^2\rangle & {\displaystyle =\frac{1}{2\pi }{\displaystyle \int \mathrm{d}}{k}_y{\mathcal{E}}_{\mathrm{\Psi }}(\omega ,\mathbf{k})k_x^2,}\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill \langle {\left|{\widehat{\zeta }}_{\mathrm{turb}}(y_{\mathrm{obs}})\right|}^2\rangle & {\displaystyle =\frac{1}{2\pi }{\displaystyle \int \mathrm{d}}{k}_y{\mathcal{E}}_{\mathrm{\Psi }}(\omega ,\mathbf{k}){(k_x^2+k_y^2)}^2.}\hfill \end{array}$$(23)

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,

$$\begin{array}{c}\hfill U(\xi )=-\overline{U}{\xi }^2,\xi \equiv y/R=sin\lambda ,\end{array}$$(24)

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 $\overline{U}=244$ m s⁻¹ 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⁻¹⁵ m⁻¹s⁻¹. Finally, the turbulent viscosity, ν_turb, is specified through the turbulent Reynolds number,

$$\begin{array}{c}\hfill {\mathrm{Re}}_{\mathrm{turb}}\equiv \frac{\overline{U}R}{{\nu }_{\mathrm{turb}}},\end{array}$$(25)

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

$$\begin{array}{c}\hfill T_n(\xi )=cos(narccos\xi ),\end{array}$$(26)

which reduces to a polynomial expression after some algebraic manipulations. These polynomials are orthogonal to each other with respect to the following inner product,

$$\begin{array}{c}\hfill {\displaystyle f|g\equiv {\displaystyle {\int }_{-1}^1}\frac{f(\xi )g(\xi )}{\sqrt{1-{\xi }^2}}\mathrm{d}\xi ,}\end{array}$$(27)

in the sense that

$$\begin{array}{c}\hfill T_i|T_j=\frac{\pi }{2}{c}_i{\delta }_{\mathit{ij}},\end{array}$$(28)

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

$$\begin{array}{c}\hfill G(\xi ,{\xi }_s)={\displaystyle \sum _{i=0}^{\infty }}{\mathcal{G}}_i({\xi }_s)T_i(\xi ),\end{array}$$(29)

where we also introduced the non-dimensionalised source position ξ_s ≡ y_s/R. The vector 𝒢(y_s) is the solution of the following linear system:

$$\begin{array}{c}\hfill \mathcal{M}\mathcal{G}({\xi }_s)=\mathcal{B}({\xi }_s),\end{array}$$(30)

where the column vector on the right-hand side comprises the projections of the Dirac distribution on the Chebyshev polynomials,

$$\begin{array}{c}\hfill {\mathcal{B}}_i\equiv \frac{2}{\pi c_i}{T}_i|\delta (\xi -{\xi }_s)=\frac{2T_i({\xi }_s)}{\pi c_i\sqrt{1-{\xi }_s}},\end{array}$$(31)

and the matrix ℳ on the left-hand side is defined by

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathit{ij}}\equiv \frac{2}{\pi c_i}{T}_i|\mathcal{L}{T}_j.\end{array}$$(32)

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 right-hand 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:

$$\begin{array}{cc}\hfill G(\xi =-1,{\xi }_s)& =G(\xi =1,{\xi }_s)=0,\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill {\frac{\partial G}{\partial \xi }|}_{\xi =-1,{\xi }_s}& ={\frac{\partial G}{\partial \xi }|}_{\xi =1,{\xi }_s}=0.\hfill \end{array}$$(34)

The Chebyshev polynomials of the first kind verify T_i(1) = 1, T_i(−1) = (−1)ⁱ, $T_i^{\prime }(1)=i^2$ and $T_i^{\prime }(-1)={(-1)}^{i+1}{i}^2$, so that enforcing these boundary conditions amounts to ensuring that the solution 𝒢(ξ_s) of Eq. (30) verifies

$$\begin{array}{cc}\hfill {\displaystyle \sum _{i=0}^{N-1}}{\mathcal{G}}_i& ={\displaystyle \sum _{i=0}^{N-1}}{(-1)}^i{\mathcal{G}}_i=0,\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill {\displaystyle \sum _{i=0}^{N-1}}{i}^2{\mathcal{G}}_i& ={\displaystyle \sum _{i=0}^{N-1}}{(-1)}^{i+1}{i}^2{\mathcal{G}}_i=0.\hfill \end{array}$$(36)

These conditions can be enforced in Eq. (30) by replacing the last four lines of the matrix ℳ by

$$\begin{array}{cc}\hfill {\mathcal{M}}_{N-4,j}& =1,{\mathcal{M}}_{N-3,j}={(-1)}^j,\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill {\mathcal{M}}_{N-2,j}& =j^2,{\mathcal{M}}_{N-1,j}={(-1)}^{j+1}{j}^2,\hfill \end{array}$$(38)

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 high-frequency 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

$$\begin{array}{c}\hfill {\displaystyle {\varphi }_{\mathit{ij}}(\omega ,\mathbf{k})\equiv {\displaystyle \int \mathrm{d}}\tau {\mathrm{d}}^2\mathbf{x}\langle u_{i,\mathrm{turb}}(T,\mathbf{X})u_{j,\mathrm{turb}}(T+\tau ,\mathbf{X}+\mathbf{x})\rangle {\mathrm{e}}^{j(\omega \tau -\mathbf{k}·\mathbf{x})}.}\end{array}$$(39)

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

$$\begin{array}{c}\hfill {\varphi }_{\mathit{ij}}=(k^2{\delta }_{\mathit{ij}}-k_i{k}_j){\mathcal{E}}_{\mathrm{\Psi }}.\end{array}$$(40)

The turbulent velocity spectrum is usually expressed as (e.g. Lesieur 2008)

$$\begin{array}{c}\hfill {\varphi }_{\mathit{ij}}=\mathrm{E}(\mathbf{k},\omega )({\delta }_{\mathit{ij}}-\frac{k_i{k}_j}{k^2}).\end{array}$$(41)

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

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{\Psi }}(\omega ,\mathbf{k})=\frac{E(\omega ,\mathbf{k})}{k^2},\end{array}$$(42)

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

$$\begin{array}{c}\hfill {\displaystyle E(\mathbf{k})\equiv {\displaystyle {\int }_{-\infty }^{+\infty }}\mathrm{d}\omega E(\mathbf{k},\omega )}\end{array}$$(43)

and then define the temporal part of the spectrum as

$$\begin{array}{c}\hfill {\chi }_{\mathbf{k}}(\omega )\equiv \frac{E(\mathbf{k},\omega )}{E(\mathbf{k})}.\end{array}$$(44)

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,

$$\begin{array}{c}\hfill E(\mathbf{k})=\frac{E_{\mathrm{iso}}(k)}{2\pi k}.\end{array}$$(45)

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

$$\begin{array}{c}\hfill {\chi }_{\mathbf{k}}(\omega )=\frac{1}{{\omega }_k}\chi (\stackrel{\sim }{\omega }),\stackrel{\sim }{\omega }\equiv \omega /{\omega }_k,\end{array}$$(46)

where we have introduced the reduced frequency $\stackrel{\sim }{\omega }$, 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

$$\begin{array}{c}\hfill u_k\equiv {\left({\displaystyle {\displaystyle {\int }_{k<k^{\prime }<2k}}\mathrm{d}{k}^{\prime }{E}_{\mathrm{iso}}(k^{\prime })}\right)}^{1/2}.\end{array}$$(47)

All in all, we can write

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{\Psi }}=\frac{E_{\mathrm{iso}}(k)u_k\chi (\stackrel{\sim }{\omega })}{2\pi k^2},\end{array}$$(48)

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:

$$\begin{array}{c}\hfill E_{\mathrm{iso}}(k)=\{\begin{array}{cc}{C}_1{\left(\frac{k}{k_{\mathrm{ref}}}\right)}^{{\alpha }_1}\hfill & \mathrm{if}k\le k_1,\hfill \\ C_2{\left(\frac{k}{k_{\mathrm{ref}}}\right)}^{{\alpha }_2}\hfill & \mathrm{if}{k}_1<k<k_2,\hfill \\ C_3{\left(\frac{k}{k_{\mathrm{ref}}}\right)}^{{\alpha }_3}\hfill & \mathrm{if}{k}_2\le k,\hfill \end{array}\end{array}$$(49)

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₁ and k₂, 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 α₁ = 0.27, α₂ = −1.88 and α₃ = −7.60; in particular, the slope α₂ 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₁R = 140 and k₂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 $\stackrel{\sim }{\omega }$, 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:

$$\begin{array}{c}\hfill \chi (\stackrel{\sim }{\omega })=\frac{A}{\sqrt{2\pi {\sigma }^2}}{e}^{-{\stackrel{\sim }{\omega }}^2/(2{\sigma }^2)}\end{array}$$(50)

and

$$\begin{array}{c}\hfill \chi (\stackrel{\sim }{\omega })=\frac{A}{\pi \sigma }\frac{1}{1+{(\stackrel{\sim }{\omega }/\sigma )}^2}.\end{array}$$(51)

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

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_xR, 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

$$\begin{array}{c}\hfill {\mathcal{L}}_1{\mathrm{\Psi }}_n={\omega }_n\widehat{\mathrm{\Delta }}{\mathrm{\Psi }}_n,\end{array}$$(52)

where

$$\begin{array}{c}\hfill {\mathcal{L}}_1\equiv k_{x}U\widehat{\mathrm{\Delta }}+k_x(\beta -U^{″})+j{\nu }_{\mathrm{turb}}{\widehat{\mathrm{\Delta }}}^2.\end{array}$$(53)

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.

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

$$\begin{array}{c}\hfill {\displaystyle {\displaystyle {\int }_{-1}^1}\mathrm{d}\xi {\left({\mathcal{L}}^{†}g\right)}^{\ast }f={\displaystyle {\int }_{-1}^1}\mathrm{d}\xi g^{\ast }\left(\mathcal{L}f\right).}\end{array}$$(54)

Performing integration by parts we find that the adjoint eigenfrequencies ${\omega }_n^{†}$ and eigenfunctions ${\mathrm{\Psi }}_n^{†}$ are the solution of the following generalised boundary eigenvalue problem (Salwen & Grosch 1981)

$$\begin{array}{c}\hfill {\mathcal{L}}_1^{†}{\mathrm{\Psi }}_n^{†}={\omega }_n^{†}\widehat{\mathrm{\Delta }}{\mathrm{\Psi }}_n^{†},\end{array}$$(55)

where

$$\begin{array}{c}\hfill {\mathcal{L}}_1^{†}\equiv k_{x}U\widehat{\mathrm{\Delta }}+k_x(\beta +2U^{\prime }\frac{\mathrm{d}}{\mathrm{d}\xi })-j{\nu }_{\mathrm{turb}}{\widehat{\mathrm{\Delta }}}^2.\end{array}$$(56)

It is easy to show that ${\omega }_n^{†}={\omega }_n^{\ast }$. Furthermore, the eigenfunctions and adjoint eigenfunctions form a biorthogonal set, in the sense that

$$\begin{array}{c}\hfill {\displaystyle {\displaystyle {\int }_{-1}^1}\mathrm{d}\xi (\frac{\mathrm{d}{\mathrm{\Psi }}_n^{\ast }}{\mathrm{d}y}\frac{\mathrm{d}{\mathrm{\Psi }}_m^{†}}{\mathrm{d}y}+k_x^2{\mathrm{\Psi }}_n^{\ast }{\mathrm{\Psi }}_m^{†})\propto {\delta }_{\mathit{nm}}.}\end{array}$$(57)

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 ω₀ 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

$$\begin{array}{c}\hfill A={\left({\displaystyle {\displaystyle {\int }_{{\omega }_a}^{{\omega }_b}}\mathrm{d}\omega P(\omega )}\right)}^{1/2},\end{array}$$(58)

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 = ω₀ ∓ Γ/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$\sqrt{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, ω₀, 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

$$\begin{array}{c}\hfill {\omega }_R=-\frac{\beta -U^{″}}{k_x}.\end{array}$$(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 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² 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

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_R={\nu }_{\mathrm{turb}}{k}_x^2.\end{array}$$(60)

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² s⁻¹. 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² s⁻¹ 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 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₀ can only drive modes with azimuthal wavenumbers k_x = m/R < 2k₀ (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.

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 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 focus 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.2 that, in this high-m regime, the amplitudes of the equatorial Rossby modes are fairly independent of m, while their linewidths increase as m². Because the peaks become wider and wider for the same total power, their spectral height decreases as 1/m², 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.

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:

$$\begin{array}{c}\hfill U(y_c)=\frac{\omega }{k_x}=\frac{\omega R}{m},\end{array}$$(61)

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 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 low-m to high-m 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, 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 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 mid-latitude inertial modes, with latitudinal velocities ranging between ∼0.1 and ∼1.5 m.s⁻¹, 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 frequency–latitude plane have a very different qualitative behaviour depending 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⁻¹. 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², 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 high-latitude 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.

Acknowledgments

We are very grateful to Zhi-Chao Liang for providing the Rossby-mode observations and Christian Baumgartner for the surface vorticity data. This work is supported in part by the ERC Synergy Grant WHOLESUN 810218.

References

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

$$\begin{array}{cc}& \frac{\partial u_x}{\partial t}+U\frac{\partial u_x}{\partial x}+u_y{U}^{\prime }+u_x\frac{\partial u_x}{\partial x}+u_y\frac{\partial u_x}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+f{u}_y,\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& \frac{\partial u_y}{\partial t}+U\frac{\partial u_y}{\partial x}+u_x\frac{\partial u_y}{\partial x}+u_y\frac{\partial u_y}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial y}-f{u}_x,\hfill \end{array}$$(A.2)

where U′=dU/dy. Assuming that the residual flow is incompressible, there exists a scalar field Ψ, the stream function, such that

$$\begin{array}{c}\hfill \mathbf{u}=\mathbf{\nabla }\wedge \left(\mathrm{\Psi }{\mathbf{e}}_z\right),\end{array}$$(A.3)

where e_z is the radial unit vector. Explicitly,

$$\begin{array}{cc}& u_x=\frac{\partial \mathrm{\Psi }}{\partial y},\hfill \end{array}$$(A.4)

$$\begin{array}{cc}& u_y=-\frac{\partial \mathrm{\Psi }}{\partial x}.\hfill \end{array}$$(A.5)

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:

$$\begin{array}{cc}& \frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial t}+U\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial x}-U^{″}\frac{\partial \mathrm{\Psi }}{\partial x}+\frac{\partial \mathrm{\Psi }}{\partial y}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial x}-\frac{\partial \mathrm{\Psi }}{\partial x}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial y}\hfill \\ \hfill & =\frac{1}{{\rho }^2}(\frac{\partial \rho }{\partial y}\frac{\partial p}{\partial x}-\frac{\partial \rho }{\partial x}\frac{\partial p}{\partial y})-\beta \frac{\partial \mathrm{\Psi }}{\partial x},\hfill \end{array}$$(A.6)

where $\mathrm{\Delta }\equiv {\partial }_x^2+{\partial }_y^2$ 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,

$$\begin{array}{c}\hfill p=p(\rho ).\end{array}$$(A.7)

Then the gas pressure gradient and density gradient are aligned, and we have

$$\begin{array}{c}\hfill \frac{\partial \rho }{\partial y}\frac{\partial p}{\partial x}-\frac{\partial \rho }{\partial x}\frac{\partial p}{\partial y}=(\mathbf{\nabla }p\wedge \mathbf{\nabla }\rho )·{\mathbf{e}}_z=0.\end{array}$$(A.8)

We thus obtain a purely mechanical equation,

$$\begin{array}{c}\hfill (\frac{\partial }{\partial t}+U\frac{\partial }{\partial x})\mathrm{\Delta }\mathrm{\Psi }+(\beta -U^{″})\frac{\partial \mathrm{\Psi }}{\partial x}+\frac{\partial \mathrm{\Psi }}{\partial y}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial x}-\frac{\partial \mathrm{\Psi }}{\partial x}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial y}=0.\end{array}$$(A.9)

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

$$\begin{array}{cc}& \delta (\frac{\partial \mathrm{\Psi }}{\partial y}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial x}-\frac{\partial \mathrm{\Psi }}{\partial x}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial y})\equiv \frac{\partial \mathrm{\Psi }}{\partial y}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial x}-\frac{\partial \mathrm{\Psi }}{\partial x}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial y}\hfill \\ \hfill & -{\langle \frac{\partial \mathrm{\Psi }}{\partial y}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial x}-\frac{\partial \mathrm{\Psi }}{\partial x}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial y}\rangle }_h,\hfill \end{array}$$(A.10)

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

$$\begin{array}{c}\hfill {\langle \frac{\partial \mathrm{\Psi }}{\partial y}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial x}-\frac{\partial \mathrm{\Psi }}{\partial x}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial y}\rangle }_h={\langle (\mathbf{u}·\mathbf{\nabla })\mathrm{\Delta }\mathrm{\Psi }\rangle }_h=\mathbf{\nabla }·{⟨\mathbf{u}\mathrm{\Delta }\mathrm{\Psi }⟩}_h,\end{array}$$(A.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

$$\begin{array}{c}\hfill ⟨\mathbf{u}\mathrm{\Delta }\mathrm{\Psi }⟩=-{\nu }_{\mathrm{turb}}\mathbf{\nabla }{\langle \mathrm{\Delta }\mathrm{\Psi }\rangle }_h.\end{array}$$(A.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

$$\begin{array}{cc}& (\frac{\partial }{\partial t}+U\frac{\partial }{\partial x})\mathrm{\Delta }\mathrm{\Psi }+(\beta -U^{″})\frac{\partial \mathrm{\Psi }}{\partial x}-\mathbf{\nabla }·({\nu }_{\mathrm{turb}}\mathbf{\nabla }\mathrm{\Delta }{⟨\mathrm{\Psi }⟩}_h)\hfill \\ \hfill & +\delta (\frac{\partial \mathrm{\Psi }}{\partial y}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial x}-\frac{\partial \mathrm{\Psi }}{\partial x}\frac{\partial \mathrm{\Delta }\mathrm{\Psi }}{\partial y})=0.\hfill \end{array}$$(A.13)

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:

$$\begin{array}{c}\hfill \mathrm{\Psi }={\mathrm{\Psi }}_{\mathrm{osc}}+{\mathrm{\Psi }}_{\mathrm{turb}}.\end{array}$$(A.14)

In the region of excitation, we assume that

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathrm{osc}}\ll {\mathrm{\Psi }}_{\mathrm{turb}}.\end{array}$$(A.15)

This is justified in the bulk of the convective region, where small-scale 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⁻¹ (Stein & Nordlund 1998), while the typical amplitudes of inertial modes are of the order of several m s⁻¹ 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

$$\begin{array}{c}\hfill {⟨\mathrm{\Psi }⟩}_h={\mathrm{\Psi }}_{\mathrm{osc}}.\end{array}$$(A.16)

Keeping only the first-order contributions in Ψ_osc, but all orders in Ψ_turb, we obtain, in the region of wave excitation,

$$\begin{array}{cc}& (\frac{\partial }{\partial t}+U\frac{\partial }{\partial x})\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{osc}}+(\beta -U^{″})\frac{\partial {\mathrm{\Psi }}_{\mathrm{osc}}}{\partial x}-{\nu }_{\mathrm{turb}}{\mathrm{\Delta }}^2{\mathrm{\Psi }}_{\mathrm{osc}}\hfill \\ \hfill & +\delta (\frac{\partial {\mathrm{\Psi }}_{\mathrm{osc}}}{\partial y}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}+\frac{\partial {\mathrm{\Psi }}_{\mathrm{turb}}}{\partial y}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{osc}}}{\partial x}\hfill \\ \hfill & -\frac{\partial {\mathrm{\Psi }}_{\mathrm{osc}}}{\partial x}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}}{\partial y}-\frac{\partial {\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{osc}}}{\partial y})\hfill \\ \hfill & =-U\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}-(\beta -U^{″})\frac{\partial {\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}\hfill \\ \hfill & +\delta (\frac{\partial {\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}}{\partial y}-\frac{\partial {\mathrm{\Psi }}_{\mathrm{turb}}}{\partial y}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}),\hfill \end{array}$$(A.17)

where we have gathered all inhomogeneous terms on the right-hand 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

$$\begin{array}{cc}& (\frac{\partial }{\partial t}+U\frac{\partial }{\partial x})\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{osc}}+(\beta -U^{″})\frac{\partial {\mathrm{\Psi }}_{\mathrm{osc}}}{\partial x}-{\nu }_{\mathrm{turb}}{\mathrm{\Delta }}^2{\mathrm{\Psi }}_{\mathrm{osc}}\hfill \\ \hfill & =\delta (\frac{\partial {\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}}{\partial y}-\frac{\partial {\mathrm{\Psi }}_{\mathrm{turb}}}{\partial y}\frac{\partial \mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}}{\partial x}).\hfill \end{array}$$(A.18)

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)

$$\begin{array}{c}\hfill {\displaystyle \mathcal{I}(y_s)\equiv {\displaystyle \int }\mathrm{d}Y\langle \widehat{S}(y_s){\widehat{S}}^{\ast }(y_s+Y)\rangle ,}\end{array}$$(B.1)

where we recall that $\widehat{.}$ 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

$$\begin{array}{c}\hfill S(t,x,y)=\delta \left({ϵ}_{\mathit{ijz}}{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}\right),\end{array}$$(B.2)

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

$$\begin{array}{cc}& {\displaystyle \mathcal{I}(y_s)=\frac{1}{T_{\mathrm{obs}}{X}_{\mathrm{obs}}}{\displaystyle \int }\mathrm{d}Y\mathrm{d}t\mathrm{d}{t}^{\prime }\mathrm{d}x\mathrm{d}{x}^{\prime }{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}}\hfill \\ \hfill & \langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t,x,y_s}{{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t^{\prime },x^{\prime },y_s+Y}\rangle e^{j(\omega (t-t^{\prime })-k_x(x-x^{\prime }))}.\hfill \end{array}$$(B.3)

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

$$\begin{array}{c}\hfill \langle abcd\rangle =\langle ab\rangle \langle cd\rangle +\langle ac\rangle \langle bd\rangle +\langle ad\rangle \langle bc\rangle .\end{array}$$(B.4)

Then Eq. B.3 becomes

$$\begin{array}{c}\hfill \mathcal{I}(y_s)={\mathcal{I}}_a(y_s)+{\mathcal{I}}_b(y_s)+{\mathcal{I}}_c(y_s),\end{array}$$(B.5)

where

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_a(y_s)\equiv \frac{1}{T_{\mathrm{obs}}{X}_{\mathrm{obs}}}{\displaystyle \int }\mathrm{d}Y\mathrm{d}t\mathrm{d}{t}^{\prime }\mathrm{d}x\mathrm{d}{x}^{\prime }{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}}\hfill \\ \hfill & \langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t,x,y_s}\rangle \langle {{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t^{\prime },x^{\prime },y_s+Y}\rangle \hfill \\ \hfill & e^{j(\omega (t-t^{\prime })-k_x(x-x^{\prime }))},\hfill \end{array}$$(B.6)

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_b(y_s)\equiv \frac{1}{T_{\mathrm{obs}}{X}_{\mathrm{obs}}}{\displaystyle \int }\mathrm{d}Y\mathrm{d}t\mathrm{d}{t}^{\prime }\mathrm{d}x\mathrm{d}{x}^{\prime }{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}}\hfill \\ \hfill & \langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t,x,y_s}{{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t^{\prime },x^{\prime },y_s+Y}\rangle \langle {{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t,x,y_s}{{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t^{\prime },x^{\prime },y_s+Y}\rangle \hfill \\ \hfill & e^{j(\omega (t-t^{\prime })-k_x(x-x^{\prime }))},\hfill \end{array}$$(B.7)

and

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_c(y_s)\equiv \frac{1}{T_{\mathrm{obs}}{X}_{\mathrm{obs}}}{\displaystyle \int }\mathrm{d}Y\mathrm{d}t\mathrm{d}{t}^{\prime }\mathrm{d}x\mathrm{d}{x}^{\prime }{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}}\hfill \\ \hfill & \langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t,x,y_s}{{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t^{\prime },x^{\prime },y_s+Y}\rangle \langle {{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t,x,y_s}{{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}|}_{t^{\prime },x^{\prime },y_s+Y}\rangle \hfill \\ \hfill & e^{j(\omega (t-t^{\prime })-k_x(x-x^{\prime }))}.\hfill \end{array}$$(B.8)

The first integral ℐ_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,

$$\begin{array}{c}\hfill T\equiv t,\tau \equiv t^{\prime }-t,X\equiv x,\delta x\equiv x^{\prime }-x,\end{array}$$(B.9)

so that

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_b(y_s)=\frac{1}{T_{\mathrm{obs}}{X}_{\mathrm{obs}}}{\displaystyle \int }\mathrm{d}Y\mathrm{d}T\mathrm{d}\tau \mathrm{d}X\mathrm{d}\delta x{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}}\hfill \\ \hfill & \langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}|}_{T,X,y_s}{{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}|}_{T+\tau ,X+\delta x,y_s+Y}\rangle \hfill \\ \hfill & \langle {{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{T,X,y_s}{{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{T+\tau ,X+\delta x,y_s+Y}\rangle e^{j(k_x\delta x-\omega \tau )},\hfill \end{array}$$(B.10)

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_c(y_s)=\frac{1}{T_{\mathrm{obs}}{X}_{\mathrm{obs}}}{\displaystyle \int }\mathrm{d}Y\mathrm{d}T\mathrm{d}\tau \mathrm{d}X\mathrm{d}\delta x{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}}\hfill \\ \hfill & \langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}|}_{T,X,y_s}{{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{T+\tau ,X+\delta x,y_s+Y}\rangle \hfill \\ \hfill & \langle {{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{T,X,y_s}{{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}|}_{T+\tau ,X+\delta x,y_s+Y}\rangle e^{j(k_x\delta x-\omega \tau )}.\hfill \end{array}$$(B.11)

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

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_b={\displaystyle \int }\mathrm{d}\tau \mathrm{d}\delta x\mathrm{d}Y{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}\langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}|}_{0,0,0}{{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}|}_{\tau ,\delta x,Y}\rangle }\hfill \\ \hfill & \langle {{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{0,0,0}{{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{\tau ,\delta x,Y}\rangle e^{j(k_x\delta x-\omega \tau )}\hfill \end{array}$$(B.12)

and

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_c={\displaystyle \int }\mathrm{d}\tau \mathrm{d}\delta x\mathrm{d}Y{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}\langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}|}_{0,0,0}{{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{\tau ,\delta x,Y}\rangle }\hfill \\ \hfill & \langle {{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{0,0,0}{{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}|}_{\tau ,\delta x,Y}\rangle e^{j(k_x\delta x-\omega \tau )}.\hfill \end{array}$$(B.13)

If we introduce the following functions,

$$\begin{array}{cc}& f_{b,ik}(\tau ,\delta \mathbf{x})\equiv \langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}|}_{0,\mathbf{0}}{{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}|}_{\tau ,\delta \mathbf{x}}\rangle ,\hfill \\ \hfill & g_{b,jl}(\tau ,\delta \mathbf{x})\equiv \langle {{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{0,\mathbf{0}}{{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{\tau ,\delta \mathbf{x}}\rangle ,\hfill \\ \hfill & f_{c,il}(\tau ,\delta \mathbf{x})\equiv \langle {{\partial }_i{\mathrm{\Psi }}_{\mathrm{turb}}|}_{0,\mathbf{0}}{{\partial }_l\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{\tau ,\delta \mathbf{x}}\rangle ,\hfill \\ \hfill & g_{c,jk}(\tau ,\delta \mathbf{x})\equiv \langle {{\partial }_j\mathrm{\Delta }{\mathrm{\Psi }}_{\mathrm{turb}}|}_{0,\mathbf{0}}{{\partial }_k{\mathrm{\Psi }}_{\mathrm{turb}}|}_{\tau ,\delta \mathbf{x}}\rangle ,\hfill \end{array}$$(B.14)

where δx ≡ δxe_x + Ye_y, then these two integrals can be rewritten more compactly as

$$\begin{array}{cc}& {\mathcal{I}}_b=\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}\stackrel{\sim }{f_{b,ik}{g}_{b,jl}}(-\omega ,-k_x,0),\hfill \end{array}$$(B.15)

$$\begin{array}{cc}& {\mathcal{I}}_c=\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}\stackrel{\sim }{f_{c,il}{g}_{c,jk}}(-\omega ,-k_x,0),\hfill \end{array}$$(B.16)

where the notation $\stackrel{\sim }{.}$ denotes the Fourier transform in (t, x, y), and is defined by

$$\begin{array}{c}\hfill {\displaystyle \stackrel{\sim }{f}(\omega ,k_x,k_y)\equiv \frac{1}{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}}{\displaystyle \int }\mathrm{d}t\mathrm{d}x\mathrm{d}yf(t,x,y)e^{j(\omega t-k_{x}x-k_{y}y)}.}\end{array}$$(B.17)

Then, expanding the Fourier transform of the products as convolution integrals, and exploiting the fact that $\stackrel{\sim }{g}(-\omega ,-\mathbf{k})={\stackrel{\sim }{g}}^{*}(\omega ,\mathbf{k})$ we obtain

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_b=\frac{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}{{(2\pi )}^3}{\displaystyle \int }\mathrm{d}{\omega }^{\prime }{\mathrm{d}}^2{\mathbf{k}}^{\prime }{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}\stackrel{\sim }{f_{b,ik}}({\omega }^{\prime },{\mathbf{k}}^{\prime })}\hfill \\ \hfill & {\stackrel{\sim }{g_{b,jl}}}^{\ast }({\omega }^{\prime }+\omega ,{\mathbf{k}}^{\prime }+\mathbf{k}),\hfill \end{array}$$(B.18)

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_c=\frac{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}{{(2\pi )}^3}{\displaystyle \int }\mathrm{d}{\omega }^{\prime }{\mathrm{d}}^2{\mathbf{k}}^{\prime }{ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}\stackrel{\sim }{f_{c,il}}({\omega }^{\prime },{\mathbf{k}}^{\prime })}\hfill \\ \hfill & {\stackrel{\sim }{g_{c,jk}}}^{\ast }({\omega }^{\prime }+\omega ,{\mathbf{k}}^{\prime }+\mathbf{k}),\hfill \end{array}$$(B.19)

where k ≡ k_xe_x.

Next, we need to express the quantities $\stackrel{\sim }{f_{b,ik}}$, $\stackrel{\sim }{g_{b,jl}}$, $\stackrel{\sim }{f_{c,il}}$ and $\stackrel{\sim }{g_{c,jk}}$ appearing in these integrals. Taking the Fourier transforms of each of Eq. B.14, we find

$$\begin{array}{cc}& {\displaystyle \stackrel{\sim }{f_{b,ik}}({\omega }^{\prime },{\mathbf{k}}^{\prime })=\frac{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}}{{(2\pi )}^3}{\displaystyle \int }\mathrm{d}{\omega }^{″}{\mathrm{d}}^2\mathbf{k}″(-k_i^{″}{k}_k^{\prime })}\hfill \\ \hfill & \langle {\stackrel{\sim }{\mathrm{\Psi }}}_{\mathrm{turb}}({\omega }^{″},\mathbf{k}″){\stackrel{\sim }{\mathrm{\Psi }}}_{\mathrm{turb}}({\omega }^{\prime },{\mathbf{k}}^{\prime })\rangle ,\hfill \end{array}$$(B.20)

$$\begin{array}{cc}& {\displaystyle \stackrel{\sim }{g_{b,jl}}({\omega }^{\prime },{\mathbf{k}}^{\prime })=\frac{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}}{{(2\pi )}^3}{\displaystyle \int }\mathrm{d}{\omega }^{″}{\mathrm{d}}^2\mathbf{k}″(-k_j^{″}{k}^{″2}{k}_l^{\prime }{k}^{\prime 2})}\hfill \\ \hfill & \langle {\stackrel{\sim }{\mathrm{\Psi }}}_{\mathrm{turb}}({\omega }^{″},\mathbf{k}″){\stackrel{\sim }{\mathrm{\Psi }}}_{\mathrm{turb}}({\omega }^{\prime },{\mathbf{k}}^{\prime })\rangle ,\hfill \end{array}$$(B.21)

$$\begin{array}{cc}& {\displaystyle \stackrel{\sim }{f_{c,il}}({\omega }^{\prime },{\mathbf{k}}^{\prime })=\frac{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}}{{(2\pi )}^3}{\displaystyle \int }\mathrm{d}{\omega }^{″}{\mathrm{d}}^2\mathbf{k}″(+k_i^{″}{k}_l^{\prime }{k}^{\prime 2})}\hfill \\ \hfill & \langle {\stackrel{\sim }{\mathrm{\Psi }}}_{\mathrm{turb}}({\omega }^{″},\mathbf{k}″){\stackrel{\sim }{\mathrm{\Psi }}}_{\mathrm{turb}}({\omega }^{\prime },{\mathbf{k}}^{\prime })\rangle ,\hfill \end{array}$$(B.22)

and

$$\begin{array}{cc}& {\displaystyle \stackrel{\sim }{g_{c,jk}}({\omega }^{\prime },{\mathbf{k}}^{\prime })=\frac{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}}{{(2\pi )}^3}{\displaystyle \int }\mathrm{d}{\omega }^{″}{\mathrm{d}}^2\mathbf{k}″(+k_j^{″}{k}^{″2}{k}_k^{\prime })}\hfill \\ \hfill & \langle {\stackrel{\sim }{\mathrm{\Psi }}}_{\mathrm{turb}}({\omega }^{″},\mathbf{k}″){\stackrel{\sim }{\mathrm{\Psi }}}_{\mathrm{turb}}({\omega }^{\prime },{\mathbf{k}}^{\prime })\rangle .\hfill \end{array}$$(B.23)

Because we have assumed a homogeneous and stationary turbulence, the correlation products can be rewritten in terms of Dirac distributions, yielding non-zero contributions only if ω″= − ω′ and k″ = k′. We obtain

$$\begin{array}{cc}& \stackrel{\sim }{f_{b,ik}}({\omega }^{\prime },{\mathbf{k}}^{\prime })=\frac{1}{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}}{k}_i^{\prime }{k}_k^{\prime }{\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime },\mathbf{k}\prime ),\hfill \end{array}$$(B.24)

$$\begin{array}{cc}& \stackrel{\sim }{g_{b,jl}}({\omega }^{\prime },{\mathbf{k}}^{\prime })=\frac{1}{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}}{k}_j^{\prime }{k}_l^{\prime }{k}^{\prime 4}{\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime },\mathbf{k}\prime ),\hfill \end{array}$$(B.25)

$$\begin{array}{cc}& \stackrel{\sim }{f_{c,il}}({\omega }^{\prime },{\mathbf{k}}^{\prime })=-\frac{1}{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}}{k}_i^{\prime }{k}_l^{\prime }{k}^{\prime 2}{\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime },\mathbf{k}\prime ),\hfill \end{array}$$(B.26)

$$\begin{array}{cc}& \stackrel{\sim }{g_{c,jk}}({\omega }^{\prime },{\mathbf{k}}^{\prime })=-\frac{1}{\sqrt{T_{\mathrm{obs}}{X}_{\mathrm{obs}}{Y}_{\mathrm{obs}}}}{k}_j^{\prime }{k}_k^{\prime }{k}^{\prime 2}{\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime },\mathbf{k}\prime ).\hfill \end{array}$$(B.27)

The function ℰ_Ψ denotes the stream function turbulent spectrum, defined by

$$\begin{array}{c}\hfill {\displaystyle {\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime },\mathbf{k}\prime )\equiv {\displaystyle \int \mathrm{d}}\tau {\mathrm{d}}^2\mathbf{x}\langle {\mathrm{\Psi }}_{\mathrm{turb}}(T,\mathbf{X}){\mathrm{\Psi }}_{\mathrm{turb}}(T+\tau ,\mathbf{X}+\mathbf{x})\rangle e^{j({\omega }^{\prime }\tau -\mathbf{k}\prime ·\mathbf{x})}.}\end{array}$$(B.28)

Plugging these into Eq. B.18 and Eq. B.19, we find

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_b=\frac{1}{{(2\pi )}^3}{\displaystyle \int \mathrm{d}}{\omega }^{\prime }{\mathrm{d}}^2\mathbf{k}\prime {ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}{k}_i^{\prime }(k_j^{\prime }+k_j)k_k^{\prime }(k_l^{\prime }+k_l)}\hfill \\ \hfill & {(\mathbf{k}\prime +\mathbf{k})}^4{\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime },\mathbf{k}\prime ){{\mathcal{E}}_{\mathrm{\Psi }}}^{\ast }({\omega }^{\prime }+\omega ,\mathbf{k}\prime +\mathbf{k}),\hfill \end{array}$$(B.29)

and

$$\begin{array}{cc}& {\displaystyle {\mathcal{I}}_c=\frac{1}{{(2\pi )}^3}{\displaystyle \int \mathrm{d}}{\omega }^{\prime }{\mathrm{d}}^2\mathbf{k}\prime {ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}{k}_i^{\prime }(k_j^{\prime }+k_j)(k_k^{\prime }+k_k)k_l^{\prime }{k}^{\prime 2}}\hfill \\ \hfill & {(\mathbf{k}\prime +\mathbf{k})}^2{\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime },\mathbf{k}\prime ){{\mathcal{E}}_{\mathrm{\Psi }}}^{\ast }({\omega }^{\prime }+\omega ,\mathbf{k}\prime +\mathbf{k}).\hfill \end{array}$$(B.30)

Flipping the indices k and l in the second integral, and exploiting the fact that ϵ_klz = −ϵ_lkz, the sum of these two integrals becomes

$$\begin{array}{cc}& {\displaystyle \mathcal{I}=\frac{1}{{(2\pi )}^3}{\displaystyle \int \mathrm{d}}{\omega }^{\prime }{\mathrm{d}}^2\mathbf{k}\prime {ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}{k}_i^{\prime }(k_j^{\prime }+k_j)k_k^{\prime }(k_l^{\prime }+k_l)}\hfill \\ \hfill & {(\mathbf{k}\prime +\mathbf{k})}^2({(\mathbf{k}\prime +\mathbf{k})}^2-{\mathbf{k}}^{\prime 2}){\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime },\mathbf{k}\prime ){{\mathcal{E}}_{\mathrm{\Psi }}}^{\ast }({\omega }^{\prime }+\omega ,\mathbf{k}\prime +\mathbf{k}).\hfill \end{array}$$(B.31)

Finally, we explicitly expanded the index contractions. A basic property of the Levi-Civita symbol is

$$\begin{array}{c}\hfill {ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}={\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}{\delta }_{\mathit{zz}}+{\delta }_{\mathit{il}}{\delta }_{\mathit{jz}}{\delta }_{\mathit{kz}}+{\delta }_{\mathit{iz}}{\delta }_{\mathit{jk}}{\delta }_{\mathit{lz}}-{\delta }_{\mathit{ik}}{\delta }_{\mathit{jz}}{\delta }_{\mathit{lz}}-{\delta }_{\mathit{iz}}{\delta }_{\mathit{jl}}{\delta }_{\mathit{kz}}-{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}}{\delta }_{\mathit{zz}}.\end{array}$$(B.32)

Seeing as neither k′ nor k has a component along the z-axis, we have

$$\begin{array}{cc}\hfill {ϵ}_{\mathit{ijz}}{ϵ}_{\mathit{klz}}{k}_i^{\prime }(k_j^{\prime }+k_j)k_k^{\prime }(k_l^{\prime }+k_l)& =({\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}-{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}})k_i^{\prime }(k_j^{\prime }+k_j)k_k^{\prime }(k_l^{\prime }+k_l)\hfill \\ \hfill & =k^{\prime 2}{(\mathbf{k}\prime +\mathbf{k})}^2-{(\mathbf{k}\prime ·(\mathbf{k}\prime +\mathbf{k}))}^2\hfill \\ \hfill & =k^{\prime 2}{k}^2-{(\mathbf{k}\prime ·\mathbf{k})}^2.\hfill \end{array}$$(B.33)

Thus we obtain

$$\begin{array}{cc}& {\displaystyle \mathcal{I}=\frac{1}{{(2\pi )}^3}{\displaystyle \int \mathrm{d}}{\omega }^{\prime }{\mathrm{d}}^2\mathbf{k}\prime [k^{\prime 2}{k}^2-{(\mathbf{k}\prime ·\mathbf{k})}^2]{(\mathbf{k}\prime +\mathbf{k})}^2}\hfill \\ \hfill & ({(\mathbf{k}\prime +\mathbf{k})}^2-k^{\prime 2}){\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime },\mathbf{k}\prime ){{\mathcal{E}}_{\mathrm{\Psi }}}^{\ast }({\omega }^{\prime }+\omega ,\mathbf{k}\prime +\mathbf{k}).\hfill \end{array}$$(B.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 ℐ as a function of the stream function turbulent spectrum, ℰ_Ψ,

$$\begin{array}{cc}& {\displaystyle \mathcal{I}=\frac{1}{4{\pi }^3}{\displaystyle \int \mathrm{d}}{\omega }^{\prime }{\mathrm{d}}^2\mathbf{k}\prime k_x^3{k}_x^{\prime }{k}_y^{\prime 2}{|\mathbf{k}\prime +\frac{\mathbf{k}}{2}|}^2}\hfill \\ \hfill & {\mathcal{E}}_{\mathrm{\Psi }}({\omega }^{\prime }-\frac{\omega }{2},\mathbf{k}\prime -\frac{\mathbf{k}}{2}){{\mathcal{E}}_{\mathrm{\Psi }}}^{\ast }({\omega }^{\prime }+\frac{\omega }{2},\mathbf{k}\prime +\frac{\mathbf{k}}{2}).\hfill \end{array}$$(B.35)

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

$$\begin{array}{c}\hfill T_n(\xi )=cos(narccos\xi ).\end{array}$$(C.1)

Using the definition of the linear operator, we can write

$$\begin{array}{cc}\hfill {\mathcal{M}}_{\mathit{ij}}=& -(\omega k_x^2+k_x(\beta -U^{″})+j{\nu }_{\mathrm{turb}}{k}_x^4){\mathcal{M}}_{\mathit{ij}}^{(1)}\hfill \end{array}$$(C.2)

$$\begin{array}{cc}& +(\omega +2j{\nu }_{\mathrm{turb}}{k}_x^2){\mathcal{M}}_{\mathit{ij}}^{(2)}\hfill \end{array}$$(C.3)

$$\begin{array}{cc}& -j{\nu }_{\mathrm{turb}}{\mathcal{M}}_{\mathit{ij}}^{(3)}+k_x^3{\mathcal{M}}_{\mathit{ij}}^{(4)}-k_x{\mathcal{M}}_{\mathit{ij}}^{(5)},\hfill \end{array}$$(C.4)

where

$$\begin{array}{cc}& {\mathcal{M}}_{\mathit{ij}}^{(1)}\equiv \frac{2}{\pi c_i}{T}_i|T_j={\delta }_{\mathit{ij}},\hfill \end{array}$$(C.5)

$$\begin{array}{cc}& {\mathcal{M}}_{\mathit{ij}}^{(2)}\equiv \frac{2}{\pi c_i}{T}_i|T_j^{″},\hfill \end{array}$$(C.6)

$$\begin{array}{cc}& {\mathcal{M}}_{\mathit{ij}}^{(3)}\equiv \frac{2}{\pi c_i}{T}_i|T_j^{″″},\hfill \end{array}$$(C.7)

$$\begin{array}{cc}& {\mathcal{M}}_{\mathit{ij}}^{(4)}\equiv \frac{2}{\pi c_i}{T}_i|U{T}_j,\hfill \end{array}$$(C.8)

$$\begin{array}{cc}& {\mathcal{M}}_{\mathit{ij}}^{(5)}\equiv \frac{2}{\pi c_i}{T}_i|U{T}_j^{″},\hfill \end{array}$$(C.9)

and $T_j^{″}$ and T_j⁗ denote the second and fourth derivatives with respect to ξ, respectively. We derive ${\mathcal{M}}_{\mathit{ij}}^{(2)}$ and ${\mathcal{M}}_{\mathit{ij}}^{(3)}$ in Sect. C.1, and ${\mathcal{M}}_{\mathit{ij}}^{(4)}$ and ${\mathcal{M}}_{\mathit{ij}}^{(5)}$ in Sect. C.2.

C.1. Projection of the derivatives

Given a function

$$\begin{array}{c}\hfill f(\xi )={\displaystyle \sum }{\lambda }_i{T}_i(\xi ),\end{array}$$(C.10)

we can write, for any positive integer, p,

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^{p}f}{\mathrm{d}{y}^p}={\displaystyle \sum }_{i=0}^{+\infty }{\lambda }_i^{(p)}{T}_i,{\lambda }_i^{(0)}={\lambda }_i.\end{array}$$(C.11)

In order to compute the coefficients ${\mathcal{M}}_{\mathit{ij}}^{(2)}$ and ${\mathcal{M}}_{\mathit{ij}}^{(3)}$, we needed to find the recurrence relation between the ${\lambda }_i^{(p)}$ and the ${\lambda }_i^{(p-1)}$, that is, a recurrence relation for all the derivatives of f.

For any p ≥ 1 we can write

$$\begin{array}{c}\hfill {\displaystyle \sum _{i=0}^{+\infty }}{\lambda }_i^{(p)}{T}_i=\frac{\mathrm{d}}{\mathrm{d}y}{\displaystyle \sum _{i=0}^{+\infty }}{\lambda }_i^{(p-1)}{T}_i.\end{array}$$(C.12)

On the other hand, from Eq. C.1, it is easily seen that

$$\begin{array}{c}\hfill T_i=\frac{c_i}{2(i+1)}{T}_{i+1}\prime -\frac{c_i^{\prime }}{2(i-1)}{T}_{i-1}\prime ,\end{array}$$(C.13)

where we recall that c_i = 1 + δ_i0, and we have also introduced a new factor $c_i^{\prime }$ (defined by $c_0^{\prime }=c_1^{\prime }=0$ and $c_{i\ge 2}^{\prime }=1$). Therefore, we can also write

$$\begin{array}{cc}\hfill {\displaystyle \sum _{i=0}^{+\infty }}{\lambda }_i^{(p)}{T}_i& =\frac{\mathrm{d}}{\mathrm{d}y}{\displaystyle \sum _{i=0}^{+\infty }}{\lambda }_i^{(p)}(\frac{c_i}{2(i+1)}{T}_{i+1}-\frac{c_i^{\prime }}{2(i-1)}{T}_{i-1})\hfill \\ \hfill & =\frac{\mathrm{d}}{\mathrm{d}y}{\displaystyle \sum _{i=1}^{+\infty }}(\frac{c_{i-1}}{2i}{\lambda }_{i-1}^{(p)}-\frac{1}{2i}{\lambda }_{i+1}^{(p)})T_i.\hfill \end{array}$$(C.14)

Equating the T_i coefficients in Eqs. C.12 and C.14, we find

$$\begin{array}{c}\hfill c_{i-1}{\lambda }_{i-1}^{(p)}-{\lambda }_{i+1}^{(p)}=2i{\lambda }_i^{(p-1)},\end{array}$$(C.15)

which constitutes a recurrence relation in terms of p. This recurrence relation, combined with the condition that we must have $\underset{i\to +\infty }{\lim }{\lambda }_i^{(p)}=0$ for any value of p, is easily solved to yield

$$\begin{array}{c}\hfill c_i{\lambda }_i^{(p)}=2{\displaystyle \sum _{\begin{array}{c}j=i+1\\ i+j\equiv 1[2]\end{array}}^{+\infty }}j{\lambda }_j^{(p-1)}.\end{array}$$(C.16)

where the notation ≡a[b] means ‘equal to a modulo b’.

For example, the coefficients of the first derivative straightforwardly read

$$\begin{array}{c}\hfill c_i{\lambda }_i^{(1)}=2{\displaystyle \sum _{\begin{array}{c}j=i+1\\ i+j\equiv 1[2]\end{array}}^{+\infty }}j{\lambda }_j.\end{array}$$(C.17)

The coefficients of the second derivative can be computed in the following way

$$\begin{array}{cc}\hfill c_i{\lambda }_i^{(2)}& =2{\displaystyle \sum _{\begin{array}{c}j=i+1\\ i+j\equiv 1[2]\end{array}}^{+\infty }}j{\lambda }_j^{(1)}=4{\displaystyle \sum _{\begin{array}{c}j=i+1\\ i+j\equiv 1[2]\end{array}}^{+\infty }}j{\displaystyle \sum _{\begin{array}{c}k=j+1\\ k+j\equiv 1[2]\end{array}}^{+\infty }}k{\lambda }_k\hfill \\ & =4{\displaystyle \sum _{\begin{array}{c}k=i+2\\ k+i\equiv 0[2]\end{array}}^{+\infty }}k{\lambda }_k{\displaystyle \sum _{\begin{array}{c}j=i+1\\ j+k\equiv 1[2]\end{array}}^{k-1}}j={\displaystyle \sum _{\begin{array}{c}k=i+2\\ k+i\equiv 0[2]\end{array}}^{+\infty }}k(k^2-i^2){\lambda }_k.\hfill \end{array}$$(C.18)

Similarly, we find

$$\begin{array}{cc}& c_i{\lambda }_i^{(3)}=\frac{1}{4}{\displaystyle \sum _{\begin{array}{c}k=i+3\\ k+i\equiv 1[2]\end{array}}^{+\infty }}k[k^2(k^2-2)-2k^2{i}^2+i^2(i^2-2)+1]{\lambda }_k,\hfill \end{array}$$(C.19)

$$\begin{array}{cc}& c_i{\lambda }_i^{(4)}=\frac{1}{24}{\displaystyle \sum _{\begin{array}{c}k=i+4\\ k+i\equiv 0[2]\end{array}}^{+\infty }}k[k^2{(k^2-4)}^2-3k^4{i}^2+3k^2{i}^4-i^2{(i^2-4)}^2]{\lambda }_k.\hfill \end{array}$$(C.20)

From these expressions it directly follows that

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathit{ij}}^{(2)}=\{\begin{array}{cc}\frac{1}{c_i}j(j^2-i^2)\hfill & \mathrm{if}j⩾i+2\mathrm{and}i+j\equiv 0[2]\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\end{array}$$(C.21)

and

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathit{ij}}^{(3)}=\{\begin{array}{c}\frac{1}{24c_i}j[j^2{(j^2-4)}^2-3j^4{i}^2+3j^2{i}^4-i^2{(i^2-4)}^2]\hfill \\ \mathrm{if}j⩾i+4\mathrm{and}i+j\equiv 0[2]\hfill & \\ 0\mathrm{otherwise}\hfill & .\hfill \end{array}\end{array}$$(C.22)

C.2. Projection of the products

In order to compute the coefficients ${\mathcal{M}}_{\mathit{ij}}^{(4)}$ and ${\mathcal{M}}_{\mathit{ij}}^{(5)}$, 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

$$\begin{array}{cc}& f(\xi )={\displaystyle \sum _{i=0}^{+\infty }}{\lambda }_i{T}_i(\xi ),\hfill \end{array}$$(C.23)

$$\begin{array}{cc}& g(\xi )={\displaystyle \sum _{i=0}^{+\infty }}{\mu }_i{T}_i(\xi ).\hfill \end{array}$$(C.24)

It will be useful, in the following, to introduce the following functions

$$\begin{array}{c}\hfill {\stackrel{\sim }{T}}_n(\xi )\equiv e^{jnarccos\xi },\end{array}$$(C.25)

so that we have

$$\begin{array}{c}\hfill 2T_i={\stackrel{\sim }{T}}_i+{\stackrel{\sim }{T}}_{-i}.\end{array}$$(C.26)

It directly follows that

$$\begin{array}{cc}& 2f={\displaystyle \sum _{i=-\infty }^{+\infty }}{c}_{|i|}{\lambda }_{|i|}{\stackrel{\sim }{T}}_i,\hfill \end{array}$$(C.27)

$$\begin{array}{cc}& 2g={\displaystyle \sum _{i=-\infty }^{+\infty }}{c}_{|i|}{\mu }_{|i|}{\stackrel{\sim }{T}}_i.\hfill \end{array}$$(C.28)

Forming the product of these two expressions, and using the identity ${\stackrel{\sim }{T}}_i{\stackrel{\sim }{T}}_j={\stackrel{\sim }{T}}_{i+j}$, we obtain

$$\begin{array}{c}\hfill 4fg={\displaystyle \sum _{i=-\infty }^{+\infty }\sum _{j=-\infty }^{+\infty }}{c}_{|i|}{c}_{|j|}{\lambda }_{|i|}{\mu }_{|j|}{\stackrel{\sim }{T}}_{i+j},\end{array}$$(C.29)

which can be rewritten as

$$\begin{array}{c}\hfill 4fg={\displaystyle \sum _{i=-\infty }^{+\infty }}\left({\displaystyle \sum _{j=-\infty }^{+\infty }}{c}_{|j|}{c}_{|i-j|}{\lambda }_{|j|}{\mu }_{|i-j|}\right){\stackrel{\sim }{T}}_i.\end{array}$$(C.30)

Finally, using Eq. C.26 to revert back to an expansion on the T_i, one finally finds

$$\begin{array}{c}\hfill fg={\displaystyle \sum _{i=0}^{+\infty }}\left(\frac{1}{2c_i}{\displaystyle \sum _{j=-\infty }^{+\infty }}{c}_{|j|}{c}_{|i-j|}{\lambda }_{|j|}{\mu }_{|i-j|}\right)T_i.\end{array}$$(C.31)

The differential rotation profile in the expressions of ${\mathcal{M}}_{\mathit{ij}}^{(4)}$ and ${\mathcal{M}}_{\mathit{ij}}^{(5)}$ is given by $U(\xi )=-\overline{U}{\xi }^2$, which is straightforwardly expanded as

$$\begin{array}{c}\hfill U=-\frac{1}{2}\overline{U}(T_0+T_2),\end{array}$$(C.32)

so that

$$\begin{array}{c}\hfill fU=-\frac{\overline{U}}{4}{\displaystyle \sum _{i=0}^{+\infty }}\frac{1}{c_i}({\lambda }_{i+2}+2c_i{\lambda }_i+c_{|i-2|}{\lambda }_{|i-2|})T_i.\end{array}$$(C.33)

From this expression directly stem the following

$$\begin{array}{cc}& {\mathcal{M}}_{\mathit{ij}}^{(4)}=-\frac{\overline{U}}{4c_i}({\delta }_{i+2,j}+2c_i{\delta }_{\mathit{ij}}+c_{|i-2|}{\delta }_{|i-2|,j}),\hfill \end{array}$$(C.34)

$$\begin{array}{cc}& {\mathcal{M}}_{\mathit{ij}}^{(5)}=-\frac{\overline{U}}{4c_i}({\mathcal{M}}_{i+2,j}^{(2)}+2c_i{\mathcal{M}}_{\mathit{ij}}^{(2)}+c_{|i-2|}{\mathcal{M}}_{|i-2|,j}^{(2)}).\hfill \end{array}$$(C.35)

All Figures