Issue 
A&A
Volume 656, December 2021



Article Number  A95  
Number of page(s)  18  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202141483  
Published online  08 December 2021 
Coupling between turbulence and solarlike oscillations: A combined Lagrangian PDF/SPH approach
I. The stochastic wave equation
LESIA, Observatoire de Paris, PSL Research University, CNRS, Université Pierre et Marie Curie, Université Paris Diderot, 92195 Meudon, France
email: jordan.philidet@obspm.fr
Received:
7
June
2021
Accepted:
3
September
2021
Context. The development of spaceborne missions such as CoRoT and Kepler now provides us with numerous and precise asteroseismic measurements that allow us to put better constraints on our theoretical knowledge of the physics of stellar interiors. In order to utilise the full potential of these measurements, however, we need a better theoretical understanding of the coupling between stellar oscillations and turbulent convection.
Aims. The aim of this series of papers is to build a new formalism specifically tailored to study the impact of turbulence on the global modes of oscillation in solarlike stars. In building this formalism, we circumvent some fundamental limitations inherent to the more traditional approaches, in particular the need for separate equations for turbulence and oscillations, and the reduction of the turbulent cascade to a unique length and timescale. In this first paper we derive a linear wave equation that directly and consistently contains the turbulence as an input to the model, and therefore naturally contains the information on the coupling between the turbulence and the modes through the stochasticity of the equations.
Methods. We use a Lagrangian stochastic model of turbulence based on probability density function methods to describe the evolution of the properties of individual fluid particles through stochastic differential equations. We then transcribe these stochastic differential equations from a Lagrangian frame to a Eulerian frame more adapted to the analysis of stellar oscillations. We combine this method with smoothed particle hydrodynamics, where all the mean fields appearing in the Lagrangian stochastic model are estimated directly from the set of fluid particles themselves, through the use of a weighting kernel function allowing to filter the particles present in any given vicinity. The resulting stochastic differential equations on Eulerian variables are then linearised. As a first step the gas is considered to follow a polytropic relation, and the turbulence is assumed anelastic.
Results. We obtain a stochastic linear wave equation governing the time evolution of the relevant wave variables, while at the same time containing the effect of turbulence. The wave equation generalises the classical, unperturbed propagation of acoustic waves in a stratified medium (which reduces to the exact deterministic wave equation in the absence of turbulence) to a form that, by construction, accounts for the impact of turbulence on the mode in a consistent way. The effect of turbulence consists of a nonhomogeneous forcing term, responsible for the stochastic driving of the mode, and a stochastic perturbation to the homogeneous part of the wave equation, responsible for both the damping of the mode and the modal surface effects.
Conclusions. The stochastic wave equation obtained here represents our baseline framework to properly infer properties of turbulenceoscillation coupling, and can therefore be used to constrain the properties of the turbulence itself with the help of asteroseismic observations. This will be the subject of the rest of the papers in this series.
Key words: methods: analytical / stars: oscillations / stars: solartype / turbulence
© J. Philidet et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Solarlike 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 pmode frequencies (Dziembowski et al. 1988; ChristensenDalsgaard et al. 1996; Rosenthal et al. 1999). The variety of physical processes responsible for the impact of turbulent convection on pmode 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; ChristensenDalsgaard 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Γ_{1} (GGM) or reducedΓ_{1} (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 timedependent mixinglength 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 mixinglength 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 cutoff 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 solartype 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 solarlike oscillations. Solarlike pmodes 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 timedependent nonadiabatic turbulent convection model able to include the oscillations. Subsequently, modelling attempts have focused on the use of mixinglength 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 mixinglength 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 timedependent mixinglength 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 higherorder 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 turbulenceinduced 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 FokkerPlanck equation (Gardiner 1994). Because the FokkerPlanck 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 turbulenceinduced 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 turbulenceoscillation 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 secondorder 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 secondorder level is replaced by equations on secondorder moments, where thirdorder 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 FokkerPlanck 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, gridbased 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.
Glossary of the notations used in this paper.
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^{1} 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
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)
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
where δ_{ij} is the Kronecker symbol and the notation 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 phasespace 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 righthand 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, jcomponent 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)
where C_{0} 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_{0} 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_{0} = 2.1 (Haworth & Pope 1986).
For the drift term we adopt the general expression given by the generalised Langevin model (Pope 1983)
where , , and are the Reynolds average of the fluid density, gas pressure, and gravitational acceleration, respectively; G_{ij} is a secondorder tensor that has the dimension of an inverse time, to which we refer as the drift tensor; and is the Favre average of the fluid velocity, with the massaverage (or Favre average) of any quantity ϕ being defined as
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)
where 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
The mean fields , , , 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^{2} 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 timeindependent average position is denoted by x. Its instantaneous position at time t is written as an explicit function of x and t
where ξ is the particle displacement around its mean position^{3}, the mean position being interpreted as an Eulerian variable.
For any given Eulerian quantity ϕ, we define its Lagrangian counterpart as
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
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 yettobedetermined 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)
where D/Dt ≡ ∂_{t} + u ⋅ ∇ denotes the particle time derivative, and the operator ⟨D⟩_{L} is defined by
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
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
so that
Putting together Eqs. (10), (17), and (21) we obtain
By construction, η(x, t)≡dW/dt is a multivariate Gaussian process whose values at two distinct locations are completely uncorrelated, and which verifies^{4}
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 turbulenceinduced 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 timeindependent 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 u_{L} (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 V_{L} represents a thin tube of fluid initially lying along the horizontal direction. The fluid displacement and velocity ξ and u_{L} 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
It must be noted that while the velocity appearing on the righthand 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 lefthand 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 righthand 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 , gas pressure , and velocity ), as well as the turbulent dissipation ϵ appearing in the diffusion tensor, and the Reynolds stress tensor and shear tensor 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 largeeddy simulation (LES) (or direct numerical simulation) to simulate the largescale 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, is the Reynoldsaveraged 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 particlelevel 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 compactsupport vicinity of x.
Thus, a kernel function K(r) is introduced, which serves as a weighting function to implement the particleposition 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 compactsupport 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)
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^{3}) is defined by the normalisation condition^{5}, 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
where Δm^{(i)} is the mass carried by the particle i, and ρ^{(i)} is the mass density characterising particle i. As such, the quantity Δm^{(i)}/ρ^{(i)} appearing under the sum corresponds to the lumped volume of fluid that the particle represents.
Setting Q^{⋆} = ρ, Q^{⋆} = ρu^{⋆}, and alternatively in Eq. (28), we find respectively the Reynoldsaveraged density, the massaveraged velocity, and the Reynolds stress tensor
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 ≡ ρ_{0}(y)d^{3}y, where ρ_{0} 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
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 . For the sake of simplicity, we consider a polytropic relation between the gas pressure and density in the form
where p_{0} is the equilibrium gas pressure (defined in the same way as ρ_{0}), 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 γ = Γ_{1}, where Γ_{1} is the equilibrium first adiabatic exponent.
The second point concerns the turbulent dissipation rate ϵ, or equivalently the turbulent frequency ω_{t} defined through
where k is the turbulent kinetic energy. Physically, ω_{t} can be interpreted as the inverse of the characteristic lifetime associated with the energybearing 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 massaveraged dissipation rate , which is very similar to the approach followed in twoequation models of turbulence, such as the k − ϵ model (Jones & Launder 1972), or else by adding 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 timeindependent 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 , the massaveraged velocity , and the Reynolds stress tensor 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 ρ_{0}(x), gas pressure p_{0}(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
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 ReynoldsAveraged NavierStokes (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 solarlike oscillators, the typical turbulent velocities have much higher amplitudes than the oscillatory velocities; the former are of the order of a few km.s^{−1}, while the latter are of the order of a few tens of cm.s^{−1}. This allows us to treat u_{osc} as a firstorder perturbation compared to u_{t}, and any second or higherorder 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_{0}/dr)^{−1} 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 firstorder perturbation compared to all length scales relevant to the problem, and any second or higherorder occurrence of ξ_{osc} will be discarded.
(H3). We adopt the anelastic approximation for turbulence, in the sense that we consider ρ_{t} ≪ ρ_{0}, where ρ_{t} is the turbulent fluctuation of density, and ρ_{0} 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 ∇ ⋅ (ρ_{0}ξ_{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 backreaction 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 backreaction 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 Reynoldsaveraged 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_{0}, 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 . Plugging Eq. (37) into Eq. (32), we find
Then, using (H2) to linearise in terms of the displacement, we find
The first term on the righthand side corresponds to the kernel estimate of ρ_{0} at x. By construction, kernel estimation is a representation of ensemble averaging, but ρ_{0} is already an equilibrium quantity, and therefore is equal to its own ensemble average. Furthermore, the last term on the righthand side can be discarded on account of hypothesis (H3). Performing an integration by part makes the quantity ∇ ⋅ (ρ_{0}ξ_{t}) appear under the integral sign. Therefore, we eventually find
The quantity ρ_{1} represents the Eulerian modal fluctuations of density, but of important note is the fact that at no point did we explicitly decompose into an equilibrium value ρ_{0} 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
where
where is the equilibrium sound speed squared, and we have introduced the xcentred kernel function K^{x}(y)≡K(y − x). The righthand 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 nonstochastic zerothorder terms in the linearised equations vanish under the unperturbed hydrostatic equilibrium condition, and then that the linear forcing (i.e. the stochastic zerothorder 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 righthand 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. 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, contains all the terms that are linear in ξ_{osc} and u_{osc} and explicitly depend on ξ_{t}, u_{t}, or η. Finally, L_{0} 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 lefthand side of Eq. (43), together with its righthand 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_{0} 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 solarlike 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 nonlinear Lagrangian, turbulent fluctuations of entropy, which is widely recognised as another source of stochastic driving for solarlike pmodes, 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 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 nonlinear advection term in the momentum equation, as represented by the first two terms on its righthand side, or from the joint effect of turbulent dissipation, buoyancy, and pressurerateofstrain 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 thirdorder contributions alike. As a whole, the term 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 . 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
Ultimately, this leads to the following simplification of Eq. (44):
Hence from Eqs. (42) and (43), which in this limit read
we obtain the following homogeneous, secondorder wave equation
with
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 backreaction 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 ∇ ⋅ (ρ_{0}ξ_{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 nonlocal. As we mention in Sect. 3.1, its domain of validity is the highradialorder 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 nonadiabatic 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 nonadiabatic effects in the model. Essentially, adopting a nonadiabatic 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 turbulenceinduced 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 energycontaining 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 nonzero ⟨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 solarlike 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 turbulenceinduced 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 selfconsistent way possible, the terms arising from the turbulenceoscillation 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 turbulentinduced 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 oscillationindependent inputs to the model.

All aspects of turbulenceoscillation 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 mixinglength hypothesis, which is crucial as this hypothesis is both almost inescapable in current convection modelling, and very invalid close to the radiativeconvective 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 nonradial 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 turbulenceoscillation coupling, and in particular explicitly and simultaneously yield the excitation rates of the modes, their lifetimes, as well as their turbulenceinduced frequency corrections. Finally, as we mentioned in Sect. 4, extending this work to a nonadiabatic 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.
In the vocabulary of asteroseismology, the term ‘adiabatic’ can sometimes be used to express the absence of energy transfer between the oscillations and the background. We insist that this is not the case here. This term is meant to apply to the thermodynamic transformations undergone by the flow, not to the oscillations. In particular, the background can still inject energy into the modes or take energy from them, allowing the modes to be driven and damped.
This statement may seem odd, as Lagrangian variables are actually often used in the analysis of stellar oscillations. However, in this study the term Lagrangian refers to a frame of reference attached to the total velocity of the flow (including both the turbulent velocity and the oscillation velocity), while the usual sense is rather meant to describe a frame of reference attached to the oscillations alone, and actually only ever refers to a pseudoLagrangian frame.
The stochastic process W(t) is not defined as an ordinary function of time, and therefore its derivative η cannot be defined the classical way; in fact, W is nowhere differentiable, as can be seen from Eq. (4). However, η can be defined formally, with its statistical properties given in the sense of distributions rather than ordinary functions. These definitions are at the heart of Ito stochastic calculus (Gardiner 1994, Chap 4).
The reason the value of c given here is different from the value c = 4/(5h) given in Welton & Pope (1997) is that they considered the 1D case, whereas we consider the 3D case.
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
 Anand, M. S., Mongia, H. C., & Pope, S. B. 1989, A PDF method for turbulent recirculating flows, 672 [Google Scholar]
 Andrews, D. G., & McIntyre, M. E. 1978, J. Fluid Mech., 89, 609 [NASA ADS] [CrossRef] [Google Scholar]
 Bakosi, J., & Ristorcelli, J. R. 2011, J. Turbul., 12, 19 [Google Scholar]
 Ball, W. H., & Gizon, L. 2014, A&A, 568, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Balmforth, N. J. 1992a, MNRAS, 255, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Balmforth, N. J. 1992b, MNRAS, 255, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Balmforth, N. J. 1992c, MNRAS, 255, 639 [NASA ADS] [Google Scholar]
 Baudin, F., Barban, C., Belkacem, K., et al. 2011, A&A, 529, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., Kupka, F., & Baudin, F. 2006, A&A, 460, 183 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., & Dupret, M. A. 2008, A&A, 478, 163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Dupret, M. A., & Noels, A. 2010, A&A, 510, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Dupret, M. A., Baudin, F., et al. 2012, A&A, 540, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Kupka, F., Samadi, R., & GrimmStrele, H. 2019, A&A, 625, A20 [EDP Sciences] [Google Scholar]
 Belkacem, K., Kupka, F., Philidet, J., & Samadi, R. 2021, A&A, 646, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bühler, O. 2009, Waves and Mean Flows [CrossRef] [Google Scholar]
 Chaplin, W. J., Houdek, G., Elsworth, Y., et al. 2005, MNRAS, 360, 859 [NASA ADS] [CrossRef] [Google Scholar]
 Chaplin, W. J., Houdek, G., Karoff, C., Elsworth, Y., & New, R. 2009, A&A, 500, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ChristensenDalsgaard, J. 2012, Astron. Nachr., 333, 914 [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [Google Scholar]
 Das, S. K., & Durbin, P. A. 2005, Phys. Fluids, 17, 025109 [CrossRef] [Google Scholar]
 Delarue, B. J., & Pope, S. B. 1997, Phys. Fluids, 9, 2704 [CrossRef] [Google Scholar]
 Dupret, M. A., Barban, C., Goupil, M. J., et al. 2006, in Proceedings of SOHO 18/GONG 2006/HELAS I, Beyond the spherical Sun, eds. K. Fletcher, & M. Thompson, ESA Spec. Publ., 624, 97 [Google Scholar]
 Dziembowski, W. A., Paterno, L., & Ventura, R. 1988, A&A, 200, 213 [NASA ADS] [Google Scholar]
 Gabriel, M., Scuflaire, R., Noels, A., & Boury, A. 1975, A&A, 40, 33 [Google Scholar]
 Gardiner, C. W. 1994, Handbook of stochastic methods for physics, chemistry and the natural sciences [Google Scholar]
 Goldreich, P., & Keeley, D. A. 1977a, ApJ, 211, 934 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Keeley, D. A. 1977b, ApJ, 212, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Kumar, P. 1991, ApJ, 374, 366 [NASA ADS] [CrossRef] [Google Scholar]
 Gough, D. O. 1977, ApJ, 214, 196 [NASA ADS] [CrossRef] [Google Scholar]
 Grigahcène, A., Dupret, M. A., Gabriel, M., Garrido, R., & Scuflaire, R. 2005, A&A, 434, 1055 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haworth, D. C., & El Tahry, S. H. 1991, AIAA J., 29, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Haworth, D. C., & Pope, S. B. 1986, Phys. Fluids, 29, 387 [CrossRef] [Google Scholar]
 Heinz, S., & Buckingham, A. 2004, Appl. Mech. Rev., 57, B28 [NASA ADS] [CrossRef] [Google Scholar]
 Houdek, G. 1996, PhD Thesis [Google Scholar]
 Houdek, G., & Dupret, M.A. 2015, Liv. Rev. Sol. Phys., 12, 8 [Google Scholar]
 Houdek, G., Balmforth, N. J., ChristensenDalsgaard, J., & Gough, D. O. 1999, A&A, 351, 582 [NASA ADS] [Google Scholar]
 Houdek, G., Trampedach, R., Aarslev, M. J., & ChristensenDalsgaard, J. 2017, MNRAS, 464, L124 [Google Scholar]
 Houdek, G., Lund, M. N., Trampedach, R., et al. 2019, MNRAS, 487, 595 [Google Scholar]
 Hsu, A. T., Tsai, Y. L. P., & Raju, M. S. 1994, AIAA J., 32, 1407 [CrossRef] [Google Scholar]
 Jones, W., & Launder, B. 1972, Int. J. Heat Mass Transfer, 15, 301 [CrossRef] [Google Scholar]
 Kjeldsen, H., Bedding, T. R., & ChristensenDalsgaard, J. 2008, ApJ, 683, L175 [Google Scholar]
 Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301 [Google Scholar]
 Kolmogorov, A. N. 1962, J. Fluid Mech., 13, 82 [Google Scholar]
 Liu, M. B., & Liu, G. R. 2010, Arch. Comput. Methods Eng., 17, 25 [CrossRef] [Google Scholar]
 Monaghan, J. J. 1992, ARA&A, 30, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, Å., & Stein, R. F. 2001, ApJ, 546, 576 [Google Scholar]
 Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 2 [Google Scholar]
 O’Brien, E. E. 1980, in The Probability Density Function (pdf) Approach to Reacting Turbulent Flows, eds. P. A. Libby, & F. A. Williams, 44, 185 [CrossRef] [Google Scholar]
 Obukhov, A. M. 1959, Adv. Geophys., 6, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Pope, S. B. 1981, Phys. Fluids, 24, 588 [CrossRef] [Google Scholar]
 Pope, S. B. 1983, Phys.f Fluids, 26, 3448 [NASA ADS] [CrossRef] [Google Scholar]
 Pope, S. B. 1985, Prog. Energy Combust. Sci., 11, 119 [Google Scholar]
 Pope, S. B. 1991, Phys. Fluids A, 3, 1947 [NASA ADS] [CrossRef] [Google Scholar]
 Pope, S. B. 1994, Ann. Rev. Fluid Mech., 26, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Pope, S. B. 2000, Turbulent Flows [Google Scholar]
 Pope, S. B., & Chen, Y. L. 1990, Phys. Fluids A, 2, 1437 [NASA ADS] [CrossRef] [Google Scholar]
 Roekaerts, D. 1991, Appl. Sci. Res., 48, 271 [CrossRef] [Google Scholar]
 Rosenthal, C. S., ChristensenDalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689 [NASA ADS] [Google Scholar]
 Samadi, R., & Goupil, M. J. 2001, A&A, 370, 136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Goupil, M. J., Alecian, E., et al. 2005, JApA, 26, 171 [NASA ADS] [Google Scholar]
 Samadi, R., Kupka, F., Goupil, M. J., Lebreton, Y., & van’t VeerMenneret, C. 2006, A&A, 445, 233 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Belkacem, K., & Sonoi, T. 2015, EAS Publ. Ser., 7374, 111 [CrossRef] [EDP Sciences] [Google Scholar]
 Sonoi, T., Samadi, R., Belkacem, K., et al. 2015, A&A, 583, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sonoi, T., Belkacem, K., Dupret, M. A., et al. 2017, A&A, 600, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Soward, A. M. 1972. Philos. Trans. R. Soc. London Ser. A, 272, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2010, ARA&A, 48, 391 [NASA ADS] [CrossRef] [Google Scholar]
 Stratonovich, R. L. 1965, Topics in the Theory of Random Noise, Vol. I and II (New York: Gordon and Breach) [Google Scholar]
 Unno, W. 1967, PASJ, 19, 140 [NASA ADS] [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars [Google Scholar]
 Van Slooten, P. R., Jayesh, & Pope S. B. 1998, Phys. Fluids, 10, 246 [NASA ADS] [CrossRef] [Google Scholar]
 Welton, W. C. 1998, J. Comput. Phys., 139, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Welton, W. C., & Pope, S. B. 1997, J. Comput. Phys., 134, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Xiong, D. R., Cheng, Q. L., & Deng, L. 2000, MNRAS, 319, 1079 [CrossRef] [Google Scholar]
 Zhou, Y., Asplund, M., Collet, R., & Joyce, M. 2020, MNRAS, 495, 4904 [CrossRef] [Google Scholar]
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 secondorder 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
and
where δ_{ij} denotes the Kronecker symbol, and we have introduced the pseudoLagrangian particle derivative .
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 lefthand 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 NavierStokes equation, which reads
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 lefthand side), but also the production term (the first three terms on the righthand side). However, the last three terms are not modelled exactly, and are (from left to right) the buoyancy contribution, the pressurerateofstrain 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
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:
Furthermore, the drift tensor is usually decomposed into an isotropic and anisotropic part, according to
where is the turbulent kinetic energy. This decomposition ensures that, in the special case of incompressible, homogeneous, isotropic turbulence, if we take , the evolution of the Reynolds stress tensor reduces to the exact, analytical solution.
Equation (A.5) can be rearranged to yield
Equation (A.8) allows us to interpret the collective effect of buoyancy and pressurerateofstrain correlation on the evolution of the Reynolds stresses. First, all the terms on the righthand 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 compressioninduced production term (i.e. the first three terms on the righthand side of Eq. (A.4)).
Furthermore, it is readily seen that the last term on the righthand side of Eq. (A.8) tends to isotropise the Reynolds stress tensor since for isotropic turbulence we would precisely have . The rate at which this term makes the Reynolds stress decay towards isotropy is equal to (1 + 3C_{0}/2)ω_{t}, where ω_{t} is the turbulent dissipation rate defined by Eq. (36). On the other hand, the other two terms on the righthand 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 LagrangiantoEulerian 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
where ϕ is an arbitrary quantity. The usual chain rules for derivation then yield
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
Plugging Eq. (B.5) into Eq. (B.3), we find
where D/Dt ≡ ∂_{t} + u_{i}∂_{i}. In turn, plugging Eq. (B.4), this transforms into
thus yielding the required identity
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^{6} position, and that the function ξ(x, t) actually denotes the fluctuating fluid displacement around this mean position x. Then by construction, we have
Point x now corresponding to a mean position, the velocity V at which it is displaced must itself be a mean quantity, so that
Let us now apply the mean operator ⟨.⟩ to Eq. (B.5); since V_{j} can be pulled out of the mean, we obtain
or in other words
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
The mean position x having no explicit time dependence, we have (∂x_{i}/∂t)_{L} = 0, (∂x_{i}/∂x_{j})_{L} = δ_{ij}, ∂X_{i}/∂t = ∂ξ_{i}/∂t, and ∂X_{i}/∂x_{j} = δ_{ij} + ∂ξ_{i}/∂x_{j}. The above equation then becomes
thus immediately yielding Eq. (16).
Secondly, if ϕ = u_{i}, then ϕ_{L} = u_{i, L}, and Eq. (B.8) becomes
But Eq. (B.4) allows us to write
where we recall that
Plugging these into Eq. (B.15), we find
Isolating the first term on the righthand side yields Eq. (17).
Appendix C: The insignificance of the backreaction of the oscillations on the turbulence
Let us formally writ e the governing equations of the flow in the following abstract form
where U represents the flow variables, ℒ is a linear operator, and ℬ a bilinear operator containing the advection terms. In the limit of small amplitudes, which are relevant for solarlike oscillations, the wave variables can be expanded as
where a is small ordering parameter. Plugging Eq. (C.2) into Eq. (C.1) and isolating the various orders in a, we obtain the following hierarchy of equations
where the first equation governs the basic flow, the second equation governs the waves, and the third equation governs the backreaction of the waves on the basic flow. In particular, Eq. (C.5) takes the form of a forced linear wave, where the linear part ℒ′≡ℒ + ℬ(U_{0}, .)+ℬ(., U_{0}) is identical to the linear part in the actual wave equation (C.4), and the forcing term is given by the righthand side of Eq. (C.5). Because ℒ′ is common to both Eqs. (C.4) and (C.5), if we denote the angular frequency of the wave as ω, we can write the homogeneous solution of Eq. (C.5) as
and the total solution (including the forcing, inhomogeneous part) formally reads
The backreaction of the waves on the turbulence is therefore driven by its resonance with the nonlinear oscillationinduced advection. In the limit t → +∞, this formal solution yields
where ‘TF’ denotes the Fourier transform. But U_{1} refers to the wave, so that its Fourier spectrum only has power around the angular frequency ω of the wave. In turn, this means that the quadratic operator ℬ applied to the velocity U_{1} has a Fourier spectrum whose power is concentrated around ω = 0 (i.e. the continuous component) as well as 2ω (twice the frequency of the waves). By contrast, it contains little to no power around the actual frequency ω of the oscillation, which justifies that the impact of the backreaction U_{2}(t) on the mean flow U_{0}(t) may be neglected.
Appendix D: Derivation of the linear wave equation
In this appendix we linearise the system comprised of Eqs. (25), (26), (32), (33), (34), (35), and (36), using the hypotheses outlined in Section 3.1. We start, in Section D.1, by linearising all the ensemble averages described in the SPH formalism (i.e. Eqs. (32), (33), (34), (35), and (36)). In Section D.2, we then plug these linearised ensemble averages to derive the linearised version of Eqs. (25) and (26). Finally, in Section D.3, we discuss which terms should be retained in the inhomogeneous forcing term of the resulting wave equation. For more clarity in the notations, we dropped all dependence on the space variable x, the space variable y used inside the integrals, and time t. It must be understood that all the quantities outside the integrals depend on x and t, and all quantities inside depend on y and t.
D.1. Linearising the mean fields
A general remark can be made beforehand concerning all ensemble averages described in the SPH formalism: the occurrence of ξ_{t} vanishes completely from their linearised version by virtue of hypothesis (H3). We have already shown, in the main body of this paper, that this is the case for the mean density , but this is also the case for the mean velocity and Reynolds stress tensor. They can both formally be written as
where Q is a function of velocity only (Q = u for the mean velocity, and for the Reynolds stress tensor), and we have introduced the xcentred kernel function K^{x}(y)≡K(y − x). Because Q only depends on the velocity variable u, and not on the displacement variable ξ, the only occurrence of ξ_{t} in the linearisation of comes from the term
Performing an integration by part yields
where the surface term vanishes because of the compact support of the kernel function K^{x}. By virtue of hypothesis (H3), the quantity ∇ ⋅ (ρ_{0}ξ_{t}) is negligible, and therefore this contribution can be safely discarded.
As we have just shown, this is true of the mean density, mean velocity, and Reynolds stress tensor. In turn, this is also true of the gas pressure (because it is given as a function of the mean density), as well as the turbulent kinetic energy k and the turbulent dissipation rate ϵ (because they are both given as a function of the Reynolds stress tensor). Therefore, ξ_{t} can indeed be neglected in the linearised version of every single ensemble average appearing in Eqs. (25) and (26).
D.1.1 Mean density
Using hypotheses (H2) and (H3), Eq. (32) can be linearised as
The first term on the righthand side corresponds to the kernel estimator of the equilibrium density ρ_{0} at x, and therefore represents the ensemble average of ρ_{0} at x. Since ρ_{0} is already an equilibrium quantity, it is equal to its own ensemble average, and this term reduces to ρ_{0}(x) itself. Finally,
with
D.1.2. Mean gas pressure
The fluctuating mean density ρ_{1} is much lower than the equilibrium density ρ_{0} on account of hypothesis (H2). Therefore, Eq. (35) can be linearised, immediately yielding
with
where is the equilibrium sound speed squared.
D.1.3. Mean velocity
Using hypotheses (H1), (H2), and (H3), Eq. (33) can be linearised as
where ρ_{1} is given by Eq. (D.6). This expression can be simplified by remarking that
because kernel averages represent ensemble averages. Since mass is locally conserved by the turbulent velocity field (i.e. the upflows carry as much mass upwards as the turbulent downdrafts carry downwards), it immediately follows that , and therefore the first two terms in Eq. (D.9) vanish. Rearranging the remaining terms, we obtain
where and
D.1.4. Mean shear tensor
We also need to express the linearised shear tensor because this quantity appears in the drift tensor G_{ij}. Differentiating Eq. (D.12) with respect to x_{i}, and noting that ∇_{x}(K^{x}(y)) = −∇_{y}(K^{x}(y)) (because we considered an isotropic kernel function), we obtain
where and
D.1.5. Reynolds stress tensor
The linearised Reynolds stress tensor is obtain from Eq. (34), using hypotheses (H1), (H2), and (H3). We find
where ρ_{1} is given by Eq. (D.6) and by Eq. (D.12). In the last two terms, can be pulled from the integral. The kernel estimator is a representation of ensemble averages, and is already an ensemble average. Once this quantity is pulled out, we recognise the same integral defined by Eq. (D.10), meaning that these terms vanish. Additionally, the third, fourth, and fifth terms can be conveniently merged together, and ρ_{1} can be replaced by its explicit expression (D.6), so that we finally obtain
where
and
Additionally, we immediately deduce the linearisation of the turbulent kinetic energy, which corresponds to half the trace of the Reynolds stress tensor
Likewise, from Eq. (36), we find the linearised turbulent dissipation
D.1.6. Drift tensor
The drift tensor G_{ij} being dependent on the mean flow, it also needs to be linearised. We recall that, in its most general form, it can be written as an arbitrary function of the Reynolds stress tensor, the shear tensor and the turbulent dissipation
and therefore its linearisation reads
with
and
In Eq. (D.23), is given by Eq. (D.17) and ϵ_{0} by Eq. (D.20). In Eq. (D.24), the derivatives , , and ∂G_{ij}/∂ϵ only depend on the functional form of the drift tensor; is given by Eq. (D.18); is given by Eq. (D.14); and ϵ_{1} is given by Eq. (D.20).
D.2. Linearising the displacement and motion equations
Hypotheses (H1) and (H2), in addition to the linearised mean fields computed in the previous section, allow us to write the linearised version of Eqs. (25) and (26) as
and
where ρ_{0} and p_{0} are the equilibrium density and gas pressure, ρ_{1} is given by Eq. (D.6), p_{1} by Eq. (D.8), G_{ij, 0} by Eq. (D.23), G_{ij, 1} by Eq. (D.24), by Eq. (D.12), and k_{0} and k_{1} by (D.19).
Furthermore, we split the righthand side of Eq. (D.26) three ways: we gather all the terms that do not depend on the oscillatory variables u_{osc} and ξ_{osc} in a quantity L_{0}, all the terms that depend on ξ_{osc} and/or u_{osc} but not on any of the turbulent fields ξ_{t} or u_{t} in a quantity , and all the terms that depend on both the oscillatory variables and the turbulent fields in a quantity . This leads us to the following linear equations
where
where is the equilibrium sound speed squared, and the quantities , , and k_{1} are given by Eqs. (D.18), (D.14), and Eq. (D.19) respectively.
D.3. The forcing term
In Eqs. (D.27) and (D.28) the righthand side represent inhomogeneous stochastic forcing terms. For the moment we have kept all zeroth order terms (i.e. all the terms that are independent of the wave variables ξ_{osc} and u_{osc}) on these righthand sides, but L_{0} can be rearranged into a more compact form, and some terms will in fact prove negligible. Firstly, let us rewrite the first term on the righthand side of Eq. (D.31). Using hypothesis (H4), we write the continuity equation for u_{t} without any contribution from the oscillatory component
where ρ is the sum of the equilibrium value ρ_{0} and the turbulent fluctuations of the density ρ_{t}. Building on hypothesis (H3), we neglect ρ_{t} in Eq. (D.32), so that
finally allowing us to write
Secondly, it does not come as a surprise that the nonstochastic part of Eq. (D.31) corresponds to the hydrostatic equilibrium condition. If radiative pressure is neglected, we have
so that
It thus becomes clear that the forcing term contains the usual contribution from the fluctuations of the turbulent pressure, a contribution that is linear in u_{t}, and a contribution that is linear in η, and therefore completely uncorrelated in space. Following the discussion from Samadi & Goupil (2001), we argue that all linear contributions are negligible. The contribution of a linear term to the excitation rate of the modes has an efficiency that is based on the resonance between the lifetime of the largescale energybearing eddies and the period of the modes, which the authors showed was negligible. Naturally, the same argument can be used to neglect the third term as well since it has no coherence in either space or time. The nonlinear term, on the other hand, is able to couple different length scales together, and therefore leads to a nonnegligible contribution to the excitation rate. Finally, after having filtered out those terms we deemed negligible, we obtain
The righthand side of Eq. (D.27) can be treated similarly: the term u_{t} being linear in the turbulent fields, its contribution to mode driving can be neglected, thus only leaving the second term.
All Tables
All Figures
Fig. 1. Visualisation of the fluid displacement ξ (in green), Lagrangian velocity u_{L} (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 V_{L} represents a thin tube of fluid initially lying along the horizontal direction. The fluid displacement and velocity ξ and u_{L} 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 the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.