© J. Philidet et al. 2022

1. Introduction

Stellar oscillations are a powerful tool for probing the interior of stars. In solar-type stars, these oscillations have a finite lifetime, meaning that it is possible to measure not only the frequency of the resonant modes, but also their line profile in the Fourier domain. While the frequencies give invaluable information about the equilibrium structure of the star, the line profiles carry the signature of the energetic aspects of the modes, which are tightly related to the physics of the turbulence occurring at the top of the convective envelope. As a result, solar-like oscillations offer a unique opportunity to constrain the properties of stellar turbulent convection, which remains to this day one of the major challenges in stellar physics (e.g., Kupka & Muthsam 2017).

The coupling between solar-like oscillations and turbulent convection manifests itself in the seismic observables in several ways. First, it is responsible for the discrepancy between the observed p-mode frequencies and those theoretically computed from oscillation calculations – commonly referred to as the ‘surface effect’ (Dziembowski et al. 1988; Christensen-Dalsgaard et al. 1996; Rosenthal et al. 1999). This discrepancy is the most important hurdle on the road to a full exploitation of the asteroseismic diagnosis potential of solar-like p-mode frequencies, and it prevents us from accurately inferring the internal equilibrium structure of these stars or their global parameters. It is therefore of utmost importance to correct the surface effect, a task to which numerous studies have been devoted, using either theoretical prescriptions (e.g., Gabriel et al. 1975; Balmforth 1992b; Houdek 1996; Rosenthal et al. 1999; Grigahcène et al. 2005) or empirical formulations (e.g., Kjeldsen et al. 2008; Christensen-Dalsgaard 2012; Ball & Gizon 2014; Sonoi et al. 2015).

Secondly, the turbulent convective motions in the superadiabatic region are responsible for both the stochastic excitation and the linear damping of solar-like oscillations. As a result, the signature of turbulent convection is also carried by the observed amplitudes and linewidths of the modes, and the latter can be used to constrain the former. This task calls for the development of theoretical prescriptions for both the stochastic driving (e.g., Goldreich & Keeley 1977a,b; Balmforth 1992a,c; Samadi & Goupil 2001; Chaplin et al. 2005; Samadi et al. 2005, 2006; Belkacem et al. 2006, 2008, 2010) and the linear damping (e.g., Goldreich & Kumar 1991; Balmforth 1992a; Grigahcène et al. 2005; Dupret et al. 2006; Belkacem et al. 2012) of solar-like p-modes.

But because the coupling between turbulent convection and solar-like oscillations occurs predominantly in the superadiabatic region, where heat transfers, convective motions, and oscillations all share the same typical spatial scales and timescales, establishing a theoretical prescription for this coupling constitutes a challenging task. As a result, simplifying assumptions are very often necessary to obtain tractable relations. Many previous theoretical efforts are based on time-dependent mixing-length theories (e.g., Gabriel et al. 1975; Houdek 1996; Grigahcène et al. 2005; Sonoi et al. 2017; Houdek et al. 2017, 2019), which amounts to reducing the entire turbulent flow to a single typical spatial scale and timescale, in direct contradiction with the picture of the turbulent cascade in high Reynolds number flows. This constitutes a major flaw of these formalisms, as the turbulent pressure and the turbulent dissipation both derive directly from this turbulent cascade and are deemed to play a major role in turbulence–oscillation coupling. Furthermore, while the mixing-length hypothesis is relevant in the bulk of the convective envelope, it ceases to be valid closer to the surface of the star, as highlighted by 3D hydrodynamic simulations (see Nordlund et al. 2009, for a review). Another problem with these approaches is the necessity to separate the equations of hydrodynamics into distinct sets of equations for the turbulent convective motions and the oscillations, respectively, either in the form of a separation between typical wavenumbers (e.g., Grigahcène et al. 2005) or by performing horizontal averaging in a simulation (e.g., Nordlund & Stein 2001). While other approaches have been investigated as an alternative, such as Reynolds-stress models (e.g., Xiong et al. 2000) or the direct extraction of modes from 3D large-eddy simulations (e.g., Belkacem et al. 2019; Zhou et al. 2020), no ideal solution seems to stand out. Therefore, it remains imperative to seek a more adapted theoretical framework in which to describe this coupling.

In this series of papers, we investigate a novel modelling approach designed to address the issue of turbulence–oscillation coupling in the context of solar-like oscillators, based on Lagrangian stochastic models of turbulence. In the first paper of this series (Philidet et al. 2021), we used a Lagrangian stochastic model to obtain a linear stochastic wave equation that is representative of solar-like p-modes and simultaneously contains the effect of the turbulence on the modes. This linear stochastic wave equation constitutes our baseline theoretical framework for quantifying turbulence–oscillation coupling. In this second paper, we exploit it to directly relate the asteroseismic observables – namely, the amplitudes and linewidths of the modes, as well as the surface effect incurred by their eigenfrequencies – to the statistical properties of the underlying turbulence. The present formalism differs from previous approaches in several ways: (1) the turbulence is included in its most general form, which allows for any prescription of turbulent convection to be tested against observations; (2) the treatment of turbulent dissipation is more realistic than in mixing-length theory, which allows its impact on turbulence–oscillation coupling to be quantified; (3) all seismic observables (surface effect, mode amplitudes, and linewidths) are quantified simultaneously, thus strengthening the constraints we can put on the properties of turbulent convection.

This paper is structured as follows. In Sect. 2 we recall the linear stochastic wave equation we obtained in the first paper of this series (Philidet et al. 2021), and we give a succinct physical interpretation of the various terms that it exhibits. In Sect. 3 we detail how this stochastic wave equation can be reduced to coupled equations for the complex amplitude of all the individual normal modes of the system, through the application of the simplified amplitude equation formalism, and we explicitly derive this amplitude equation in the case of solar-like p-modes. In doing so, we keep the simplifying assumption that the entropy fluctuations of the gas can be neglected. In Sect. 4 we consider the single-mode case, which, for main-sequence solar-like oscillators, is likely to constitute the dominant contribution, even in the presence of mode coupling. This yields simultaneous analytical expressions for the excitation rate, the damping rate, and the modal surface effect associated with any given mode, as an explicit function of the statistical properties of the underlying turbulent velocity field in their most general form, and allows us to separate the contributions to these quantities. To illustrate the usefulness of the present formalism, we apply it to a simplified case in Sect. 4.4, where the relative importance of the physical contributions to both mode damping and the modal surface effect can be quantified. Finally, we summarise our findings and draw our conclusions in Sect. 5.

2. The stochastic wave equation

In this section we summarise the results obtained in Philidet et al. (2021, hereafter referred to as Paper I). Paper I established a stochastic linear wave equation for solar-like acoustic modes, which consistently encompass the effect of the turbulence, and therefore naturally contains the information on the coupling between turbulence and oscillations. As such, this linear wave equation is stochastic by nature (i.e. it contains an intrinsically random part), and the stochastic part represents the effect of the random turbulent fluctuations on the behaviour of solar-like p-modes. To derive the linear stochastic wave equation, we adopt the three following steps.

First, we consider a Lagrangian stochastic model of turbulence, whereby the flow is represented by a large set of individual fluid particles. As a first step, these fluid particles are only characterised by their position and velocity. This Lagrangian stochastic model is then modified to yield stochastic differential equations for Eulerian variables, which are more suited to the description of oscillations. This operation leads us to two stochastic differential equations governing the evolution of the fluid displacement ξ(x, t) and velocity u(x, t) (see Eqs. (25) and (26) of Paper I).

Second, these two equations rely on the estimate of the instantaneous ensemble average of several fluid quantities – namely the density, pressure, velocity and Reynolds-stress tensor. These instantaneous ensemble averages are estimated directly from the set of fluid particles, in a similar fashion to the smoothed particle hydrodynamics (SPH) formalism. In particular, this step requires the introduction of a kernel function K(r), which serves as a spatial weighting function designed to filter the neighbouring fluid particles for the estimation of the local instantaneous ensemble averages. This step leads us to explicit expressions for these averages, as a function of ξ(x, t) and u(x, t) only (see Eqs. (32)–(35) in Paper I).

Third, we then treat the turbulent part of the displacement variable ξ_t (resp. u_t for the velocity variable) as an external input, while the residual ξ_osc ≡ ξ − ξ_t (resp. u_osc ≡ u − u_t) represent the oscillatory component of these variables, including their coupling with turbulence. We note that the equations themselves are not split into turbulence equations and oscillation equations, so that we retain one system of equations for both at the same time, instead of two separate sets. The full system of equations – namely the fluid equations obtained in step (1) and the mean field expressions obtained in step (2) – is then linearised in terms of the two wave variables ξ_osc and u_osc. This is done by adopting a certain number of hypotheses and approximations, all of which are explicitly itemised in Paper I (see hypotheses H1 to H5 in Sect. 3.1 therein, as well as H6 to H8 in Sect. 4 therein). Of particular importance is the adiabatic approximation, which we adopted as a first step, and whereby we considered any fluid parcel to conserve its entropy during its evolution. This freed us from having to include an entropy or energy equation in the system, and we instead adopted a polytropic relation between gas pressure and density. In the following, this is what we mean by ‘adiabatic’. By contrast, we do not consider the oscillations to be adiabatic in the mechanical sense of the term, meaning that they can still exchange energy with the background, thus allowing the contributions of both the driving and the damping of the modes to be included (see discussion in Sect. 4).

Ultimately, this allows us to establish the following linear stochastic wave equation for solar-like acoustic modes in the presence of turbulent convection (see Eqs. (42) and (43) in Paper I):

(1)

(2)

The vectors L₀, , and are given by (see Eqs. (44)–(46) in Paper I)

(3)

(4)

and

(5)

where we have defined

(6)

The indices ‘0’ and ‘1’ refer to equilibrium values and fluctuations around the equilibrium values, respectively. The notation ∂_i denotes the ith component of the gradient operator, and we have used Einstein summation convention. The dimensionless quantity C₀ denotes the Kolmogorov constant, for which a commonly accepted experimental value is C₀ = 2.1 (e.g., Haworth & Pope 1986). We have introduced the centred kernel function K^x(y)≡K(y − x), and the perturbations of the Reynolds-stress tensor and mean shear tensor are given by (see Eqs. (47) and (48) in Paper I)

(7)

and

(8)

The quantities ρ₀, p₀, and γ are the equilibrium density, equilibrium gas pressure, and polytropic exponent, respectively, and is the equilibrium sound speed squared. Finally, by construction, η_i(x, t) is defined as a multivariate Gaussian process whose values at two distinct locations are completely uncorrelated, and which verifies

(9)

(10)

where δ_ij is the Kronecker symbol, δ is the Dirac distribution, and .̄ denotes the ensemble average.

Equations (1) and (2) put together govern the time evolution of the wave variables ξ_osc and u_osc as an explicit function of either the wave variables themselves or the turbulent fields ξ_t and u_t. More specifically, Eq. (2) is split three ways. contains all the terms that are linear in the wave variables, but independent of the turbulent fields (they are therefore deterministic). contains all the terms that are linear in the wave variables, and also depend on the turbulent fields (they are therefore stochastic). Finally, L₀ contains all the terms that are independent of the wave variables, and encompass all the inhomogeneous forcing terms.

Formally, the oscillations can be described in a Hilbert space of infinite dimension, where the ket |z⟩, in Dirac notation, is defined as¹

(11)

where 𝒱 is the volume of the star, and ω refers to the angular frequency of the mode under consideration. Then, Eqs. (1) and (2) can be rewritten in the following way:

(12)

where the deterministic, time-independent linear operator ℒ^d is given by

(13)

the stochastic, time-dependent linear operator ℒ^s(t) is given by

(14)

the stochastic vector |θ(t)⟩ is given by

(15)

and we recall that the vectors , , and L₀ are given respectively by Eqs. (3)–(5). We note that, since the space dependence of the wave variables is already contained in the infinite dimension of the Hilbert space where the oscillations are described, |z⟩ only depends on time, such that d/dt is a total time derivative.

The impact of turbulent convection on the frequency and the damping rate of the acoustic modes is contained within the linear stochastic operator ℒ^s(t), whereas the driving of the acoustic modes is contained within the additive noise vector |θ(t)⟩. Because of their stochastic nature, a large portion of these two quantities is actually incoherent or statistically uncorrelated with the acoustic modes themselves, and therefore incapable of efficiently impacting the oscillations. It is therefore necessary to filter in the part of these stochastic perturbations that does lead to a significant impact on the modes. A naive way to do this would be to simply take the ensemble average of Eq. (12) (or, equivalently, of Eqs. (1) and (2)). However, this would not work, because even the zero-mean stochastic part of the turbulent perturbations, being auto-correlated over timescales comparable to the period of the modes, can impact the behaviour of the oscillations. It is therefore necessary to incorporate the effective impact of the finite memory time of the stochastic perturbations on the oscillations. This can be done in the framework of the simplified amplitude equation formalism (Stratonovich et al. 1965; Buchler & Goupil 1984; Buchler et al. 1993).

3. The simplified amplitude equations

The goal of this section is to use the simplified amplitude equation formalism to properly account for the effect of convective noise on the excitation and damping rates of the normal modes of oscillation. To do so, we proceeded in a two-step approach.

First, the linear stochastic wave equation is transformed into a finite set of stochastic differential equations governing the evolution of the complex amplitudes of the normal modes of oscillation. This is done by treating the stochastic part of the wave equation as perturbations to an otherwise deterministic wave equation. This leads to a set of coupled amplitude equations (Buchler & Goupil 1984). We present this first step in Sect. 3.1.

The amplitude equations contain a stochastic part that, as mentioned above, has a finite memory time. In a second step, these are substituted for simplified amplitude equations, which are also stochastic differential equations governing the evolution of the complex amplitude of the modes, but whose stochastic part, unlike the original amplitude equations, has zero memory time, and therefore has no long-term impact on the oscillations. This step involves two procedures. First the original stochastic amplitude equations are substituted for an equation for the probability density function (PDF) of the complex amplitude of the modes, in the form of a Fokker-Planck equation. We present this Fokker-Planck equation in Sect. 3.2. Then another set of stochastic differential equations is constructed, under the constraints that (1) it must reduce to the exact same Fokker-Planck equation, while (2) only involving Markov processes (i.e. stochastic processes with zero memory time). We derive these simplified amplitude equations in Sect. 3.3.

3.1. The full amplitude equations

The stochastic amplitude equations are obtained in two steps (Stratonovich et al. 1965): first deterministic amplitude equations are derived as if there were no stochastic perturbations; then the stochastic part is treated in a perturbative framework, and stochastic corrections are added to the deterministic amplitude equations after the fact. The non-perturbed wave equation is obtained by setting both ℒ^s(t) and |θ(t)⟩ to zero in Eq. (12), which yields

(16)

The linear operator ℒ^d can be diagonalised, and any vector |z⟩ can be decomposed on a basis of eigenvectors for ℒ^d. Because this linear operator is real, the associated eigenvalues are either real or come in pairs of complex conjugates. Each such pair of eigenvectors with complex-conjugated eigenvalues correspond to a single mode of oscillation of the system; we denote them as |Ψ^μ⟩ and |Ψ^†μ⟩ (where the exponent notation μ is used to index the different modes) and their associated eigenvalues as κ_μ ± jω_μ, where j denotes the imaginary unit.

Without loss of generality, we consider that the eigenvectors are normalised to unity,

(17)

where ⟨.|.⟩ refers to the scalar product in the |z⟩ Hilbert space. This scalar product needs to be specified. A necessary condition is that distinct eigenvectors |Ψ^μ⟩ of the non-perturbed wave equation (Eq. (16)) – that is, eigenvectors of the linear operator ℒ^d – must be orthogonal to one another in the sense of this scalar product. For adiabatic, non-radial oscillations, this is verified by the following scalar product (Unno et al. 1989):

(18)

where ψ_ξ(x) (resp. ψ_u(x)) is the part of |ψ⟩ associated with the oscillatory displacement (resp. velocity) at location x. For any eigenvector |Ψ⟩ of the problem, these two quantities are simply the normalised displacement and velocity eigenfunctions. Because of the normalisation condition (Eq. (17)), they are related to the non-normalised eigenfunctions ξ_osc and u_osc through (see Sect. B.4 for more details)

(19)

where ℐ_μ is the inertia of the mode, defined by

(20)

The general solution of Eq. (16) is (Buchler et al. 1993)

(21)

where ‘c.c’ denotes the complex conjugate. In general, |z⟩ would be written as an arbitrary linear combination of |Ψ^μ⟩ and |Ψ^†μ⟩; however, if |z⟩ is initially real, Eq. (16) shows that it will remain real at all later times (because ℒ^d itself is real), so the |Ψ^†μ⟩ component of |z⟩ is necessarily the complex conjugate of its |Ψ^μ⟩ component. As for a_μ(t), it denotes the slowly varying complex amplitude of the mode μ. In the general case, plugging the solution given by Eq. (21) into Eq. (16), one finds the following, very simple amplitude equation for a_μ(t):

(22)

As expected, while the imaginary part ω_μ of the eigenvalue associated with the mode represents its oscillatory angular frequency, its real part κ_μ represents the rate at which its amplitude varies – either a growing rate or a damping rate, depending on the intrinsic stability of the mode in the absence of turbulence. However, in the present adiabatic case, since the operator ℒ^d is Hermitian, we actually have κ_μ = 0. The quantity κ_μ therefore contains the contribution of all non-adiabatic effects to the damping of the modes. As it contributes to defining the variation timescale of the mode amplitudes, we keep this term; however, we do not consider it more explicitly, as our present goal is to estimate the additional contribution coming from the interaction with turbulent convection.

Turning now to the full, stochastically perturbed wave equation (Eq. (12)), the addition of the two terms ℒ^s(t)|z⟩ and |θ(t)⟩ leads to the following modifications to the amplitude equations (Buchler et al. 1993):

(23)

where

(24)

The form of the stochastic linear operator ℒ_s(t) and stochastic vector |θ(t)⟩ being completely specified by Eqs. (14) and (15), respectively, and the scalar product being defined by Eq. (18), the stochastic processes c_i(t) (i = 1, 2, 3) constitute ‘knowns’ of the system. These are given by

(25)

(26)

where is obtained by replacing ξ_osc and u_osc by the eigenfunctions and in Eq. (4).

It will be more practical, in the following, to separate the complex amplitude of the modes into their amplitude A_μ(t) and phase Φ_μ(t), both being real functions of time, such that

(27)

In place of a single complex stochastic equation for the evolution of a_μ(t), we obtain two real stochastic equations for the evolution of A_μ(t) and Φ_μ(t), which take the following form:

(28)

and

(29)

where Re and Im denote the real and imaginary parts, respectively. We warn the reader that the function G_μ(A_ν, Φ_ν, t) defined by the right-hand side of Eq. (28) must not be confused with the drift tensor G_ij appearing in Eqs. (3) and (4).

3.2. The Fokker-Planck equation for mode amplitude

A fundamental hypothesis we need to make concerning ℒ^s(t) and |θ(t)⟩ is that their correlation timescale is very small compared to the timescale over which the amplitude and the phase of the mode typically vary. In the context of solar-like oscillations, the correlation timescale of these stochastic perturbations correspond to the turnover time of the turbulent eddies, which, close to the surface of the star, is similar to the period of the p-modes (i.e. ∼5 min). On the other hand, the amplitude and phase of the modes vary over typical timescales that correspond to their life time, which is indeed much longer than their period (τ ∼ 3 h for the shortest-lived solar modes). As such, this hypothesis is largely verified in solar-like oscillators. Consequently, all the turbulent perturbations in the wave equation (Eq. (12)) can be approximated by Markov processes (in the sense that their memory time, while finite, is much smaller than the evolution timescale of the amplitude of the modes), and subsequently, so can the stochastic processes A_μ(t) and Φ_μ(t). It is then possible to replace the set of 2N (where N is the total number of modes) stochastic differential equations on these quantities (i.e. Eqs. (28) and (29)) with an equivalent, single Fokker-Planck equation governing the evolution of their joint PDF w(A_μ, Φ_μ, t), but whose coefficients are carefully defined to incorporate in a rigorous manner the effect of the finite memory time of the processes A_μ(t) and Φ_μ(t). The Fokker-Planck equation takes the general form (Stratonovich et al. 1965)

(30)

where 𝒢_μ and ℋ_μ represent the probability fluxes in (A_μ, Φ_μ) space, and 𝒟_μν, ℰ_μν and ℱ_μν are the elements of the 2N × 2N probability diffusion matrix. The explicit computation of the probability fluxes and diffusion coefficients is detailed in Appendix A, and these coefficients read

(31)

(32)

(33)

(34)

(35)

where ‘sym.’ refers to a swapping between indices μ and ν, and we have introduced the following autocorrelation spectra:

(36)

(37)

(38)

(39)

(40)

(41)

(42)

It can already be seen, without having to specify the stochastic processes c_i(t), that their autocorrelation spectra are independent of the two stochastic variables A and Φ. Furthermore, under the assumption that the turbulence characterising the convective motions at the top of the envelope of solar-like oscillators is stationary, the autocorrelation functions appearing in the definition of these coefficients are only dependent on the time increment τ, and not on the absolute time t. As such, the α_i are also independent of time t. All things considered, they are therefore simply complex constants.

3.3. The equivalent simplified amplitude equations

As we have mentioned above, either the stochastic differential equations (Eqs. (28) and (29)) together or the Fokker-Planck equation (Eq. (30)), both of which are equivalent to one another, can be used to model the time evolution of the real amplitude and phase of the modes. However, both are equally impractical to use, albeit for different reasons. Numerically integrating the stochastic equations proves extremely expensive, because a large range of very different timescales must be resolved. Indeed, the total integration time must far exceed the typical timescale of the slowly varying mode amplitude. But at the same time, the rapidly varying phase ωt appearing in Eqs. (28) and (29) must be accurately resolved. Finally, the whole range of memory timescales associated with the processes c_i(t) must also be resolved, which is problematic, as it corresponds to the range of timescales in the turbulent cascade, and is therefore very wide in high Reynolds number flows. Consequently, the numerical integration of Eqs. (28) and (29) requires an unreasonably small time step compared to the total integration time.

By contrast, the Fokker-Planck equation (Eq. (30)) does not have this timescale problem. Indeed, in computing its probability fluxes and diffusion coefficients, we have filtered out all rapidly oscillating features (see Appendix A for more details), and as such, the numerical integration of the Fokker-Planck equation does not require these very short timescales to be resolved. However, it poses other difficulties inherent to the integration of Fokker-Planck equations in general. Indeed, because the PDF w is a function of the entire parameter space (i.e. not only of time, but also of the stochastic variables A_μ and Φ_μ), its numerical integration would require the discretisation of all three variables. For that reason, it is usually very impractical to directly integrate the Fokker-Planck equation in time.

Fortunately, while a given stochastic differential equation is equivalent to a single Fokker-Planck equation, the opposite is not true. Indeed, there exists an infinite number of stochastic models for the mode amplitudes A_μ and phases Φ_μ that possess the same Fokker-Planck equation (Eq. (30)) and therefore contain the exact same statistical information as Eqs. (28) and (29) while having a much simpler form. They are the simplified amplitude equations (Stratonovich et al. 1965; see also Buchler et al. 1993)

(43)

(44)

where 𝒟^1/2 and ℱ^1/2 are the square-root of the positive definite matrices² 𝒟 and ℱ, and dW_Aν and dW_Φν are the increment over dt of independent Wiener processes, that is, Gaussian processes with zero mean and autocorrelation ⟨dW_AμdW_Aν⟩=⟨dW_ΦμdW_Φν⟩=δ_μνdt. The symbol ∘ is the differential notation of the Stratonovich stochastic integral. We note that, in order for the equivalence between the Fokker-Planck equation (Eq. (30)) and the stochastic differential equations (Eqs. (43) and (44)) to be justified, a sufficient condition is that the probability fluxes 𝒢_ν and ℋ_ν, as well as the elements of the diffusion matrix 𝒟_μν and ℱ_μν, only depend on time t and on the variables A_μ and Φ_μ themselves. In fact, in the present case, the situation is even simpler, since these coefficients actually only depend on the real amplitudes A_ν.

These equations are much more practical to handle than either the exact amplitude equations derived in Sect. 3.1 or the Fokker-Planck equation derived in Sect. 3.2, in the sense that they allow us to circumvent the problems outlined above. Indeed, similarly to the Fokker-Planck equation, all rapidly oscillating terms have been averaged out, which means that the short timescales ∼ω⁻¹ need not be resolved. Furthermore, the stochastic part of these equations now has zero memory – because the effect of the finite width of the turbulent cascade timescale range has been properly and rigorously incorporated in the coefficients of the Fokker-Planck equation – which also drastically reduces the range of timescales that we need to resolve. In addition, Eqs. (43) and (44) are much easier to integrate numerically than the corresponding Fokker-Planck equation, as only the time variable needs to be discretised. The only cost is that a large number of independent realisations must be integrated in order to reconstruct the moments of the amplitude and phase of the modes.

4. Single-mode case

The final result of Sect. 3 consists in the coupled simplified amplitude equations (Eqs. (43) and (44)). One important aspect of these equations is the fact that they are coupled: even in the limit of small amplitudes, where the modes behave linearly, the non-linear nature of the turbulent convection with which they are coupled makes it possible to transfer energy from mode to mode, and for their respective phases to evolve in a dependent manner. This was already noted by previous studies concerning mode scattering in solar-like oscillators (e.g., Goldreich & Murray 1994) or non-linear mode interactions (e.g., Kumar & Goldreich 1989). This can potentially have an impact on the energy distribution across the p-mode spectrum of the star, and may also impact their frequency and their lifetime. That being said, as a first application of the present formalism, we temporarily discard the question of mode coupling and instead focus on the single-mode case. The coupled-mode case, to which we briefly return at the end of Sect. 5, will be the subject of a later paper in this series.

4.1. Simplified amplitude equations in the single-mode case

In the single-mode case, we obtain

(45)

(46)

where the probability fluxes and diffusion matrix elements become

(47)

(48)

(49)

(50)

and the α_i are defined by

(51)

(52)

(53)

(54)

ω now being the angular frequency of the single mode under consideration. Plugging Eqs. (47)–(50) into Eqs. (45) and (46), we obtain

(55)

(56)

Alternatively, we can merge these two equations into a single stochastic differential equation on the complex amplitude a = Aexp(jΦ) of the mode, which yields³

(57)

While Eqs. (55) and (56) – or, equivalently, Eq. (57) – feature all of the above autocorrelation spectra, α_i, the first terms on their respective right-hand sides (which, as will be discussed in more detail in Sects. 4.2 and 4.3, governs the mean evolution of the mode amplitudes) only depend on α₁ and α₃. These two autocorrelation spectra can be computed by plugging the explicit expressions of (Eq. (4)) and L₀ (Eq. (5)) into the expressions for c₁(t) and c₃(t) (Eqs. (25) and (26)), and then the latter into Eqs. (51) and (54). Since c₁(t) and c₃(t) depend on the turbulent velocity field, u_t, their autocorrelation spectrum can be described in terms of the autocorrelation spectrum of the turbulent velocity itself. We introduce

(58)

where

(59)

Then we can define the second-order turbulent velocity spectrum as

(60)

as well as the following fourth-order spectra

(61)

(62)

(63)

(64)

The derivation of the autocorrelation spectra α₁ and α₃ is detailed in Appendix B, and we eventually obtain

(65)

and

(66)

where we recall that the operator ℑ is defined by Eq. (6), k is the (space-dependent) wave vector of the mode, and the functions F⁽¹⁾, F⁽²⁾, F^(3a) and F^(3b) are given by⁴

(67)

(68)

(69)

and

(70)

These four functions can be broken down into several contributions, thus allowing us to separate the effect of each physical process to the turbulence–oscillation coupling. In the present case, one can identify the turbulent pressure contribution, represented by the advection term in the stochastic velocity equation (see Eq. (26) in Paper I), and which is responsible for the entirety of α₁, as well as the first two terms on the right-hand side of Eq. (67). On the other hand, all the other contributions to α₃ stem from the G_ij term in the stochastic velocity equation, which we recall encompasses the collective effect of the buoyancy force, the fluctuating gas pressure force, and the turbulent dissipation (see Paper I, and more specifically Appendix A therein, for more details). We also note that, Eqs. (65) and (66) being written in the form of spatial integrals spanning across the entire volume of the star, this formalism allows us to determine which regions of the star are most responsible for the coupling between the turbulent motions and the oscillatory motions.

4.2. Mode driving

Multiplying Eq. (55) by 2A and taking the ensemble average, we obtain the following evolution equation for the mean⁵ energy of the mode E_m(t)≡⟨A(t)²⟩:

(71)

We note that, because W_A(t) is a Wiener processes, the diffusion part of Eq. (55) (i.e. the second, stochastic term on its right-hand side) is necessarily of zero mean, even when the diffusion coefficient depends explicitly on the stochastic variables themselves (e.g., Gardiner 1994). The right-hand side of Eq. (71) contains a linear contribution that corresponds to the linear damping of the mode and an additive term, 𝒫 ≡ Re(α₁), which corresponds to the excitation rate of the mode (i.e the amount of energy injected into the mode per unit time). We focus on the latter for the moment. It can be seen that 𝒫 only contains contributions from the turbulence-induced stochastic perturbation of the wave equation (Eq. (12)), in the sense that if c₁(t) was identically zero, 𝒫 would also be zero. This is in accordance with the widely acknowledged picture of solar-like oscillations being stochastically excited by highly turbulent motions of the plasma at the top of the convection zone.

Using Eq. (54), the excitation rate 𝒫 of the mode becomes

(72)

where k and ω denote the wave vector and angular frequency of the mode. This expression is fully similar to the formulation obtained by previous models for the excitation of solar-like oscillations (see, for instance, Samadi & Goupil 2001; Chaplin et al. 2005), which gives further support to the consistency and validity of the method presented here. First, it is inversely proportional to the inertia of the mode: the larger the mass flow entailed by an oscillating mode, the harder it is to get the flow to actually move. Secondly, it appears that the efficiency of mode driving by turbulence is directly related to the spectrum of the fourth-order moment of turbulent velocity. This is in accordance with widely accepted results of past studies. Finally, the integrand is weighted by a quantity that takes the form of the product of two different components of the wave vector and two different components of the modal fluid displacement. Qualitatively, it can be seen that this is closely related to the square of the mode compressibility |∇ ⋅ ξ_osc|. In other words, the turbulence drives the mode much more efficiently in regions where the compressibility of the mode is high. This is also not surprising, as the turbulent pressure must actually be able to transfer mechanical work to the mode in order to give it energy, and mechanical work is only transferable if the flow undergoes successive compression and dilatation phases.

4.3. Mode damping and the modal surface effect

We now turn to the first term on the right-hand side of Eq. (71). Solar-like oscillations being intrinsically stable, it gives the damping rate of the mode:

(73)

The total damping rate, η_tot, is directly related to the observed linewidth at half maximum Γ_obs of the resonant peak corresponding to the mode in the Fourier power spectrum, through

(74)

In the Sun, Γ_obs typically has the order of a few μHz, while the eigenfrequency of the modes is much larger (∼3000 μHz at the maximum spectral height). Eq. (73) shows that the total damping rate of the mode – and therefore its observed linewidth – consists of two different contributions. The first one corresponds to the non-turbulent, deterministic contribution κ. It represents the effect of the non-adiabatic energy exchanges between the p-modes and the background in which they develop. Explicitly modelling κ requires the inclusion of non-adiabatic effects in the model, through an energy equation, which is outside the scope of this study (see the discussion in Sect. 5 however). The second one represents the contribution from turbulence, and encompasses a whole array of individual physical contributions from turbulent pressure, turbulent dissipation or convective enthalpy flux for instance. In the following, we only investigate this second contribution, and which we denote as η_turb.

As a side note, we remark that taking dE_m/dt = 0 in Eq. (71) (i.e. assuming that the transient has died away and that the mode has reached a stationary state) yields the following expression for the stationary average energy of the mode:

(75)

This is in perfect accordance with Eq. (5) of Goldreich & Keeley (1977a) for instance, and reflects the well-established picture of the equilibrium energy of the mode resulting from a balance between driving and damping processes.

In parallel, taking the ensemble average of Eq. (56), we obtain the following simple evolution equation for the mean phase Φ_m(t)≡⟨Φ(t)⟩:

(76)

which is easily integrated to yield

(77)

where Φ_m, 0 is an arbitrary initial average phase. In turn, this yields the following expression for the average global phase of the mode ϕ(t)≡ωt + Φ_m(t):

(78)

This amounts to a systematic shift, δω, in the angular frequency of the mode compared to the angular frequency in the absence of turbulence. This shift actually represents what is commonly referred to as the modal – or ‘intrinsic’ – part of the surface effect (Balmforth 1992b), that is, the contribution of turbulent convection not to the equilibrium structure but to the propagation of the waves themselves.

All in all, the turbulence-induced damping rate and the modal surface effect turn out to correspond respectively to the real and imaginary parts of the same complex autocorrelation spectrum α₃, according to

(79)

(80)

This is not quite surprising, as the damping rate and oscillatory frequency of a mode are themselves simply two sides of the complex eigenvalue associated with the mode – more specifically, its real and imaginary parts. This is further illustrated by the fact that, in the first term on the right-hand side of Eq. (57), the linear contribution (i.e. the one proportional to a) involves directly the complex quantity α₃. This is a very interesting output of this model as mode excitation on the one hand and mode damping and the surface effect on the other hand are usually treated separately. Furthermore, while mixing-length formalisms allow mode damping rates and the surface effect to be investigated simultaneously, they must nevertheless be the subject of separate fitting procedures, owing to the large number of free parameters involved (see for example Houdek et al. 2017). By contrast, the model presented here allows us to treat all three aspects – mode excitation, in addition to mode damping and the surface effect – in the same consistent framework, and without having to resort to separate procedures.

Using the expression for α₃ given by Eq. (66), we find

(81)

(82)

where we have defined

(83)

and we recall that the wave vector, k, and the angular frequency, ω, are those of the mode under consideration.

Some general results can be drawn from Eqs. (81) and (82). First it can be seen that both quantities go as ℐ⁻²: similarly to the excitation rate 𝒫, the higher the mass flow pertaining to the oscillation, the harder it is to take energy from it. This is due to the fact that the wave variable ξ_osc appears to the fourth power in Eqs. (81) and (82). Since the inertia of the mode depends on the square of ξ_osc, the normalisation condition on the eigenfunction naturally involves the square of the inertia (see Sect. B.4 for more details). This is similar to the work integral formulation of the damping rate (e.g., Samadi et al. 2015), where the wave variables appear to the second power, and therefore the damping rate goes as ℐ⁻¹ instead. We note that the reason why the integrals on the right-hand side of Eqs. (81) and (82) involve the fourth power of the wave variables is that we compute the mean effect of the damping on the mode, which corresponds to the autocorrelation of the instantaneous damping: since the latter depends on the square of the wave variables, the autocorrelation depends on their fourth power.

Secondly, while the driving source is a fourth-order quantity in terms of the turbulent velocity (see Eq. (72) for the excitation rate 𝒫), the damping rate and the modal surface effect, by contrast, have both a second-order and a fourth-order dependence in terms of the turbulent velocity. These are explicitly broken down in Eq. (83).

Finally, it is possible to separate the different physical contributions to the turbulent damping rate η_turb and the modal surface effect δω. To that end, we decomposed the function (defined by Eq. (67)) into three different terms, in the following way:

(84)

where

(85)

(86)

(87)

The first term represents the effect of turbulent pressure; the second term represents the joint effect of the pressure-rate-of-strain correlation and the dissipation of turbulent kinetic energy into heat, both of which are collectively modelled by means of the drift tensor, G_ij; and the third term stems from the variations in the drift tensor, G_ij, under the influence of both the turbulence and the oscillations. By contrast, the functions , , and , defined by Eqs. (68)–(70), respectively, entirely originate from the modulation of the drift tensor. Then Eqs. (81) and (82) can be rewritten in the following way:

(88)

(89)

where the various physical contributions to the total damping rate are given by

(90)

(91)

(92)

(93)

those of the modal surface effect are given by

(94)

(95)

(96)

(97)

and we have defined

(98)

(99)

(100)

(101)

All the functions F⁽¹⁾, F⁽²⁾ and F⁽³⁾ are implicitly evaluated at X, and all the turbulent spectra ϕ⁽²⁾ and ϕ⁽⁴⁾ are implicitly evaluated at (2k, 2ω), even though we dropped the notation for the sake of clarity. The contributions labelled ‘pt’, ‘G’, and ‘δG’ have the same meaning as in Eq. (84), and those labelled ‘cross’ are cross terms between these different contributions. This decomposition shows that the first two contributions only depend on the second-order turbulent velocity spectrum, while the others depend on both the second- and fourth-order spectra.

4.4. Contributions of turbulent dissipation and redistribution versus turbulent pressure

In Sect. 4.3 we derived general expressions for the various physical contributions to mode damping (see Eqs. (90)–(93)) and the modal surface effect (see Eqs. (94)–(97)). The resulting expressions are valid for any specification of the drift tensor G_ij, and regardless of whether the mode is radial or not. While we postpone more realistic applications of these expressions to a later paper in this series, it is still worthwhile to examine their output in the following simple case. First, we assume that the drift tensor is constant, in the sense that it does not depend on the Reynolds-stress tensor, , the mean shear tensor, , or the turbulent dissipation rate, ϵ

(102)

Physically, this means that the rate at which the different components of the turbulent motions exchange energy, or are dissipated into heat, is fixed. We note that this assumption is not very realistic, but that it allows for a simple illustration of the usefulness of this formalism. However, we do not specify how the behaviour of these different components compare to each other, meaning that we do not make any assumption regarding the relative value of the different components of G_ij. Secondly, we consider a radial mode, such that both the displacement eigenfunction and the wave vector of the mode are radial:

(103)

(104)

Then Eqs. (90)–(92) reduce to

(105)

(106)

(107)

and Eqs. (94)–(96) become

(108)

(109)

(110)

We note that, to keep the discussion simple, we omit the crossed contributions η_cross and δω_cross in the following. The only remaining contributions to both the damping rate of the modes and their modal surface effect are the contribution of turbulent pressure, and the contribution of the drift tensor G_ij, which we recall here jointly models the pressure-rate-of-strain correlations and the turbulent dissipation.

These simplified expressions allow us to roughly estimate which process contributes predominantly to both the damping rate and the modal surface effect, depending on the frequency regime. Under the assumption that the region of the star where the damping of the modes and the surface effects occur is a thin subsurface layer, then we can write

(111)

(112)

Explicitly expanding the index summation in the denominator, we obtain

(113)

where the index ‘h’ refers to an arbitrary horizontal direction. Then

(114)

(115)

where the anisotropy factors are defined by

(116)

(117)

(118)

The anisotropy factors Φ_u, η and Φ_u, δω concern the turbulent velocity field. By contrast, the anisotropy factor Φ_G concerns the drift tensor G_ij, and physically represents the ratio between the rate at which vertical and horizontal turbulent motions exchange energy, and the rate at which the vertical motions are dissipated into heat.

Regardless of the relative weight of the different components of G_ij, their order of magnitude is always ∼ϵ/k = ω_t (e.g., Pope 1994), where we recall that ϵ is the turbulent dissipation rate, k is the turbulent kinetic energy, and the ratio between the two defines the turbulent frequency ω_t already introduced above. We note that in the present model, ω_t is assumed constant, but not necessarily uniform. Physically, it represents the lifetime of the energy-bearing eddies in the turbulent cascade, and corresponds to the autocorrelation timescale of granulation. Typically, close to the surface of the Sun, this timescale is τ_conv ∼ 300−500 s (e.g., Samadi et al. 2015), in which case ω_t ∼ 0.01−0.02 rad.s⁻¹. In particular, in the Simplified Langevin Model, every diagonal component of the drift tensor is equal to (e.g., Pope 2000)

(119)

where C₀ = 2.1 is the Kolmogorov constant. Then, we finally obtain

(120)

In the following, we consider indifferently either one of the two ratios R_η or R_δω, which we simply denote as R. The turbulent velocity anisotropy factor will indifferently refer to Φ_u, η or Φ_u, δω, which we simply denote as Φ_u. It can be thought of as representing the ratio between the squared horizontal and vertical turbulent velocities, and can therefore be related to the perhaps more familiar anisotropy factor Φ introduced by Gough (1977) through

(121)

In Fig. 1, we plot the ratio R between the two contributions to either mode damping or the modal surface effect, as a function of mode frequency. While the turbulent frequency ω_t is fixed (we set it to ω_t = 0.03 rad.s⁻¹), we vary the two anisotropy factors Φ_u and Φ_G (see caption). As expected from Eq. (120), the turbulent pressure contribution tends to dominate for high-frequency modes, while the joint effect of the pressure-rate-of-strain correlation and of the turbulent dissipation tends to prevail for low-frequency modes. However, the threshold between these two regimes drastically depends on Φ_u and Φ_G. More specifically, for low values of Φ_G, this threshold is on the low-frequency edge of the solar damping plateau, around 2 mHz, and does not depend on Φ_u. For higher values of Φ_G – and in particular for Φ_G = 1, which corresponds to a vertical-horizontal turbulent energy redistribution rate similar to the vertical turbulent dissipation rate –, this threshold is shifted to higher frequencies, and has a much stronger dependence on Φ_u.

To better illustrate the effect of varying the model parameters, we plot in Fig. 2 the same ratio R, but this time as a function of the drift tensor anisotropy factor Φ_G (left panel) and turbulent frequency ω_t (right panel), for a fixed mode frequency ν = ν_max = 3 mHz. These figures show that, if the turbulent velocity anisotropy factor is fixed to the value observed in large-eddy simulations (i.e. Φ_u = 0.5; (e.g., Samadi et al. 2003)), then the turbulent pressure contribution dominates for Φ_G ≲ 1 and ν_t ≡ ω_t/(2π)≲5 mHz, respectively. On the other hand, for a higher anisotropy factor Φ_G or a higher turbulent frequency ν_t, the turbulent dissipation plays a much more important role in both mode damping and the modal surface effect. While we only considered a simplified case for illustration purposes, this still tends to show how important it is to account for turbulent dissipation in a realistic way when modelling the lifetime of the modes, or when trying to correct the surface effect.

Naturally, this picture is a simplistic one, as it stems primarily from hypothesis (H7) in Paper I: the reduction of the turbulent cascade to a single timescale, . In a more realistic picture, one can assume that within the continuous range of turbulent timescales, some eddies act on the modes through turbulent pressure, and others through turbulent dissipation and pressure-rate-of-strain correlation, depending on the criteria given by Eq. (120) with being replaced with their lifetime. Nevertheless, the results showcased in Figs. 1 and 2 illustrate how the formalism presented in this paper allows us to break down the relative importance of the various physical contributions to mode damping and the modal surface effect.

5. Conclusions

This series of papers aims to investigate Lagrangian stochastic models of turbulence for the study of the coupling between solar-like oscillations and turbulent convection, as an alternative to more traditional methods based on mixing-length theory for example. In Paper I (Philidet et al. 2021) we established a linear stochastic wave equation designed to govern the oscillations while encompassing the full effect of the turbulence thereon. In this second paper we exploit this wave equation to obtain stochastic equations that govern the coupled evolution of the complex amplitudes of all the modes of the system simultaneously. As a first application of these equations, we consider the single-mode case, which yields explicit and simultaneous expressions for the excitation rate (Eq. (72)), the damping rate (Eq. (81)), and the modal – or ‘intrinsic’ – surface effect (Eq. (82)) associated with any radial or non-radial solar-like p-mode, directly as a function of the statistical properties of the underlying turbulent velocity field.

The expressions for the excitation rate, damping rate, and modal surface effect provide a useful insight into the physical origin of the energetic aspects of the solar-like oscillations (at least those stemming from mechanical work), as well as the discrepancy between observed and theoretical p-mode frequencies. The main conclusions are as follows.

1. The seismic observables – amplitudes, linewidths, and the modal surface effect – all simultaneously arise in a coherent manner from the same formalism. Furthermore, they are explicitly related to the statistical properties of the underlying turbulent convection.

2. The excitation rate of the mode (Eq. (72)) reduces to the same theoretical formulation obtained by previous studies on the subject (e.g., Samadi & Goupil 2001; Chaplin et al. 2005), which further supports the validity of this formalism.

3. We also recover the fact that the damping rate and modal surface effect (Eqs. (81) and (82)) are simply the real and imaginary part of the same complex quantity. This means that it is necessary to treat them not as separate phenomena, but as two aspects of the same physical process, which should be reflected in the way they are theoretically predicted or empirically treated.

4. Our formalism allows for the separation of the different physical contributions to both the damping rate (Eq. (88)) and the modal surface effect (Eq. (89)). In the present quasi-adiabatic formalism (i.e. under the assumption that fluid particles conserve their entropy during their motion), we can separate the effect of the turbulent pressure from the joint effect of the pressure-rate-of-strain correlation and that of the turbulent dissipation.

5. We apply this formalism to a simple case where the relative importance of these two contributions (turbulent pressure vs. pressure-rate-of-strain correlation and turbulent dissipation) can be quantified (Eq. (120)) in terms of (1) the turbulent frequency, (2) the degree of anisotropy of the Reynolds-stress tensor, and (3) the degree of anisotropy of the drift tensor, G_ij (i.e. the ratio between the turbulent energy vertical–horizontal redistribution rate and its dissipation rate). We show that the turbulent pressure tends to dominate for high-frequency modes, and that the turbulent dissipation and pressure-rate-of-strain correlation tend to prevail for low-frequency modes.

We stress that the Lagrangian stochastic model of turbulence from which we started in Paper I, and which we have continued to adopt in the present study, was deliberately oversimplified so as to provide a proof-of-concept to illustrate the validity of this formalism and its relevance for the study of turbulence–oscillation coupling in the context of solar-like oscillators. In particular, among the various simplifying assumptions that we adopted in Paper I, two of them stand out. The first is hypothesis H6, according to which the individual fluid particles conserve their entropy over time. This hypothesis is substantially drastic in the context of solar-like oscillations, as the local thermal timescale actually has the same order of magnitude as the period of the modes, especially in the superadiabatic region. Therefore, this approximation needs to be dropped to obtain realistic expressions for the driving rate, damping rate, and modal surface effect associated with solar-like oscillations. This can be done from the start by adding a stochastic differential equation for the internal energy associated with the individual fluid particles in the Lagrangian stochastic model. The second is hypothesis H7, according to which the lifetime of the turbulent eddies is assumed to be constant in time – albeit not necessarily uniform – and is parameterised by the turbulent frequency, ω_t. This is an equally drastic approximation, as the range of timescales associated with the turbulent eddies is actually very large due to the very high Reynolds number characteristic of stellar turbulent convection at the top of the convective envelope. Like the quasi-adiabatic approximation, obtaining a realistic analytical prescription for the excitation rate, damping rate, and modal surface effect will require this assumption to be discarded. This can be done by treating ω_t as a physical property of each individual fluid particle separately. In doing so, it would also be necessary to add a stochastic differential equation for the turbulent frequency associated with a given fluid particle in the Lagrangian stochastic model from the start. Adding turbulent frequency and internal energy to the set of wave variables used to describe the solar-like oscillations is therefore an essential task in the development of this formalism, and will be the subject of a future paper in this series. While it does not constitute a conceptual challenge, and while many studies in the fluid dynamics community have already been devoted to the inclusion of turbulent frequency in turbulence modelling (e.g., Pope & Chen 1990; Pope 1991), the application to solar-like oscillations is likely to involve many more theoretical developments and yield much more complex relations: it may become necessary, at some point, to proceed towards a numerical resolution of these equations instead.

It will also be necessary, in the long run, to discard the assumption that the PDF of the turbulent flow is Gaussian. Indeed, the Lagrangian stochastic model from which we started in Paper I is originally constructed under the implicit assumption that the equivalent Fokker-Planck equation yields a multivariate Gaussian PDF. By contrast, as shown by 3D hydrodynamic simulations, stellar turbulent convection is characterised by a typical structure composed of upflows and downdrafts, each separately exhibiting Gaussian turbulence, such that the total flow has a bimodal distribution. Because of the asymmetry between the upflows and the downdrafts – the latter being colder and more turbulent than the former – the total distribution is not Gaussian. A possible solution that would allow the non-Gaussian nature of turbulence to be accounted for would be to use a two-flow model, where the distribution is determined by the mean and variance of two Gaussian distributions, or equivalently by the first four moments of the whole distribution. It is possible to implement two-flow prescriptions in Lagrangian stochastic models (see for instance Rodean 1996) by using different sets of stochastic differential equations to model the evolution of the fluid particles, depending on whether their vertical velocity is positive or negative and on whether their temperature is higher or lower than the local mean temperature. Concretely, this would influence the relation between the turbulent fluctuations – especially the turbulent velocity – and the perturbations ℒ^s and |θ⟩ to the wave equation. This would constitute an important step towards making the present approach more realistic, and its results quantitatively applicable to solar-like oscillators and thus suitable for comparison with observations.

As a final, important note, we recall that the results presented from Sect. 4 onwards, and in particular the interpretation of the various terms in the averaged simplified amplitude equations in terms of mode excitation, mode damping, and the surface effect, were obtained by neglecting mode coupling, that is, by reducing the coupled amplitude and phase equations (Eqs. (43) and (44)) to two equations on the amplitude and phase of a single mode (Eqs. (55) and (56)). While such coupling is deemed negligible in solar-type, main-sequence stars, where it is usually assumed that modes can be studied in isolation from one another, this assumption is in fact not warranted by any a priori physical argument. This is a most important caveat of the findings presented in this study; in the future, it will be necessary to go beyond the single-mode case, and instead consider the full, coupled system of equations, especially when it comes to applying the present formalism to the case of a more evolved star. A later paper in this series will be dedicated to this task.

---

Acknowledgments

We thank the anonymous referee for very helpful comments that helped improve the quality of this manuscript. Part of this work was conducted with the support of the Max Planck Institute for Solar System Research.

References

Appendix A: Derivation of the Fokker-Planck equation

In this appendix, we detail the derivation of the probability fluxes 𝒢_μ and ℋ_μ and the components of the diffusion matrix 𝒟_μν, ℰ_μν and ℱ_μν appearing in the joint-amplitude-phase Fokker-Planck equation (Eq. 30). In order to account for the finite memory time of G_μ(A_ν, Φ_ν, t) and H_μ(A_ν, Φ_ν, t) (defined as the right-hand sides of Eqs. 28 and (29), respectively), Stratonovich et al. (1965) showed that one can define effective probability fluxes and diffusion coefficients in the following way:

(A.1)

(A.2)

(A.3)

(A.4)

(A.5)

where ⟨.⟩ denotes an ensemble average, and the bilinear operator {. ;  .} is defined by

(A.6)

We can rewrite each of the stochastic processes c_1μ(t), c_2μν(t) and c_3μν(t) in the following polar form:

(A.7)

(A.8)

(A.9)

where the functions C_i(t) and ϕ_i(t) are both real. Then the right-hand sides of Eqs. 28 and (29) can be rewritten

(A.10)

and

(A.11)

Their respective derivatives with respect to A_ν and Φ_ν are

(A.12)

(A.13)

(A.14)

and

(A.15)

Plugging these into Eqs. A.1–A.5, it may seem at first glance that the expansion of the bilinear operators { ... ;  ... } involves a large number of terms, all of which can be put in the following general form

(A.16)

with various amplitudes C_α/β, angular frequencies Ω_α/β and phases ϕ_α/β (we note that sines can always be written as cosines by redefining their phase). Plugging this general form into Eq. A.6, we find

(A.17)

The first term oscillates in time t with an angular frequency Ω_β + Ω_α, while the second one oscillates with an angular frequency Ω_β − Ω_α. However, since we are interested in the long-term variations in the complex amplitude of the modes, which occur on timescales much larger than , all terms oscillating at frequencies comparable to any ω_λ can be averaged out. This means that Eq. A.17 only yields non-vanishing contributions if either Ω_α = Ω_β or Ω_α = −Ω_β. By redefining the phase inside the cosine, we can always restrict ourselves to the first case. Otherwise stated, this means that the only terms that can ‘interact’ through the operator { ... ;  ... } to yield a non-oscillating contribution are those that share the same frequency in the first place. This can happen either because the two frequencies are rigorously identical or because there is an accidental resonance of the type ∑_in_iω_i = 0, where the sum concerns two or more modes, and n_i are integers. In the following we refer to the first kind of contribution as ‘cross terms’ and to the second kind as ‘resonances’. We treat each kind separately in the following two sections.

A.1. Cross terms

The contributions from the cross terms in the curly bracket operators appearing on the right-hand sides of Eqs. A.1–A.5 are given by Eqs. A.19–A.24 below⁶. It can be seen that in each of the individual curly brackets therein, the frequencies in both cosines/sines are indeed rigorously identical.

(A.19)

(A.20)

(A.21)

(A.22)

(A.23)

(A.24)

Applying the same procedure underlined in Eq. A.17, Eqs. A.19–A.24 become

(A.25)

(A.26)

(A.27)

(A.28)

(A.29)

and

(A.30)

In turn, plugging these expressions into Eqs. A.1–A.5 yields the final expressions for the probability flux and diffusion matrix elements, as given in the main body of this paper (Eqs. 31–35).

A.2. Resonances

To illustrate the emergence of these resonance terms in the coefficients of the Fokker-Planck equation, we consider the curly bracket operator acting on the second term on the right-hand side of ∂G_μ/∂A_ν and the third term on the right-hand side of G_ν. These two terms oscillate at frequencies ω_μ − ω_ν and ω_ν − ω_λ, respectively. As we saw before, the resulting curly bracket only yields non-zero contributions to the long-term evolution of the mode amplitudes and phases if these two frequencies are either equal or opposite. Obviously, for λ = μ the two frequencies are the exact opposite of each other: from this stems one of the cross terms accounted for in the previous section (more specifically, the first term on the right-hand side of Eq. A.19). But there may also be another mode λ for which the two above frequencies coincide (i.e. ω_μ + ω_λ − 2ω_ν = 0). To be perfectly accurate, this resonance condition only needs to be verified to within the sum of the inverse lifetimes of all the modes involved (i.e. about 10 μHz close to the frequency of maximum spectral height ν_max for the Sun). Following the lead of Kumar & Goldreich (1989) for instance, we argue that given the high density of modes per unit frequency in the solar p-mode spectrum, the sums over frequencies can be treated as integrals, and the resonance condition enforced with a Dirac distribution. The resulting contribution then reads

(A.31)

It can be seen that the major difference between this contribution and those computed in the previous section is the additional factor exp^{j(2Φ_ν − Φ_λ − Φ_μ)} inside the integral (this factor being rigorously equal to unity in each of the cross terms above). Under the hypothesis that the phase differences between the modes are fairly independent from one pair of modes to another, it is reasonable to assume that all these resonance terms cancel each other out, thus yielding a total net contribution that can be neglected compared to the cross terms computed in the previous section.

Appendix B: Derivation of the coefficients α_i

The evolution equations on the mean mode energy and phase, given by Eqs. 71 and 76, only depend on the constant, complex values of the coefficients α₁ and α₃. We recall here, for clarity, that they are defined as the complex autocorrelation spectra of the stochastic processes c₁(t) and c₃(t) evaluated at angular frequencies ω and 2ω, respectively, where ω is the angular eigenfrequency of the mode, such that (see Eqs. 51 and 54)

(B.1)

(B.2)

In turn, these stochastic processes are given by (see Eqs. 25 and 26)

(B.3)

(B.4)

where we recall that the operator ℑ is given by Eq. 6, Ψ_ξ and Ψ_u are the normalised displacement and velocity eigenfunctions respectively, and we have (see Eqs. 4 and 5)

(B.5)

and

(B.6)

and the perturbations of the Reynolds-stress tensor and mean shear tensor are given by (see Eqs. 7 and 8)

(B.7)

and

(B.8)

The apparition of the factor ω⁻¹ in conjunction with every occurrence of Ψ_ξ stems from the fact that the latter is defined in terms of ωξ_osc rather than simply ξ_osc.

B.1. Contribution of the turbulent displacement field

In Paper I we used a certain number of approximations to derive our formalism. One of these approximations – which we labelled (H3) – consisted in adopting the anelastic approximation for turbulence, in the sense that we considered ρ_t ≪ ρ₀, where ρ_t is the turbulent fluctuation of density, and ρ₀ is the equilibrium density. Using the continuity equation, this amounts to neglecting the quantity ∇ ⋅ (ρ₀ξ_t). As we will now see, this allows us to discard the first term on the right-hand side of Eq. B.3 and the second term on the right-hand side of Eq. B.4. Considering the former and performing an integration by part, we can put it in the following form:

(B.9)

The surface integral (i.e. the first term on the right-hand side) systematically vanishes, on account of the product ρ₀Ψ_ξ vanishing on the surface of the star. Furthermore, the typical length scale over which ξ_t varies is much smaller than the wavelength of the mode. This allows us to pull Ψ_ξ, i out of the gradient in the last term, and we recognise ∇ ⋅ (ρ₀ξ_t), which we neglected on account of our hypothesis (H3). The same procedure can be applied to the second term on the right-hand side of Eq. B.4. As a result, in Eqs. B.3 and B.4, the only source of stochasticity comes from the turbulent velocity field u_t and the Wiener process η, with no contribution from the turbulent displacement field ξ_t.

B.2. Derivation of α₁

Plugging Eq. B.6 into Eq. B.3, neglecting the first term on the right-hand side (see Section B.1 for a justification), and performing an integration by parts, the quantity c₁(t) can be rewritten in the following way:

(B.10)

where the surface contribution of the integration by part vanishes, because it involves ρ₀Ψ_u at the outer boundary of the star, where the oscillation is evanescent. This allows for a simple physical interpretation of this integral: in broad strokes, the quantity ∂_j(ρ₀Ψ_u, i)/ρ₀ corresponds to the local compressibility of the mode, while the quantity [ρ₀u_t, iu_t, j]′ corresponds to the instantaneous fluctuating turbulent pressure. The product between the two can therefore be interpreted as a local instantaneous ‘pdV’ mechanical work exerted by the turbulent pressure on the mode; integrating over the entire stellar volume yields the total instantaneous work c₁; and forming the autocorrelation function of c₁ yields the effective work exerted on the mode over a large number of cycles, which corresponds to the excitation rate of the mode.

We remark that only the fluctuation of the turbulent pressure around its equilibrium state appears in this integral. Expanding it into the difference between the total turbulent pressure and its equilibrium value leads to four different integrals. However, it is readily seen that only one of them – the one not involving the equilibrium turbulent pressure at all – will have a non-zero ω-component in its Fourier transform in time, and therefore will contribute to α₁. As such, [u_t, iu_t, j]′ can be replaced by u_t, iu_t, j. Then, forming the autocorrelation function of c₁, we obtain

(B.11)

In the scope of the Jeffreys-Wentzel-Kramers-Brillouin (JWKB) approximation, the velocity eigenfunction can be locally approximated by the following expression

(B.12)

where Ψ_u, 0(x) is the slowly varying amplitude in space of the velocity eigenfunction, and k(x) is the space-dependent wave vector of the mode. Both Ψ_u, 0 and k are slowly varying functions of space (meaning that they vary on length scales much larger than |k|⁻¹). As such, we have

(B.13)

so

(B.14)

Then, we implement the following change of variables:

(B.15)

(B.16)

where X and δx represent the slow and fast space variables, respectively. The implementation of this change of variables in Eq. B.14 allows us to completely decouple the slowly varying contributions of Ψ_u, 0 and k on the one hand, and the rapidly varying contributions of exp(jk ⋅ x) and the turbulent velocity two-point correlation products on the other. We obtain⁷

(B.17)

Taking advantage of the scale decoupling, we can separate the integrals over X and δx,

(B.18)

Finally, using Eq. B.1, and taking advantage of the definition of the operator ℑ (see Eq. 6), the autocorrelation spectrum α₁ can thus be expressed as

(B.19)

where is the fourth-order correlation spectrum of the turbulent velocity u_t, defined as

(B.20)

and the operator 𝒞_ω, k is defined by

(B.21)

and

(B.22)

Naturally, also depends on the slow variable X and time t, even though they do not appear explicitly in Eq. B.19.

B.3. Derivation of α₃

Plugging Eq. B.5 into Eq. B.4 (from which the second term on the right-hand side was discarded, as per the argument developed in Section B.1), we obtain

(B.23)

where the perturbation of the Reynolds-stress tensor and mean shear tensor are given by Eqs. B.7 and B.8, respectively.

We consider the first four terms on the right-hand side of Eq. B.23. They can be rearranged with integrations by part to yield

(B.24)

(B.25)

The last term can be further simplified by permuting the integral over x and the integral defining the operator ℑ (see Eq. 6), which yields

(B.26)

where we have used the isotropy of the kernel function K to write K^x(y)=K^y(x). It can be seen that the integral over x corresponds to the kernel estimate at point y of the quantity Ψ_u, iG_ij, 0/(ωρ₀), which only involves quantities that are not stochastic. As such, this kernel estimation equals the actual value of this quantity at y, and c_3d(t) reduces to

(B.27)

As for the fifth, sixth, and seventh terms on the right-hand side of Eq. B.23, we can expand them in the following way. For the sake of clarity, we define the following symbol:

(B.28)

Then the fifth term on the right-hand side of Eq. B.23 can be split into

(B.29)

(B.30)

(B.31)

and

(B.32)

the sixth term can be split into

(B.33)

(B.34)

(B.35)

and

(B.36)

and the seventh term can be split into

(B.37)

(B.38)

(B.39)

Finally, in the last term on the right-hand side of Eq. B.23, the quantity under the integral is proportional to the stochastic process η_i(x, t), which, by construction, is δ-correlated in both space and time. In particular, its correlation length scale is infinitesimally small compared to that of the turbulent velocity u_t(x, t). But as we saw in Section B.2, it is precisely the spatial coherence of the stochastic perturbations to the wave equation that explains its ability to impact the complex amplitude of the modes. As such, this part will not actually contribute to the final expression of α₃, or to the stochastic amplitude equations in any way, and will be discarded in the following.

Formally, c₃(t) can be written as a sum of contributions that are first-, second-, or third-order in terms of the turbulent velocity, u_t,

(B.40)

where

(B.41)

(B.42)

and

(B.43)

Forming the autocorrelation product of c₃(t) from Eq. B.40, we see that the expansion involves correlation products of the turbulent velocity field u_t of various orders, ranging from 2 to 6. In the following, we make the assumptions, often used in the context of Gaussian turbulence, that 1) these moments can be cut at fourth order, and 2) the contribution of the third-order moment can be neglected compared to that of the second- or fourth-order. With this approximation in mind, we then proceed to adopt the JWKB approximation in the same form as in Section B.2 (see Eq. B.12), which yields

(B.44)

where the subscript τ means that the turbulent velocity field is evaluated at time t + τ, and we have introduced

(B.45)

(B.46)

(B.47)

(B.48)

(B.49)

Equation B.44 can be drastically simplified by remarking that any integral involving the product of a function f with the kernel function K correspond, by construction, to the ensemble average of said function f (see Paper I for more details). Similarly, if the integral is weighted by the gradient of the kernel function, then it corresponds to the ensemble average of the gradient of f (with a minus sign). But every quantity appearing in Eq. B.44 is either already an ensemble average or an equilibrium quantity, or else the normalised velocity eigenfunction Ψ_u. None of these are stochastic quantities, which means they are equal to their own ensemble average. This allows us to perform the kind of simplification illustrated by Eq. 50 of Paper I, and we eventually find

(B.50)

where we defined F⁽¹⁾ ≡ F^(1a) + F^(1b).

Finally, plugging this into Eq. B.2, we find the following expression:

(B.51)

where we recall that the wave vector k and the angular frequency, ω, are those of the mode under consideration, and we have defined the following spectra on the same template as Eq. B.21:

(B.52)

(B.53)

(B.54)

(B.55)

(B.56)

We note, as we did in Section B.2, that F⁽¹⁾, F⁽²⁾, F^(3a), and F^(3b) are also functions of X, and that , , , , and also depend on both X and t, even though these dependences do not appear explicitly in Eq. B.51.

B.4. Mode normalisation and final form of α_i

The last remaining modification to the explicit expressions of α₁ and α₃ concerns the normalisation condition on the mode |Ψ⟩. Indeed, we need to relate Ψ_u, 0(x) to the actual modal velocity fluctuation u_osc(x), as it can be obtained through an oscillation code for instance. By construction of the ket |Ψ⟩, we have this very simple proportionality relation:

(B.57)

where the proportionality factor is given by the condition that |Ψ⟩ must be normalised to unity, such that

(B.58)

Plugging Eqs. 11 and 18 into Eq. B.58, this becomes

(B.59)

and since Ψ_u = jΨ_ξ,

(B.60)

Finally, plugging Eq. B.57, we find

(B.61)

which we recognise as the inertia of the mode. We then find the relation between Ψ_u, 0 and u_osc by plugging Eq. B.12 into Eq. B.57. In turn, plugging Eq. B.57 into Eqs. B.19 and B.51, we find the final expressions for α₁ and α₃ reproduced in the main body of the paper (Eqs. 65 and 66).

All Figures