© J. Philidet et al. 2021

1. Introduction

Solar-like oscillations are coupled with turbulent convection in a complex manner, especially in the highly turbulent subsurface layers of the star (see Samadi et al. 2015; Houdek & Dupret 2015, for a review). This coupling impacts the behaviour of the modes in several major ways. One of the most prominent effects concerns mode frequencies, and explains in a large part the systematic discrepancy between the theoretical and observed p-mode frequencies (Dziembowski et al. 1988; Christensen-Dalsgaard et al. 1996; Rosenthal et al. 1999). The variety of physical processes responsible for the impact of turbulent convection on p-mode frequencies is collectively referred to as ‘surface effects’. These surface effects constitute a major obstacle preventing us from using the full potential of modal frequencies for an accurate probing of stellar interiors or for a precise inference of stellar global parameters.

Many efforts have thus been devoted to the correction of surface effects, either from theoretical modelling (e.g., Gabriel et al. 1975; Balmforth 1992b; Houdek 1996; Rosenthal et al. 1999; Grigahcène et al. 2005) or through empirical formulae (e.g., Kjeldsen et al. 2008; Christensen-Dalsgaard 2012; Ball & Gizon 2014; Sonoi et al. 2015). Some aspects, however, are very complicated to model, and existing models make use of assumptions that can barely be justified, if at all. For instance, turbulent pressure modulations are usually described in the Gas-Γ₁ (GGM) or reduced-Γ₁ (RGM) approximations (Rosenthal et al. 1999), which amounts to neglecting the effects of turbulent dissipation and buoyancy on the mode (Belkacem et al. 2021). Another problem is the use of time-dependent mixing-length formalisms (Unno 1967; Gough 1977) to account for modal surface effects (e.g., Gabriel et al. 1975; Houdek 1996; Grigahcène et al. 2005; Sonoi et al. 2017; Houdek et al. 2017, 2019). While useful for the bulk of the convective region, the mixing-length hypothesis is no longer valid in the superadiabatic region just beneath the surface of the star, as shown by 3D hydrodynamic simulations of stellar atmospheres (see Nordlund et al. 2009, for a review). Finally, such formalisms require that the oscillations be separated from the convective motions, thus yielding separate equations. This is done either by assuming a cut-off in wavelength space, with oscillations having much shorter wavelengths than turbulent convection (Grigahcène et al. 2005), or by using 3D hydrodynamic simulations and separating the oscillations from convection though horizontal averaging (Nordlund & Stein 2001). The necessity to separate the equations for oscillations and convection is fundamentally problematic as there is no rigorous way to disentangle the two components, mainly because, in solar-type stars, they have the same characteristic lengths and timescales (Samadi et al. 2015). This is even truer if one wishes to model their mutual coupling.

Turbulent convection also has a crucial impact on the energetic aspects of solar-like oscillations. Solar-like p-modes are stochastically excited and damped by turbulent convection at the top of the convective zone. As such, understanding the energetic processes pertaining to the oscillations leads to better constraints on the highly turbulent layers located beneath the surface of these stars. Many theoretical modelling efforts were deployed on the subject of mode 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), as well as mode damping (e.g., Goldreich & Kumar 1991; Balmforth 1992a; Grigahcène et al. 2005; Dupret et al. 2006; Belkacem et al. 2012). The fact that these energetic processes take place in the superadiabatic region, however, makes any predictive model extremely complicated to design as this requires a time-dependent non-adiabatic turbulent convection model able to include the oscillations. Subsequently, modelling attempts have focused on the use of mixing-length formalisms to account for mode damping. This approach, however, presents the considerable disadvantage of reducing the turbulent cascade to a single length scale, and is therefore unable to correctly account for the contribution of turbulent dissipation or turbulent pressure to mode damping. Alternative approaches have been followed in an attempt to go beyond the mixing-length hypothesis, either through a Reynolds stress model (Xiong et al. 2000) or through the use of 3D hydrodynamic simulations of stellar atmospheres to directly measure mode linewidths (Belkacem et al. 2019; Zhou et al. 2020).

These traditional approaches therefore show some fundamental limitations, which prevent them from being able to fully describe the interaction between turbulent convection and oscillations, whether it be to explain the surface effects on mode frequency or the energetic aspects of global modes of oscillation regarding their driving and damping physical processes. Among these limitations, we can include the following.

These approaches require that the turbulent convection and the oscillations be separated into two distinct sets of equations from the start. This is usually justified either by a separation of spatial scales or timescales, or else by performing some averaging process designed to separate an average component from a fluctuating component. The necessity of artificially separating these two intertwined phenomena from the outset is problematic when it comes to modelling their coupling.

Most of these approaches are based on a time-dependent mixing-length formalism, which oversimplifies the behaviour of convection in the superadiabatic region. In addition to poorly describing the structure of the convective motions close to the surface, it reduces the turbulent cascade to a single characteristic length scale, thus only offering a crude understanding of turbulent dissipation, a phenomenon deemed crucial to turbulence–oscillation coupling.

In these approaches, designing a closure relation for the model equations is a complicated process. In particular, it is very difficult to properly relate the chosen closure to the underlying physical assumptions. This is illustrated, for instance, by the wealth of free parameters that need to be adjusted in approaches based on mixing length theory (MLT) or in Reynolds stress approaches where higher-order moments need to be closed at the mean flow level. Approaches based on 3D simulations are not spared, as is illustrated by the need to rely on assumptions like the GGM or RGM, which are not clearly physically grounded.

The multiple free parameters needed in these approaches, and the fact that they are not easily constrained physically, also presents the distinct disadvantage of robbing these models from their predictive power. This becomes problematic when, for example, mode damping rates are used in scaling relations for seismic diagnosis purposes (Houdek et al. 1999; Chaplin et al. 2009; Baudin et al. 2011; Belkacem et al. 2012). The exponent in these scaling relations is difficult to determine, and varies substantially across the Hertzsprung–Russell diagram. Being able to predict the damping rates of stars with different global parameters would go a long way towards a more effective use of this quantity in such scaling relations.

Model parameters for the surface effects on the one hand, and mode damping rates on the other are usually constrained by completely separate adjustment procedures. This is also problematic, as these two quantities are closely related, and are actually just two sides of the same phenomenon: the real and imaginary part of the turbulence-induced shift in the complex eigenfrequency of the modes.

These limitations form substantial hurdles towards a correct turbulence–oscillation coupling model, and circumventing them requires going beyond the methods presented above. Therefore, this series of papers follows a completely different approach. More precisely, the fundamental motivation behind this work is to provide a method that (1) does not initially rely on a separation between convection equations and oscillation equations, but instead encompasses both components at the same time, and therefore naturally contains their coupling; (2) avoids the reduction of length scales in the problem to a unique scale, but instead accounts for the full description of the turbulent cascade; (3) simultaneously describes all effects of turbulent convection on mode properties, including the surface effects and the energetic aspects pertaining to mode driving and damping, in a single consistent framework; and (4) includes the properties of turbulence in such a way that they can be easily related to the observed properties of the modes.

In this paper we therefore build a formalism for the modelling of turbulence–oscillation coupling, which is based on probability density function (PDF) models of turbulence (e.g., O’Brien 1980; Pope 1985, 1991, 1994; Pope & Chen 1990; Van Slooten & Jayesh 1998). The quantity we model is the PDF associated with the random flow variables, whose evolution follows a transport equation that takes the form of a Fokker-Planck equation (Gardiner 1994). Because the Fokker-Planck equation is impractical to handle both analytically and numerically, the PDF is usually represented by a set of fluid particles constituting the flow. The properties of the particles evolve according to stochastic differential equations, and are then used to reconstruct any given statistics of the flow. This is at the heart of Lagrangian stochastic models of turbulence, which have been used extensively by the fluid dynamics community, first for incompressible flows (e.g., Pope 1981; Anand et al. 1989; Haworth & El Tahry 1991; Roekaerts 1991), and then for compressible flows as well (e.g., Hsu et al. 1994; Delarue & Pope 1997; Welton & Pope 1997; Welton 1998; Das & Durbin 2005; Bakosi & Ristorcelli 2011). In this series of papers we present a general way of using such a Lagrangian stochastic model of turbulence to derive a linear stochastic wave equation applicable to the stellar context. The wave equation is designed to govern the physics of the modes, while simultaneously and consistently encompassing the impact of turbulent convection thereon. This first paper describes how the linear stochastic wave equation is obtained. A subsequent paper will present how this wave equation can be used to simultaneously model the turbulence-induced surface effects, as well as the stochastic driving and damping of the modes by turbulent convection.

This paper is structured as follows. In Sect. 2 we introduce the stochastic model of turbulence that we will use throughout this study in terms of Lagrangian variables. Then we carry out a variable transformation to obtain stochastic equations on Eulerian quantities, which are more suitable for stellar oscillation analysis. In Sect. 3 we linearise the Eulerian stochastic equations to obtain a linear stochastic wave equation, and then discuss how it relates to other more familiar forms of the wave equation found in the literature and obtained through more traditional methods. Finally, in Sect. 4 we return to the various simplifications and approximations adopted in the present derivation, and what they entail as regards the resulting properties of turbulence-oscillation coupling. Conclusions are drawn in Sect. 5.

2. Stochastic model of turbulence

In MLT formalisms the modelled quantities pertain to the mean flow (e.g., mean density, velocity, entropy), and the second-order moments appearing in the mean equations must be expressed in terms of the mean flow in order to close the system. In Reynolds stress formalisms the closure at second-order level is replaced by equations on second-order moments, where third-order moments must be similarly closed. These moments are all defined as ensemble averages of stochastic processes such as flow velocity and entropy. The core idea behind PDF models is to replace these numerous equations on various statistical moments of turbulent quantities by a single equation on the PDF of these quantities, in the form of a Fokker-Planck equation. These models present several advantages, which are of special interest given the issues raised in the previous section. By nature, the modelled PDF contains all the required statistical information on the flow, which includes both the turbulent convection and the oscillating modes. As such, this type of model is perfectly suited for the study of turbulence–oscillation coupling. In addition, all the usual quantities can be obtained from the PDF.

However, the direct modelling of the PDF, using its time evolution equation, can quickly become very cumbersome. The reason is that the PDF is a function not only of space and time, but also of each of the turbulent quantities used to represent the flow (starting with the three components of the velocity and the entropy). This makes the PDF equation computationally heavy to integrate, and quite impractical to handle analytically. This is why PDF models are often implemented in a Lagrangian particle framework, where the flow is no longer represented by a set of Eulerian, grid-based fluid quantities, but rather by a set of individual fluid particles whose properties (including their position) are tracked over time. Using Monte Carlo methods, the flow PDF can be reconstructed directly from the set of particles, so that the set contains the exact same statistical information as the PDF itself. In order to represent the turbulent nature of the flow, and to model the PDF accurately, particle properties must evolve according to stochastic differential equations rather than ordinary ones. Therefore, PDF models of turbulence go hand in hand with the implementation of Lagrangian stochastic methods, primarily because it makes their numerical integration much easier and more tractable. In this section we introduce the Lagrangian stochastic model, and we present how it can be rearranged to yield stochastic differential equations for Eulerian quantities instead.

We note that this paper aims to show that the method we present is relevant to the study of turbulence–oscillation coupling, and therefore serves as a proof of concept for this approach. As such, we do not claim to use the most realistic turbulence model possible, but rather we wish to limit the level of complexity so that the basics of this method may be understood in the most efficient way. We leave the use of a more realistic turbulence model for a later paper.

2.1. Lagrangian description: The generalised Langevin model

We consider the simplified case of an adiabatic¹ flow, in the sense that we do not include an energy equation in the system, and instead adopt a polytropic relation between the pressure and the density of the gas. In terms of Eulerian transport equations that would mean only considering the density, the mean velocity, and the Reynolds stress tensor as relevant fluid quantities, with the mean pressure being given, for instance, by the ideal gas law. In the framework of a Lagrangian stochastic model, however, that means that the only fluid particle properties whose evolution we need to put into equations are their position and velocity.

The equation for the particle position is derived by stating that it must evolve according to its own velocity. It reads

$$\begin{array}{c}\hfill \mathrm{d}{x}_i^{\star }=u_i^{\star }\mathrm{d}t,\end{array}$$(1)

where x^⋆ and u^⋆ are the position and velocity of the fluid particle, which only depend on the time variable (as well as the initial state). In general, in the following the notation ^⋆ will denote a stochastic variable. In order to account for the turbulent nature of the flow, the equation on velocity must take the form of a stochastic differential equation (SDE), instead of an ordinary one. In its most general form, an SDE takes the form (Gardiner 1994, Chap. 3)

$$\begin{array}{c}\hfill \mathrm{d}{u}_i^{\star }=a_i({\mathbf{x}}^{\star },{\mathbf{u}}^{\star },t)\mathrm{d}t+b_{\mathit{ij}}({\mathbf{x}}^{\star },{\mathbf{u}}^{\star },t)\mathrm{d}{W}_j,\end{array}$$(2)

where we use the Einstein convention on repeated indices, a_i and b_ij are functions of the particle properties (and time), and W(t) is an isotropic Wiener process. The last is a stochastic process (i.e. a random variable whose statistical properties depend on time) whose PDF at any given time t is Gaussian, and which verifies

$$\begin{array}{cc}& \overline{\mathbf{W}(\mathbf{t})}=\mathbf{0},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill & \overline{W_i(t^{\prime })W_j(t)}=(t^{\prime }-t){\delta }_{\mathit{ij}},\hfill \end{array}$$(4)

where δ_ij is the Kronecker symbol and the notation $\overline{.}$ refers to an ensemble average. We note that this is not a simplifying assumption regarding the stochastic part of the SDE, but rather a very general property, which is necessary for the resulting particle trajectory in phase-space to be continuous in time (Gardiner 1994). In terms of dimension, the drift vector a_i is an acceleration, while W is the square root of a time, and the diffusion tensor b_ij is a velocity divided by the square root of a time.

On the right-hand side of Eq. (2) the first term corresponds to the deterministic part of the force exerted on the fluid particle, while the randomness of the equation is only brought about by the second term. Physically, the stochastic part of Eq. (2) stems from the fluctuating components of both the pressure and viscous stress forces, which in turn are brought about by the underlying highly fluctuating turbulent velocity field. An illuminating analogy to consider is Brownian motion, which can also be described by means of Eq. (2), and where the stochastic part describes the random collision undergone by the colloidal particle from the water molecules. In the vocabulary of stochastic processes the function a_i(x^⋆, u^⋆, t) is the ith component of the drift vector, while b_ij(x^⋆, u^⋆, t) is the i, j-component of the diffusion tensor. In order to close the system, an explicit expression is needed for these two coefficients.

The specification of the drift and diffusion terms in Eq. (2) is the subject of an abundant amount of literature on turbulence modelling (see Heinz & Buckingham 2004, for a review). It has long been recognised that, in order to be consistent with the Kolmogorov hypotheses, both original (Kolmogorov 1941) and refined (Kolmogorov 1962), the diffusion coefficient has to take the form (Obukhov 1959)

$$\begin{array}{c}\hfill b_{\mathit{ij}}({\mathbf{x}}^{\star },t)=\sqrt{C_0ϵ({\mathbf{x}}^{\star },t)}{\delta }_{\mathit{ij}},\end{array}$$(5)

where C₀ is a dimensionless constant and ϵ is the local dissipation rate of turbulent kinetic energy into heat. This is especially verified in the high Reynolds number limit (which is relevant in the stellar context), where C₀ then actually corresponds to the Kolmogorov constant. This constant is not universal per se; however, it tends asymptotically to a universal value for very high Reynolds numbers, in which case its value is fairly well constrained. An accepted experimental value is C₀ = 2.1 (Haworth & Pope 1986).

For the drift term we adopt the general expression given by the generalised Langevin model (Pope 1983)

$$\begin{array}{c}\hfill a_i({\mathbf{x}}^{\star },{\mathbf{u}}^{\star },t)=-{\displaystyle \frac{1}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}+\overline{g_i}+G_{\mathit{ij}}(u_j^{\star }-\stackrel{\sim }{u_j}),\end{array}$$(6)

where $\overline{\rho }$, $\overline{p}$, and $\overline{\mathbf{g}}$ are the Reynolds average of the fluid density, gas pressure, and gravitational acceleration, respectively; G_ij is a second-order tensor that has the dimension of an inverse time, to which we refer as the drift tensor; and $\stackrel{\sim }{\mathbf{u}}$ is the Favre average of the fluid velocity, with the mass-average (or Favre average) of any quantity ϕ being defined as

$$\begin{array}{c}\hfill \stackrel{\sim }{\varphi }\equiv {\displaystyle \frac{\overline{\rho \varphi }}{\overline{\rho }}}.\end{array}$$(7)

All these Reynolds or Favre averages are local and instantaneous quantities, and therefore depend on both time and space. In Eq. (6) they are evaluated at time t and at the position x^⋆ of the particle.

The various terms in Eq. (6) can be interpreted in the following way. The first two terms are the mean pressure gradient and the gravitational force exerted on the particle, and correspond to the mean force in the momentum equation, the only ones that remain in the absence of turbulence; we note that rotation and magnetic fields are not accounted for in this model. On the other hand, the last term ensures that, were the turbulent sources to disappear, the particle velocity would decay towards the local mean velocity, thus ensuring that the Reynolds stresses are dissipated. More precisely, the drift tensor can be thought of as the rate at which the various Reynolds stresses decay towards zero. In this paper we need not specify the form of the drift tensor, only to say that in the standard approach it is written as a function of the Reynolds stresses, the mean velocity gradients, and the turbulent dissipation (Haworth & Pope 1986)

$$\begin{array}{c}\hfill G_{\mathit{ij}}=f(\widetilde{u_i^{″}{u}_j^{″}},{\partial }_i\stackrel{\sim }{u_j},ϵ),\end{array}$$(8)

where ${\mathbf{u}}^{″}\equiv \mathbf{u}-\stackrel{\sim }{\mathbf{u}}$ denotes the fluctuation of the turbulent velocity around its local Favre average. In particular, G_ij only depends on the mean fields and not on the particle properties themselves.

Putting together Eqs. (1), (2), (5), and (6), we obtain

$$\begin{array}{cc}& \mathrm{d}{x}_i^{\star }=u_i^{\star }\mathrm{d}t,\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill & \mathrm{d}{u}_i^{\star }=[-{\displaystyle \frac{1}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}+g_i+G_{\mathit{ij}}(u_j^{\star }-\stackrel{\sim }{u_j})]\mathrm{d}t+\sqrt{C_0ϵ}\mathrm{d}{W}_i.\hfill \end{array}$$(10)

The mean fields $\overline{\rho }$, $\overline{p}$, $\stackrel{\sim }{\mathbf{u}}$, G_ij, and ϵ still need to be closed; we return to this matter in Sect. 2.3.

The stochastic Eqs. (9) and (10) contain more information than the corresponding average equations on the mean velocity and Reynolds stress tensor, the same way the PDF of a distribution carries more statistical information than its first few moments. We do not make use of these corresponding mean equations in the following; nevertheless, we provide them explicitly in Appendix A, to which the reader can refer for a better grasp on the origin of the SDE used in this study.

2.2. From Lagrangian to Eulerian variables

2.2.1. The Lagrangian mean trajectory formalism

By construction, the turbulence model given by Eqs. (9) and (10) is a Lagrangian model as it pertains to the properties of fluid particles followed along their trajectories. By contrast, we would like to obtain equations on stochastic variables pertaining to the stochastic properties of the flow at a fixed point. This would allow us to ultimately obtain a wave equation where the wave variables can be easily related to the known properties of the modes, something for which a purely Lagrangian² description is extremely impractical.

A very general approach to this transcription from Lagrangian to Eulerian variables is the Lagrangian mean trajectories formalism (Soward 1972; Andrews & McIntyre 1978). In the following, we give the general ideas and the main steps of the derivation; more detailed calculations are provided in Appendix B, to which we will refer each time an important step is reached. Let us consider a fluid particle whose time-independent average position is denoted by x. Its instantaneous position at time t is written as an explicit function of x and t

$$\begin{array}{c}\hfill \mathbf{X}(\mathbf{x},t)=\mathbf{x}+\mathit{\xi }(\mathbf{x},t),\end{array}$$(11)

where ξ is the particle displacement around its mean position³, the mean position being interpreted as an Eulerian variable.

For any given Eulerian quantity ϕ, we define its Lagrangian counterpart as

$$\begin{array}{c}\hfill {\varphi }_L(\mathbf{x},t)\equiv \varphi (\mathbf{X}(\mathbf{x},t),t).\end{array}$$(12)

In particular, we denote by u_L the velocity field evaluated at X, in other words the Lagrangian velocity, and by u_i, L the ith component of this velocity. Similarly, for any Eulerian averaging process ⟨.⟩, we define the corresponding Lagrangian mean ⟨.⟩_L as

$$\begin{array}{c}\hfill {⟨\varphi ⟩}_L\equiv ⟨\varphi (\mathbf{X}(\mathbf{x},t),t)⟩.\end{array}$$(13)

For the time being, we do not yet specify the averaging process ⟨.⟩ as this formalism is very general and can be used regardless of how the means are defined. It is important to note here that the mean position x of the particle is defined in terms of this yet-to-be-determined averaging process. In the following we simply refer to ⟨.⟩ as the ‘Eulerian mean’, but let us keep in mind that it does not necessarily correspond to an ensemble average.

With the above notations and definitions, the following identity can be derived (Andrews & McIntyre 1978)

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{D\varphi }{\mathit{Dt}}}\right)}_L={⟨D⟩}_L\left({\varphi }_L\right),\end{array}$$(14)

where D/Dt ≡ ∂_t + u ⋅ ∇ denotes the particle time derivative, and the operator ⟨D⟩_L is defined by

$$\begin{array}{c}\hfill {⟨D⟩}_L\equiv {\partial }_t+{⟨\mathbf{u}⟩}_L·\mathbf{\nabla },\end{array}$$(15)

and ⟨u⟩_L is the Lagrangian mean of the flow velocity. For a detailed derivation of this identity, we refer the reader to Appendix B.1. Because the Lagrangian and Eulerian frames are in motion with respect to one another, the index _L does not commute with either the space or time derivative. For instance, ∂(ϕ_L)/∂t corresponds to the time derivative of the quantity ϕ as seen from the point of view of a fluid parcel (i.e. in the Lagrangian frame), while (∂ϕ/∂t)_L is the time derivative of the quantity ϕ as seen from an Eulerian point of view, and then evaluated at a given Lagrangian coordinate, after the fact. Essentially, Eq. (14) describes how the material time derivative commutes with the passage from Lagrangian to Eulerian variables, and will therefore be useful for transcribing our Lagrangian model into a Eulerian one.

Applying Eq. (14) on position and velocity respectively yields

$$\begin{array}{cc}& {\displaystyle \frac{\partial {\xi }_i}{\partial t}}=u_{i,L}-{⟨u_i⟩}_L-{⟨u_j⟩}_L{\partial }_j{\xi }_i,\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial \left(u_{i,L}\right)}{\partial t}}={\left({\displaystyle \frac{\partial u_i}{\partial t}}\right)}_L+[u_{j,L}{\delta }_{\mathit{jk}}-{⟨u_j⟩}_L{\delta }_{\mathit{jk}}-{⟨u_j⟩}_L{\partial }_j{\xi }_k]{\left({\displaystyle \frac{\partial u_i}{\partial x_k}}\right)}_L.\hfill \end{array}$$(17)

The derivation of these two equations is given in detail in Appendix B.2. We note that, for the moment, the displacement ξ and velocity u are flow variables, which is why they are not denoted with a ^⋆. We now relate these flow quantities to the position x^⋆ and velocity u^⋆ of the fluid particles. Since ξ and u_L correspond to the displacement and velocity of the particle whose mean position is x, then for any fixed x we have

$$\begin{array}{cc}& {\mathbf{x}}^{\star }(t)=\mathbf{x}+\mathit{\xi }(\mathbf{x},t),\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill & {\mathbf{u}}^{\star }(t)={\mathbf{u}}_L(\mathbf{x},t),\hfill \end{array}$$(19)

so that

$$\begin{array}{cc}& {\displaystyle \frac{\mathrm{d}{\mathbf{x}}^{\star }}{\mathrm{d}t}}={\displaystyle \frac{\partial \mathit{\xi }}{\partial t}},\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\mathrm{d}{\mathbf{u}}^{\star }}{\mathrm{d}t}}={\displaystyle \frac{\partial \left({\mathbf{u}}_L\right)}{\partial t}}.\hfill \end{array}$$(21)

Putting together Eqs. (10), (17), and (21) we obtain

$$\begin{array}{cc}& {\left({\displaystyle \frac{\partial u_i}{\partial t}}\right)}_L+[u_{j,L}-{⟨u_j⟩}_L-{⟨u_k⟩}_L{\left({\displaystyle \frac{\partial {\xi }_j}{\partial x_k}}\right)}_L]{\left({\displaystyle \frac{\partial u_i}{\partial x_j}}\right)}_L\hfill \\ \hfill & =-{\displaystyle \frac{1}{{\overline{\rho }}_L}}{\left({\displaystyle \frac{\partial \overline{p}}{\partial x_i}}\right)}_L+g_{i,L}+G_{ij,L}(u_{j,L}-{\stackrel{\sim }{u_j}}_{,L})+\sqrt{C_0{ϵ}_L}{\eta }_{i,L}.\hfill \end{array}$$(22)

By construction, η(x, t)≡dW/dt is a multi-variate Gaussian process whose values at two distinct locations are completely uncorrelated, and which verifies⁴

$$\begin{array}{cc}& \overline{\mathit{\eta }(\mathbf{x},t)}=\mathbf{0},\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill & \overline{{\eta }_i(\mathbf{x},t^{\prime }){\eta }_j(\mathbf{x},t)}=\delta (t^{\prime }-t){\delta }_{\mathit{ij}},\hfill \end{array}$$(24)

where δ(t) is the Dirac distribution.

Finally, we note that in this form, all the quantities present in Eq. (22) are evaluated at the instantaneous position X. Insofar as the transformation x ↦ X is invertible (i.e. for any x and t there exists y such that X(y, t) = x), we can drop the notation _L from this equation, except in ⟨u⟩_L, thus yielding a stochastic equation for the evolution of the Eulerian velocity u.

2.2.2. Specification of the averaging process

For the moment, we still have not specified the nature of the mean which defines both the mean particle position x and the Lagrangian mean velocity ⟨u⟩_L. The specification of this averaging process is crucial, and the fact that this formalism applies to any averaging process is of the utmost importance. Usually, the Lagrangian mean trajectory formalism is used to transform the exact Eulerian hydrodynamics equations into equations on carefully defined Lagrangian means, which happen to be much more suitable to the study of hydrodynamic waves (Andrews & McIntyre 1978). In that context it is customary to consider that ⟨.⟩ actually denotes an ensemble average, so that the mean values contain the oscillations as well as the background equilibrium, whereas the fluctuations contain the turbulent fields.

In the present context, however, this is not the picture after which we are. Instead, we want the fluctuating part to contain the waves in addition to the turbulent fluctuations, whereas the means should only contain the background equilibrium. Only then can we obtain a wave equation directly containing turbulence-induced fluctuations, and therefore the turbulence–oscillation coupling. As such, we will define ⟨.⟩ as a time average over timescales that are very long compared to the typical turbulent timescale, the period of the oscillations, and their lifetime. This ensures that the mean values only contain the time-independent equilibrium, and the fluctuating part does indeed contain both the waves and the turbulent fields.

The fact that ⟨.⟩ denotes a time average also considerably simplifies Eqs. (16) and (22). The Lagrangian mean velocity ⟨u⟩_L is constructed in such a way that when the fluid velocity at X is u_L, then the mean position x is displaced with the velocity ⟨u⟩_L (Andrews & McIntyre 1978). A perhaps more illustrative way of interpreting the quantities ξ, u_L, and ⟨u⟩_L is given in Fig. 1, in the case where ⟨.⟩ is defined as a spatial average in a given direction x. If we isolate a thin tube of fluid lying along this axis, then ξ corresponds to the local deformation of the tube, u_L corresponds to the instantaneous velocity of the local portion of tube, and ⟨u⟩_L corresponds to the velocity of the centre of mass of the tube.

Fig. 1. Visualisation of the fluid displacement ξ (in green), Lagrangian velocity uL (in red), and Lagrangian mean velocity ⟨u⟩L (in blue), in the case where the mean ⟨.⟩ is defined as an average over a given axis (spanning horizontally in the figure). The grey volume VL represents a thin tube of fluid initially lying along the horizontal direction. The fluid displacement and velocity ξ and uL pertain to the deformation and local velocity of the tube, while ⟨u⟩L refers to the velocity of the centre of mass of the tube, represented by the dashed horizontal line. This illustration is inspired by Andrews & McIntyre (1978) (see bottom panel of their Fig. 1).

In our case where ‘average’ refers to time average, ⟨u⟩_L refers to the movement of the centre of mass of a given parcel of fluid in the absence of turbulent convection and oscillations, or in other words, to the background fluid velocity. We note that the effects of rotation, if we were to take them into account, would be encompassed in ⟨u⟩_L. In the absence of rotation, however, we have ⟨u⟩_L = 0, and Eqs. (16) and (22) then reduce to

$$\begin{array}{cc}& {\displaystyle \frac{\partial {\xi }_i}{\partial t}}=u_i(\mathbf{x}+\mathit{\xi },t),\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial u_i}{\partial t}}+u_j{\displaystyle \frac{\partial u_i}{\partial x_j}}=-{\displaystyle \frac{1}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}+g_i+G_{\mathit{ij}}(u_j-\stackrel{\sim }{u_j})+\sqrt{C_0ϵ}{\eta }_i.\hfill \end{array}$$(26)

It must be noted that while the velocity appearing on the right-hand side of Eq. (25) is evaluated at the Lagrangian position x + ξ(x, t), every quantity in Eq. (26), by contrast, is evaluated at the Eulerian position x. Furthermore, it should be noted that everything depends on time, even when not explicitly specified. The key difference between Eqs. (9) and (10), on the one hand, and (25) and (26) on the other, is that the variables whose evolution is described are no longer the Lagrangian quantities x^⋆(t) and u^⋆(t) pertaining to a set of fluid particles, but the Eulerian quantities ξ(x, t) and u(x, t) pertaining to a set of Eulerian fixed positions x. This description will allow for a much more practical derivation of the wave equation in Sect. 3.1.

In particular, the momentum equation now makes the contribution of the turbulent pressure explicitly appear, in the form of the advection term on the left-hand side of Eq. (26). The contribution of the turbulent fluctuations of the gas pressure and the turbulent dissipation, on the other hand, are still contained in the last two terms on the right-hand side. Furthermore, the procedure described in this section allowed us to rigorously separate the effect of the equilibrium background (defined by a time average over very long timescales) from the joint contributions of the turbulence and of the oscillations (both of which are contained in the fluctuations around the equilibrium background), and thus will allow us to study their mutual coupling in a consistent framework.

2.3. Evaluating the mean fields: Smoothed particle hydrodynamics

The stochastic model on fluid displacement ξ and velocity u, comprised of Eqs. (25) and (26), is still not in closed form, as it contains several mean quantities (mean density $\overline{\rho }$, gas pressure $\overline{p}$, and velocity $\stackrel{\sim }{\mathbf{u}}$), as well as the turbulent dissipation ϵ appearing in the diffusion tensor, and the Reynolds stress tensor $\widetilde{u_{i^{″}}{u}_{j^{″}}}$ and shear tensor ${\partial }_i\stackrel{\sim }{u_j}$ on which the drift tensor G_ij depends. These equations must therefore be supplemented with a model for the mean fields.

One possibility is to make use of a large-eddy simulation (LES) (or direct numerical simulation) to simulate the large-scale flow using the exact equations of hydrodynamics, and then use the mean fields yielded by this simulation as external inputs into the stochastic model. However, the mean fields appearing in Eq. (26) are instantaneous averages (for example, $\overline{\rho }$ is the Reynolds-averaged density at a given time t), not time averages. As a result, the ergodic principle cannot be used to extract the means from a LES. The only way to do this would be to consider that horizontal averages in a 3D LES provide an accurate estimate of instantaneous ensemble averages, and would only work in the scope of a 1D model. Furthermore, this procedure would defeat the purpose of what we are trying to achieve; since the mean fields contain the information on the oscillations, without containing the turbulence, treating them as external inputs would effectively amount to modelling the turbulence and the oscillations in a separate manner, which is what we are trying to avoid.

An alternative method makes use of the particle representation we adopt in Sect. 2.1. The set of fluid particles used to represent the flow contains all the required statistical information, so that the mean fields can actually be estimated directly from the set of fluid particles themselves (Welton & Pope 1997). This is the core idea behind particle methods, and particularly smoothed particle hydrodynamics (SPH). The reader can refer to Liu & Liu (2010) or Monaghan (1992) for a comprehensive review on the subject, or to Springel (2010) for the use of SPH in the astrophysical context, but we give an outline of this method in the following.

Ideally, we would like to estimate all local means at a given Eulerian position x by averaging the corresponding particle-level quantity over all fluid particles conditioned on their being located at x. However, implementing this last condition exactly does not yield the required result: for any given position x, any individual fluid particle has exactly zero probability of finding itself at this exact location. Therefore, it is necessary to relax the condition on particle position, and instead of computing means over particles exactly located at x, we compute them over particles within a given compact-support vicinity of x.

Thus, a kernel function K(r) is introduced, which serves as a weighting function to implement the particle-position condition in the estimation of the means. The exact form of K is not important here, but we mention some of its properties, namely that it is a compact-support function, it is normalised to unity, and it is isotropic. The first two properties are mandatory, as the first one ensures that distant particles cannot impact local means, and the second that the estimation of the means is unbiased. The third property makes the subsequent calculations much easier to carry out. A good example is the kernel function used by Welton & Pope (1997)

$$\begin{array}{c}\hfill \begin{array}{ccc}& K(\mathbf{r})=c(1+3{\displaystyle \frac{|\mathbf{r}|}{h}}){(1-{\displaystyle \frac{|\mathbf{r}|}{h}})}^3\hfill & \mathrm{if}|\mathbf{r}|<h,\hfill \\ \hfill & =0\hfill & \mathrm{if}|\mathbf{r}|>h,\hfill \end{array}\end{array}$$(27)

where r is the position of the particle with respect to the centre of the kernel (where the mean is estimated), c = 105/(16πh³) is defined by the normalisation condition⁵, and h is the size of the kernel compact support. This expression ensures that the kernel function and its first two derivatives are continuous at the surface of its support.

The SPH formalism is best formulated if we temporarily return to the representation of the flow as a large set of N particles, whose position and velocity we denote by x^⋆(i) and u^⋆(i), respectively, where i is the index used to identify each particle. For any quantity Q pertaining to the flow representation, if there is an equivalent quantity Q^⋆ in the particle representation, we can estimate the mean value of Q at any Eulerian position x and time t through the following kernel estimator

$$\begin{array}{c}\hfill \overline{Q(\mathbf{x},t)}={\displaystyle \sum _{i=1}^N{\displaystyle \frac{\mathrm{\Delta }{m}^{(i)}}{{\rho }^{(i)}}}{Q}^{\star (i)}(t)K({\mathbf{x}}^{\star (i)}(t)-\mathbf{x}),}\end{array}$$(28)

where Δm⁽ⁱ⁾ is the mass carried by the particle i, and ρ⁽ⁱ⁾ is the mass density characterising particle i. As such, the quantity Δm⁽ⁱ⁾/ρ⁽ⁱ⁾ appearing under the sum corresponds to the lumped volume of fluid that the particle represents.

Setting Q^⋆ = ρ, Q^⋆ = ρu^⋆, and $Q^{\star }=\rho (u_i^{\star }-\stackrel{\sim }{u_i})(u_j^{\star }-\stackrel{\sim }{u_j})$ alternatively in Eq. (28), we find respectively the Reynolds-averaged density, the mass-averaged velocity, and the Reynolds stress tensor

$$\begin{array}{cc}& \overline{\rho }(\mathbf{x},t)={\displaystyle \sum _{i=1}^N\mathrm{\Delta }{m}^{(i)}K({\mathbf{x}}^{\star (i)}(t)-\mathbf{x}),}\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill & \stackrel{\sim }{\mathbf{u}}(\mathbf{x},t)={\displaystyle \frac{1}{\overline{\rho }(\mathbf{x},t)}}{\displaystyle \sum _{i=1}^N\mathrm{\Delta }{m}^{(i)}{\mathbf{u}}^{\star (i)}(t)K({\mathbf{x}}^{\star (i)}(t)-\mathbf{x}),}\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill & \widetilde{u_i^{″}{u}_j^{″}}(\mathbf{x},t)={\displaystyle \frac{1}{\overline{\rho }(\mathbf{x},t)}}{\displaystyle \sum _{i=1}^N\mathrm{\Delta }{m}^{(i)}(u_i^{\star (i)}(t)-\stackrel{\sim }{u_i}(\mathbf{x},t))}\hfill \\ \hfill & \times (u_j^{\star (i)}(t)-\stackrel{\sim }{u_j}(\mathbf{x},t))K({\mathbf{x}}^{\star (i)}(t)-\mathbf{x}).\hfill \end{array}$$(31)

In particular, we note that in the SPH formalism, the local mean density is computed by counting the particles present in the vicinity. This means that the continuity condition is automatically met in the particle representation, thus lowering the order of the set of equations needed to describe the flow.

We then rewrite Eqs. (29)–(31) in the representation chosen in Sect. 2.2, specifically in terms of ξ(x) and u(x) rather than x^⋆(i) and u^⋆(i). To do so, we perform the change of variables given by Eqs. (18) and (19). Furthermore, the sum over infinitesimally small masses can be replaced by a continuous integral over dm ≡ ρ₀(y)d³y, where ρ₀ is the equilibrium fluid density (which can be thought of as an average of the local fluid density over very long timescales so as to only contain the background value). Finally, in this new representation, the SPH formalism yields

$$\begin{array}{cc}& {\displaystyle \overline{\rho }(\mathbf{x},t)={\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y})K(\mathbf{y}+\mathit{\xi }(\mathbf{y},t)-\mathbf{x}),}\hfill \end{array}$$(32)

$$\begin{array}{cc}\hfill & {\displaystyle \stackrel{\sim }{\mathbf{u}}(\mathbf{x},t)={\displaystyle \frac{1}{\overline{\rho }(\mathbf{x},t)}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y})\mathbf{u}(\mathbf{y}+\mathit{\xi }(\mathbf{y},t),t)K(\mathbf{y}+\mathit{\xi }(\mathbf{y},t)-\mathbf{x}),}\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill & {\displaystyle \widetilde{u_i^{″}{u}_j^{″}}(\mathbf{x},t)={\displaystyle \frac{1}{\overline{\rho }(\mathbf{x},t)}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y})(u_i(\mathbf{y}+\mathit{\xi }(\mathbf{y},t),t)-\stackrel{\sim }{u_i}(\mathbf{y},t))}\hfill \\ \hfill & \times (u_j(\mathbf{y}+\mathit{\xi }(\mathbf{y},t),t)-\stackrel{\sim }{u_j}(\mathbf{y},t))K(\mathbf{y}+\mathit{\xi }(\mathbf{y},t)-\mathbf{x}).\hfill \end{array}$$(34)

While these integrals span across the entire volume of the star, we note that they actually only involve the compact support vicinity of x defined by the kernel function K.

Two more points need to be addressed here. The first concerns the mean pressure $\overline{p}$. For the sake of simplicity, we consider a polytropic relation between the gas pressure and density in the form

$$\begin{array}{c}\hfill ln\left({\displaystyle \frac{\overline{p}(\mathbf{x},t)}{p_0(\mathbf{x})}}\right)=\gamma (\mathbf{x})ln\left({\displaystyle \frac{\overline{\rho }(\mathbf{x},t)}{{\rho }_0(\mathbf{x})}}\right),\end{array}$$(35)

where p₀ is the equilibrium gas pressure (defined in the same way as ρ₀), and we allow the polytropic exponent γ to depend on space. We note that we do not consider the possibility that the oscillations may entail fluctuations in the polytropic index itself. We also note that we can recover the isentropic case at any point by setting γ = Γ₁, where Γ₁ is the equilibrium first adiabatic exponent.

The second point concerns the turbulent dissipation rate ϵ, or equivalently the turbulent frequency ω_t defined through

$$\begin{array}{c}\hfill {\omega }_t\equiv {\displaystyle \frac{ϵ}{k}}={\displaystyle \frac{2ϵ}{\widetilde{u_i^{″}{u}_i^{″}}}},\end{array}$$(36)

where k is the turbulent kinetic energy. Physically, ω_t can be interpreted as the inverse of the characteristic lifetime associated with the energy-bearing eddies. The turbulent kinetic energy k is given in closed form by the velocity part of the model (here it is given by half the trace of Eq. (34)), and we still need to model ω_t. Usually, this is done either by adding a model equation for the mass-averaged dissipation rate $\stackrel{\sim }{{\omega }_t}$, which is very similar to the approach followed in two-equation models of turbulence, such as the k − ϵ model (Jones & Launder 1972), or else by adding ${\omega }_t^{\star }$ to the particle properties in the Lagrangian stochastic model, such as in the refined Langevin model (Pope & Chen 1990). However, in the present work, and in the scope of the generalised Langevin model, we regard ω_t as a time-independent equilibrium quantity, which can still, however, depend on x. Physically, this amounts to assuming that all eddies have the same typical lifetime, regardless of their size, but that it can depend on the depth at which they are located. In the long run, it will be necessary to go beyond this drastic assumption.

To recap Sect. 2, the model equations are Eqs. (25) and (26), which are stochastic differential equations governing the evolution of the fluid displacement ξ(x, t) and the Eulerian velocity u(x, t) for any given Eulerian position x. The mean density $\overline{\rho }$, the mass-averaged velocity $\stackrel{\sim }{\mathbf{u}}$, and the Reynolds stress tensor $\widetilde{u_{i^{″}}{u}_{j^{″}}}$ are given by Eqs. (32)–(34), respectively, as explicit functions of ξ(x, t) and u(x, t) only. The mean pressure is given by Eq. (35) as a function of mean density, and the turbulent dissipation rate ϵ is given by Eq. (36). Therefore, all the quantities appearing in the model equations are written explicitly as functions of the modelled variables ξ(x, t) and u(x, t) themselves: the model is in closed form. The only inputs of the model are (1) the equilibrium density ρ₀(x), gas pressure p₀(x), and polytropic exponent γ(x), which can be extracted from an equilibrium model of the star; (2) the functional form of the drift tensor G_ij (see Eq. (8)), which can be constrained using direct numerical simulations or experimental measurements (Pope 1994); and (3) the equilibrium turbulent frequency ω_t(x), which can be constrained using a 3D hydrodynamic simulation of the atmosphere of the star, or on the contrary serve as a control parameter for turbulence, which can be varied for a parametric study.

3. The stochastic wave equation

We now set out to linearise the closed set of equations derived in Sect. 2 to obtain a linear stochastic wave equation that was designed to govern the physics of the mode while simultaneously encompassing the effect of turbulence on the mode. We then discuss the properties of this wave equation, and how it relates to other forms of the wave equation obtained in previous studies.

3.1. Linearisation of the stochastic model

The system to linearise is comprised of Eqs. (25), (26), (32)–(36). The only variables in these equations are the fluid displacement ξ(x, t) and the Eulerian velocity u(x, t). We define

$$\begin{array}{cc}& {\mathit{\xi }}_{\mathrm{osc}}(\mathbf{x},t)\equiv \mathit{\xi }(\mathbf{x},t)-{\mathit{\xi }}_t(\mathbf{x},t),\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill & {\mathbf{u}}_{\mathrm{osc}}(\mathbf{x},t)\equiv \mathbf{u}(\mathbf{x},t)-{\mathbf{u}}_t(\mathbf{x},t),\hfill \end{array}$$(38)

where ξ_t and u_t are the fluid displacement and velocity that would be obtained if there were no oscillations, and that represent the turbulent component of the fluid displacement and velocity, and ξ_osc and u_osc are the oscillatory displacement and velocity. We note that while the variables are split into a turbulent and an oscillatory component, the system of equations itself is not split into equations for turbulence and equations for oscillations; instead, there is still only one system of equations containing both components. This is in contrast, for instance, with MLT or Reynolds-Averaged Navier-Stokes (RANS) approaches where the equations are averaged from the start, thus implicitly separating the two components whose coupling we wish to study.

Obtaining a linear wave equation requires the adoption of a certain number of hypotheses regarding the fluid variables, which we itemise here.

(H1). We consider |u_osc|≪|u_t|. This ordering is justified by the fact that, at the top of the convective envelope of solar-like oscillators, the typical turbulent velocities have much higher amplitudes than the oscillatory velocities; the former are of the order of a few km.s⁻¹, while the latter are of the order of a few tens of cm.s⁻¹. This allows us to treat u_osc as a first-order perturbation compared to u_t, and any second- or higher-order occurrence of u_osc will be discarded.

(H2). We consider |ξ_osc|≪h, H_p, where we recall that h is the size of the averaging kernel function K, and H_p ≡ −(dlnp₀/dr)⁻¹ is the pressure scale height. In other words, the modal fluid displacement is much smaller than the stratification length scale, and the width of the kernel function must be sufficiently large. The first hypothesis is justified by the fact that, in the Sun for instance, the modal displacement is of the order of a few tens of meters, while H_p is of the order of a few hundreds of kilometers. The second hypothesis, on the other hand, constitutes a constraint on h. This allows us to treat ξ_osc as a first-order perturbation compared to all length scales relevant to the problem, and any second- or higher-order occurrence of ξ_osc will be discarded.

(H3). We adopt the anelastic approximation for turbulence, in the sense that we consider ρ_t ≪ ρ₀, where ρ_t is the turbulent fluctuation of density, and ρ₀ the equilibrium density. This is the most severe approximation we make in this section. Nevertheless, the anelastic approximation is widely used in analytical models of turbulent convection in these regions, on the grounds that the flow is subsonic (with turbulent Mach numbers peaking at around 0.3 in the superadiabatic region), as shown by 3D hydrodynamic simulations of the atmosphere of these stars (Nordlund et al. 2009). Using the continuity equation, this amounts to neglecting the quantity ∇ ⋅ (ρ₀ξ_t). As will become apparent in the following, this allows us to discard all ξ_t-dependent contributions in the linearisation of the ensemble averages in the SPH formalism.

(H4). We consider that the turbulent velocity field u_t is the same as it would be without the presence of an oscillating velocity u_osc; in other words, we neglect the back-reaction of the oscillations on the turbulent motions of the gas. We justify this approximation in Appendix C on the basis of a discussion that can be found in Bühler (2009, see their Sect. 5.1.1). We note that the back-reaction being neglected here concerns both the equilibrium part and the stochastic part (i.e. both the equilibrium structure of the star and the turbulent velocity field). This assumption allows us to consider u_t as an input to the model, whose statistical properties (average, covariance, autocorrelation function) are considered completely known.

(H5). We consider that the gravitational potential is not perturbed by the turbulent motions of the gas or by its oscillatory motions. These are actually two separate approximations. The first is justified by the fact that the Reynolds-averaged mass flow through any given horizontal layer due to turbulence is zero, meaning the total mass present beneath this layer is always the same. The second corresponds to the Cowling approximation, and is justified for modes that feature a large number of radial nodes. These two approximations put together allow us to replace the gravitational acceleration g by its equilibrium value g₀, which only depends on the hydrostatic equilibrium of the star.

We insist on the fact that these approximations, with the exception of (H5), only concern the fluid displacement and velocity. By contrast, no specific approximation is adopted concerning the mean fields; a linearised form of these mean fields naturally arises from the SPH formalism and the hypotheses (H1) through (H4). As an example, let us consider the mean density $\overline{\rho }$. Plugging Eq. (37) into Eq. (32), we find

$$\begin{array}{c}\hfill {\displaystyle \overline{\rho }(\mathbf{x},t)={\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y})K(\mathbf{y}+{\mathit{\xi }}_t(\mathbf{y},t)+{\mathit{\xi }}_{\mathrm{osc}}(\mathbf{y},t)-\mathbf{x}).}\end{array}$$(39)

Then, using (H2) to linearise in terms of the displacement, we find

$$\begin{array}{cc}& {\displaystyle \overline{\rho }(\mathbf{x},t)={\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y})K(\mathbf{y}-\mathbf{x})+{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y}){({\mathit{\xi }}_{\mathrm{osc}}(\mathbf{y},t)·\mathbf{\nabla })K|}_{\mathbf{y}-\mathbf{x}}}\hfill \\ \hfill & {\displaystyle +{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y}){({\mathit{\xi }}_t(\mathbf{y},t)·\mathbf{\nabla })K|}_{\mathbf{y}-\mathbf{x}}.}\hfill \end{array}$$(40)

The first term on the right-hand side corresponds to the kernel estimate of ρ₀ at x. By construction, kernel estimation is a representation of ensemble averaging, but ρ₀ is already an equilibrium quantity, and therefore is equal to its own ensemble average. Furthermore, the last term on the right-hand side can be discarded on account of hypothesis (H3). Performing an integration by part makes the quantity ∇ ⋅ (ρ₀ξ_t) appear under the integral sign. Therefore, we eventually find

$$\begin{array}{c}\hfill \overline{\rho }(\mathbf{x},t)={\rho }_0(\mathbf{x})+\underset{\equiv {\rho }_1(\mathbf{x},t)}{\underbrace{{\displaystyle {\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y}){({\mathit{\xi }}_{\mathrm{osc}}(\mathbf{y},t)·\mathbf{\nabla })K|}_{\mathbf{y}-\mathbf{x}}}}}.\end{array}$$(41)

The quantity ρ₁ represents the Eulerian modal fluctuations of density, but of important note is the fact that at no point did we explicitly decompose $\overline{\rho }$ into an equilibrium value ρ₀ and a residual, oscillatory part; instead, the decomposition (41) arises naturally from the SPH formalism, and hypotheses (H2) and (H3).

The linear wave equation is derived in detail in Appendix D using the hypotheses listed above. Ultimately, we obtain

$$\begin{array}{cc}& {\displaystyle \frac{\partial {\mathit{\xi }}_{\mathrm{osc}}}{\partial t}}-{\mathbf{u}}_{\mathrm{osc}}-({\mathit{\xi }}_{\mathrm{osc}}·\mathbf{\nabla }){\mathbf{u}}_t-({\mathit{\xi }}_{\mathit{t}}·\mathbf{\nabla }){\mathbf{u}}_{\mathrm{osc}}=({\mathit{\xi }}_{\mathit{t}}·\mathbf{\nabla }){\mathbf{u}}_t,\hfill \end{array}$$(42)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial {\mathbf{u}}_{\mathrm{osc}}}{\partial t}}-{\mathbf{L}}_1^d-{\mathbf{L}}_1^s={\mathbf{L}}_0,\hfill \end{array}$$(43)

where

$$\begin{array}{cc}& {\displaystyle L_{1,i}^d=[{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\partial p_0}{\partial x_i}}-{\displaystyle \frac{\partial c_0^2}{\partial x_i}}]{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y}){\left({\xi }_{\mathrm{osc},\mathrm{j}}{\partial }_j{K}^{\mathbf{x}}\right)|}_{\mathbf{y},\mathbf{t}}}\hfill \\ \hfill & {\displaystyle +{\displaystyle \frac{c_0^2}{{\rho }_0}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y}){\left({\xi }_{\mathrm{osc},\mathrm{j}}{\partial }_j{\partial }_i{K}^{\mathbf{x}}\right)|}_{\mathbf{y},\mathbf{t}}}\hfill \\ \hfill & +G_{ij,0}\left({\displaystyle u_{\mathrm{osc},\mathrm{j}}-{\displaystyle \frac{1}{{\rho }_0(\mathbf{x})}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y})u_{\mathrm{osc},\mathrm{j}}(\mathbf{y})K^{\mathbf{x}}(\mathbf{y})}\right),\hfill \end{array}$$(44)

$$\begin{array}{cc}& L_{1,i}^s=-u_{\mathrm{osc},\mathrm{j}}{\partial }_j{u}_{t,i}-u_{t,j}{\partial }_j{u}_{\mathrm{osc},\mathrm{i}}\hfill \\ \hfill & {\displaystyle -G_{ij,0}{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y}){\left({\xi }_{\mathrm{osc},\mathrm{k}}{\partial }_k\left(u_{t,j}{K}^{\mathbf{x}}\right)\right)|}_{\mathbf{y},\mathbf{t}}}\hfill \\ \hfill & +({\displaystyle \frac{\partial G_{\mathit{ij}}}{\partial \widetilde{u_k^{″}{u}_l^{″}}}}{\widetilde{u_k^{″}{u}_l^{″}}}_1+{\displaystyle \frac{\partial G_{\mathit{ij}}}{\partial ({\partial }_k\stackrel{\sim }{u_l})}}{\partial }_k\stackrel{\sim }{u_l}+{\displaystyle \frac{\partial G_{\mathit{ij}}}{\partial ϵ}}{\omega }_t{k}_1)u_{t,j}\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{C_0{\omega }_t}{k_0}}}{k}_1{\eta }_i,\hfill \end{array}$$(45)

$$\begin{array}{c}\hfill L_{0,i}=-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\partial ({\rho }_0{u}_{t,i}{u}_{t,j}-{\rho }_0\overline{u_{t,i}{u}_{t,j}})}{\partial x_j}},\end{array}$$(46)

$$\begin{array}{cc}& {\displaystyle {\widetilde{u_i^{″}{u}_j^{″}}}_1=-{\displaystyle \frac{{\widetilde{u_i^{″}{u}_j^{″}}}_0}{{\rho }_0}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0{\xi }_{\mathrm{osc},\mathrm{k}}{\partial }_k{K}^{\mathbf{x}}}\hfill \\ \hfill & {\displaystyle +{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0{\xi }_{\mathrm{osc},\mathrm{k}}{\partial }_k\left(u_{t,i}{u}_{t,j}{K}^{\mathbf{x}}\right)}\hfill \\ \hfill & {\displaystyle +{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0{u}_{t,i}{u}_{\mathrm{osc},\mathrm{j}}{K}^{\mathbf{x}}}\hfill \\ \hfill & {\displaystyle +{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0{u}_{t,j}{u}_{\mathrm{osc},\mathrm{i}}{K}^{\mathbf{x}},}\hfill \end{array}$$(47)

$$\begin{array}{cc}& {\displaystyle {({\partial }_i\stackrel{\sim }{u_j})}_1=-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0{u}_{\mathrm{osc},\mathrm{j}}{\partial }_i{K}^{\mathbf{x}}-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0{\xi }_{\mathrm{osc},\mathrm{k}}{\partial }_k{u}_{t,j}{\partial }_i{K}^{\mathbf{x}}}\hfill \\ \hfill & {\displaystyle -{\displaystyle \frac{1}{{\rho }_0^2}}{\displaystyle \frac{\partial {\rho }_0}{\partial x_i}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0{u}_{\mathrm{osc},\mathrm{j}}{K}^{\mathbf{x}}}\hfill \\ \hfill & {\displaystyle -{\displaystyle \frac{1}{{\rho }_0^2}}{\displaystyle \frac{\partial {\rho }_0}{\partial x_i}}{\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0{\xi }_{\mathrm{osc},\mathrm{k}}{\partial }_k\left(u_{t,j}{K}^{\mathbf{x}}\right),}\hfill \end{array}$$(48)

$$\begin{array}{c}\hfill k_1={\displaystyle \frac{1}{2}}\widetilde{u_i^{″}{u}_i^{″}},\end{array}$$(49)

where $c_0^2\equiv p_0\gamma /{\rho }_0$ is the equilibrium sound speed squared, and we have introduced the x-centred kernel function K^x(y)≡K(y − x). The right-hand sides of Eqs. (42) and (43) constitute inhomogeneous forcing terms (see Sect. 3.2 for more details), and it was therefore possible to filter out certain negligible contributions. We refer the reader to the details given in Appendix D.3, where we essentially argue first that the non-stochastic zeroth-order terms in the linearised equations vanish under the unperturbed hydrostatic equilibrium condition, and then that the linear forcing (i.e. the stochastic zeroth-order terms that are linear in the stochastic processes ξ_t, u_t or η) is negligible. The last point is related to the fact that the turbulent spectrum has most of its power in wavevectors and angular frequencies far removed from those characteristic of the modes, and therefore unable to provide with efficient mode driving.

Formally, Eqs. (42) and (43) take the form of a linear stochastic, inhomogeneous wave equation in a completely closed form, in the sense that the various terms on their right-hand side are written as explicit functions of the wave variables ξ_osc and u_osc themselves or the turbulent fields ξ_t and u_t, whose statistical properties are considered known (see hypothesis H4). In writing Eq. (43), we split the velocity equation into three components. ${\mathbf{L}}_1^d$ contains all the terms that are linear in ξ_osc and u_osc, but do not explicitly contain either stochastic processes ξ_t, u_t, or η. It represents the deterministic contribution to the homogeneous part of the wave equation, and corresponds to the classical propagation of acoustic waves, without any impact from the turbulence. On the other hand, ${\mathbf{L}}_1^s$ contains all the terms that are linear in ξ_osc and u_osc and explicitly depend on ξ_t, u_t, or η. Finally, L₀ contains all the terms that are independent of ξ_osc and u_osc. The reason for this specific decomposition will become apparent in a moment, when we discuss the physical role played by each term.

3.2. Effects of turbulence on the wave equation

The last term on the left-hand side of Eq. (43), together with its right-hand side, contain the contribution of turbulence to the wave equation, which arises from the action of the turbulent fields on the oscillations. We can see that one effect of turbulence is to add the inhomogeneous part L₀ to the wave equation. This part acts as a forcing term, and Eq. (46) shows that it corresponds to the fluctuations of the turbulent pressure around its ensemble average. This is in perfect accordance with the widely accepted picture that the stochastic excitation of the global modes of oscillation in solar-like stars is due mainly to quadrupolar turbulent acoustic emission (Samadi & Goupil 2001). Furthermore, we note that we only kept the contributions to mode excitation that are not linear in the turbulent velocity field as linear contributions turn out to be negligible (see Appendix D.3 for a more developed discussion). We note that the non-linear Lagrangian, turbulent fluctuations of entropy, which is widely recognised as another source of stochastic driving for solar-like p-modes, does not arise from the above formalism. The only reason is because we considered a polytropic equation of state from the start, and as such neglected to model entropy fluctuations in both the oscillations and the turbulence.

The second effect of the turbulent fields on the oscillations is to modify the linear part, that is to say the propagation of the waves. This stochastic correction corresponds to the term ${\mathbf{L}}_1^s$ defined by Eq. (45). This term models two effects that are usually studied as distinct phenomena, but are actually intertwined and cannot be considered separately: a shift in the eigenfrequency of the resonant modes of the system (commonly referred to as the modal or ‘intrinsic’ part of the surface effects), and the absorption, or damping, of the energy of the waves as they travel through the turbulent medium. Equation (45) shows that these phenomena arise either from the non-linear advection term in the momentum equation, as represented by the first two terms on its right-hand side, or from the joint effect of turbulent dissipation, buoyancy, and pressure-rate-of-strain correlations, as jointly represented by all the other terms. It is apparent, in particular, that while the former is linear in u_t, the latter has a more complicated multipolar decomposition in terms of u_t, with first-, second-, and third-order contributions alike. As a whole, the term ${\mathbf{L}}_1^s$ in the velocity equation plays the same role, for instance, as 𝒟(v_osc) in Samadi & Goupil (2001) (see their Eq. (26)).

3.3. Limiting case: The standard wave equation

We now explore the limiting case where there is no turbulence, in which case the only term that remains in Eq. (43) is ${\mathbf{L}}_1^d$. In the absence of turbulence the integrals appearing in Eq. (44) are drastically simplified because, in this limit, the wave variables ξ_osc and u_osc are equal to their own ensemble average, that is to say to their own kernel estimate. This allows us to write, for instance

$$\begin{array}{c}\hfill {\displaystyle {\displaystyle \int }{\mathrm{d}}^3\mathbf{y}{\rho }_0(\mathbf{y}){\xi }_{\mathrm{osc},\mathrm{i}}(\mathbf{y},t){{\displaystyle \frac{\partial K}{\partial x_j}}|}_{\mathbf{y}-\mathbf{x}}=-{\displaystyle \frac{\partial {\rho }_0{\xi }_{\mathrm{osc},\mathrm{i}}}{\partial x_j}}.}\end{array}$$(50)

Ultimately, this leads to the following simplification of Eq. (44):

$$\begin{array}{c}\hfill L_{1,i}^d=[{\displaystyle \frac{\partial c_0^2}{\partial x_i}}-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\partial p_0}{\partial x_i}}]{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\partial {\rho }_0{\xi }_{\mathrm{osc},\mathrm{j}}}{\partial x_j}}+{\displaystyle \frac{c_0^2}{{\rho }_0}}{\displaystyle \frac{{\partial }^2{\rho }_0{\xi }_{\mathrm{osc},\mathrm{j}}}{\partial x_i\partial x_j}}.\end{array}$$(51)

Hence from Eqs. (42) and (43), which in this limit read

$$\begin{array}{cc}& {\displaystyle \frac{\partial {\mathit{\xi }}_{\mathrm{osc}}}{\partial t}}-{\mathbf{u}}_{\mathrm{osc}}=\mathbf{0},\hfill \end{array}$$(52)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial {\mathbf{u}}_{\mathrm{osc}}}{\partial t}}-{\mathbf{L}}_1^d=\mathbf{0},\hfill \end{array}$$(53)

we obtain the following homogeneous, second-order wave equation

$$\begin{array}{c}\hfill ({\displaystyle \frac{{\partial }^2}{\partial t^2}}-\mathcal{L}){\mathbf{u}}_{\mathrm{osc}}=\mathbf{0},\end{array}$$(54)

with

$$\begin{array}{c}\hfill \mathcal{L}\left({\mathbf{u}}_{\mathrm{osc}}\right)={\displaystyle \frac{1}{{\rho }_0}}[\left(\mathbf{\nabla }{c}_0^2(\mathbf{\nabla }·{\rho }_0{\mathbf{u}}_{\mathrm{osc}})\right)-{\displaystyle \frac{1}{{\rho }_0}}\mathbf{\nabla }{p}_0(\mathbf{\nabla }·{\rho }_0{\mathbf{u}}_{\mathrm{osc}})].\end{array}$$(55)

We recover the equation for free acoustic oscillations in a stratified medium in its exact form, provided the Cowling approximation is adopted (see hypothesis H5). It corresponds exactly to the homogeneous part of the wave equation derived, for instance by Goldreich & Keeley (1977b) (see their Eq. (16)); or, equivalently, to the equation presented in Unno et al. (1989) (see their Eqs. (14.2) and (14.3)), although these are written in terms of displacement and pressure fluctuation rather than displacement and velocity (see also Samadi & Goupil 2001, their Eq. (16)). The only exception is the absence of the term depending on the entropy gradient (which does not appear here because the gas is polytropic).

4. Discussion

The linear stochastic wave equation comprised of Eqs. (42) and (43) was obtained in the scope of a certain number of hypotheses and approximations, whose validity we now discuss. We can split these hypotheses into two categories: those pertaining to the establishment of the stochastic differential equations, and those pertaining to the linearisation of these equations.

All the hypotheses we adopted in the linearisation are itemised in Sect. 3.1. Hypotheses (H1) through (H5) are actually of two different natures. Hypotheses (H1) on the smallness of u_osc, (H2) on the smallness of ξ_osc, and (H4) on the absence of back-reaction of the oscillations on the turbulence define the framework in which we performed the linearisation, and are therefore necessary assumptions. On the other hand, hypotheses (H3) on the neglect of ∇ ⋅ (ρ₀ξ_t) and (H5) on the neglect of the perturbed gravitational potential are simplifying assumptions that are not necessary, strictly speaking, but help simplify the formalism considerably. Hypothesis (H5) is a common assumption in the analysis of stellar oscillations: without it, because gravity is an unscreened force acting on long distances, the resulting equations would be highly non-local. As we mention in Sect. 3.1, its domain of validity is the high-radial-order modes of oscillation, but it is usually adopted throughout the entire oscillation spectrum. Hypothesis (H3), on the other hand, may require some more discussion. As we briefly mention above, it corresponds to the anelastic approximation, and amounts to neglecting the turbulent fluctuations of the fluid density ρ_t. Taking these fluctuations into account would require having knowledge of the statistical properties of ρ_t, the same way we consider the properties of u_t known. However, current models of compressible turbulence are not yet able to fully account for ρ_t without any underlying simplifying assumptions, such as the Boussinesq approximation or the anelastic approximation. It is difficult to assess how sensible to this assumption the results obtained for the behaviour of turbulent convection are, and a fortiori its coupling with oscillations, but for lack of a more realistic treatment of turbulence compressibility, we nevertheless chose to adopt hypothesis (H3).

We have also adopted a number of approximations in order to establish the closed system of equations (Eqs. (25), (26), (32)–(36)) in Sect. 2. All of them consist in simplifying assumptions, that we adopt not because they are necessary to build the formalism, but because the aim of this paper is to make the basics of this method as clear as possible, rather than adopting the most realistic turbulence model possible. As such, we do not attempt to give a physical justification for the following hypotheses, but instead discuss how they affect the final stochastic wave equation, and how one would go about circumventing these simplifications.

(H6). We consider the flow to be adiabatic, in the sense that the only fluid particle properties that need to be described in the Lagrangian stochastic model are the position and velocity of the particles. In the scope of this hypothesis, the energy equation is replaced with a relation between the mean density and pressure that we chose to be polytropic, without specifying the associated polytropic exponent γ, which means that the non-adiabatic effects pertaining to the oscillations are not contained in the formalism presented in this paper. This includes the perturbation of the convective flux and the radiative flux by the oscillations, which are in reality susceptible to affect the damping rate of the modes as well as the surface effects. Avoiding hypothesis (H6) would allow for the inclusion of all non-adiabatic effects in the model. Essentially, adopting a non-adiabatic framework would require an additional SDE for the internal energy of the fluid particles (or any other alternative thermodynamic variable), thus leading to the introduction of an additional thermodynamic wave variable e_osc, to be linearised around the turbulence-induced energy fluctuations e_t. This would then increase the order of the system of equations, and would require the statistical properties of the additional turbulent field e_t, including its correlation with u_t, to be known.

(H7). We consider that the turbulent frequency ω_t, defined by Eq. (36) as the ratio of the dissipation rate ϵ to the turbulent kinetic energy k, takes a constant value. The turbulent frequency represents the rate at which k would decay towards zero if there were no production of turbulence whatsoever, and can be interpreted as the inverse lifetime of the energy-containing turbulent eddies. In essence, this amounts to assuming the existence of a single timescale associated with the entire turbulent cascade, which is at odds with even the simplest picture of turbulence. Avoiding hypothesis (H7) would allow a much more realistic modelling of the turbulent dissipation and its perturbation by the oscillations, which is likely to play an important role in both mode damping and surface effects. This would require including the turbulent frequency as a fluid particle property, with its own SDE. As for velocity or internal energy (see hypothesis (H6) above), this would lead to the introduction of ω_t, osc as an additional wave variable, to be linearised around a new turbulent field ω_t, t, whose statistical properties would have to be input in the model.

(H8). We consider that the time average of the flow velocity over very long timescale, in other words the velocity associated with the equilibrium background, is zero. This amounts to neglecting rotation, whether it be global or differential. Taking rotation into account would require either a non-zero ⟨u⟩_L field to be included, or else a Coriolis inertial force to be added in the velocity SDE.

In summary, hypotheses (H1), (H2), and (H4) are fundamental in building the formalism, and cannot be avoided, but they are also firmly and physically grounded. Hypotheses (H3) and (H5) are simplifying assumptions that are not strictly necessary, nor as clearly valid, but which are unavoidable given the current state of our capabilities. Finally, hypotheses (H6), (H7), and (H8) are also simplifying assumptions, and are very much invalid; however, we adopted them here to provide a simple framework serving as a proof of concept for the formalism presented in this paper. In particular, hypotheses (H6) and (H7) must be discarded as soon as possible if a realistic model of turbulence is to be adopted. This is left to a future work in this series.

5. Conclusion

In this series of papers we investigate Lagrangian stochastic models of turbulence as a rigorous way of modelling the various phenomena arising from the interaction between the highly turbulent motions of the gas at the top of the convection zone in solar-like stars and the global acoustic modes of oscillation developing in these stars. These include the stochastic excitation of the modes, their stochastic damping, and the turbulence-induced shift in their frequency called surface effects.

In this first paper we presented a very simple polytropic Lagrangian stochastic turbulence model, serving as a proof of concept for the novel method presented here, and we showed how it can be used to derive stochastic differential equations (SDEs) governing the evolution of Eulerian fluid variables relevant to the study of oscillations. We then linearised these SDEs to obtain a linear stochastic wave equation containing, in the most self-consistent way possible, the terms arising from the turbulence-oscillation coupling. This wave equation correctly reduces to the classical propagation of free acoustic waves in a stratified medium in the limit where turbulence is neglected. It also exactly models the stochastic forcing term due to turbulent acoustic emission, arising from coherent fluctuations in the turbulent pressure. In addition, the resulting stochastic wave equation contains the turbulent-induced correction to the linear operator governing the propagation of the waves, thus allowing for the modelling of both mode damping and modal surface effects. The method presented here offers multiple, key advantages:

• At no point does it require separating the equations of the flow into a turbulent equation and an oscillation equation, thus allowing the turbulent contribution to naturally and consistently arise in the wave equation. Instead, we leave the statistical properties of the turbulence as known oscillation-independent inputs to the model.

• All aspects of turbulence-oscillation interaction are modelled simultaneously, within the same stochastic wave equation, thus shedding a more consistent light on these intertwined phenomena.

• This method completely circumvents the need to adopt the mixing-length hypothesis, which is crucial as this hypothesis is both almost inescapable in current convection modelling, and very invalid close to the radiative-convective transition zone. The reason we do not need to adopt this assumption stems from the fact that the starting model is at particle level, where equations are much easier to close.

• The parameters appearing in Lagrangian stochastic models are much more easily linked to the underlying physical assumptions, and therefore easier to constrain, with the help of 3D hydrodynamic simulations. They are also more firmly physically grounded.

• In addition, this formalism applies to radial and non-radial oscillations alike.

However, this paper only constitutes a first step. In the following paper in this series we will show how such a stochastic wave equation can be used to yield a set of stochastic differential equations governing the temporal evolution of the complex amplitude of the modes. These simplified amplitude equations (Stratonovich 1965) are much more practical for the study of turbulence-oscillation coupling, and in particular explicitly and simultaneously yield the excitation rates of the modes, their lifetimes, as well as their turbulence-induced frequency corrections. Finally, as we mentioned in Sect. 4, extending this work to a non-adiabatic model (i.e. discarding hypothesis H6), and with a more realistic treatment of eddy lifetimes (i.e. discarding hypothesis H7), constitutes an essential and unavoidable step to apply this formalism to the actual stellar case, and will be the subject of a subsequent paper.

Acknowledgments

The authors wish to thank the anonymous referee for his/her insightful comments, which helped improve the clarity and quality of this manuscript.

References

Appendix A: The equivalent Reynolds stress model

The procedure leading from a given Lagrangian stochastic model to the equivalent Reynolds stress model (i.e. the corresponding transport equations for the first- and second-order moments of the flow velocity) can be found, for instance, in Pope (2000, Chap. 12). For the generalised Langevin model considered in this paper, it yields

$$\begin{array}{c}\hfill {\displaystyle \frac{D\overline{\rho }}{\mathit{Dt}}}+\overline{\rho }{\displaystyle \frac{\partial \stackrel{\sim }{u_i}}{\partial x_i}}=0,\end{array}$$(A.1)

$$\begin{array}{c}\hfill {\displaystyle \frac{D\stackrel{\sim }{u_i}}{\mathit{Dt}}}+{\displaystyle \frac{1}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{\rho }\widetilde{u_i^{″}{u}_j^{″}}}{\partial x_j}}=-{\displaystyle \frac{1}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}+g_i,\end{array}$$(A.2)

and

$$\begin{array}{cc}& {\displaystyle \frac{D\widetilde{u_i^{″}{u}_j^{″}}}{\mathit{Dt}}}+{\displaystyle \frac{1}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{\rho }\widetilde{u_i^{″}{u}_j^{″}{u}_k^{″}}}{\partial x_k}}=-\widetilde{u_i^{″}{u}_j^{″}}{\displaystyle \frac{\partial \stackrel{\sim }{u_k}}{\partial x_k}}-\widetilde{u_i^{″}{u}_k^{″}}{\displaystyle \frac{\partial \stackrel{\sim }{u_j}}{\partial x_k}}-\widetilde{u_k^{″}{u}_j^{″}}{\displaystyle \frac{\partial \stackrel{\sim }{u_i}}{\partial x_k}}\hfill \\ \hfill & +G_{\mathit{ik}}\widetilde{u_j^{″}{u}_k^{″}}+G_{\mathit{jk}}\widetilde{u_i^{″}{u}_k^{″}}+C_0ϵ{\delta }_{\mathit{ij}},\hfill \end{array}$$(A.3)

where δ_ij denotes the Kronecker symbol, and we have introduced the pseudo-Lagrangian particle derivative $D/Dt\equiv {\partial }_t+\stackrel{\sim }{u_i}{\partial }_i$.

Equation (A.1) yields the continuity equation in its exact form, without having to include an evolution equation for density at particle level. This is due to the fact that particle positions are advanced through time using their own individual velocities; since each particle carries its own unchanging mass, then by construction there can be no local mass loss or gain.

Equation (A.2) also yields the mean momentum equation in its exact form, primarily because the mean force in the stochastic model is already included in its exact form from the start. We note, however, that the transport term (i.e. the second term on the left-hand side of Eq. A.2) is also modelled exactly, even though it is not explicitly included in any way in Eqs. (9) and (10). This is, once again, because of the Lagrangian nature of the stochastic model, and is incidentally one of its most interesting features: all advection terms are implicitely and exactly modelled because trajectories integrated through Eqs. (9) and (10) coincide with actual fluid particle trajectories.

Equation (A.3) differs slightly from the exact Reynolds stress equation derived directly from the Navier-Stokes equation, which reads

$$\begin{array}{cc}& {\displaystyle \frac{D\widetilde{u_i^{″}{u}_j^{″}}}{\mathit{Dt}}}+{\displaystyle \frac{1}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{\rho }\widetilde{u_i^{″}{u}_j^{″}{u}_k^{″}}}{\partial x_k}}=-\widetilde{u_i^{″}{u}_j^{″}}{\displaystyle \frac{\partial \stackrel{\sim }{u_k}}{\partial x_k}}-\widetilde{u_i^{″}{u}_k^{″}}{\displaystyle \frac{\partial \stackrel{\sim }{u_j}}{\partial x_k}}-\widetilde{u_k^{″}{u}_j^{″}}{\displaystyle \frac{\partial \stackrel{\sim }{u_i}}{\partial x_k}}\hfill \\ \hfill & \mathrm{sym}(-{\displaystyle \frac{\overline{u_i^{″}}}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{p}}{\partial x_j}}-{\displaystyle \frac{1}{\overline{\rho }}}\overline{u_i^{″}{\displaystyle \frac{\partial p^{\prime }}{\partial x_j}}}-{ϵ}_{\mathit{ij}}),\hfill \end{array}$$(A.4)

where ‘sym’ refers to the symmetric part of the tensor inside the brackets. Several contributions are still modelled in their exact form, including the transport term (the second term on the left-hand side), but also the production term (the first three terms on the right-hand side). However, the last three terms are not modelled exactly, and are (from left to right) the buoyancy contribution, the pressure-rate-of-strain tensor, and the dissipation tensor. Comparing Eqs. (A.3) and (A.4), it can be seen that these contributions are collectively modelled by the last two terms in Eq. (10), which correspond to the fluctuating part of the force acting upon the fluid particle. More specifically, we obtain

$$\begin{array}{c}\hfill \mathrm{sym}(-{\displaystyle \frac{\overline{u_i^{″}}}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}-{\displaystyle \frac{1}{\overline{\rho }}}\overline{u_i^{″}{\displaystyle \frac{\partial p^{\prime }}{\partial x_i}}}-{ϵ}_{\mathit{ij}})=G_{\mathit{ik}}\widetilde{u_j^{″}{u}_k^{″}}+G_{\mathit{jk}}\widetilde{u_i^{″}{u}_k^{″}}+C_0ϵ{\delta }_{\mathit{ij}}.\end{array}$$(A.5)

This equation is more readily interpreted if we remember that, in the high Reynolds number limit, the dissipation tensor is isotropic, which allows for the definition of the scalar dissipation ϵ appearing in the stochastic model:

$$\begin{array}{c}\hfill {ϵ}_{\mathit{ij}}\equiv {\displaystyle \frac{2}{3}}ϵ{\delta }_{\mathit{ij}}.\end{array}$$(A.6)

Furthermore, the drift tensor is usually decomposed into an isotropic and anisotropic part, according to

$$\begin{array}{c}\hfill G_{\mathit{ij}}=-({\displaystyle \frac{1}{2}}+{\displaystyle \frac{3}{4}}{C}_0){\displaystyle \frac{ϵ}{k}}{\delta }_{\mathit{ij}}+G_{\mathit{ij}}^a,\end{array}$$(A.7)

where $k\equiv \widetilde{u_i^{″}{u}_i^{″}}/2$ is the turbulent kinetic energy. This decomposition ensures that, in the special case of incompressible, homogeneous, isotropic turbulence, if we take $G_{ij}^a=0$, the evolution of the Reynolds stress tensor reduces to the exact, analytical solution.

Equation (A.5) can be rearranged to yield

$$\begin{array}{cc}& \mathrm{sym}(-{\displaystyle \frac{\overline{u_i^{″}}}{\overline{\rho }}}{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}-{\displaystyle \frac{1}{\overline{\rho }}}\overline{u_i^{″}{\displaystyle \frac{\partial p^{\prime }}{\partial x_i}}})=G_{\mathit{ik}}^a\widetilde{u_j^{″}{u}_k^{″}}+G_{\mathit{jk}}^a\widetilde{u_i^{″}{u}_k^{″}}\hfill \\ \hfill & -(1+{\displaystyle \frac{3}{2}}{C}_0){\displaystyle \frac{ϵ}{k}}(\widetilde{u_i^{″}{u}_j^{″}}-{\displaystyle \frac{2}{3}}k{\delta }_{\mathit{ij}}).\hfill \end{array}$$(A.8)

Equation (A.8) allows us to interpret the collective effect of buoyancy and pressure-rate-of-strain correlation on the evolution of the Reynolds stresses. First, all the terms on the right-hand side are traceless, which means that they only have a redistributive role; they redistribute energy among the different components of the Reynolds stress tensor, without ever resulting in a net loss or gain of energy. By contrast, it is the scalar dissipation ϵ which is responsible for the decay of kinetic turbulent energy, an effect that is only counterbalanced by the shear- and compression-induced production term (i.e. the first three terms on the right-hand side of Eq. (A.4)).

Furthermore, it is readily seen that the last term on the right-hand side of Eq. (A.8) tends to isotropise the Reynolds stress tensor since for isotropic turbulence we would precisely have $\widetilde{u_i^{″}{u}_j^{″}}=2k{\delta }_{\mathit{ij}}/3$. The rate at which this term makes the Reynolds stress decay towards isotropy is equal to (1 + 3C₀/2)ω_t, where ω_t is the turbulent dissipation rate defined by Eq. (36). On the other hand, the other two terms on the right-hand side of Eq. (A.8) create anisotropy in the Reynolds stress tensor, and we can intuitively understand that the anisotropy of the stationary Reynolds stress results from a balance between these two effects.

Appendix B: A detailed derivation for the Lagrangian-to-Eulerian change of variables

The goal of this appendix is to provide a detailed derivation of the various steps in the procedure described in Section 2.2. We first derive the general identity (14), which is valid for any fluid quantity; we then apply this general identity to the displacement and velocity variables.

B.1. Derivation of identity (14)

Let us consider, for the moment, that the function ξ(x, t) is an arbitrary function of space and time, which we do not specify at first. We recall the following notations

$$\begin{array}{cc}& \mathbf{X}(\mathbf{x},t)\equiv \mathbf{x}+\mathit{\xi }(\mathbf{x},t),\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill & {\varphi }_L(\mathbf{x},t)\equiv \varphi (\mathbf{X}(\mathbf{x},t),t),\hfill \end{array}$$(B.2)

where ϕ is an arbitrary quantity. The usual chain rules for derivation then yield

$$\begin{array}{cc}& {\displaystyle \frac{\partial ({\varphi }_L)}{\partial t}}={\left({\displaystyle \frac{\partial \varphi }{\partial t}}\right)}_L+{\displaystyle \frac{\partial X_i}{\partial t}}{\left({\displaystyle \frac{\partial \varphi }{\partial x_i}}\right)}_L,\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial ({\varphi }_L)}{\partial x_i}}={\displaystyle \frac{\partial X_j}{\partial x_i}}{\left({\displaystyle \frac{\partial \varphi }{\partial x_j}}\right)}_L.\hfill \end{array}$$(B.4)

If the function x ↦ X(x, t) (the time t being fixed) is bijective, then for any velocity field u(x, t), there necessarily exists an associated field V(x, t) such that, were point x to move with velocity V, point X would then move with the actual fluid velocity u_L. Otherwise stated, V corresponds to the advective velocity in the material derivative of X, so that

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial X_i}{\partial t}}+V_j{\displaystyle \frac{\partial X_i}{\partial x_j}}=u_{i,L}.\end{array}$$(B.5)

Plugging Eq. (B.5) into Eq. (B.3), we find

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial ({\varphi }_L)}{\partial t}}& ={\left({\displaystyle \frac{\partial \varphi }{\partial t}}\right)}_L+(u_{i,L}-V_j{\displaystyle \frac{\partial X_i}{\partial x_j}}){\left({\displaystyle \frac{\partial \varphi }{\partial x_i}}\right)}_L\hfill \\ \hfill & ={\left({\displaystyle \frac{D\varphi }{\mathit{Dt}}}\right)}_L-V_j{\displaystyle \frac{\partial X_i}{\partial x_j}}{\left({\displaystyle \frac{\partial \varphi }{\partial x_i}}\right)}_L,\hfill \end{array}$$(B.6)

where D/Dt ≡ ∂_t + u_i∂_i. In turn, plugging Eq. (B.4), this transforms into

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial ({\varphi }_L)}{\partial t}}={\left({\displaystyle \frac{D\varphi }{\mathit{Dt}}}\right)}_L-V_j{\displaystyle \frac{\partial ({\varphi }_L)}{\partial x_j}},\end{array}$$(B.7)

thus yielding the required identity

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{D\varphi }{\mathit{Dt}}}\right)}_L={⟨D⟩}_L({\varphi }_L),\end{array}$$(B.8)

where ⟨D⟩_L ≡ ∂_t + V_i∂_i.

This is all valid regardless of the definition of the function ξ(x, t). However, let us now consider that x does actually correspond to a mean⁶ position, and that the function ξ(x, t) actually denotes the fluctuating fluid displacement around this mean position x. Then by construction, we have

$$\begin{array}{c}\hfill ⟨{\xi }_i⟩=0.\end{array}$$(B.9)

Point x now corresponding to a mean position, the velocity V at which it is displaced must itself be a mean quantity, so that

$$\begin{array}{c}\hfill ⟨\mathbf{V}⟩=\mathbf{V}.\end{array}$$(B.10)

Let us now apply the mean operator ⟨.⟩ to Eq. (B.5); since V_j can be pulled out of the mean, we obtain

$$\begin{array}{cc}\hfill {⟨u_i⟩}_L& ={\displaystyle \frac{\partial ⟨X_i⟩}{\partial t}}+V_j{\displaystyle \frac{\partial ⟨X_i⟩}{\partial x_j}}\hfill \\ \hfill & ={\displaystyle \frac{\partial x_i}{\partial t}}+{\displaystyle \frac{\partial ⟨{\xi }_i⟩}{\partial t}}+V_j({\displaystyle \frac{\partial x_i}{\partial x_j}}+{\displaystyle \frac{\partial ⟨{\xi }_i⟩}{\partial x_j}})\hfill \\ \hfill & =0+0+V_j({\delta }_{\mathit{ij}}+0),\hfill \end{array}$$(B.11)

or in other words

$$\begin{array}{c}\hfill \mathbf{V}={⟨\mathbf{u}⟩}_L.\end{array}$$(B.12)

To summarise, the identity (B.8) is always verified, but in general, the velocity V appearing in the definition of the operator ⟨D⟩_L is not easily specified. Only when the displacement function ξ is judiciously defined as a fluctuating particle displacement does the velocity V reduce to the mean Lagrangian velocity ⟨u⟩_L.

B.2. Derivation of Eqs. (16) and (17)

Let us apply Eq. (B.8) to ϕ = x_i and ϕ = u_i alternatively. First, if ϕ = x_i, then ϕ_L = X_i, and Eq. (B.8) becomes

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{\partial x_i}{\partial t}}\right)}_L+u_{j,L}{\left({\displaystyle \frac{\partial x_i}{\partial x_j}}\right)}_L={\displaystyle \frac{\partial X_i}{\partial t}}+{⟨u_j⟩}_L{\displaystyle \frac{\partial X_i}{\partial x_j}}.\end{array}$$(B.13)