Open Access
Issue
A&A
Volume 673, May 2023
Article Number A124
Number of page(s) 19
Section The Sun and the Heliosphere
DOI https://doi.org/10.1051/0004-6361/202245666
Published online 23 May 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model.

Open access funding provided by Max Planck Society.

1. Introduction

Multiple types of waves can propagate in the interior of a star. In the case of a 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

x R ϕ , $$ \begin{aligned}&x \equiv R\phi , \end{aligned} $$(1)

y R sin λ , $$ \begin{aligned}&{ y} \equiv R\sin \lambda , \end{aligned} $$(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 (vx and vy). The equations of motion in the rotating frame and in the inviscid limit become

v x t + ( v · ) v x = 1 ρ p x + f v y , $$ \begin{aligned}&\frac{\partial { v}_x}{\partial t} + (\mathbf v \cdot \boldsymbol{\nabla }) { v}_x = -\dfrac{1}{\rho } \dfrac{\partial p}{\partial x} + f { v}_{ y}, \end{aligned} $$(3)

v y t + ( v · ) v y = 1 ρ p y f v x , $$ \begin{aligned}&\dfrac{\partial { v}_{ y}}{\partial t} + (\mathbf v \cdot \boldsymbol{\nabla }) { v}_{ y} = -\dfrac{1}{\rho } \dfrac{\partial p}{\partial { y}} - f { v}_x, \end{aligned} $$(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) ex, 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

u = ( Ψ e z ) , $$ \begin{aligned} \mathbf u = \boldsymbol{\nabla } \wedge \left( \Psi \mathbf e_z \right), \end{aligned} $$(5)

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

( t + U x ) Δ Ψ osc + ( β U ) Ψ osc x ν turb Δ 2 Ψ osc = δ ( Ψ turb x Δ Ψ turb y Ψ turb y Δ Ψ turb x ) , $$ \begin{aligned}&\left( \dfrac{\partial }{\partial t} + U\dfrac{\partial }{\partial x}\right) \Delta {\Psi _{\rm osc}} + (\beta - U^{\prime \prime }) \dfrac{\partial {\Psi _{\rm osc}}}{\partial x} - {\nu _{\rm turb}} \Delta ^2 {\Psi _{\rm osc}} \nonumber \\&\quad = \delta \left( \dfrac{\partial {\Psi _{\rm turb}}}{\partial x}\dfrac{\partial \Delta {\Psi _{\rm turb}}}{\partial { y}} - \dfrac{\partial {\Psi _{\rm turb}}}{\partial { y}}\dfrac{\partial \Delta {\Psi _{\rm turb}}}{\partial x} \right), \end{aligned} $$(6)

where U″ is the second derivative of U(y), Δ= x 2 + y 2 $ \Delta = \partial_x^2 + \partial_{\it 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  −  ⟨qh 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 − jkxx), where j denotes the imaginary unit. We used the following convention for the Fourier transform,

f ^ ( ω , k x , y ) 1 T obs X obs d t d x f ( t , x , y ) e j ( ω t k x x ) , $$ \begin{aligned} \widehat{f}(\omega , k_x, { y}) \equiv \dfrac{1}{\sqrt{T_{\rm obs} X_{\rm obs}}} \displaystyle \int \mathrm{d} t \mathrm{d} x ~ f(t, x, { y}) e^{{j}(\omega t - k_x x)}, \end{aligned} $$(7)

where Tobs and Xobs are respectively the t and x windows over which the integral defining the Fourier transform is computed. In the Fourier domain, Eq. (6) becomes

L Ψ ^ osc = S ^ . $$ \begin{aligned} \mathcal{L} \ {\widehat{\Psi }_{\rm osc}} = \widehat{S}. \end{aligned} $$(8)

The notation S ^ $ \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

L ( ω k x U ) Δ ^ k x ( β U ) j ν turb Δ ^ 2 , $$ \begin{aligned} \mathcal{L} \equiv (\omega - k_x U) \widehat{\Delta } - k_x (\beta - U^{\prime \prime }) - {j} {\nu _{\rm turb}} \widehat{\Delta }^2, \end{aligned} $$(9)

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

The solution Ψ ^ osc $ {\widehat{\Psi}_{\mathrm{osc}}} $ to Eq. (8) can be expressed in terms of the Green function G(y, ys), which is defined as the solution of the following differential equation (and with the same boundary conditions as the full linear problem),

L G ( y , y s ) = δ ( y y s ) , $$ \begin{aligned} \mathcal{L} \ G({ y},{ y}_s) = \delta ({ y} - { y}_s), \end{aligned} $$(10)

where δ is the Dirac distribution, so that

Ψ ^ osc ( y ) = R R d y s G ( y , y s ) S ^ ( y s ) . $$ \begin{aligned} {\widehat{\Psi }_{\rm osc}}({ y}) = \displaystyle \int _{-R}^{R} \mathrm{d} { y}_s ~ G({ y},{ y}_s) ~ \widehat{S}({ y}_s). \end{aligned} $$(11)

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

u ^ x , osc ( y ) = R R d y s G y S ^ ( y s ) , $$ \begin{aligned} {\widehat{u}_{x,\mathrm{osc}}}({ y})&= \displaystyle \int _{-R}^R \mathrm{d} { y}_s ~ \dfrac{\partial G}{\partial { y}} ~ \widehat{S}({ y}_s), \end{aligned} $$(12)

u ^ y , osc ( y ) = j k x R R d y s G ( y , y s ) S ^ ( y s ) , $$ \begin{aligned} {\widehat{u}_{{ y},\mathrm{osc}}}({ y})&= -{j} k_x \displaystyle \int _{-R}^R \mathrm{d} { }{ y}_s ~ G({ y},{ y}_s) ~ \widehat{S}({ y}_s), \end{aligned} $$(13)

ζ ^ osc ( y ) = R R d y s ( k x 2 G 2 G y 2 ) S ^ ( y s ) . $$ \begin{aligned} {\widehat{\zeta }_{\rm osc}}({ y})&= \displaystyle \int _{-R}^R \mathrm{d} { y}_s ~ \left(k_x^2 G - \dfrac{\partial ^2 G}{\partial { y}^2} \right)~ \widehat{S}({ y}_s). \end{aligned} $$(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

| u ^ x , osc ( y obs ) | 2 = R R d y s | G y obs | 2 I ( y s ) , $$ \begin{aligned} \left\langle \left| {\widehat{u}_{x,\mathrm{osc}}}({{ y}_{\rm obs}}) \right|^2 \right\rangle&= \displaystyle \int _{-R}^R \mathrm{d} { y}_s ~ \left| \dfrac{\partial G}{\partial {{ y}_{\rm obs}}} \right|^2 \mathcal{I} ({ y}_s), \end{aligned} $$(15)

| u ^ y , osc ( y obs ) | 2 = R R d y s | k x G ( y obs , y s ) | 2 I ( y s ) , $$ \begin{aligned} \left\langle \left| {\widehat{u}_{{ y},\mathrm{osc}}}({{ y}_{\rm obs}}) \right|^2 \right\rangle&= \displaystyle \int _{-R}^R \mathrm{d} { y}_s ~ \left| k_x G({{ y}_{\rm obs}},{ y}_s) \right|^2 \mathcal{I} ({ y}_s), \end{aligned} $$(16)

| ζ ^ osc ( y obs ) | 2 = R R d y s | ( k x 2 G 2 G y obs 2 ) | 2 I ( y s ) , $$ \begin{aligned} \left\langle \left| {\widehat{\zeta }_{\rm osc}}({{ y}_{\rm obs}}) \right|^2 \right\rangle&= \displaystyle \int _{-R}^R \mathrm{d} { y}_s ~ \left| \left(k_x^2 G - \dfrac{\partial ^2 G}{\partial {{ y}_{\rm obs}}^2} \right) \right|^2 \mathcal{I} ({ y}_s), \end{aligned} $$(17)

where yobs is the latitudinal coordinate at which the power spectrum is evaluated, and ⟨.⟩ denotes an ensemble average. The function ℐ(ys) denotes the source covariance, and is defined by

I ( y s ) d Y S ^ ( y s ) S ^ ( y s + Y ) . $$ \begin{aligned} \mathcal{I} ({ y}_s) \equiv \displaystyle \int \mathrm{d} Y ~ \left\langle \widehat{S}({ y}_s) ~ \widehat{S}^*({ y}_s + Y) \right\rangle . \end{aligned} $$(18)

We recall that the source term S ^ $ \widehat{S} $ depends on ω and kx, 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 ys. Eventually, we find (see Eq. (B.35))

I = 1 4 π 3 d ω d 2 k k x 3 k x k y 2 | k + k / 2 | 2 × E Ψ ( ω ω 2 , k k 2 ) E Ψ ( ω + ω 2 , k + k 2 ) . $$ \begin{aligned}&\mathcal{I} = \dfrac{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} \left| \mathbf{k }^{\prime } + \mathbf k /2 \right|^2 \nonumber \\&\quad \times {\mathcal{E} _\Psi }\left(\omega^\prime - \dfrac{\omega }{2}, \mathbf{k }^{\prime } - \dfrac{\mathbf{k}}{2} \right) {\mathcal{E} _\Psi }^*\left(\omega^\prime + \dfrac{\omega }{2}, \mathbf{k }^{\prime } + \dfrac{\mathbf{k}}{2} \right). \end{aligned} $$(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))

E Ψ ( ω , k ) d τ d 2 x Ψ turb ( T , X ) Ψ turb ( T + τ , X + x ) e j ( ω τ k · x ) . $$ \begin{aligned} {\mathcal{E} _\Psi }(\omega , \mathbf k ) \equiv \displaystyle \int \mathrm{d} \tau \mathrm{d} ^2\mathbf x \left\langle {\Psi _{\rm turb}}(T, \mathbf X ) {\Psi _{\rm turb}}(T+\tau , \mathbf X +\mathbf x ) \right\rangle ~ e^{{j} (\omega \tau - \mathbf k \cdot \mathbf x )}. \end{aligned} $$(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

| u ^ x , turb ( y obs ) | 2 = 1 2 π d k y E Ψ ( ω , k ) k y 2 , $$ \begin{aligned} \left\langle \left| {\widehat{u}_{x,\mathrm{turb}}}({{ y}_{\rm obs}}) \right|^2 \right\rangle&= \dfrac{1}{2\pi } \displaystyle \int \mathrm{d} k_{ y} ~ {\mathcal{E} _\Psi }(\omega , \mathbf k )~ k_{ y}^2, \end{aligned} $$(21)

| u ^ y , turb ( y obs ) | 2 = 1 2 π d k y E Ψ ( ω , k ) k x 2 , $$ \begin{aligned} \left\langle \left| {\widehat{u}_{{ y},\mathrm{turb}}}({{ y}_{\rm obs}}) \right|^2 \right\rangle&= \dfrac{1}{2\pi } \displaystyle \int \mathrm{d} k_{ y} ~ {\mathcal{E} _\Psi }(\omega , \mathbf k )~ k_x^2, \end{aligned} $$(22)

| ζ ^ turb ( y obs ) | 2 = 1 2 π d k y E Ψ ( ω , k ) ( k x 2 + k y 2 ) 2 . $$ \begin{aligned} \left\langle \left| {\widehat{\zeta }_{\rm turb}}({{ y}_{\rm obs}}) \right|^2 \right\rangle&= \dfrac{1}{2\pi } \displaystyle \int \mathrm{d} k_{ y} ~ {\mathcal{E} _\Psi }(\omega , \mathbf k )~ \left(k_x^2+k_{ y}^2 \right)^2. \end{aligned} $$(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, ys) associated with the linear operator ℒ, and (ii) the turbulent stream function spectrum ℰΨ.

2.2. The Green function

The angular frequency, ω, and azimuthal wavenumber, kx, 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,

U ( ξ ) = U ¯ ξ 2 , ξ y / R = sin λ , $$ \begin{aligned} U(\xi ) = -\overline{U} \xi ^2, \qquad \xi \equiv { y} / R = \sin \lambda , \end{aligned} $$(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 U ¯ = 244 $ \overline{U} = 244 $ m s−1 as Gizon et al. (2020), which approximates the solar differential rotation at low latitudes. With this value of Ωeq, the Coriolis parameter becomes β = 8.18 × 10−15 m−1s−1. Finally, the turbulent viscosity, νturb, is specified through the turbulent Reynolds number,

Re turb U ¯ R ν turb , $$ \begin{aligned} {\mathrm{Re} _{\rm turb}} \equiv \dfrac{\overline{U} R}{{\nu _{\rm turb}}}, \end{aligned} $$(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

T n ( ξ ) = cos ( n arccos ξ ) , $$ \begin{aligned} T_n(\xi ) = \cos (n\arccos \xi ), \end{aligned} $$(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,

f | g 1 1 f ( ξ ) g ( ξ ) 1 ξ 2 d ξ , $$ \begin{aligned} {f|g} \equiv \displaystyle \int _{-1}^{1} \dfrac{f(\xi )g(\xi )}{\sqrt{1-\xi ^2}} ~ \mathrm{d} \xi , \end{aligned} $$(27)

in the sense that

T i | T j = π 2 c i δ ij , $$ \begin{aligned} {T_i|T_j} = \dfrac{\pi }{2} c_i \delta _{ij}, \end{aligned} $$(28)

where δij is the Kronecker delta and ci = 1 + δi0. We denote the column vector containing the decomposition of the Green function on the Chebyshev basis as 𝒢(ξs), so that

G ( ξ , ξ s ) = i = 0 G i ( ξ s ) T i ( ξ ) , $$ \begin{aligned} G(\xi ,\xi _s) = \sum _{i=0}^\infty \mathcal{G} _i(\xi _s) T_i(\xi ), \end{aligned} $$(29)

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

M G ( ξ s ) = B ( ξ s ) , $$ \begin{aligned} \mathcal{M} \mathcal{G} (\xi _s) = \mathcal{B} (\xi _s), \end{aligned} $$(30)

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

B i 2 π c i T i | δ ( ξ ξ s ) = 2 T i ( ξ s ) π c i 1 ξ s , $$ \begin{aligned} \mathcal{B} _i \equiv \dfrac{2}{\pi c_i} {T_i|\delta (\xi - \xi _s)} = \dfrac{2T_i(\xi _s)}{\pi c_i \sqrt{1-\xi _s}}, \end{aligned} $$(31)

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

M ij 2 π c i T i | L T j . $$ \begin{aligned} \mathcal{M} _{ij} \equiv \dfrac{2}{\pi c_i} {T_i|\mathcal{L} T_j}. \end{aligned} $$(32)

We note that the factor 2/(πci) 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 kx, 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:

G ( ξ = 1 , ξ s ) = G ( ξ = 1 , ξ s ) = 0 , $$ \begin{aligned} G\left(\xi =-1, \xi _s\right)&= G\left(\xi =1, \xi _s\right) = 0, \end{aligned} $$(33)

G ξ | ξ = 1 , ξ s = G ξ | ξ = 1 , ξ s = 0 . $$ \begin{aligned} \left.\dfrac{\partial G}{\partial \xi } \right|_{\xi =-1,\xi _s}&= \left.\dfrac{\partial G}{\partial \xi } \right|_{\xi =1,\xi _s} = 0. \end{aligned} $$(34)

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

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

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

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

M N 4 , j = 1 , M N 3 , j = ( 1 ) j , $$ \begin{aligned} \mathcal{M} _{N-4,j}&= 1~, \qquad \mathcal{M} _{N-3,j}\,\, = (-1)^j,\end{aligned} $$(37)

M N 2 , j = j 2 , M N 1 , j = ( 1 ) j + 1 j 2 , $$ \begin{aligned} \mathcal{M} _{N-2,j}&= j^2~, \qquad \mathcal{M} _{N-1,j} = (-1)^{j+1} j^2, \end{aligned} $$(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 uturb rather than the stream function

ϕ ij ( ω , k ) d τ d 2 x u i , turb ( T , X ) u j , turb ( T + τ , X + x ) e j ( ω τ k · x ) . $$ \begin{aligned} {\phi _{ij}}(\omega , \mathbf k ) \equiv \displaystyle \int \mathrm{d} \tau \mathrm{d} ^2\mathbf x \left\langle u_{i,\mathrm{turb} }(T, \mathbf X ) u_{j,\mathrm{turb} }(T+\tau , \mathbf X +\mathbf x ) \right\rangle \mathrm{e} ^{{j} (\omega \tau - \mathbf k \cdot \mathbf x )}. \end{aligned} $$(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

ϕ ij = ( k 2 δ ij k i k j ) E Ψ . $$ \begin{aligned} {\phi _{ij}} = \left(k^2 \delta _{ij} - k_i k_j \right) {\mathcal{E} _\Psi }. \end{aligned} $$(40)

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

ϕ ij = E ( k , ω ) ( δ ij k i k j k 2 ) . $$ \begin{aligned} {\phi _{ij}} = \mathrm{E} (\mathbf k , \omega ) \left(\delta _{ij} - \dfrac{k_i k_j}{k^2} \right). \end{aligned} $$(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

E Ψ ( ω , k ) = E ( ω , k ) k 2 , $$ \begin{aligned} {\mathcal{E} _\Psi }(\omega , \mathbf k ) = \dfrac{{E}(\omega , \mathbf k )}{k^2}, \end{aligned} $$(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

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

and then define the temporal part of the spectrum as

χ k ( ω ) E ( k , ω ) E ( k ) . $$ \begin{aligned} \chi _\mathbf{k }(\omega ) \equiv \dfrac{{E}(\mathbf k ,\omega )}{E(\mathbf k )}. \end{aligned} $$(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,

E ( k ) = E iso ( k ) 2 π k . $$ \begin{aligned} E(\mathbf k ) = \dfrac{{E_{\rm iso}}(k)}{2\pi k}. \end{aligned} $$(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 ∼ uk/λk (or, equivalently, angular frequencies around ωk ∼ kuk), where uk 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

χ k ( ω ) = 1 ω k χ ( ω ) , ω ω / ω k , $$ \begin{aligned} \chi _\mathbf{k }(\omega ) = \dfrac{1}{\omega _k} \chi (\widetilde{\omega })~, \qquad \widetilde{\omega } \equiv \omega / \omega _k, \end{aligned} $$(46)

where we have introduced the reduced frequency ω $ \widetilde{\omega} $, and χ is now a function of frequency that does not depend on the eddy size k. Stein (1967) proved that the typical velocity uk associated with these eddies is determined by the spatial spectrum through

u k ( k < k < 2 k d k E iso ( k ) ) 1 / 2 . $$ \begin{aligned} u_k \equiv \left( \displaystyle \int _{k < k^\prime < 2k} \mathrm{d} k^\prime ~ {E_{\rm iso}}(k^\prime ) \right)^{1/2}. \end{aligned} $$(47)

All in all, we can write

E Ψ = E iso ( k ) u k χ ( ω ) 2 π k 2 , $$ \begin{aligned} {\mathcal{E} _\Psi } = \dfrac{{E_{\rm iso}}(k) ~ u_k ~ \chi (\widetilde{\omega })}{2\pi k^2}, \end{aligned} $$(48)

where uk is given by Eq. (47). Describing the whole turbulent spectrum requires only two ingredients, namely (i) the ω-independent spatial spectrum Eiso 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 Eiso(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:

E iso ( k ) = { C 1 ( k k ref ) α 1 if k k 1 , C 2 ( k k ref ) α 2 if k 1 < k < k 2 , C 3 ( k k ref ) α 3 if k 2 k , $$ \begin{aligned} {E_{\rm iso}}(k) = {\left\{ \begin{array}{ll} C_1 \left(\dfrac{k}{k_{\rm ref}}\right)^{\alpha _1}&\mathrm{if} ~~ k \le k_1, \\ C_2 \left(\dfrac{k}{k_{\rm ref}}\right)^{\alpha _2}&\mathrm{if} ~~ k_1 < k < k_2, \\ C_3 \left(\dfrac{k}{k_{\rm ref}}\right)^{\alpha _3}&\mathrm{if} ~~ k_2 \le k, \end{array}\right.} \end{aligned} $$(49)

thumbnail Fig. 1.

Spatial turbulent spectrum and turbulent frequency. Left: Spatial turbulent spectrum, Eiso(k), extracted from solar observations (solid blue line). The spectrum comprises roughly three sections, each with a different power law in kα. The dashed orange line shows the best fit of the data to a three-piece-wise power law. The exponents and the boundaries between the three regimes are indicated on the plot. Right: Turbulent frequency ωk ≡ kuk as a function of wavenumber, k, where uk is the typical velocity of the turbulent eddies of inverse size k, given by Eq. (47).

where the choice of kref is completely arbitrary and can be absorbed in the factors Ci. The spectrum is therefore parameterised by (i) the three exponents α1, 2, 3, (ii) the two wavenumber cutoffs k1 and k2, and (iii) the total turbulent kinetic energy per unit mass ∫Eiso(k)dk. The best fit to the observational data is also shown in the left panel of Fig. 1. We find exponents α1 = 0.27, α2 = −1.88 and α3 = −7.60; in particular, the slope α2 of the middle section is quite close to the −5/3 power law theoretically predicted for the inertial range under the Kolmogorov hypotheses. As for the wavenumber cutoffs – also indicated in the left panel of Fig. 1 – we find k1R = 140 and k2R = 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 = kuk as a function of k, where uk 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 ω $ \widetilde{\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:

χ ( ω ) = A 2 π σ 2 e ω 2 / ( 2 σ 2 ) $$ \begin{aligned} \chi (\widetilde{\omega }) = \dfrac{A}{\sqrt{2\pi \sigma ^2}} e^{-\widetilde{\omega }^2 / (2 \sigma ^2)} \end{aligned} $$(50)

thumbnail Fig. 2.

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

and

χ ( ω ) = A π σ 1 1 + ( ω / σ ) 2 . $$ \begin{aligned} \chi (\widetilde{\omega }) = \dfrac{A}{\pi \sigma } \dfrac{1}{1 + \left( \widetilde{\omega } / \sigma \right)^2}. \end{aligned} $$(51)

In each case, the dimensionless factor σ is introduced to account for the uncertainty on the relation ωk = kuk. 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 uy. This is obtained by setting yobs = 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 kx and angular frequency ω. It is then straightforward to transform it into a power spectral density per unit ω only, for any given azimuthal order m ≡ kxR, 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

L 1 Ψ n = ω n Δ ^ Ψ n , $$ \begin{aligned} \mathcal{L} _1 \Psi _n = \omega _n \widehat{\Delta }\Psi _n, \end{aligned} $$(52)

where

L 1 k x U Δ ^ + k x ( β U ) + j ν turb Δ ^ 2 . $$ \begin{aligned} \mathcal{L} _1 \equiv k_x U \widehat{\Delta } + k_x(\beta - U^{\prime \prime }) + {j}{\nu _{\rm turb}}\widehat{\Delta }^2. \end{aligned} $$(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. 35, 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.

thumbnail Fig. 3.

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

thumbnail Fig. 4.

Same as Fig. 3, but for m = 5.

thumbnail Fig. 5.

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

We also need to compute the complex eigenfrequencies and eigenfunctions of the associated adjoint system ℒΨ = 0, where the adjoint linear operator ℒ is defined in such a way that for any function f and g satisfying the boundary conditions of the problem, we have

1 1 d ξ ( L g ) f = 1 1 d ξ g ( L f ) . $$ \begin{aligned} \displaystyle \int _{-1}^1 \mathrm{d} \xi ~ \left(\mathcal{L} ^\dagger g\right)^*f = \displaystyle \int _{-1}^1 \mathrm{d} \xi ~ g^*\left(\mathcal{L} f\right). \end{aligned} $$(54)

Performing integration by parts we find that the adjoint eigenfrequencies ω n $ \omega_n^\dagger $ and eigenfunctions Ψ n $ \Psi_n^\dagger $ are the solution of the following generalised boundary eigenvalue problem (Salwen & Grosch 1981)

L 1 Ψ n = ω n Δ ^ Ψ n , $$ \begin{aligned} \mathcal{L} _1^\dagger \Psi _n^\dagger = \omega _n^\dagger \widehat{\Delta } \Psi _n^\dagger , \end{aligned} $$(55)

where

L 1 k x U Δ ^ + k x ( β + 2 U d d ξ ) j ν turb Δ ^ 2 . $$ \begin{aligned} \mathcal{L} _1^\dagger \equiv k_x U \widehat{\Delta } + k_x \left(\beta + 2 U^\prime \dfrac{\mathrm{d} }{\mathrm{d} \xi } \right) - {j}{\nu _{\rm turb}}\widehat{\Delta }^2. \end{aligned} $$(56)

It is easy to show that ω n = ω n $ \omega_n^\dagger = \omega_n^\ast $. Furthermore, the eigenfunctions and adjoint eigenfunctions form a biorthogonal set, in the sense that

1 1 d ξ ( d Ψ n d y d Ψ m d y + k x 2 Ψ n Ψ m ) δ nm . $$ \begin{aligned} \displaystyle \int _{-1}^1 \mathrm{d} \xi ~ \left( \dfrac{\mathrm{d} \Psi _n^*}{\mathrm{d} y} \dfrac{\mathrm{d} \Psi _m^\dagger }{\mathrm{d} y} + k_x^2 \Psi _n^*\Psi _m^\dagger \right) \propto \delta _{nm}. \end{aligned} $$(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 uy spectrum associated with each mode separately.

The results are shown in the bottom panels of Figs. 35, for azimuthal orders m = 3, 5, and 10, respectively. The total spectrum is represented by the solid black line, while the coloured dashed lines represent the individual contributions of the most prominent modes in the discrete spectrum. We note that for the most part, and as is expected of a resonant mode whose driving source has a broadband frequency dependence, the spectrum contribution of each individual mode takes the form of a Lorentzian profile. For low values of m (more specifically for m ≤ 5), the spectrum is clearly dominated by two main peaks, and the frequencies of these peaks makes them clearly identifiable, as they fall very close to the real part of discrete eigenmodes of the homogeneous system, as shown in the top panels.

As m increases, the peaks become wider (because the imaginary part of the eigenfrequencies increase in modulus), and for m ≥ 6 the spectrum becomes dominated by one mode only. From the eigenspectrum shown in the top panels it can be seen that the dominant mode corresponds to the equatorial Rossby mode. Interestingly, we find that several modes have an amplitude that, if taken individually, should lead to a visible peak in the spectrum, but remain invisible in the total spectrum, where all modes are accounted for at once. This seems to suggest that the reason these modes cannot be extracted from the spectrum is not the inefficiency of the excitation process but rather the fact that they form a mutually destructive interference pattern with each other. This is possible because the modes all share the same driving source. This interference phenomenon can only occur between modes that share the same azimuthal order m, because the problem is axisymmetric, and therefore different m are completely decoupled. Furthermore, the modes must share similar eigenfrequencies (more specifically, the frequency difference must be of the order or smaller than the inverse of their lifetime).

3.1.2. Frequencies, amplitudes, and linewidths

For the resonant peaks that are sufficiently separated from each other in the synthetic spectrum, it is possible to directly infer, from the model, not only their frequency, but also their amplitude and linewidths. We define the angular frequency ω0 of each peak as the location of their local maximum, and their full linewidth at half maximum Γ as the angular frequency range where the power spectral density is above half the height of the peak. The amplitude is defined as

A = ( ω a ω b d ω P ( ω ) ) 1 / 2 , $$ \begin{aligned} A = \left( \displaystyle \int _{\omega _a}^{\omega _b} \mathrm{d} \omega ~ P(\omega ) \right)^{1/2}, \end{aligned} $$(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 = ω0 ∓ Γ/2. In the case of a Lorentzian profile, this range encloses exactly half the energy of the mode, so we only have to apply a factor of 2 $ \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, ω0, is systematically negative.

The frequencies of the modes are showcased, as a function of m, in Fig. 6. The coloured diamonds are identical on each panel: blue and green symbols represent each of the upper mode branches in the spectrum (see Sect. 3.1.1 for a description of these categories), while red symbols represent equatorial Rossby modes. On the background of each panel is superimposed an image of the spectrum in terms of different physical quantities (latitudinal velocity on top, azimuthal velocity in the middle, and vorticity at the bottom), in the m-ω plane, where each vertical slice is normalised so that the maximum is unity. One can distinguish the same transition, already mentioned above, between low azimuthal orders, for which several clearly identifiable resonant peaks can be resolved in the spectrum, and high azimuthal orders, where the distinction is no longer as clear. The theoretical Rossby mode dispersion relation is also plotted in each panel: in Cartesian coordinates, and in the presence of a background azimuthal jet U, it is given by

ω R = β U k x . $$ \begin{aligned} \omega _R = -\dfrac{\beta - U^{\prime \prime }}{k_x}. \end{aligned} $$(59)

thumbnail Fig. 6.

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

The high m frequencies match the theoretical dispersion relation quite well. The agreement, however, is not as close for lower values of m. We also note that while the low-m equatorial Rossby modes are quite distinguishable in the uy spectrum or the vorticity equatorial spectrum, they do not show in the equatorial ux spectrum, due to the fact that their ux 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 Returb. 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 Returb = 300 case, the same m2 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

Γ R = ν turb k x 2 . $$ \begin{aligned} \Gamma _R = {\nu _{\rm turb}} k_x^2. \end{aligned} $$(60)

thumbnail Fig. 7.

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

We note that the value of the turbulent Reynolds number that seems to give the best agreement between the theoretical dispersion relation and the observations is Returb = 300, which corresponds to a turbulent viscosity νturb = 570 km2 s−1. This value is in accordance with the surface value inferred by Gizon et al. (2020). However, it is significantly larger than the upper limit of 100 km2 s−1 inferred by Gizon et al. (2021), which was obtained under the assumption that the turbulent viscosity is constant over the entire convection zone.

Finally, we report the amplitude of the synthetic Rossby modes in the right panel of Fig. 7, in terms of latitudinal velocity uy, for several models corresponding to different values of Returb. 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 Returb = 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 Returb) causes both an earlier start of the plateau and a lower value thereof. For instance, a value Returb = 300 causes the equatorial Rossby modes to reach an amplitude of 1.5 m s−1 after m ∼ 10, while a value Returb = 1000 causes them to reach an amplitude of 2.5 m s−1 after m ∼ 15.

The fact that the amplitudes predicted by our simple 2D model compare well with the solar observations deserves some discussion. This good agreement suggests that neither the radial eigenfunctions nor the source function depend strongly on radius r. Regarding the eigenfunctions, linear eigenmode computations in 3D by Bekki et al. (2022a) indicate that the dependence on r is less pronounced for the smallest m than for the larger m. For example, the velocity eigenfunction of the equatorial Rossby mode scales like the slowly varying function rm for m = 3 and 4. Regarding the source of excitation, a turbulent vorticity spectrum that peaks at a scale k0 can only drive modes with azimuthal wavenumbers kx = m/R < 2k0 (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.

thumbnail Fig. 8.

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

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 uy 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 uy power spectrum over this latitude range (that is, we average between yobs = ±sin(15° ) instead of simply taking yobs = 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.

thumbnail Fig. 9.

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

We find that for m = 1 and m = 2, the amplitude of the modes is lower than the turbulent noise level. On this, our model prediction agrees with observational evidence; Liang et al. (2019) had indeed reported that no significant peak could be detected for such low values of the azimuthal order.

For m = 3 to 5, we recall that several peaks could be distinguished in our synthetic spectrum. However, as illustrated for m = 3 in the top panel of Fig. 9, it turns out that only one of these resonant peaks has enough amplitude to rise above the turbulent noise level. This is also in agreement with the observations.

For m ≥ 5, we recall that only one peak is predicted to dominate in the equatorial spectrum, corresponding to the equatorial Rossby mode. At the start of this m ≥ 5 regime, this mode is clearly visible above the noise level (as illustrated for m = 10 in the bottom panel of Fig. 9). However, we find that the visibility of the mode becomes more and more questionable as m increases. This is in agreement with the solar observations reported by Liang et al. (2019), who could not detect any significant peak in the equatorial spectrum for m ≥ 16.

Our model allows some light to be shed on this 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 m2. Because the peaks become wider and wider for the same total power, their spectral height decreases as 1/m2, 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 uy 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 yc(ω) given by the following implicit relation:

U ( y c ) = ω k x = ω R m , $$ \begin{aligned} U({ y}_c) = \dfrac{\omega }{k_x} = \dfrac{\omega R}{m}, \end{aligned} $$(61)

thumbnail Fig. 10.

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

where U is given by Eq. (24).

The shape of the power spectrum clearly transitions from a 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 ux 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 uy 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 uy spectrum. This is simply a symptom of the qualitative difference between the uy and ux 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.

thumbnail Fig. 11.

Same as Fig. 10, but for the power spectrum of the azimuthal velocity, ux.

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 uy or ux eigenfunctions. However, the low-m to high-m transition described above can still be observed in the vorticity power spectra.

thumbnail Fig. 12.

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

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−1, similar to those reported, for instance, by Liang et al. (2019). The amplitude of the linearly stable inertial modes observed in the equatorial region of the Sun is therefore consistent with a stochastic excitation by turbulent convection. However, we did not treat the case of the unstable, high-latitude inertial modes.

We also show that the power spectra in the 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−1. This confirms constraints previously obtained in the literature for the convection-induced turbulent viscosity. We also show that the linewidths of the equatorial Rossby modes vary as m2, 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

  1. Allison, M. 1990, Icarus, 83, 282 [NASA ADS] [CrossRef] [Google Scholar]
  2. Augustson, K. C., Mathis, S., & Astoul, A. 2020, ApJ, 903, 90 [NASA ADS] [CrossRef] [Google Scholar]
  3. Balmforth, N. J. 1992, MNRAS, 255, 639 [NASA ADS] [Google Scholar]
  4. Bekki, Y., Cameron, R. H., & Gizon, L. 2022a, A&A, 666, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Bekki, Y., Cameron, R. H., & Gizon, L. 2022b, A&A, 662, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Belkacem, K., Samadi, R., Goupil, M. J., et al. 2010, A&A, 522, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Boussinesq, J. 1877, Essai sur la théorie des eaux courantes, Mémoires présentés par divers savants à l’Académie des sciences de l’Institut national de France (Imprimerie Nationale) [Google Scholar]
  8. Chaplin, W. J., Houdek, G., Elsworth, Y., et al. 2005, MNRAS, 360, 859 [NASA ADS] [CrossRef] [Google Scholar]
  9. Christensen-Dalsgaard, J. 2002, Rev. Modern Phys., 74, 1073 [NASA ADS] [CrossRef] [Google Scholar]
  10. Deubner, F. L. 1975, A&A, 44, 371 [Google Scholar]
  11. Dikpati, M., Gilman, P. A., Guerrero, G. A., et al. 2022, ApJ, 931, 117 [NASA ADS] [CrossRef] [Google Scholar]
  12. Fournier, D., Gizon, L., & Hyest, L. 2022, A&A, 664, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Gizon, L., Fournier, D., & Albekioni, M. 2020, A&A, 642, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Gizon, L., Cameron, R. H., Bekki, Y., et al. 2021, A&A, 652, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Goldreich, P., & Keeley, D. A. 1977, ApJ, 212, 243 [NASA ADS] [CrossRef] [Google Scholar]
  16. Hanson, C. S., Hanasoge, S., & Sreenivasan, K. R. 2022, Nat. Astron., 6, 708 [NASA ADS] [CrossRef] [Google Scholar]
  17. Kraichnan, R. H. 1965, Phys. Fluids, 8, 575 [NASA ADS] [CrossRef] [Google Scholar]
  18. Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 581, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Leighton, R. B., Noyes, R. W., & Simon, G. W. 1962, ApJ, 135, 474 [NASA ADS] [CrossRef] [Google Scholar]
  20. Lesieur, M. 2008, Turbulence in Fluids (Springer) [CrossRef] [Google Scholar]
  21. Liang, Z.-C., Gizon, L., Birch, A. C., & Duvall, T. L. 2019, A&A, 626, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Lighthill, M. J. 1967, in Aerodynamic Phenomena in Stellar Atmospheres, ed. R. N. Thomas, 28, 429 [NASA ADS] [Google Scholar]
  23. Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat. Astron., 2, 568 [Google Scholar]
  24. Mack, L. M. 1976, J. Fluid Mech., 73, 497 [NASA ADS] [CrossRef] [Google Scholar]
  25. Mathis, S., Neiner, C., & Tran Minh, N. 2014, A&A, 565, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Miesch, M. S., Brun, A. S., DeRosa, M. L., & Toomre, J. 2008, ApJ, 673, 557 [NASA ADS] [CrossRef] [Google Scholar]
  27. Millionshchikov, M. 1941, Dokl. Akad. Nauk SSSR, 32, 611 [Google Scholar]
  28. Musielak, Z. E., Rosner, R., Stein, R. F., & Ulmschneider, P. 1994, ApJ, 423, 474 [Google Scholar]
  29. Orszag, S. A. 1971, J. Fluid Mech., 50, 689 [NASA ADS] [CrossRef] [Google Scholar]
  30. Papaloizou, J., & Pringle, J. E. 1978, MNRAS, 182, 423 [NASA ADS] [CrossRef] [Google Scholar]
  31. Rossby, C.-G. 1939, J. Marine Res., 2, 38 [CrossRef] [Google Scholar]
  32. Salwen, H., & Grosch, C. E. 1981, J. Fluid Mech., 104, 445 [NASA ADS] [CrossRef] [Google Scholar]
  33. Samadi, R., & Goupil, M. J. 2001, A&A, 370, 136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Samadi, R., Georgobiani, D., Trampedach, R., et al. 2007, A&A, 463, 297 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Sánchez-Lavega, A., Río-Gaztelurrutia, T., Hueso, R., et al. 2014, Geophys. Rev. Lett., 41, 1425 [CrossRef] [Google Scholar]
  36. Stein, R. F. 1967, Sol. Phys., 2, 385 [NASA ADS] [CrossRef] [Google Scholar]
  37. Stein, R. F., & Nordlund, Å. 1998, ApJ, 499, 914 [Google Scholar]
  38. Tennekes, H., & Lumley, J. 1978, in A First Course in Turbulence (MIT Press) [Google Scholar]
  39. Triana, S. A., Guerrero, G., Barik, A., & Rekier, J. 2022, ApJ, 934, L4 [NASA ADS] [CrossRef] [Google Scholar]
  40. Xiong, D.-R. 1989, A&A, 209, 126 [NASA ADS] [Google Scholar]

Appendix A: Vorticity wave equation in the equatorial β plane

A.1. Transport equation for the stream function

The total flow velocity v is decomposed into a background zonal flow U ≡ U(y) ex, which represents differential rotation, and a residual flow u ≡ ux(x, y)ex + uy(x, y)ey, which contains both the waves and the convective noise. The equations of motion become

u x t + U u x x + u y U + u x u x x + u y u x y = 1 ρ p x + f u y , $$ \begin{aligned}&\dfrac{\partial u_x}{\partial t} + U\dfrac{\partial u_x}{\partial x} + u_{ y} U^{\prime } + u_x \dfrac{\partial u_x}{\partial x} + u_{ y} \dfrac{\partial u_x}{\partial y} = -\dfrac{1}{\rho } \dfrac{\partial p}{\partial x} + f u_{ y}~, \end{aligned} $$(A.1)

u y t + U u y x + u x u y x + u y u y y = 1 ρ p y f u x , $$ \begin{aligned}&\dfrac{\partial u_{ y}}{\partial t} + U\dfrac{\partial u_{ y}}{\partial x} + u_x \dfrac{\partial u_{ y}}{\partial x} + u_{ y} \dfrac{\partial u_{ y}}{\partial y} = -\dfrac{1}{\rho } \dfrac{\partial p}{\partial y} - f u_x~ , \end{aligned} $$(A.2)

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

u = ( Ψ e z ) , $$ \begin{aligned} \mathbf u = \boldsymbol{\nabla } \wedge \left( \Psi \mathbf e _z \right) , \end{aligned} $$(A.3)

where ez is the radial unit vector. Explicitly,

u x = Ψ y , $$ \begin{aligned}&u_x = \dfrac{\partial \Psi }{\partial y}~, \end{aligned} $$(A.4)

u y = Ψ x . $$ \begin{aligned}&u_{ y} = -\dfrac{\partial \Psi }{\partial x}~. \end{aligned} $$(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:

Δ Ψ t + U Δ Ψ x U Ψ x + Ψ y Δ Ψ x Ψ x Δ Ψ y = 1 ρ 2 ( ρ y p x ρ x p y ) β Ψ x , $$ \begin{aligned}&\dfrac{\partial \Delta \Psi }{\partial t} + U\dfrac{\partial \Delta \Psi }{\partial x} - U^{\prime \prime } \dfrac{\partial \Psi }{\partial x} + \dfrac{\partial \Psi }{\partial y} \dfrac{\partial \Delta \Psi }{\partial x} - \dfrac{\partial \Psi }{\partial x} \dfrac{\partial \Delta \Psi }{\partial y} \nonumber \\&\quad = \dfrac{1}{\rho ^2}\left( \dfrac{\partial \rho }{\partial y} \dfrac{\partial p}{\partial x} - \dfrac{\partial \rho }{\partial x} \dfrac{\partial p}{\partial y} \right) - \beta \dfrac{\partial \Psi }{\partial x}~, \end{aligned} $$(A.6)

where Δ x 2 + y 2 $ \Delta \equiv \partial_x^2 + \partial_{\it 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,

p = p ( ρ ) . $$ \begin{aligned} p = p(\rho )~. \end{aligned} $$(A.7)

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

ρ y p x ρ x p y = ( p ρ ) · e z = 0 . $$ \begin{aligned} \dfrac{\partial \rho }{\partial y} \dfrac{\partial p}{\partial x} - \dfrac{\partial \rho }{\partial x} \dfrac{\partial p}{\partial y} = (\boldsymbol{\nabla } p \wedge \boldsymbol{\nabla }\rho ) \cdot \mathbf e _z = 0~. \end{aligned} $$(A.8)

We thus obtain a purely mechanical equation,

( t + U x ) Δ Ψ + ( β U ) Ψ x + Ψ y Δ Ψ x Ψ x Δ Ψ y = 0 . $$ \begin{aligned} \left( \dfrac{\partial }{\partial t} + U\dfrac{\partial }{\partial x}\right) \Delta \Psi + \left(\beta - U^{\prime \prime }\right) \dfrac{\partial \Psi }{\partial x} + \dfrac{\partial \Psi }{\partial y} \dfrac{\partial \Delta \Psi }{\partial x} - \dfrac{\partial \Psi }{\partial x} \dfrac{\partial \Delta \Psi }{\partial y} = 0~. \end{aligned} $$(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

δ ( Ψ y Δ Ψ x Ψ x Δ Ψ y ) Ψ y Δ Ψ x Ψ x Δ Ψ y Ψ y Δ Ψ x Ψ x Δ Ψ y h , $$ \begin{aligned}&\delta \left( \dfrac{\partial \Psi }{\partial y} \dfrac{\partial \Delta \Psi }{\partial x} - \dfrac{\partial \Psi }{\partial x} \dfrac{\partial \Delta \Psi }{\partial y} \right) \equiv \dfrac{\partial \Psi }{\partial y} \dfrac{\partial \Delta \Psi }{\partial x} - \dfrac{\partial \Psi }{\partial x} \dfrac{\partial \Delta \Psi }{\partial y} \nonumber \\&\quad - \left\langle \dfrac{\partial \Psi }{\partial y} \dfrac{\partial \Delta \Psi }{\partial x} - \dfrac{\partial \Psi }{\partial x} \dfrac{\partial \Delta \Psi }{\partial y} \right\rangle _h~, \end{aligned} $$(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

Ψ y Δ Ψ x Ψ x Δ Ψ y h = ( u · ) Δ Ψ h = · u Δ Ψ h , $$ \begin{aligned} \left\langle \dfrac{\partial \Psi }{\partial y} \dfrac{\partial \Delta \Psi }{\partial x} - \dfrac{\partial \Psi }{\partial x} \dfrac{\partial \Delta \Psi }{\partial y} \right\rangle _h = \left\langle (\mathbf u \cdot \boldsymbol{\nabla }) \Delta \Psi \right\rangle _h = \boldsymbol{\nabla } \cdot \langle \mathbf u ~ \Delta \Psi \rangle _h~, \end{aligned} $$(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

u Δ Ψ = ν turb Δ Ψ h . $$ \begin{aligned} \langle \mathbf u ~ \Delta \Psi \rangle = -{\nu _{\rm turb}} \boldsymbol{\nabla } \left\langle \Delta \Psi \right\rangle _h~. \end{aligned} $$(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

( t + U x ) Δ Ψ + ( β U ) Ψ x · ( ν turb Δ Ψ h ) + δ ( Ψ y Δ Ψ x Ψ x Δ Ψ y ) = 0 . $$ \begin{aligned}&\left( \dfrac{\partial }{\partial t} + U\dfrac{\partial }{\partial x}\right) \Delta \Psi + \left(\beta - U^{\prime \prime }\right) \dfrac{\partial \Psi }{\partial x} - \boldsymbol{\nabla } \cdot ({\nu _{\rm turb}} \boldsymbol{\nabla } \Delta \langle \Psi \rangle _h) \nonumber \\&\quad + \delta \left( \dfrac{\partial \Psi }{\partial y} \dfrac{\partial \Delta \Psi }{\partial x} - \dfrac{\partial \Psi }{\partial x} \dfrac{\partial \Delta \Psi }{\partial y} \right) = 0~ . \end{aligned} $$(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:

Ψ = Ψ osc + Ψ turb . $$ \begin{aligned} \Psi = {\Psi _{\rm osc}} + {\Psi _{\rm turb}} ~. \end{aligned} $$(A.14)

In the region of excitation, we assume that

Ψ osc Ψ turb . $$ \begin{aligned} {\Psi _{\rm osc}} \ll {\Psi _{\rm turb}}~. \end{aligned} $$(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−1 (Stein & Nordlund 1998), while the typical amplitudes of inertial modes are of the order of several m s−1 at most. However, there must be a smooth transition between the region of wave excitation and the neighbouring radiative zone, where Ψturb is negligible. Therefore, the ordering Ψosc ≪ Ψturb cannot be valid everywhere. In the following, we work in the region of wave excitation. We also note that, by definition of the horizontal average used to introduce the turbulent viscosity, we have

Ψ h = Ψ osc . $$ \begin{aligned} \langle \Psi \rangle _h = {\Psi _{\rm osc}}~. \end{aligned} $$(A.16)

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

( t + U x ) Δ Ψ osc + ( β U ) Ψ osc x ν turb Δ 2 Ψ osc + δ ( Ψ osc y Δ Ψ turb x + Ψ turb y Δ Ψ osc x Ψ osc x Δ Ψ turb y Ψ turb x Δ Ψ osc y ) = U Δ Ψ turb x ( β U ) Ψ turb x + δ ( Ψ turb x Δ Ψ turb y Ψ turb y Δ Ψ turb x ) , $$ \begin{aligned}&\left( \dfrac{\partial }{\partial t} + U\dfrac{\partial }{\partial x}\right) \Delta {\Psi _{\rm osc}} + (\beta - U^{\prime \prime }) \dfrac{\partial {\Psi _{\rm osc}}}{\partial x} - {\nu _{\rm turb}} \Delta ^2 {\Psi _{\rm osc}} \nonumber \\&\quad + \delta \left( \dfrac{\partial {\Psi _{\rm osc}}}{\partial y} \dfrac{\partial \Delta {\Psi _{\rm turb}}}{\partial x} + \dfrac{\partial {\Psi _{\rm turb}}}{\partial y} \dfrac{\partial \Delta {\Psi _{\rm osc}}}{\partial x} \right.\nonumber \\&\quad \left. - \dfrac{\partial {\Psi _{\rm osc}}}{\partial x} \dfrac{\partial \Delta {\Psi _{\rm turb}}}{\partial y} - \dfrac{\partial {\Psi _{\rm turb}}}{\partial x} \dfrac{\partial \Delta {\Psi _{\rm osc}}}{\partial y} \right) \nonumber \\&\quad = -U\dfrac{\partial \Delta {\Psi _{\rm turb}}}{\partial x} - (\beta - U^{\prime \prime })\dfrac{\partial {\Psi _{\rm turb}}}{\partial x} \nonumber \\&\quad + \delta \left( \dfrac{\partial {\Psi _{\rm turb}}}{\partial x}\dfrac{\partial \Delta {\Psi _{\rm turb}}}{\partial y} - \dfrac{\partial {\Psi _{\rm turb}}}{\partial y}\dfrac{\partial \Delta {\Psi _{\rm turb}}}{\partial x} \right)~, \end{aligned} $$(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

( t + U x ) Δ Ψ osc + ( β U ) Ψ osc x ν turb Δ 2 Ψ osc = δ ( Ψ turb x Δ Ψ turb y Ψ turb y Δ Ψ turb x ) . $$ \begin{aligned}&\left( \dfrac{\partial }{\partial t} + U\dfrac{\partial }{\partial x}\right) \Delta {\Psi _{\rm osc}} + (\beta - U^{\prime \prime }) \dfrac{\partial {\Psi _{\rm osc}}}{\partial x} - {\nu _{\rm turb}} \Delta ^2 {\Psi _{\rm osc}} \nonumber \\&\quad = \delta \left( \dfrac{\partial {\Psi _{\rm turb}}}{\partial x}\dfrac{\partial \Delta {\Psi _{\rm turb}}}{\partial y} - \dfrac{\partial {\Psi _{\rm turb}}}{\partial y}\dfrac{\partial \Delta {\Psi _{\rm turb}}}{\partial x} \right)~. \end{aligned} $$(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)

I ( y s ) d Y S ^ ( y s ) S ^ ( y s + Y ) , $$ \begin{aligned} \mathcal{I} ({ y}_s) \equiv \displaystyle \int \mathrm{d} Y ~ \left\langle \widehat{S}({ y}_s) ~ \widehat{S}^*({ y}_s + Y) \right\rangle ~, \end{aligned} $$(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

S ( t , x , y ) = δ ( ϵ ijz i Ψ turb j Δ Ψ turb ) , $$ \begin{aligned} S(t, x, y) = \delta \left( \epsilon _{ijz} \partial _i {\Psi _{\rm turb}} \partial _j \Delta {\Psi _{\rm turb}} \right)~, \end{aligned} $$(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

I ( y s ) = 1 T obs X obs d Y d t d t d x d x ϵ ijz ϵ klz i Ψ turb j Δ Ψ turb | t , x , y s k Ψ turb l Δ Ψ turb | t , x , y s + Y e j ( ω ( t t ) k x ( x x ) ) . $$ \begin{aligned}&\mathcal{I} ({ y}_s) = \dfrac{1}{T_{\rm obs} X_{\rm obs}} \displaystyle \int \mathrm{d} Y \mathrm{d} t \mathrm{d} t^{\prime } \mathrm{d} x \mathrm{d} x^{\prime } ~ \epsilon _{ijz} \epsilon _{klz} \nonumber \\&\quad \left\langle \left. \partial _i {\Psi _{\rm turb}} \partial _j \Delta {\Psi _{\rm turb}} \right|_{t, x, { y}_s} ~ \left. \partial _k {\Psi _{\rm turb}} \partial _l \Delta {\Psi _{\rm turb}} \right|_{t^{\prime }, x^{\prime }, { y}_s + Y} \right\rangle e^{{j} \left( \omega (t-t^{\prime }) - k_x (x-x^{\prime }) \right)} ~. \end{aligned} $$(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

a b c d = a b c d + a c b d + a d b c . $$ \begin{aligned} \left\langle abcd \right\rangle = \left\langle ab \right\rangle \left\langle cd \right\rangle + \left\langle ac \right\rangle \left\langle bd \right\rangle + \left\langle ad \right\rangle \left\langle bc \right\rangle ~. \end{aligned} $$(B.4)

Then Eq. B.3 becomes

I ( y s ) = I a ( y s ) + I b ( y s ) + I c ( y s ) , $$ \begin{aligned} \mathcal{I} ({ y}_s) = \mathcal{I} _a({ y}_s) + \mathcal{I} _b({ y}_s) + \mathcal{I} _c({ y}_s)~, \end{aligned} $$(B.5)

where

I a ( y s ) 1 T obs X obs d Y d t d t d x d x ϵ ijz ϵ klz i Ψ turb j Δ Ψ turb | t , x , y s k Ψ turb l Δ Ψ turb | t , x , y s + Y e j ( ω ( t t ) k x ( x x ) ) , $$ \begin{aligned}&\mathcal{I} _a({ y}_s) \equiv \dfrac{1}{T_{\rm obs} X_{\rm obs}} \displaystyle \int \mathrm{d} Y \mathrm{d} t \mathrm{d} t^{\prime } \mathrm{d} x \mathrm{d} x^{\prime } ~ \epsilon _{ijz} \epsilon _{klz} \nonumber \\&\quad \left\langle \left. \partial _i {\Psi _{\rm turb}} \partial _j \Delta {\Psi _{\rm turb}} \right|_{t, x, { y}_s} \right\rangle ~ \left\langle \left. \partial _k {\Psi _{\rm turb}} \partial _l \Delta {\Psi _{\rm turb}} \right|_{t^{\prime }, x^{\prime }, { y}_s + Y} \right\rangle \nonumber \\&\quad e^{{j} \left( \omega (t-t^{\prime }) - k_x (x-x^{\prime }) \right)}~, \end{aligned} $$(B.6)

I b ( y s ) 1 T obs X obs d Y d t d t d x d x ϵ ijz ϵ klz i Ψ turb | t , x , y s k Ψ turb | t , x , y s + Y j Δ Ψ turb | t , x , y s l Δ Ψ turb | t , x , y s + Y e j ( ω ( t t ) k x ( x x ) ) , $$ \begin{aligned}&\mathcal{I} _b({ y}_s) \equiv \dfrac{1}{T_{\rm obs} X_{\rm obs}} \displaystyle \int \mathrm{d} Y \mathrm{d} t \mathrm{d} t^{\prime } \mathrm{d} x \mathrm{d} x^{\prime } ~ \epsilon _{ijz} \epsilon _{klz} \nonumber \\&\quad \left\langle \left. \partial _i {\Psi _{\rm turb}} \right|_{t, x, { y}_s} \left. \partial _k {\Psi _{\rm turb}} \right|_{t^{\prime }, x^{\prime }, { y}_s + Y} \right\rangle ~ \left\langle \left. \partial _j \Delta {\Psi _{\rm turb}} \right|_{t, x, { y}_s} \left. \partial _l \Delta {\Psi _{\rm turb}} \right|_{t^{\prime }, x^{\prime }, { y}_s + Y} \right\rangle \nonumber \\&\quad e^{{j} \left( \omega (t-t^{\prime }) - k_x (x-x^{\prime }) \right)}~, \end{aligned} $$(B.7)

and

I c ( y s ) 1 T obs X obs d Y d t d t d x d x ϵ ijz ϵ klz i Ψ turb | t , x , y s l Δ Ψ turb | t , x , y s + Y j Δ Ψ turb | t , x , y s k Ψ turb | t , x , y s + Y e j ( ω ( t t ) k x ( x x ) ) . $$ \begin{aligned}&\mathcal{I} _c({ y}_s) \equiv \dfrac{1}{T_{\rm obs} X_{\rm obs}} \displaystyle \int \mathrm{d} Y \mathrm{d} t \mathrm{d} t^{\prime } \mathrm{d} x \mathrm{d} x^{\prime } ~ \epsilon _{ijz} \epsilon _{klz} \nonumber \\&\quad \left\langle \left. \partial _i {\Psi _{\rm turb}} \right|_{t, x, { y}_s} \left. \partial _l \Delta {\Psi _{\rm turb}} \right|_{t^{\prime }, x^{\prime }, { y}_s + Y} \right\rangle ~ \left\langle \left. \partial _j \Delta {\Psi _{\rm turb}} \right|_{t, x, { y}_s} \left. \partial _k {\Psi _{\rm turb}} \right|_{t^{\prime }, x^{\prime }, { y}_s + Y} \right\rangle \nonumber \\&\quad e^{{j} \left( \omega (t-t^{\prime }) - k_x (x-x^{\prime }) \right)}~. \end{aligned} $$(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,

T t , τ t t , X x , δ x x x , $$ \begin{aligned} T \equiv t~, \qquad \tau \equiv t^{\prime } - t ~, \qquad X \equiv x ~, \qquad \delta x \equiv x^{\prime } - x~, \end{aligned} $$(B.9)

so that

I b ( y s ) = 1 T obs X obs d Y d T d τ d X d δ x ϵ ijz ϵ klz i Ψ turb | T , X , y s k Ψ turb | T + τ , X + δ x , y s + Y j Δ Ψ turb | T , X , y s l Δ Ψ turb | T + τ , X + δ x , y s + Y e j ( k x δ x ω τ ) , $$ \begin{aligned}&\mathcal{I} _b({ y}_s) = \dfrac{1}{T_{\rm obs} X_{\rm obs}} \displaystyle \int \mathrm{d} Y \mathrm{d} T \mathrm{d} \tau \mathrm{d} X \mathrm{d} \delta x ~ \epsilon _{ijz} \epsilon _{klz} \nonumber \\&\quad \left\langle \left. \partial _i {\Psi _{\rm turb}} \right|_{T, X, { y}_s} \left. \partial _k {\Psi _{\rm turb}} \right|_{T+\tau , X+\delta x, { y}_s+Y} \right\rangle \nonumber \\&\quad \left\langle \left. \partial _j \Delta {\Psi _{\rm turb}} \right|_{T, X, { y}_s} \left. \partial _l \Delta {\Psi _{\rm turb}} \right|_{T+\tau , X+\delta x, { y}_s+Y} \right\rangle ~ e^{{j} ( k_x \delta x - \omega \tau )}~, \end{aligned} $$(B.10)

I c ( y s ) = 1 T obs X obs d Y d T d τ d X d δ x ϵ ijz ϵ klz i Ψ turb | T , X , y s l Δ Ψ turb | T + τ , X + δ x , y s + Y j Δ Ψ turb | T , X , y s k Ψ turb | T + τ , X + δ x , y s + Y e j ( k x δ x ω τ ) . $$ \begin{aligned}&\mathcal{I} _c({ y}_s) = \dfrac{1}{T_{\rm obs} X_{\rm obs}} \displaystyle \int \mathrm{d} Y \mathrm{d} T \mathrm{d} \tau \mathrm{d} X \mathrm{d} \delta x ~ \epsilon _{ijz} \epsilon _{klz} \nonumber \\&\quad \left\langle \left. \partial _i {\Psi _{\rm turb}} \right|_{T, X, { y}_s} \left. \partial _l \Delta {\Psi _{\rm turb}} \right|_{T+\tau , X+\delta x, { y}_s+Y} \right\rangle \nonumber \\&\quad \left\langle \left. \partial _j \Delta {\Psi _{\rm turb}} \right|_{T, X, { y}_s} \left. \partial _k {\Psi _{\rm turb}} \right|_{T+\tau , X+\delta x, { y}_s+Y} \right\rangle ~ e^{{j} ( k_x \delta x - \omega \tau )} ~. \end{aligned} $$(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 ys), in which case these two integrals simplify to

I b = d τ d δ x d Y ϵ ijz ϵ klz i Ψ turb | 0 , 0 , 0 k Ψ turb | τ , δ x , Y j Δ Ψ turb | 0 , 0 , 0 l Δ Ψ turb | τ , δ x , Y e j ( k x δ x ω τ ) $$ \begin{aligned}&\mathcal{I} _b = \displaystyle \int \mathrm{d} \tau \mathrm{d} \delta x \mathrm{d} Y ~ \epsilon _{ijz} \epsilon _{klz} \left\langle \left. \partial _i {\Psi _{\rm turb}} \right|_{0, 0, 0} \left. \partial _k {\Psi _{\rm turb}} \right|_{\tau , \delta x, Y} \right\rangle \nonumber \\&\quad \left\langle \left. \partial _j \Delta {\Psi _{\rm turb}} \right|_{0, 0, 0} \left. \partial _l \Delta {\Psi _{\rm turb}} \right|_{\tau , \delta x, Y} \right\rangle ~ e^{{j} ( k_x \delta x - \omega \tau )} \end{aligned} $$(B.12)

and

I c = d τ d δ x d Y ϵ ijz ϵ klz i Ψ turb | 0 , 0 , 0 l Δ Ψ turb | τ , δ x , Y j Δ Ψ turb | 0 , 0 , 0 k Ψ turb | τ , δ x , Y e j ( k x δ x ω τ ) . $$ \begin{aligned}&\mathcal{I} _c = \displaystyle \int \mathrm{d} \tau \mathrm{d} \delta x \mathrm{d} Y ~ \epsilon _{ijz} \epsilon _{klz} \left\langle \left. \partial _i {\Psi _{\rm turb}} \right|_{0, 0, 0} \left. \partial _l \Delta {\Psi _{\rm turb}} \right|_{\tau , \delta x, Y} \right\rangle \nonumber \\&\quad \left\langle \left. \partial _j \Delta {\Psi _{\rm turb}} \right|_{0, 0, 0} \left. \partial _k {\Psi _{\rm turb}} \right|_{\tau , \delta x, Y} \right\rangle ~ e^{{j} ( k_x \delta x - \omega \tau )}~. \end{aligned} $$(B.13)

If we introduce the following functions,

f b , i k ( τ , δ x ) i Ψ turb | 0 , 0 k Ψ turb | τ , δ x , g b , j l ( τ , δ x ) j Δ Ψ turb | 0 , 0 l Δ Ψ turb | τ , δ x , f c , i l ( τ , δ x ) i Ψ turb | 0 , 0 l Δ Ψ turb | τ , δ x , g c , j k ( τ , δ x ) j Δ Ψ turb | 0 , 0 k Ψ turb | τ , δ x , $$ \begin{aligned}&f_{b, ik}(\tau , \delta \mathbf x ) \equiv \left\langle \left. \partial _i {\Psi _{\rm turb}} \right|_{0, \mathbf 0 } \left. \partial _k {\Psi _{\rm turb}} \right|_{\tau , \delta \mathbf x } \right\rangle ~, \nonumber \\&g_{b, jl}(\tau , \delta \mathbf x ) \equiv \left\langle \left. \partial _j \Delta {\Psi _{\rm turb}} \right|_{0, \mathbf 0 } \left. \partial _l \Delta {\Psi _{\rm turb}} \right|_{\tau , \delta \mathbf x } \right\rangle ~, \nonumber \\&f_{c, il}(\tau , \delta \mathbf x ) \equiv \left\langle \left. \partial _i {\Psi _{\rm turb}} \right|_{0, \mathbf 0 } \left. \partial _l \Delta {\Psi _{\rm turb}} \right|_{\tau , \delta \mathbf x } \right\rangle ~, \nonumber \\&g_{c, jk}(\tau , \delta \mathbf x ) \equiv \left\langle \left. \partial _j \Delta {\Psi _{\rm turb}} \right|_{0, \mathbf 0 } \left. \partial _k {\Psi _{\rm turb}} \right|_{\tau , \delta \mathbf x } \right\rangle ~, \end{aligned} $$(B.14)

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

I b = T obs X obs Y obs ϵ ijz ϵ klz f b , i k g b , j l ( ω , k x , 0 ) , $$ \begin{aligned}&\mathcal{I} _b = \sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}} \epsilon _{ijz} \epsilon _{klz} \widetilde{f_{b, ik} g_{b, jl}}(-\omega , -k_x, 0)~, \end{aligned} $$(B.15)

I c = T obs X obs Y obs ϵ ijz ϵ klz f c , i l g c , j k ( ω , k x , 0 ) , $$ \begin{aligned}&\mathcal{I} _c = \sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}} \epsilon _{ijz} \epsilon _{klz} \widetilde{f_{c, il} g_{c, jk}}(-\omega , -k_x, 0)~, \end{aligned} $$(B.16)

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

f ( ω , k x , k y ) 1 T obs X obs Y obs d t d x d y f ( t , x , y ) e j ( ω t k x x k y y ) . $$ \begin{aligned} \widetilde{f}(\omega , k_x, k_{ y}) \equiv \dfrac{1}{\sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}} \displaystyle \int \mathrm{d} t \mathrm{d} x \mathrm{d} y ~ f(t, x, y) e^{{j}(\omega t - k_x x - k_{ y} y)}~. \end{aligned} $$(B.17)

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

I b = T obs X obs Y obs ( 2 π ) 3 d ω d 2 k ϵ ijz ϵ klz f b , i k ( ω , k ) g b , j l ( ω + ω , k + k ) , $$ \begin{aligned}&\mathcal{I} _b = \dfrac{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime } \mathrm{d} ^2 \mathbf{k }^{\prime } \epsilon _{ijz} \epsilon _{klz} \widetilde{f_{b, ik}}(\omega ^{\prime }, \mathbf{k }^{\prime }) \nonumber \\&\quad \widetilde{g_{b, jl}}^*(\omega ^{\prime } + \omega , \mathbf{k }^{\prime } + \mathbf k )~, \end{aligned} $$(B.18)

I c = T obs X obs Y obs ( 2 π ) 3 d ω d 2 k ϵ ijz ϵ klz f c , i l ( ω , k ) g c , j k ( ω + ω , k + k ) , $$ \begin{aligned}&\mathcal{I} _c = \dfrac{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime } \mathrm{d} ^2 \mathbf{k }^{\prime } \epsilon _{ijz} \epsilon _{klz} \widetilde{f_{c, il}}(\omega ^{\prime }, \mathbf{k }^{\prime }) \nonumber \\&\quad \widetilde{g_{c, jk}}^*(\omega ^{\prime } + \omega , \mathbf{k }^{\prime } + \mathbf k )~, \end{aligned} $$(B.19)

where k ≡ kxex.

Next, we need to express the quantities f b , i k $ \widetilde{f_{b, ik}} $, g b , j l $ \widetilde{g_{b, jl}} $, f c , i l $ \widetilde{f_{c, il}} $ and g c , j k $ \widetilde{g_{c, jk}} $ appearing in these integrals. Taking the Fourier transforms of each of Eq. B.14, we find

f b , i k ( ω , k ) = T obs X obs Y obs ( 2 π ) 3 d ω d 2 k ( k i k k ) Ψ turb ( ω , k ) Ψ turb ( ω , k ) , $$ \begin{aligned}&\widetilde{f_{b, ik}}(\omega ^{\prime }, \mathbf{k }^{\prime }) = \dfrac{\sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime \prime } \mathrm{d} ^2\mathbf k {\prime \prime }~ (-k_i^{\prime \prime } k_k^{\prime }) \nonumber \\&\quad \left\langle {\widetilde{\Psi }_{\rm turb}}(\omega ^{\prime \prime }, \mathbf k {\prime \prime }) ~ {\widetilde{\Psi }_{\rm turb}}(\omega ^{\prime }, \mathbf{k }^{\prime }) \right\rangle ~, \end{aligned} $$(B.20)

g b , j l ( ω , k ) = T obs X obs Y obs ( 2 π ) 3 d ω d 2 k ( k j k 2 k l k 2 ) Ψ turb ( ω , k ) Ψ turb ( ω , k ) , $$ \begin{aligned}&\widetilde{g_{b, jl}}(\omega ^{\prime }, \mathbf{k }^{\prime }) = \dfrac{\sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime \prime } \mathrm{d} ^2\mathbf k {\prime \prime }~ (-k_j^{\prime \prime } k^{\prime \prime 2} k_l^{\prime } k^{\prime 2}) \nonumber \\&\quad \left\langle {\widetilde{\Psi }_{\rm turb}}(\omega ^{\prime \prime }, \mathbf k {\prime \prime }) ~ {\widetilde{\Psi }_{\rm turb}}(\omega ^{\prime }, \mathbf{k }^{\prime }) \right\rangle ~, \end{aligned} $$(B.21)

f c , i l ( ω , k ) = T obs X obs Y obs ( 2 π ) 3 d ω d 2 k ( + k i k l k 2 ) Ψ turb ( ω , k ) Ψ turb ( ω , k ) , $$ \begin{aligned}&\widetilde{f_{c, il}}(\omega ^{\prime }, \mathbf{k }^{\prime }) = \dfrac{\sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime \prime } \mathrm{d} ^2\mathbf k {\prime \prime }~ (+k_i^{\prime \prime } k_l^{\prime } k^{\prime 2}) \nonumber \\&\quad \left\langle {\widetilde{\Psi }_{\rm turb}}(\omega ^{\prime \prime }, \mathbf k {\prime \prime }) ~ {\widetilde{\Psi }_{\rm turb}}(\omega ^{\prime }, \mathbf{k }^{\prime }) \right\rangle ~, \end{aligned} $$(B.22)

and

g c , j k ( ω , k ) = T obs X obs Y obs ( 2 π ) 3 d ω d 2 k ( + k j k 2 k k ) Ψ turb ( ω , k ) Ψ turb ( ω , k ) . $$ \begin{aligned}&\widetilde{g_{c, jk}}(\omega ^{\prime }, \mathbf{k }^{\prime }) = \dfrac{\sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime \prime } \mathrm{d} ^2\mathbf k {\prime \prime }~ (+k_j^{\prime \prime } k^{\prime \prime 2} k_k^{\prime }) \nonumber \\&\quad \left\langle {\widetilde{\Psi }_{\rm turb}}(\omega ^{\prime \prime }, \mathbf k {\prime \prime }) ~ {\widetilde{\Psi }_{\rm turb}}(\omega ^{\prime }, \mathbf{k }^{\prime }) \right\rangle ~. \end{aligned} $$(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

f b , i k ( ω , k ) = 1 T obs X obs Y obs k i k k E Ψ ( ω , k ) , $$ \begin{aligned}&\widetilde{f_{b, ik}}(\omega ^{\prime }, \mathbf{k }^{\prime }) = \dfrac{1}{\sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}} k_i^{\prime } k_k^{\prime } {\mathcal{E} _\Psi }(\omega ^{\prime }, \mathbf k \prime )~, \end{aligned} $$(B.24)

g b , j l ( ω , k ) = 1 T obs X obs Y obs k j k l k 4 E Ψ ( ω , k ) , $$ \begin{aligned}&\widetilde{g_{b, jl}}(\omega ^{\prime }, \mathbf{k }^{\prime }) = \dfrac{1}{\sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}} k_j^{\prime } k_l^{\prime } k^{\prime 4} {\mathcal{E} _\Psi }(\omega ^{\prime }, \mathbf k \prime )~, \end{aligned} $$(B.25)

f c , i l ( ω , k ) = 1 T obs X obs Y obs k i k l k 2 E Ψ ( ω , k ) , $$ \begin{aligned}&\widetilde{f_{c, il}}(\omega ^{\prime }, \mathbf{k }^{\prime }) = -\dfrac{1}{\sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}} k_i^{\prime } k_l^{\prime } k^{\prime 2} {\mathcal{E} _\Psi }(\omega ^{\prime }, \mathbf k \prime )~, \end{aligned} $$(B.26)

g c , j k ( ω , k ) = 1 T obs X obs Y obs k j k k k 2 E Ψ ( ω , k ) . $$ \begin{aligned}&\widetilde{g_{c, jk}}(\omega ^{\prime }, \mathbf{k }^{\prime }) = -\dfrac{1}{\sqrt{T_{\rm obs} X_{\rm obs} Y_{\rm obs}}} k_j^{\prime } k_k^{\prime } k^{\prime 2} {\mathcal{E} _\Psi }(\omega ^{\prime }, \mathbf k \prime )~. \end{aligned} $$(B.27)

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

E Ψ ( ω , k ) d τ d 2 x Ψ turb ( T , X ) Ψ turb ( T + τ , X + x ) e j ( ω τ k · x ) . $$ \begin{aligned} {\mathcal{E} _\Psi }(\omega ^{\prime }, \mathbf k \prime ) \equiv \displaystyle \int \mathrm{d} \tau \mathrm{d} ^2\mathbf x \left\langle {\Psi _{\rm turb}}(T, \mathbf X ) {\Psi _{\rm turb}}(T+\tau , \mathbf X +\mathbf x ) \right\rangle ~ e^{{j} (\omega ^{\prime }\tau - \mathbf k \prime \cdot \mathbf x )}~. \end{aligned} $$(B.28)

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

I b = 1 ( 2 π ) 3 d ω d 2 k ϵ ijz ϵ klz k i ( k j + k j ) k k ( k l + k l ) ( k + k ) 4 E Ψ ( ω , k ) E Ψ ( ω + ω , k + k ) , $$ \begin{aligned}&\mathcal{I} _b = \dfrac{1}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime } \mathrm{d} ^2\mathbf k \prime \epsilon _{ijz} \epsilon _{klz} k_i^{\prime } (k_j^{\prime }+k_j) k_k^{\prime } (k_l^{\prime }+k_l) \nonumber \\&\quad (\mathbf k \prime + \mathbf k )^4 ~ {\mathcal{E} _\Psi }(\omega ^{\prime }, \mathbf k \prime ) {\mathcal{E} _\Psi }^*(\omega ^{\prime } + \omega , \mathbf k \prime + \mathbf k )~, \end{aligned} $$(B.29)

and

I c = 1 ( 2 π ) 3 d ω d 2 k ϵ ijz ϵ klz k i ( k j + k j ) ( k k + k k ) k l k 2 ( k + k ) 2 E Ψ ( ω , k ) E Ψ ( ω + ω , k + k ) . $$ \begin{aligned}&\mathcal{I} _c = \dfrac{1}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime } \mathrm{d} ^2\mathbf k \prime \epsilon _{ijz} \epsilon _{klz} k_i^{\prime } (k_j^{\prime } + k_j)(k_k^{\prime } + k_k) k_l^{\prime } k^{\prime 2} \nonumber \\&\quad (\mathbf k \prime + \mathbf k )^2 {\mathcal{E} _\Psi }(\omega ^{\prime }, \mathbf k \prime ) {\mathcal{E} _\Psi }^*(\omega ^{\prime } + \omega , \mathbf k \prime + \mathbf k )~. \end{aligned} $$(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

I = 1 ( 2 π ) 3 d ω d 2 k ϵ ijz ϵ klz k i ( k j + k j ) k k ( k l + k l ) ( k + k ) 2 ( ( k + k ) 2 k 2 ) E Ψ ( ω , k ) E Ψ ( ω + ω , k + k ) . $$ \begin{aligned}&\mathcal{I} = \dfrac{1}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime } \mathrm{d} ^2\mathbf k \prime \epsilon _{ijz} \epsilon _{klz} k_i^{\prime } (k_j^{\prime }+k_j) k_k^{\prime } (k_l^{\prime }+k_l) \nonumber \\&\quad (\mathbf k \prime + \mathbf k )^2 \left( (\mathbf k \prime + \mathbf k )^2 - \mathbf k ^{\prime 2} \right) ~ {\mathcal{E} _\Psi }(\omega ^{\prime }, \mathbf k \prime ) {\mathcal{E} _\Psi }^*(\omega ^{\prime } + \omega , \mathbf k \prime + \mathbf k )~. \end{aligned} $$(B.31)

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

ϵ ijz ϵ klz = δ ik δ jl δ zz + δ il δ jz δ kz + δ iz δ jk δ lz δ ik δ jz δ lz δ iz δ jl δ kz δ il δ jk δ zz . $$ \begin{aligned} \epsilon _{ijz} \epsilon _{klz} = \delta _{ik}\delta _{jl}\delta _{zz} + \delta _{il}\delta _{jz}\delta _{kz} + \delta _{iz}\delta _{jk}\delta _{lz} - \delta _{ik}\delta _{jz}\delta _{lz} - \delta _{iz}\delta _{jl}\delta _{kz} - \delta _{il}\delta _{jk}\delta _{zz}~. \end{aligned} $$(B.32)

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

ϵ ijz ϵ klz k i ( k j + k j ) k k ( k l + k l ) = ( δ ik δ jl δ il δ jk ) k i ( k j + k j ) k k ( k l + k l ) = k 2 ( k + k ) 2 ( k · ( k + k ) ) 2 = k 2 k 2 ( k · k ) 2 . $$ \begin{aligned} \epsilon _{ijz} \epsilon _{klz} k_i^{\prime } (k_j^{\prime } + k_j) k_k^{\prime } (k_l^{\prime } + k_l)&= (\delta _{ik}\delta _{jl} - \delta _{il}\delta _{jk}) k_i^{\prime } (k_j^{\prime } + k_j) k_k^{\prime } (k_l^{\prime } + k_l) \nonumber \\&= k^{\prime 2} (\mathbf k \prime + \mathbf k )^2 - (\mathbf k \prime \cdot (\mathbf k \prime + \mathbf k ))^2 \nonumber \\&= k^{\prime 2} k^2 - (\mathbf k \prime \cdot \mathbf k )^2~. \end{aligned} $$(B.33)

Thus we obtain

I = 1 ( 2 π ) 3 d ω d 2 k [ k 2 k 2 ( k · k ) 2 ] ( k + k ) 2 ( ( k + k ) 2 k 2 ) E Ψ ( ω , k ) E Ψ ( ω + ω , k + k ) . $$ \begin{aligned}&\mathcal{I} = \dfrac{1}{(2\pi )^3} \displaystyle \int \mathrm{d} \omega ^{\prime } \mathrm{d} ^2\mathbf k \prime ~ \left[ k^{\prime 2} k^2 - (\mathbf k \prime \cdot \mathbf k )^2 \right] (\mathbf k \prime + \mathbf k )^2 \nonumber \\&\quad \left( (\mathbf k \prime + \mathbf k )^2 - k^{\prime 2} \right) ~ {\mathcal{E} _\Psi }(\omega ^{\prime }, \mathbf k \prime ) {\mathcal{E} _\Psi }^*(\omega ^{\prime } + \omega , \mathbf k \prime + \mathbf k )~. \end{aligned} $$(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, ℰΨ,

I = 1 4 π 3 d ω d 2 k k x 3 k x k y 2 | k + k 2 | 2 E Ψ ( ω ω 2 , k k 2 ) E Ψ ( ω + ω 2 , k + k 2 ) . $$ \begin{aligned}&\mathcal{I} = \dfrac{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} \left| \mathbf k \prime + \dfrac{\mathbf{k}}{2} \right|^2 \nonumber \\&\quad {\mathcal{E} _\Psi }\left(\omega ^{\prime } - \dfrac{\omega }{2}, \mathbf k \prime - \dfrac{\mathbf{k}}{2} \right) {\mathcal{E} _\Psi }^*\left(\omega ^{\prime } + \dfrac{\omega }{2}, \mathbf k \prime + \dfrac{\mathbf{k}}{2} \right)~. \end{aligned} $$(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

T n ( ξ ) = cos ( n arccos ξ ) . $$ \begin{aligned} T_n(\xi ) = \cos (n \arccos \xi )~. \end{aligned} $$(C.1)

Using the definition of the linear operator, we can write

M ij = ( ω k x 2 + k x ( β U ) + j ν turb k x 4 ) M ij ( 1 ) $$ \begin{aligned} \mathcal{M} _{ij} =&- \left(\omega k_x^2 + k_x(\beta - U^{\prime \prime }) + {j} {\nu _{\rm turb}} k_x^4 \right) \mathcal{M} _{ij}^{(1)} \end{aligned} $$(C.2)

+ ( ω + 2 j ν turb k x 2 ) M ij ( 2 ) $$ \begin{aligned}&+ \left(\omega + 2 {j} {\nu _{\rm turb}} k_x^2\right) \mathcal{M} _{ij}^{(2)} \end{aligned} $$(C.3)

j ν turb M ij ( 3 ) + k x 3 M ij ( 4 ) k x M ij ( 5 ) , $$ \begin{aligned}&- {j} {\nu _{\rm turb}} \mathcal{M} _{ij}^{(3)} + k_x^3 \mathcal{M} _{ij}^{(4)} - k_x \mathcal{M} _{ij}^{(5)}~, \end{aligned} $$(C.4)

where

M ij ( 1 ) 2 π c i T i | T j = δ ij , $$ \begin{aligned}&\mathcal{M} _{ij}^{(1)} \equiv \dfrac{2}{\pi c_i} {T_i|T_j} = \delta _{ij}~, \end{aligned} $$(C.5)

M ij ( 2 ) 2 π c i T i | T j , $$ \begin{aligned}&\mathcal{M} _{ij}^{(2)} \equiv \dfrac{2}{\pi c_i} {T_i|T_j^{\prime \prime }}~, \end{aligned} $$(C.6)

M ij ( 3 ) 2 π c i T i | T j , $$ \begin{aligned}&\mathcal{M} _{ij}^{(3)} \equiv \dfrac{2}{\pi c_i} {T_i|T_j^{\prime \prime \prime \prime }}~, \end{aligned} $$(C.7)

M ij ( 4 ) 2 π c i T i | U T j , $$ \begin{aligned}&\mathcal{M} _{ij}^{(4)} \equiv \dfrac{2}{\pi c_i} {T_i|UT_j}~, \end{aligned} $$(C.8)

M ij ( 5 ) 2 π c i T i | U T j , $$ \begin{aligned}&\mathcal{M} _{ij}^{(5)} \equiv \dfrac{2}{\pi c_i} {T_i|U T_j^{\prime \prime }}~, \end{aligned} $$(C.9)

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

C.1. Projection of the derivatives

Given a function

f ( ξ ) = λ i T i ( ξ ) , $$ \begin{aligned} f (\xi ) = \sum \lambda _i T_i (\xi )~, \end{aligned} $$(C.10)

we can write, for any positive integer, p,

d p f d y p = i = 0 + λ i ( p ) T i , λ i ( 0 ) = λ i . $$ \begin{aligned} \dfrac{\mathrm{d} ^p f}{\mathrm{d} { y}^p} = \sum _{i=0}^{+\infty } \lambda _i^{(p)} T_i~, \qquad \lambda _i^{(0)} = \lambda _i~. \end{aligned} $$(C.11)

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

For any p ≥ 1 we can write

i = 0 + λ i ( p ) T i = d d y i = 0 + λ i ( p 1 ) T i . $$ \begin{aligned} \sum _{i=0}^{+\infty } \lambda _i^{(p)} T_i = \dfrac{\mathrm{d} }{\mathrm{d} y} \sum _{i=0}^{+\infty } \lambda _i^{(p-1)} T_i~. \end{aligned} $$(C.12)

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

T i = c i 2 ( i + 1 ) T i + 1 c i 2 ( i 1 ) T i 1 , $$ \begin{aligned} T_i = \dfrac{c_i}{2(i+1)}T_{i+1}\prime - \dfrac{c_i^{\prime }}{2(i-1)}T_{i-1}\prime ~, \end{aligned} $$(C.13)

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

i = 0 + λ i ( p ) T i = d d y i = 0 + λ i ( p ) ( c i 2 ( i + 1 ) T i + 1 c i 2 ( i 1 ) T i 1 ) = d d y i = 1 + ( c i 1 2 i λ i 1 ( p ) 1 2 i λ i + 1 ( p ) ) T i . $$ \begin{aligned} \sum _{i=0}^{+\infty } \lambda _i^{(p)} T_i&= \dfrac{\mathrm{d} }{\mathrm{d} y} \sum _{i=0}^{+\infty } \lambda _i^{(p)} \left( \dfrac{c_i}{2(i+1)} T_{i+1} - \dfrac{c_i^{\prime }}{2(i-1)}T_{i-1} \right) \nonumber \\&= \dfrac{\mathrm{d} }{\mathrm{d} y} \sum _{i=1}^{+\infty } \left( \dfrac{c_{i-1}}{2i} \lambda _{i-1}^{(p)} - \dfrac{1}{2i} \lambda _{i+1}^{(p)} \right) T_i~. \end{aligned} $$(C.14)

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

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

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

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

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

For example, the coefficients of the first derivative straightforwardly read

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

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

c i λ i ( 2 ) = 2 j = i + 1 i + j 1 [ 2 ] + j λ j ( 1 ) = 4 j = i + 1 i + j 1 [ 2 ] + j k = j + 1 k + j 1 [ 2 ] + k λ k = 4 k = i + 2 k + i 0 [ 2 ] + k λ k j = i + 1 j + k 1 [ 2 ] k 1 j = k = i + 2 k + i 0 [ 2 ] + k ( k 2 i 2 ) λ k . $$ \begin{aligned} c_i \lambda _i^{(2)}&= 2 \sum _{\begin{matrix} j=i+1 \\ i+j \equiv 1[2] \end{matrix}}^{+\infty } j \lambda _j^{(1)} = 4 \sum _{\begin{matrix} j=i+1 \\ i+j \equiv 1[2] \end{matrix}}^{+\infty } j \sum _{\begin{matrix} k=j+1 \\ k+j \equiv 1 [2] \end{matrix}}^{+\infty } k\lambda _k \nonumber \\&= 4\sum _{\begin{matrix} k=i+2 \\ k+i \equiv 0 [2] \end{matrix}}^{+\infty } k \lambda _k \sum _{\begin{matrix} j=i+1 \\ j+k \equiv 1 [2] \end{matrix}}^{k-1} j = \sum _{\begin{matrix} k=i+2 \\ k+i \equiv 0 [2] \end{matrix}}^{+\infty } k \left( k^2 - i^2 \right) \lambda _k~. \end{aligned} $$(C.18)

Similarly, we find

c i λ i ( 3 ) = 1 4 k = i + 3 k + i 1 [ 2 ] + k [ k 2 ( k 2 2 ) 2 k 2 i 2 + i 2 ( i 2 2 ) + 1 ] λ k , $$ \begin{aligned}&c_i \lambda _i^{(3)} = \dfrac{1}{4} \sum _{\begin{matrix} k=i+3 \\ k+i \equiv 1 [2] \end{matrix}}^{+\infty } k \left[ k^2 \left( k^2 - 2 \right) - 2 k^2 i^2 + i^2 \left( i^2 - 2 \right) + 1 \right] \lambda _k~, \end{aligned} $$(C.19)

c i λ i ( 4 ) = 1 24 k = i + 4 k + i 0 [ 2 ] + k [ k 2 ( k 2 4 ) 2 3 k 4 i 2 + 3 k 2 i 4 i 2 ( i 2 4 ) 2 ] λ k . $$ \begin{aligned}&c_i \lambda _i^{(4)} = \dfrac{1}{24} \sum _{\begin{matrix} k=i+4 \\ k+i \equiv 0 [2] \end{matrix}}^{+\infty } k \left[ k^2 \left( k^2 - 4 \right)^2 - 3 k^4 i^2 + 3 k^2 i^4 - i^2 \left( i^2 - 4 \right)^2 \right] \lambda _k~. \end{aligned} $$(C.20)

From these expressions it directly follows that

M ij ( 2 ) = { 1 c i j ( j 2 i 2 ) if j i + 2 and i + j 0 [ 2 ] 0 otherwise , $$ \begin{aligned} \mathcal{M} _{ij}^{(2)} = {\left\{ \begin{array}{ll} \dfrac{1}{c_i} j \left( j^2 - i^2 \right)&\mathrm{if} ~~ j \geqslant i+2 ~ \mathrm{and} ~ i+j \equiv 0 ~ [2] \\ 0&\mathrm{otherwise} ~, \end{array}\right.} \end{aligned} $$(C.21)

and

M ij ( 3 ) = { 1 24 c i j [ j 2 ( j 2 4 ) 2 3 j 4 i 2 + 3 j 2 i 4 i 2 ( i 2 4 ) 2 ] if j i + 4 and i + j 0 [ 2 ] 0 otherwise . $$ \begin{aligned} \mathcal{M} _{ij}^{(3)} = {\left\{ \begin{array}{ll} \dfrac{1}{24 c_i} j \left[ j^2 \left( j^2 - 4 \right)^2 - 3 j^4 i^2 + 3 j^2 i^4 - i^2 \left( i^2 - 4 \right)^2 \right] \\ \mathrm{if} ~~ j \geqslant i+4 ~ \mathrm{and} ~ i+j \equiv 0 ~ [2]&\\ 0 \mathrm{otherwise}&~. \end{array}\right.} \end{aligned} $$(C.22)

C.2. Projection of the products

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

f ( ξ ) = i = 0 + λ i T i ( ξ ) , $$ \begin{aligned}&f(\xi ) = \sum _{i=0}^{+\infty } \lambda _i T_i(\xi )~, \end{aligned} $$(C.23)

g ( ξ ) = i = 0 + μ i T i ( ξ ) . $$ \begin{aligned}&g(\xi ) = \sum _{i=0}^{+\infty } \mu _i T_i(\xi )~. \end{aligned} $$(C.24)

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

T n ( ξ ) e j n arccos ξ , $$ \begin{aligned} \widetilde{T}_n(\xi ) \equiv e^{{j} n \arccos \xi }~, \end{aligned} $$(C.25)

so that we have

2 T i = T i + T i . $$ \begin{aligned} 2T_i = \widetilde{T}_i + \widetilde{T}_{-i}~. \end{aligned} $$(C.26)

It directly follows that

2 f = i = + c | i | λ | i | T i , $$ \begin{aligned}&2f = \sum _{i=-\infty }^{+\infty } c_{|i|} \lambda _{|i|} \widetilde{T}_i~, \end{aligned} $$(C.27)

2 g = i = + c | i | μ | i | T i . $$ \begin{aligned}&2g = \sum _{i=-\infty }^{+\infty } c_{|i|} \mu _{|i|} \widetilde{T}_i~. \end{aligned} $$(C.28)

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

4 f g = i = + j = + c | i | c | j | λ | i | μ | j | T i + j , $$ \begin{aligned} 4fg = \sum _{i=-\infty }^{+\infty } \sum _{j=-\infty }^{+\infty } c_{|i|} c_{|j|} \lambda _{|i|} \mu _{|j|} \widetilde{T}_{i+j}~, \end{aligned} $$(C.29)

which can be rewritten as

4 f g = i = + ( j = + c | j | c | i j | λ | j | μ | i j | ) T i . $$ \begin{aligned} 4fg = \sum _{i=-\infty }^{+\infty } \left( \sum _{j=-\infty }^{+\infty } c_{|j|} c_{|i-j|} \lambda _{|j|} \mu _{|i-j|} \right) \widetilde{T}_i~. \end{aligned} $$(C.30)

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

f g = i = 0 + ( 1 2 c i j = + c | j | c | i j | λ | j | μ | i j | ) T i . $$ \begin{aligned} fg = \sum _{i=0}^{+\infty } \left( \dfrac{1}{2c_i} \sum _{j=-\infty }^{+\infty } c_{|j|} c_{|i-j|} \lambda _{|j|} \mu _{|i-j|} \right) T_i~. \end{aligned} $$(C.31)

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

U = 1 2 U ¯ ( T 0 + T 2 ) , $$ \begin{aligned} U = -\dfrac{1}{2}\overline{U} \left( T_0 + T_2 \right)~, \end{aligned} $$(C.32)

so that

f U = U ¯ 4 i = 0 + 1 c i ( λ i + 2 + 2 c i λ i + c | i 2 | λ | i 2 | ) T i . $$ \begin{aligned} fU = -\dfrac{\overline{U}}{4} \sum _{i=0}^{+\infty } \dfrac{1}{c_i} \left( \lambda _{i+2} + 2c_i\lambda _i + c_{|i-2|}\lambda _{|i-2|} \right) T_i~. \end{aligned} $$(C.33)

From this expression directly stem the following

M ij ( 4 ) = U ¯ 4 c i ( δ i + 2 , j + 2 c i δ ij + c | i 2 | δ | i 2 | , j ) , $$ \begin{aligned}&\mathcal{M} _{ij}^{(4)} = -\dfrac{\overline{U}}{4c_i} \left( \delta _{i+2,j} + 2 c_i \delta _{ij} + c_{|i-2|} \delta _{|i-2|,j} \right)~, \end{aligned} $$(C.34)

M ij ( 5 ) = U ¯ 4 c i ( M i + 2 , j ( 2 ) + 2 c i M ij ( 2 ) + c | i 2 | M | i 2 | , j ( 2 ) ) . $$ \begin{aligned}&\mathcal{M} _{ij}^{(5)} = -\dfrac{\overline{U}}{4c_i} \left( \mathcal{M} _{i+2,j}^{(2)} + 2 c_i \mathcal{M} _{ij}^{(2)} + c_{|i-2|} \mathcal{M} _{|i-2|,j}^{(2)} \right)~. \end{aligned} $$(C.35)

All Figures

thumbnail Fig. 1.

Spatial turbulent spectrum and turbulent frequency. Left: Spatial turbulent spectrum, Eiso(k), extracted from solar observations (solid blue line). The spectrum comprises roughly three sections, each with a different power law in kα. The dashed orange line shows the best fit of the data to a three-piece-wise power law. The exponents and the boundaries between the three regimes are indicated on the plot. Right: Turbulent frequency ωk ≡ kuk as a function of wavenumber, k, where uk is the typical velocity of the turbulent eddies of inverse size k, given by Eq. (47).

In the text
thumbnail Fig. 2.

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

In the text
thumbnail Fig. 3.

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

In the text
thumbnail Fig. 4.

Same as Fig. 3, but for m = 5.

In the text
thumbnail Fig. 5.

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

In the text
thumbnail Fig. 6.

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

In the text
thumbnail Fig. 7.

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

In the text
thumbnail Fig. 8.

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

In the text
thumbnail Fig. 9.

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

In the text
thumbnail Fig. 10.

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

In the text
thumbnail Fig. 11.

Same as Fig. 10, but for the power spectrum of the azimuthal velocity, ux.

In the text
thumbnail Fig. 12.

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

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.