# A new particle-based code for Lagrangian stochastic models applied to stellar turbulent convection

J. Philidet, K. Belkacem

► To cite this version:

J. Philidet, K. Belkacem. A new particle-based code for Lagrangian stochastic models applied to stellar turbulent convection. Astronomy & Astrophysics - A&A, 2025, 10.1051/0004-6361/202556625. hal-05383066

HAL Id: hal-05383066

https://hal.science/hal-05383066v1

Submitted on 26 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A new particle-based code for Lagrangian stochastic models applied to stellar turbulent convection

J. Philidet¹ and K. Belkacem¹

LIRA, Observatoire de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CY Cergy Paris Université, CNRS, 92190 Meudon, France

Received 29 July 2025 / Accepted 26 October 2025

## ABSTRACT

**Context.** The inclusion of convection in stellar evolution models, mostly based on Mixing Length Theories, lacks realism, especially near convective-radiative interfaces. Furthermore, the interaction of convection with oscillations is poorly understood, giving rise to surface effects that currently prevent us from accurately predicting seismic frequencies, and therefore from fully exploiting the asteroseismic data of low-mass stars.

**Aims.** We aim to develop a new formalism to model the one-point statistics of stellar convection, to implement it in a new numerical code, and to validate this implementation against benchmark cases.

**Methods.** This new formalism is based on Lagrangian Probability Density Function (PDF) methods, where a Fokker-Planck equation for the PDF of particle-based turbulent properties is integrated in time. The PDF equation is established so that the underlying transport equations for all first- and second-order moments of the turbulent flow are identical to the exact ones stemming from first principles. We then develop a Monte-Carlo implementation of this method, where the flow is represented by a large number of notional particles acting as realisations of the PDF. Notional particles interact with each other through the time- and space-dependent mean flow, which is estimated from the particle realisations through a scheme similar to Smoothed Particle Hydrodynamics.

**Results.** We establish a model for the evolution of turbulent properties along Lagrangian trajectories applicable to stellar turbulent convection, with only a minimal number of physical assumptions necessary to close the system. In particular, no closure is needed for the non-linear advection terms, which are included exactly through the Lagrangian nature of formalism. The numerical implementation of this new formalism allows us to extract time-dependent maps of the statistical properties of turbulent convection in a way which is not possible in grid-based large-eddy simulations, in particular the turbulent pressure, Reynolds stress tensor, internal energy variance and convective flux.

**Key words.** stars: interiors – convection – turbulence

## 1. Introduction

Convection is ubiquitous in stellar interiors. It causes large-scale motions that efficiently mix energy, angular momentum and chemical elements within stars. As such, it plays a key role in shaping the structure and evolution of the stellar interiors (Böhm-Vitense 1992).

For low-mass stars exhibiting stochastically excited oscillations, the modelling of the surface is particularly crucial for inferring stellar properties through asteroseismology. Indeed, the characterisation of asteroseismic targets is heavily reliant on stellar evolution models, whose shortcomings lead to biases in predicted seismic frequencies, and therefore biases in radius, mass and age determination (e.g. Lebreton et al. 2014a,b). Currently, most of the limitations come from our inability to correctly model the very turbulent superficial layers of low-mass stars harbouring a convective envelope. These surface effects (e.g. Christensen-Dalsgaard & Thompson 1997; Rosenthal et al. 1999) force us to only consider surface-independent seismic indicators, such as frequency separation ratios (Roxburgh & Vorontsov 2003; Otí Floranes et al. 2005), and therefore prevents us from exploiting all the information contained in the observed frequencies. This is particularly limiting in a time when asteroseismic observations, thanks to missions like CoRoT (Baglin et al. 2006b,a), *Kepler* (Borucki et al. 2010),

or the future mission PLATO (Rauer et al. 2025), have reached a level of precision that clearly highlights the limitations of theory.

The inclusion of convective transport in stellar models invariably relies on Mixing Length Theories (MLT), the basics of which have been devised by Böhm-Vitense (1958) (see also Kippenhahn et al. 1967). This description is founded on an overly simplified picture of convection, where convective elements are displaced upwards over a unique travel distance before depositing all their entropy in the surrounding medium. By contrast, 3D hydrodynamic simulations show the non-linear, multi-scale nature of stellar convection, not captured by Mixing-Length Theories (see Kupka & Muthsam 2017, for a comprehensive review). The diffusive hypothesis behind MLT is at odds with the advective and non-local nature of actual convection. Standard MLT, in its current implementation in stellar evolution codes, also has the greatest difficulties accounting for radiative/convective interfaces: local models of convection cannot describe the penetration of convective plumes into neighbouring stably stratified radiative regions (Zahn 1991). This has important implications: for example, the uncertainty on the age determination of stars with convective cores is dominated by our inability to quantify the increase in effective core size brought about by convective penetration. Prescriptions for angular momentum transport across interface layers are also challenging to derive, and the inadequacy of those implemented in stellar evo-
lution codes is partly responsible for the overestimation of the spin-up incurred by the core of evolved stars as they ascend the Red Giant Branch (Mosser et al. 2012).

Standard MLT is also ill-suited to describe surface convection in stars harbouring convective envelopes. Because the convective, thermal and dynamical timescales all coincide in the uppermost layers of these stars, the description of convection must account for the time-dependence of the properties of the medium from which the convective elements arise. Many attempts to incorporate a time dependence in MLT have been made (e.g. Unno 1967; Gough 1977; Balmforth 1992a; Grigahcène et al. 2005), and subsequently exploited to try to model asteroseismic surface effects (Balmforth 1992b; Grigahcène et al. 2012; Houdek et al. 2017; Sonoi et al. 2017). So far, however, the success of this approach has been limited. Several mechanisms are at play in surface effects, including for example the mechanical work exerted by turbulent pressure on the oscillations, or the wave-induced modulation of the heat exchange stemming from turbulent dissipation. Different time-dependent mixing-length formalisms have provided conflicting results regarding the nature of the dominant mechanism, even in the case of the Sun (see Belkacem & Samadi 2013; Houdek & Dupret 2015; Samadi et al. 2015, for reviews). Furthermore, the free parameters of time-dependent MLT are hard to relate to the physics of the stellar interior, because of the phenomenological nature of the formalism. As a result, agreement between the modelled and observed surface effect in the Sun can only be obtained at the price of tuning these parameters (e.g. Belkacem et al. 2012; Houdek et al. 2017, 2019). Extensions of MLT to include higher-order moments of the flow variables have been developed (Canuto 1997; Xiong et al. 1997), and subsequently applied to the stellar case to investigate the effect of convection on pulsation stability (Xiong et al. 2000), but this approach suffers essentially from the same limitations, and did not allow to solve the surface effect problem. We note however the recent promising results obtained by (Kupka et al. 2022; Ahlborn et al. 2022) for treating interfaces in the limit of efficient convection.

On the other hand, 3D Large-Eddy Simulations (LES) have been exploited with success to refine the theoretical description of stellar convection (see Kupka & Muthsam 2017). In particular, many works have been devoted to exploiting 3D LES to investigate how convection impacts oscillations. Belkacem et al. (2019) have studied the damping of acoustic modes by convection, and shown that 3D LES could help disentangle the importance of the various mechanism by which the damping can occur. Samadi & Goupil (2001); Samadi et al. (2001) have developed a formalism for the excitation of solar-like oscillations by surface convection, and used hydrodynamical simulations to show that the formalism explains the observed mode amplitudes. Conversely, Zhou et al. (2019) have used hydrodynamical simulations to assess the efficiency of acoustic mode excitation by convection throughout the Hertzsprung-Russell (HR) diagram. Also, the matter of surface effects on the frequencies have been investigated with such simulations in the case of the Sun by Schou & Birch (2020). While their approach proves very promising, it suffers from limitations as regards its use to model surface effects for a large variety of stars. For instance, the limited spatial extent of these simulations means that only a few modes are present, which complicates the study of how the different contributions to surface effects depend on frequency. The prohibitive computational cost of these simulations also prevents the creation of extensive grids of models, spanning a wide range of masses or effective temperatures. A more fundamental problem, however, is that 3D LES cannot correctly include turbulent

dissipation of energy into heat (e.g. Kupka & Muthsam 2017). Indeed, the very high Reynolds numbers characterising stellar turbulent convection means that all relevant scales cannot be resolved in 3D hydrodynamical simulations. As a result, the effect of the small sub-grid scales on the large-scale motions is included in a drastically simplified manner, either by letting the numerical viscosity of the spatial discretisation scheme dominate, or by artificially enhancing dissipative effects, for example with hyperviscosity (Borue & Orszag 1995). Either way, 3D LES fail to correctly include turbulent dissipation. This is particularly problematic for surface effects, for which it is expected to play a major role (e.g. Grigahcène et al. 2005; Belkacem et al. 2019, 2021).

In a context where hydrodynamical LES show unavoidable limitations in their ability to deal with turbulent closure, Probability Density Function (PDF) methods have emerged in the past decades as an alternative approach (e.g. Pope 2000). It consists in solving a transport equation for the one-point, one-time Eulerian joint PDF of all quantities describing the dynamical, thermodynamic, and possibly chemical state of a flow. Since their first developments (Dopazo 1975; Janicka et al. 1979; Pope 1981), PDF methods have found a considerable echo in the engineering fluid dynamics community, especially to model turbulent combustion (see Haworth 2010, for a thorough review), but also in atmospheric science (Rodean 1996). Since the seminal work of Pope (1979), the largely predominant approach in solving PDF transport equations has consisted in implementing particle methods, where the flow is viewed in a purely Lagrangian point of view, and discretised into fluid parcels. This has led to many developments to Lagrangian stochastic models, both for incompressible turbulence (Haworth & Pope 1987; Pope & Chen 1990; Pope 1994b; Dreeben & Pope 1998; Van Slooten & Jayesh 1998; Bakosi et al. 2008) and compressible turbulence (Delarue & Pope 1997, 1998; Das & Durbin 2005; Almeida & Navarro-Martinez 2021). These PDF methods offer strikingly compelling advantages for modelling stellar turbulent convection. In particular, the adopted Lagrangian point of view provides a natural way to include the advective transport of momentum and energy, exactly and without having to model it. This is a key advantage of Lagrangian PDF methods over traditional LES, because most of the difficulty in modelling stellar convection stems precisely from non-linear advective convection. While this does not entirely eliminate the need for turbulent closure, it nevertheless makes the closing procedure much easier. It also allows to include turbulent dissipation in a physical and controlled manner, by solving for the joint velocity-energy-dissipation PDF, something that LES do not allow to do.

Numerical solutions designed to implement PDF methods are predominantly based on Lagrangian-particle/Eulerian-mesh methods, where solving the particle equations requires knowledge of the means of the flow (density, velocity, etc.), and a Eulerian, mesh-based, finite-volume scheme is still needed to solve the (linear) transport equations for those means (Haworth 2010). However, such methods would not do in our case, because it would require modelling separately the large-scale mean flow and the small-scale turbulent fluctuations, and would therefore defeat one of the primary purposes of this work, namely studying the active impact of convection on waves. In this work, we thus make use of stand-alone particle methods instead, where the means of the flow are derived from the particle properties through filtered ensemble averages. Such an approach, based on kernel-estimation techniques, was developed by Welton & Pope (1997); Welton (1998), and presents the advantage of simultaneously and consistently modelling both the mean flow quanti-
ties and their turbulent fluctuations. We had previously shown (Philidet et al. 2021, 2022) how this approach can be exploited to establish a physical, semi-analytical model for the impact of convection on the main properties of the global acoustic modes of a star (its frequency, i.e. the surface effects, but also its damping rate and its excitation rate). However, this was done with an oversimplified model for the PDF transport equation.

The primary contributions of this work are two-fold. First, we establish a Lagrangian particle description adapted to stellar turbulent convection, in a way which exempts us from having to close the non-linear advection terms, from which stems most of the difficulty in devising MLT formulations. The analytical developments leading to this model are presented in Sect. 2. Second, we present a numerical implementation of this method based on a Monte-Carlo stand-alone particle algorithm, where a kernel estimation scheme is used to extract the mean flow quantities from the particles. This implementation is presented in Sect. 3, and validated against benchmark cases in Sect. 4. Finally, we showcase results pertaining to the statistics of stellar convection in Sect. 5, and we offer conclusions and perspectives in Sect. 6.

## 2. The Lagrangian stochastic model

The Monte-Carlo Lagrangian-particle method that we exploit to model stellar turbulent convection relies on the concept of “notional” particles, a large number of which can be used to describe realisations of the flow, and thus allowing to reconstruct the one-point, one-time joint PDF of all the flow variables. These notional particles evolve in time according to a set of stochastic differential equations (SDE), leading to Lagrangian stochastic models for turbulence (see Pope 1994a, 2000; Haworth 2010, for reviews and books on the subject). In this section, we establish the SDE describing turbulent convection in the stellar context, and outline the main physical assumptions involved in the derivation. Details of the derivation are given in App. A.

### 2.1. General framework

Since we wish to model the joint velocity-energy-dissipation PDF of the flow, each notional particle must be described by their position $\mathbf{x}^*$, velocity $\mathbf{u}^*$, internal energy $e^*$ and turbulent dissipation rate $\epsilon^*$ (the latter being replaced by the turbulent frequency $\omega^* \equiv \epsilon^*/k$, where $k$ is the turbulent kinetic energy). Each of these properties evolve in time according to SDEs that, in their most general form, can be expressed as (e.g. Pope 2000; Heinz 2003)

$$d\mathbf{x}^* = \mathbf{u}^* dt, \quad (1)$$

$$d\mathbf{u}^* = \mathbf{a}_u (\mathbf{x}^*, t; \mathbf{u}^*, e^*, \omega^*) dt + b_u (\mathbf{x}^*, t; \mathbf{u}^*, e^*, \omega^*) d\mathbf{W}_u, \quad (2)$$

$$d\epsilon^* = a_e (\mathbf{x}^*, t; \mathbf{u}^*, e^*, \omega^*) dt + b_e (\mathbf{x}^*, t; \mathbf{u}^*, e^*, \omega^*) dW_e, \quad (3)$$

$$d\omega^* = a_\omega (\mathbf{x}^*, t; \mathbf{u}^*, e^*, \omega^*) dt + b_\omega (\mathbf{x}^*, t; \mathbf{u}^*, e^*, \omega^*) dW_\omega, \quad (4)$$

each of which has a deterministic part (the first term on the right-hand side, $a_X$ with $X = \{u, e, \omega\}$), called the drift term, and a stochastic part (the second term on the right-hand side, $b_X$), called the diffusion term. The quantities $dW_X$ are increments over the time step $dt$ of mutually independent Wiener processes, meaning that

$$\overline{W_i(t)} = 0, \\ \overline{W_i(t)W_j(t+dt)} = \delta(dt) \delta_{ij},$$

where $\delta$ is the Dirac distribution, $\delta_{ij}$ the Kronecker symbol, and the notation $\overline{X}$ refers to the ensemble average of $X$.

### 2.2. Velocity and energy equations

The derivation of the drift and diffusion coefficients is detailed in App. A. The main idea is ensure that the exact Fokker-Planck equation derived from the equations of hydrodynamics is equivalent to the Fokker-Planck equation derived from Eqs. 1, 2, 3 and 4. This procedure is relatively standard for obtaining Lagrangian stochastic models, and has been the subject of considerable efforts by the fluid dynamics community (Haworth 2010). In this work, we adapt this method to stellar convection, where two physical ingredients are particularly important, namely (i) the compressibility of the flow, and (ii) the convective flux which plays a central role in buoyancy-driven convection. The non-linear advection terms in the momentum and energy equations are modelled exactly, without having to be included explicitly by hand. This is one of the key advantages of this class of models over MLT. Other non-linear turbulent terms, however, do require to be closed. As described in Sect. A.1, a number of assumption are required to obtained a tractable model. The main hypotheses are as follows: (1) that the fluctuations of density $\rho'$, gas pressure $p'$ and specific internal energy $e''$ are related to each other through a polytropic relation given by Eq. A.29, characterised by a polytropic index $n$, (2) that the compressibility of the turbulent velocity is uncorrelated with other turbulent fluctuations, so that $\overline{p'\partial u_i'' / \partial x_i} = \overline{e''\partial u_i'' / \partial x_i} = \overline{p'e''\partial u_i'' / \partial x_i} = 0$, and (3) that the effect of the viscous stress tensor $\sigma_{ij}$ can be neglected, except in the term $\overline{\sigma'_{ik}\partial u_i'' / \partial x_k} = \overline{\rho\epsilon/3} \delta_{ij}$. Under these assumptions, the drift and diffusion coefficients in the velocity and energy equations read (see Sect. A.1)

$$a_{u,i} = G_{0i} + G_{ij} (\overline{u_j} - \overline{u_i}) + G_{ei} (e^* - \overline{e}), \quad (5)$$

$$a_e = K_0 + K_j (\overline{u_j}^* - \overline{u_j}) + K_e (e^* - \overline{e}), \quad (6)$$

$$b_u = (C_0 \omega^* k)^{1/2}, \quad (7)$$

$$b_e = (C_1 \omega^* k_e)^{1/2}, \quad (8)$$

where

$$G_{0i} = -\frac{1}{\bar{\rho}} \frac{\partial \bar{p}}{\partial x_i} + g_i,$$

$$G_{ij} = -\left(\frac{1}{2} + \frac{3}{4}C_0\right) \bar{\omega} \delta_{ij} - \frac{1}{2}C_0 (\omega^* - \bar{\omega}) k (C^{-1})_{ij},$$

$$G_{ei} = \frac{\Gamma_1 - 1}{n-1} \frac{1}{\bar{p}} \frac{\partial \bar{p}}{\partial x_i} - \frac{1}{2} C_0 (\omega^* - \bar{\omega}) k (C^{-1})_{i0},$$

$$K_0 = -\frac{\bar{p}}{\bar{\rho}} \frac{\partial \overline{\tilde{u}_i}}{\partial x_i} - \frac{\Gamma_1 - 1}{n-1} \left( \frac{1}{\bar{\rho}} \frac{\partial \overline{\rho F_i^{(e)}}}{\partial x_i} - \frac{F_i^{(e)}}{\bar{p}} \frac{\partial \overline{p}}{\partial x_i} \right) + \kappa \frac{\partial^2 \overline{e}}{\partial x_i \partial x_i} + \overline{\omega} k, \quad (9)$$

$$K_i = -\frac{1}{2} C_1 (\omega^* - \bar{\omega}) k_e (C^{-1})_{0i},$$

$$K_e = -\left( \frac{1}{2} C_1 \bar{\omega} + (\Gamma_1 - 1) \frac{\partial \overline{u_i}}{\partial x_i} + \frac{\overline{\omega} k}{\overline{e}(n-1)} - \frac{1}{2} C_1 (\omega^* - \bar{\omega}) k_e (C^{-1})_{00} \right),$$

the notation $\overline{\tilde{X}} \equiv \overline{\rho\tilde{X}}/\bar{\rho}$ refers to the density-weighted (or Favre) ensemble average of $X$, $\bar{\rho}$ and $\overline{p}$ are the Reynolds average of
the density and gas pressure, $g$ is the gravitational acceleration (considered constant and uniform), $\Gamma_1$ is the first adiabatic exponent, $n$ is the polytropic index, $\kappa$ is the radiative diffusivity, and the coefficients of the $4 \times 4$ covariance matrix $C$ are defined by Eq. A.53 (the index 0 refers to internal energy, and each of the indices $i = 1, 2, 3$ to the three components of the velocity). The mean pressure is given by the ensemble average of the ideal gas law

$$ \bar{p} = (\Gamma_1 - 1)\bar{\rho}\bar{e}. $$

We emphasise that all averages depend on time, as well as the position of the particle. We also introduced the second order moments characterising the velocity-energy turbulent fluctuations

$$ R_{ij} \equiv \overline{u_i'' u_j''}, \quad \text{where } \mathbf{u}'' \equiv \mathbf{u} - \bar{\mathbf{u}}, $$

$$ k \equiv R_{ii}/2, $$

$$ k_e \equiv \overline{e''^2}, \quad \text{where } e'' \equiv e - \tilde{e}, $$

$$ \mathbf{F}^{(e)} \equiv \overline{e'' \mathbf{u}''}, $$

which are respectively the Reynolds stress tensor, the turbulent kinetic energy, the internal energy variance, and the internal energy flux. Repeated indices correspond to Einstein summation.

The constant $C_0$ appearing in Eq. 7 is the Kolmogorov constant, the universality of which has been widely discussed and investigated (see Heinz 2003, for a discussion). It is now understood that the value that should be adopted for $C_0$ is model-dependent, and may only be estimated by comparing model results with DNS or experimental data obtained in similar conditions. For the model considered here, the value is $C_0 = 3.5$ (Pope 1991), and is found to be rather independent of the Reynolds number. Since the constant $C_1$ plays the same role as $C_0$, but for the specific internal energy, we adopted the same value $C_1 = 3.5$.

## 2.3. Turbulent frequency equation

A stochastic model is added for the turbulent frequency of the notional particles, which allows to include the effect of internal intermittency on the turbulent time scale of velocity and energy fluctuations. While early models were constructed which led to the experimentally observed log-normal distribution of turbulent dissipation (e.g. Pope & Chen 1990; Pope 1991), these are computationally expensive because of the long tails of the distribution. An alternative has been developed by Van Slooten & Jayesh (1998), where the coefficients in Eq. 4 are given by (see App. A.2)

$$ a_{\omega} = -\Omega(\omega^* - \tilde{\omega}) - \tilde{\omega}\omega^* S_{\omega}, \quad (10) $$

$$ b_{\omega} = (2\sigma^2\omega^*\tilde{\omega}\Omega)^{1/2}, \quad (11) $$

where $\Omega$ is the conditional mean turbulent frequency (Van Slooten & Jayesh 1998)

$$ \Omega \equiv C_{\Omega} \frac{\overline{\rho\omega|\omega > \tilde{\omega}}}{\overline{\rho|\omega > \tilde{\omega}}}, \quad (12) $$

the fixed parameter $\sigma^2$ controls the variance of the $\omega^*$ distribution and $S_\omega$ is a source term for turbulent frequency.

The fundamental idea behind the introduction of this conditional frequency is to include external intermittency: if the fluid only has a certain probability of being turbulent at a given position in the flow, then the appropriate inverse mean turbulent time scale is not given by $\tilde{\omega}$ but by $\Omega$. The constant $C_\Omega$

is chosen so that $\Omega = \tilde{\omega}$ in a fully turbulent region. For homogeneous turbulence, without the source term, the stochastic model leads to a stationary distribution in the form of a $\Gamma$ function, with mean $\tilde{\omega}$ and variance $\sigma^2\tilde{\omega}^2$. In order to be consistent with DNS data for moderate Reynolds number (Yeung & Pope 1989), Van Slooten & Jayesh (1998) indicated that the adopted value for the turbulent frequency variance should be $\sigma^2 = 1/4$, in which case $C_\Omega = 0.69$. We adopt the same values in this paper.

In their approach, the inverse turbulent time scale in the velocity equation is given by $\Omega$, which is a mean quantity, instead of $\omega^*$, which is a stochastic quantity. This solves the issue about long distribution tails increasing the computational cost. Here, however, we only introduce $\Omega$ as a turbulent time scale in Eq. 4, and we keep $\omega^*$ in the velocity and energy equations (see Eqs. 7 and 8). This allows us to still include the effect of internal intermittency. By contrast with the experimentally verified log-normal distribution for $\omega$, the $\Gamma$ distribution has a much shorter tail, which alleviates at least part of the computational cost issue outlined above.

The quantity $S_\omega$ represent a source term for the turbulent frequency. In the $k-\epsilon$ standard model of turbulence, it is given by (Launder & Spalding 1974)

$$ S_{\omega} = (C_{\epsilon 2} - 1) - (C_{\epsilon 1} - 1) \frac{\mathcal{P}}{\tilde{\omega}k}, \quad (13) $$

where $\mathcal{P}$ is the rate of production of turbulent kinetic energy. The first constant term represents the decay of turbulent energy in the absence of a source, and the second term gathers all the sources of turbulent energy. The model constants are given by $C_{\epsilon 1} = 1.45$ and $C_{\epsilon 2} = 1.92$ (Launder & Sharma 1974). In incompressible turbulence, the kinetic energy production rate is given by $\mathcal{P} = 2C_\mu k S_\omega/\omega$, where $S$ is the contracted shear tensor

$$ S \equiv \frac{1}{4} \left( \frac{\partial \tilde{u}_i}{\partial x_j} + \frac{\partial \tilde{u}_j}{\partial x_i} \right) \left( \frac{\partial \tilde{u}_i}{\partial x_j} + \frac{\partial \tilde{u}_j}{\partial x_i} \right), $$

and $C_\mu = 0.09$ (Launder & Sharma 1974). In compressible turbulence, an extra source term due to buoyancy must be added, in the form of $\mathcal{P}_B = -\beta \tilde{u}_i'' T'' g_i$ (Rodi 1980), where $\beta = -(1/\rho)(\partial\rho/\partial T)_\rho$ is the coefficient of thermal expansion. With our notations, the total turbulent energy rate of production reads

$$ \mathcal{P} = C_{\mu} \frac{2k}{\tilde{\omega}} S + \frac{\Gamma_1 - 1}{n - 1} \frac{1}{\tilde{p}} \frac{\partial \tilde{p}}{\partial x_i} F_i^{(e)}, \quad (14) $$

## 3. Numerical setup

In this section, we present a two-dimensional implementation of the Lagrangian stochastic model, introduced in Sect. 2, with the objective of simulating a near-surface convective region.

### 3.1. Setup, boundary conditions, and initial conditions

We only follow two components of the fluid particles position – the component $x$ corresponding to the vertical coordinate, in the opposite direction of gravity, and the component $y$ representing one of the horizontal directions. However, all three components of the velocity are modelled, with the assumption of homogeneous turbulence in the $z$ direction. Our model is therefore 2.5-dimensional, thus increasing the realism of the simulation while keeping the computational cost manageable. The modelled variables are $x^*(t), y^*(t), u_x^*(t), u_y^*(t), u_z^*(t), e^*(t)$ and $\omega^*(t)$ for a set of $N$ fluid parcels, for a total of $7N$ variables which depend on
time *t* only. The uniform gravitational acceleration *g* and radialive diffusivity *κ* are fixed parameters.

In this Lagrangian stochastic framework, boundary conditions are enforced at particle level, every time a fluid particle exits the domain at the end of a time step. In order to conserve the total mass, we keep the number of particles constant, which means we reintroduce as many fluid particles as the ones exiting the domain. If a particle exits the domain through the top or bottom boundary, it is reintroduced with a reversed vertical velocity $u_y^*$, and its vertical position $x^*$ is mirrored with respect to the boundary (similar to a wall boundary condition). The turbulent frequency $\omega^*$, horizontal velocity $u_y^*$ and horizontal position $y^*$ of the particle are not modified. The internal energy is treated differently: in order to impose a given energy flux at the bottom and top edges of the domain, a linear mean internal energy profile $\tilde{e}(x)$ is imposed in two narrow bands at the top and bottom of the domain, spanning only a few percents of the total volume of the domain. In practice, at the end of each time step, each particle located in either of those bands is assigned an internal energy $e^* = e^{(bottom)} + \nabla_e x^*$, where $e^{(bottom)}$ is the imposed internal energy at the bottom of the simulation, and $\nabla_e = -F^{(bottom)}/\kappa$, where $F^{(bottom)}$ is the imposed energy flux. In these two narrow bands, the internal energy equation is not solved explicitly.

Periodic boundary conditions are enforced in the horizontal direction: a particle exiting through the left side of the box is reintroduced on the right side, and vice versa. The same periodic conditions are used in the averaging procedure (see Sect. 3.3): particles on the right side of the box contribute to the weighted average Eq. 19 computed on the left side, and vice-versa, in such a way that the domain folds back on itself horizontally.

For the initial state, we impose hydrostatic and radiative equilibrium, characterised by the following mean density and internal energy profiles

$$ \tilde{e}(x) = e^{(\text{bottom})} + \nabla_e x, \quad (15) $$

$$ \bar{\rho}(x) = \rho^{(\text{bottom})} \left( \frac{\tilde{e}(x)}{e_{\text{bot}}} \right)^{-\left(1+\frac{g}{\nabla_e(\Gamma_1 - 1)}\right)}, \quad (16) $$

where $\rho^{(\text{bottom})}$ is the initial mean density imposed at the bottom of the simulation. This amounts to a polytrop with index

$$ n_{\text{init}} = \left( 1 + \frac{\nabla_e(\Gamma_1 - 1)}{g} \right)^{-1} . \quad (17) $$

We impose that $n_{\text{init}}$ is identical to the polytropic index *n* introduced in the Lagrangian stochastic model, which constrains the value of $\nabla_e$.

The vertical position of each particle is initialised randomly according to the PDF $\bar{\rho}(x)$ given by Eq. 16. The particles are also assigned a random horizontal position chosen uniformly between the two horizontal boundaries of the domain. Once the position of each particle is chosen, they are assigned a random velocity according to a normal distribution with zero mean and a variance chosen arbitrarily (this arbitrary variance defines the initial Reynolds stress tensor); a random specific energy according to a normal distribution with a mean given by Eq. 15, and an arbitrary variance (which defines the initial value of $k_e$); and a random turbulent frequency according to a $\Gamma$ distribution with shape parameter $k = 4$ and a scale parameter $\theta = \omega_{\text{init}}/4$, where $\omega_{\text{init}}$ is the initial mean turbulent frequency, chosen arbitrarily. This initial turbulent frequency distribution corresponds to the solution of Eq. 4 in the stationary, homogeneous case.

## 3.2. Time integration scheme

The 7N variables of the system evolve according to Eqs. 1, 2, 3 and 4, which are stochastic differential equations of time only: one strength of this formalism is that there is no need to discretise the coordinate axes. We integrate the system forward in time using a stochastic predictor/corrector integration scheme, which is an algorithm of (weak) order 1.5 (Platen 1995). The scheme works in two steps. If we formally write the system of SDE in multivariate form $d\mathbf{X} = \mathbf{a}(\mathbf{X}, t) dt + \mathcal{B}(\mathbf{X}, t) \cdot d\mathbf{W}$, we first predict the state of the system after one time step $\Delta t$

$$ \tilde{\mathbf{X}}_{k+1} = \mathbf{X}_k + \mathbf{a}(\mathbf{X}_k, t_k)\Delta t + \mathcal{B}(\mathbf{X}_k, t_k) \cdot \Delta \mathbf{W}_k, $$

where $\Delta \mathbf{W}_k \equiv \mathbf{W}(t_{k+1}) - \mathbf{W}(t_k)$ is a set of independent random variables of mean 0 and variance $\Delta t$. Then this prediction is corrected according to (Platen 1995)

$$ \mathbf{X}_{k+1} = \mathbf{X}_k + (\alpha \mathbf{a}_\beta(\mathbf{X}_k, t_k) + (1-\alpha)\mathbf{a}_\beta(\tilde{\mathbf{X}}_{k+1}, t_{k+1}))\Delta t + (\beta \mathcal{B}(\mathbf{X}_k, t_k) + (1-\beta)\mathcal{B}(\tilde{\mathbf{X}}_{k+1}, t_{k+1}))\Delta \mathbf{W}_k. $$

The coefficients $\alpha$ and $\beta$ control the level of implicitness of the time integration scheme. They must be chosen between 0 and 1, with a value of 1 corresponding to a fully explicit scheme. In this study, we consistently chose $\alpha = \beta = 1/2$. The modified drift coefficient is defined by $\mathbf{a}_\beta \equiv \mathbf{a} - \beta \mathcal{B} \cdot (\nabla_X : \mathcal{B})$ and only differs from the actual drift coefficient if the diffusion matrix of the process explicitly depends on the stochastic process itself. In our case, this will only affect the equation on $\omega^*$.

The usual Courant–Friedrichs–Lewy (CFL) stability condition is not relevant in the present scheme, since the coordinate axes are not discretised. Here, the time step is chosen based on accuracy rather than stability considerations, and we impose that the particle Courant number of each notional particle must not exceed unity. This can be written (see Eq. (106) of Haworth 2010)

$$ \Delta t < \min_{1 \le k \le N} \left( \frac{1}{\omega^{*(k)}}, \frac{h_x}{|u_x^*|}, \frac{h_y}{|u_y^*|}, \frac{h_x^2 + h_y^2}{\kappa} \right), $$

where $h_x$ and $h_y$ are the vertical and horizontal widths of the kernel function which is used to compute the local mean flow quantities from the set of notional particles (see Sect. 3.3).

## 3.3. Ensemble average estimation using kernels

### 3.3.1. Standard kernel estimates

The Lagrangian stochastic model described by Eqs. 1, 2, 3 and 4 requires knowledge of the means (density $\bar{\rho}$, gas pressure $\bar{p}$, velocity $\tilde{\mathbf{u}}$, internal energy $\tilde{e}$, turbulent frequency $\tilde{\omega}$) and covariance (Reynolds stress tensor $R_{ij}$, energy variance $k_e$, energy flux $\mathbf{F}^{(e)}$) of the turbulent properties of the flow, at any given point of the domain. This information is extracted by averaging the appropriate particle-based stochastic variables over all neighbouring notional particles. To that end, as proposed by Welton & Pope (1997), we couple the Lagrangian stochastic model with a kernel estimation of the means (in a fashion similar to Smoothed Particle Hydrodynamics, see Liu et al. 2004, for details).

To filter local notional particles prior to averaging, we use a kernel function $K(\mathbf{r})$, where $\mathbf{r}$ is the position vector, which is centred on the location where the mean flow properties are evaluated. This function must satisfy some elementary conditions (e.g. Liu et al. 2004): it must peak at $\mathbf{r} = \mathbf{0}$, decrease outwards,
have a compact support, and be normalized so that its integral equals unity. We also choose – although it is not mandatory – a function which is symmetric, piecewise-quartic, possesses continuous first and second derivatives. Finally, with the objective of using an efficient algorithm (see below), we make the function separable so that

$$K(\mathbf{r}) = K(r_x, h_x) K(r_y, h_y),$$

where $r_x$ and $r_y$ are the vertical and horizontal component of the 2D position vector $\mathbf{r}$, and $h_x$ and $h_y$ are the vertical and horizontal extents of the kernel function. The following function (also adopted by Welton 1998) fulfil all requirements

$$K(r, h) = \begin{cases} \frac{5}{4h} \left(1 + 3\left|\frac{r}{h}\right|\right) \left(1 - \left|\frac{r}{h}\right|\right)^3, & \text{if } r < h, \\ 0, & \text{if } r > h. \end{cases} \quad (18a) \quad (18b)$$

The ensemble average of any particle-based quantity $Q^*$ at any position $\mathbf{x}$ and time $t$ can then be estimated via

$$\overline{Q}(\mathbf{x}, t) = \sum_{k=1}^{N} \mathcal{V}^{*(k)}(t) \mathcal{Q}^{*(k)}(t) K(\mathbf{x}^{*(k)}(t) - \mathbf{x}, h), \quad (19)$$

where $\mathcal{V}^{*(k)} \equiv \Delta m/\rho^{*(k)}$ is the lumped volume of fluid that particle $k$ is meant to represent. The difficulty, then, is that $\rho^{*(k)}$ is not a modelled quantity. For that reason, the kernel estimation is much more suited to yield density-averaged means

$$\overline{Q}(\mathbf{x}, t) = \frac{\overline{\rho Q}}{\overline{\rho}} = \frac{1}{\overline{\rho}} \sum_{k=1}^{N} \Delta m \mathcal{Q}^{*(k)}(t) K(\mathbf{x}^{*(k)}(t) - \mathbf{x}, h), \quad (20)$$

where the mean density is given by Eq. 19 with $Q = \rho$

$$\overline{\rho}(\mathbf{x}, t) = \sum_{k=1}^{N} \Delta m K(\mathbf{x}^{*(k)}(t) - \mathbf{x}, h).$$

The gradient of $\overline{Q}$ can also be estimated using this kernel method: because the position variable $\mathbf{x}$ only appears in the kernel function $K$, we have

$$\nabla \overline{Q} = - \sum_{k=1}^{N} \mathcal{V}^{*(k)}(t) \mathcal{Q}^{*(k)}(t) \nabla K(\mathbf{x}^{*(k)}(t) - \mathbf{x}, h),$$

and higher-order spatial derivatives may be computed in the same way.

From a practical point of view, directly computing the means and their derivatives would be very time-consuming as it would be a $O(N^2)$ algorithm. To overcome this issue, we adopt the algorithm developed by Welton (1998), which mainly consists in projecting the kernel function in Fourier space. This leads to a $O(M^2N)$ algorithm, where $M$ is the number of notional particles, and $M$ the number of Fourier modes on which the kernel is projected. In practice, we adopt a value $M = 10$, which is enough to accurately represent the kernel given by Eq. 18 up to its second derivative, while leading to a drastic decrease in computation time compared to the brute-force $O(N^2)$.

### 3.3.2. Correcting for kernel estimation inconsistency

Any estimator for a mean field $\overline{Q}(\mathbf{x}, t)$ is said to have $k$-th order consistency if it can exactly reproduce a polynomial function of the coordinates up to degree $k$. The standard formalism presented in Sect. 3.3.1 fails this condition. There are two reasons. First, the averaging kernel function $K$ is truncated by the edges of the domain: when trying to apply Eq. 19 to a position $\mathbf{x}$ close to the boundary of the domain, the lack of particles on a certain portion of the kernel support will artificially decrease $\overline{Q}$. The second reason is that an irregular particle distribution can introduce a bias not only in $\overline{Q}$, but especially in its gradient. There are several ways of correcting the standard scheme presented in Sect. 3.3.1 so as to restore its consistency. Here, we follow Korzilius et al. (2016), who provide a scheme with 1st order consistency for the means $\overline{Q}$, their gradient and their 2nd order derivatives. We describe this Corrective Smoothed Particle Method (CSPM) in App. B.

## 4. Benchmarking and validation

To validate the numerical implementation described in Sect. 3, we benchmarked the code against standard test cases used for particle-based hydrodynamics codes (Liu et al. 2004). These benchmark cases involve incompressible flows; in this code, the incompressible condition is enforced by setting the lumped volume $\mathcal{V}^{*(k)}$ of each particle to the total volume of the domain divided by the total number of particles, at each time step, regardless of particle motion.

### 4.1. Deterministic cases: plane Poiseuille and Couette flows

These benchmark cases allow us to test specifically the accuracy of the kernel estimation scheme presented in Sect. 3.3.

The plane Poiseuille flow is obtained by treating the top and bottom boundaries as two parallel plates with no-slip boundary conditions, and by exerting a constant and uniform body-force on the flow, parallel to the plates. We initially placed a square of $150 \times 150$ particles in the domain, regularly spaced. The particles are evolved according to the following set of deterministic differential equations

$$dx = u_x dt, \quad (21)$$

$$dy = u_y dt, \quad (22)$$

$$du_x = \left( -\frac{1}{\overline{\rho}} \frac{\partial \overline{p}}{\partial x} + v \Delta \overline{u}_x \right) dt, \quad (23)$$

$$du_y = \left( -\frac{1}{\overline{\rho}} \frac{\partial \overline{p}}{\partial y} + v \Delta \overline{u}_y + F \right) dt, \quad (24)$$

and the gas pressure is related to the density through $\overline{p} = \overline{\rho} c_s^2$. The exact solution can be expressed as an infinite series

$$u_y(x, t) = \frac{Fx(x - L_x)}{2\nu} + \sum_{n=0}^{+\infty} \frac{4FL_x^2}{\nu\pi^3(2n + 1)^3} \sin\left(\frac{\pi x(2n + 1)}{L_x}\right) \\ \times \exp\left(-\frac{(2n + 1)^2\pi^2\nu}{L_x^2}t\right). \quad (25)$$

Here, the code described in Sect. 3 is implemented without the stochastic terms, thus only keeping the deterministic part of the equations in the form of Eqs. 21 to 24. The domain is chosen to be $10^{-3}$ m in height and in width. The mass of each particle
<figure><img src="image_2.png" /><figcaption>Fig. 1. Incompressible Poiseuille flow benchmark case. The coloured dots show snapshots of the Poiseuille velocity profile at different times. The solid lines show the exact analytical solution given by Eq. 25. Relative errors remain under 1% throughout the simulation.</figcaption></figure>

is chosen so that the flow has a uniform density $\rho = 10^3$ kg/m$^2$. The kinematic viscosity is set to $\nu = 10^{-6}$ m$^2$/s, and the driving body acceleration to $F = 2 \times 10^{-4}$ m/s$^2$. We assumed a time step $\Delta t = 10^{-4}$ s, and we run the simulation for 5000 time steps. The size of the kernel function is set to 3% of the domain size, both horizontally and vertically. Figure 1 shows the horizontal velocity profile, as a function of vertical coordinate, at different times in the simulation. The departure from the exact solution, given by Eq. 25, remains smaller than 1% throughout the whole simulation.

The Couette flow benchmark case is similar to the Poiseuille flow, the difference being that no body-force is exerted on the flow, and the motion is driven instead by moving one of the plates at a constant velocity $v_0$. The Couette flow simulation we run has the same setup as the Poiseuille flow, except that the no-slip condition at the top boundary is replaced by a constant horizontal velocity condition. We adopted $v_0 = 2.5 \times 10^{-5}$ m/s.

The exact solution can be expressed as an infinite series

$$u_y(x, t) = \frac{v_0 x}{L_x} + \sum_{n=1}^{+\infty} (-1)^n \frac{2v_0}{n\pi} \sin\left(\frac{n\pi x}{L_x}\right) \exp\left(-\frac{n^2\pi^2 v}{L_x^2}t\right). \quad (26)$$

Figure 2 shows the Couette flow horizontal velocity profile, in the same way as Fig. 1 for the Poiseuille case. The departure from the exact solution, given by Eq. 26, remains smaller than 0.6% throughout the whole simulation.

The topology of the plane Poiseuille and plane Couette flows is mainly controlled by the viscous force. As such, the very good agreement found in these two cases showcases not only the validity of the time integration scheme implemented in the code, but also the accuracy of the corrective kernel estimation scheme up to second spatial derivative.

## 4.2. Stochastic case: incompressible Simplified Langevin Model

In order to test the accuracy of the stochastic scheme for particle evolution along with the kernel-estimation of the means, we benchmarked the code against the incompressible Simplified

<figure><img src="image_3.png" /><figcaption>Fig. 2. Incompressible Couette flow benchmark case. The symbols are the same as in Fig. 1, with the exact analytical solution being given by Eq. 26. Relative errors remain under 0.6% throughout the simulation.</figcaption></figure>

Langevin Model, where the particles evolve according to the following set of stochastic differential equations

$$dx = u_x dt, \quad dy = u_y dt, \\ du_\alpha = -\left(\frac{1}{2} + \frac{3}{4}C_0\right)\omega(u_\alpha - \tilde{u}_\alpha)dt + \sqrt{C_0\omega k} dW_\alpha,$$

where $\alpha = \{x, y, z\}$, and $\omega$ is the fixed turbulent frequency. We recall that $\tilde{u}_\alpha$ is the local density-averaged velocity, $k$ is the local turbulent kinetic energy, and $C_0$ is the Kolmogorov constant. This corresponds to the model presented in Sect. 2, in the limit where every particle's turbulent frequency $\omega$ and internal energy $e$ is kept constant (we set $\omega = 0.1$, and $e$ is irrelevant in this model) so that only their velocity and position is allowed to evolve. In this model, the turbulence is homogeneous, and the turbulent kinetic energy $k$ decays with time at a rate equals to $\omega$ (e.g. Pope 2000)

$$k(t) = k_0 \exp^{-\omega t}, \quad (27)$$

where $k_0$ is the initial turbulent kinetic energy at $t = 0$. We ran the simulation with 300,000 particles, with a time step $\Delta t = 1/(10\omega)$, for a total duration of $100/\omega$. The kernel size was set to 3% of the domain size, both in the vertical and horizontal directions.

We show in Fig. 3 the evolution of the spatially averaged turbulent kinetic energy $k$ with time in the simulation output, compared to the expected result given by Eq. 27. The relative error never exceeds 0.35%. This shows that the validity of the kernel estimation scheme extends to mean flow properties – and especially to centered second-order moments – in the presence of turbulence, and with the particles evolving according to stochastic equations rather than deterministic ones.

## 4.3. Convergence

In order to determine the appropriate number $N$ of particles needed to correctly represent the flow, as well as the appropriate kernel size $h$ and time step $\Delta t$, we ran the compressible Simplified Langevin Model case (i.e. the same model as in Sect. 4.2, but estimating mean flow density from particle positions on the
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of the spatially averaged turbulent kinetic energy $k$ with time, for the benchmark incompressible Simplified Langevin Model case. The dots show the simulation output, while the solid line shows the analytical solution given by Eq. 27. Relative errors remain under 0.35% throughout the simulation.</figcaption></figure>

fly, rather than setting a fixed value), with the same setup as in Sect. 4.2.

We first ran simulations with different kernel sizes $h$, while keeping the time step constant ($\Delta t = 0.003$), but we varied $N$ alongside $h$ in such a way that the average number of particles counting towards a kernel average is fixed. This amounts to $N \propto 1/h^2$, with a base value of $N = 200000$ for $h = L/30$, meaning that there is ~220 particles per kernel average. We then extracted the spatial average of the turbulent kinetic energy $k$, as a function of time, for each simulation. The top panel of Fig. 4 shows that the results are independent of the kernel size if the kernel size remains smaller than about 2% of the total domain size.

We did a similar comparison by varying the time step $\Delta t$, while keeping the number of particles ($N = 200,000$) and the kernel size ($h = L/30$) fixed. The middle panel of Fig. 4 shows that the results are independent of the time step for $\Delta t < 0.001$.

Finally, ran different simulations with varying numbers of particles $N$, while keeping the time step ($\Delta t = 0.003$) and the kernel size ($h = L/30$) fixed. The bottom panel of Fig. 4 shows that the results are independent of the particle number for $N > 500,000$.

Based on these results, in the rest of this study, we will adopt a kernel size of $h = L/50$, a time step of $\Delta t = 0.001$, and a particle number of $N = 500,000$.

## 5. 2D convectively unstable simulation

To showcase the capabilities of our code, we simulated a 0.45 × 2.7 Mm convectively unstable domain, using the Lagrangian stochastic model developed in Sect. 2 and the numerical scheme described in Sect. 3. We considered a uniform polytropic index $n = 9$, as well as a uniform radiative diffusivity $\kappa = 23 \text{ km}^2/\text{s}$, throughout the entire domain. The first adiabatic index is set to $\Gamma_1 = 5/3$, the gravitational acceleration to the solar surface gravity $g = 275 \text{ m/s}^2$. Based on the results outlined in Sect. 4.3, we adopted a time step of $\Delta t = 0.001$ (in dimensionless units, corresponding to a physical time step of about 0.1 s), the kernel size is

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Evolution of the spatially averaged turbulent kinetic energy $k$ with time, for the compressible Simplified Langevin Model case, for different run parameters. **Top:** comparison of different kernel sizes $h$. **Middle:** comparison of different time steps $\Delta t$. **Bottom:** comparison of different particle numbers $N$.</figcaption></figure>

2% of the domain size in both directions, and we used 500,000 particles to represent the flow.

It should be noted that in mesh-free, particle methods, the choice of the kernel size $h$ is dictated not by requirements of stability (as is the grid resolution in Eulerian methods), but by requirements of accuracy and convergence (specifically weak convergence in our case). In this matter, the important parameter is the average number of particles in each kernel estimation, which is equal to $N_{pc} = N(h/L)^2$. With the values adopted above, we
have $N_{pc} = 200$, which amounts to statistical errors in the estimation of the means of about 3% (e.g. Haworth 2010).

## 5.1. Particle properties

The primary output of the simulation consists of the physical properties of each particle (position, velocity, specific internal energy and turbulent frequency) as a function of time only. A snapshot of the simulation is shown in Fig. 5, where each dot represents a particle. We note that, for better readability, only 10,000 notional particles (5% of the total) are shown. The simulated flow clearly features a convectively unstable structure (especially visible in the vertical velocity and internal energy), with alternating hot up-flows and cold down-flows.

## 5.2. Mean flow properties

Using the kernel estimation scheme described in Sect. 3.3, we can extract mean flow properties from the particles. The novelty of the present particle-based approach lies in the fact that it yields proper ensemble averages, which are neither time averages nor horizontal averages. That way, we do not rely on the usual proxies for estimating ensemble averages of turbulent properties, and we can estimate mean flow properties that depend both on spatial coordinates and on time.

The top, middle, and bottom panels of Fig. 6 show maps of the turbulent kinetic energy $k \equiv \tilde{u}_i''\tilde{u}_i''/2$, internal energy variance $k_e \equiv e''^2$, and vertical internal energy flux $F_x^{(e)} \equiv \tilde{e}''\tilde{u}_i''$ respectively, for the same time step as in Fig. 5. The second-order statistics of the flow clearly have a different depth dependence in the upflows and downflows. In particular, the turbulent kinetic energy is higher in the downflows than in the upflows, which leads to a negative overall kinetic energy flux throughout the convectively-unstable region. This is in agreement with the known behaviour of stellar surface convection, from 3D hydrodynamic simulations (Kupka & Muthsam 2017).

## 5.3. Energy flux budget

Adding together the equation on large-scale kinetic energy – stemming from Eq. A.13 –, the equation on small-scale kinetic energy – i.e. the trace of Eq. A.15 – and the equation on mean internal energy Eq. A.14, we obtain the following total energy equation

$$ \frac{\partial}{\partial t} \left( \overline{\rho \tilde{e}} + \overline{\rho k} + \frac{1}{2} \overline{\rho \tilde{u}_i \tilde{u}_i} \right) + \frac{\partial}{\partial x_i} \left( F_i^{(\text{e})} + F_i^{(\text{kin})} + F_i^{(\text{rad})} + F_i^{(\text{grav})} \right) = 0, \quad (28) $$

where

$$ F_i^{(\text{e})} = \overline{\rho \tilde{e} u_i} + \frac{n - \Gamma_1}{n - 1} \overline{\rho e'' \tilde{u}_i''} + \overline{\rho \tilde{u}_i}, \quad (29) $$

$$ F_i^{(\text{kin})} = \overline{\rho k \tilde{u}_i} + \frac{1}{2} \overline{\rho u_i'' u_j'' u_j'' + \rho R_{ij} \tilde{u}_j} + \frac{1}{2} \overline{\rho \tilde{u}_i \tilde{u}_j \tilde{u}_j}, \quad (30) $$

$$ F_i^{(\text{rad})} = -\overline{\rho \kappa \frac{\partial \tilde{e}}{\partial x_i}}, \quad (31) $$

$$ F_i^{(\text{grav})} = \overline{\rho \tilde{u}_i \Psi}. \quad (32) $$

In the stationary state, the sum of all energy fluxes should be conserved. This is illustrated in Fig. 7, which show the vertical profiles (i.e. the time and horizontal average) of the different contributions to the energy flux budget equation Eq. 28, excluding the relaxation phase at the start of the simulation. The

<figure><img src="image_6.png" /><figcaption>Fig. 5. Snapshot of the particle properties for the convectively unstable simulation described in Sect. 5. Each dot represents one particle, and the color code refers to the individual vertical velocity $u_x^*$ (top), specific internal energy $e^*$ (middle), and turbulent frequency $\omega^*$ (bottom) of each particle. Only 10,000 particles are shown for readability, whereas the simulation contains 500,000.</figcaption></figure>

two main contributions come from the convective and radiative fluxes, with the kinetic energy flux being of much smaller amplitude, but systematically negative. Breaking down the different contributions to the convective flux, however (see middle panel), it is apparent that all contributions account for a significant portion of the total, with no one contribution standing out as dominant.

We note that the energy conservation breaks down in the two narrow bands located at the top and bottom of the simulation –
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Snapshot of the second-order moments maps for the convectively unstable simulation described in Sect. 5, for the same time step as Fig. 5. **Top:** turbulent kinetic energy $k = \tilde{u}_i^t \tilde{u}_i^t / 2$. **Middle:** internal energy variance $k_v = \tilde{e}^{\prime \prime 2}$. **Bottom:** vertical internal energy flux $F_x^{(e)} = \tilde{u}_x^{\prime \prime} \tilde{e}^{\prime \prime}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Horizontal and time average of the energy flux budget. The grey areas mark the two narrow regions at the top and bottom of the domain where the internal energy profile is imposed, and the particle-based internal energy equation is not solved. **Top:** Total energy flux budget, defined in Eq. 28. **Middle:** Convective energy flux budget, defined in Eq. 29. **Bottom:** Kinetic energy flux budget, defined in Eq. 30.</figcaption></figure>

grey bands in Fig. 7 –, where the particle-based internal energy equation given by Eq. 3 is not solved, but instead an internal energy profile is imposed (see Sect. 3.3).

The constant $C_1$ introduced in the Lagrangian stochastic model was chosen to be equal to the Kolmogorov constant $C_0$. This choice, however, is arbitrary. Since the constant $C_1$ controls the decay rate of internal energy fluctuations, changing its value can be expected to modify the second order moments of turbulent convection, and in particular the energy flux budget. In order

to investigate this, we ran a simulation entirely identical to the
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison between the energy flux budgets of the two simulations described in the text, the only modified parameters being the constant $C_1$ in the Lagrangian stochastic model (solid lines: $C_1 = 3.5$; dashed lines: $C_1 = 2.0$). The different colours represent different terms in the energy flux budget, as indicated on the plot.</figcaption></figure>

one described in rest of this section, but adopting $C_1 = 2.0 \neq C_0$. The energy flux budgets are compared in Fig. 8. Modifications to the energy flux budget remain small, with the quantity most impacted (in terms of relative values) being the turbulent kinetic energy flux.

### 5.4. Up- and down-flow asymmetry

Because of the presence of upflows and downdrafts, the vertical velocity and internal energy distributions are not expected to be Gaussian (e.g. Kupka & Robinson 2007). To measure the departure from Gaussian distributions, it is customary to define the skewness $S_w$ (resp. $S_\theta$) and kurtosis $K_w$ (resp. $K_\theta$) of the vertical velocity (resp. internal energy) distribution as

$$S_w = \overline{u_x''^3} / \overline{u_x''^2}^{3/2}, \quad K_w = \overline{u_x''^4} / \overline{u_x''^2}^2, \\ S_\theta = \overline{e'''^3} / \overline{e'''^2}^{3/2}, \quad K_\theta = \overline{e'''^4} / \overline{e'''^2}^2.$$

where $u_x'' = u_x - \tilde{u}_x$ and $e'' = e - \tilde{e}$. The velocity skewness gives a measure of the asymmetry between the flows directed upwards and downwards (positive skewness means that upflows take up less space but are stronger, negative skewness means the opposite), while the kurtosis measures the rarity of intermittent, stronger-than-normal events (large kurtosis means that such events are frequent, small kurtosis means that they are rare, and a Gaussian distribution has a kurtosis of 3). The internal energy skewness and kurtosis have a similar meaning for the asymmetry between hot and cold flows.

We sliced the domain into horizontal layers spanning 1% of the domain each, and for each layer we extracted the vertical velocity and internal energy of all notional particles located in that slice, thus reconstructing the velocity and energy distributions as a function of height. We show the skewness and kurtosis of each distribution, as a function of height, in Fig. 9. As expected, the velocity distribution is skewed towards downwards velocities, pointing to the downdrafts taking less space than the upflows, but being more turbulent (as visible also from the top panel of Fig. 6). Furthermore, while the kurtosis of the energy distribution

<figure><img src="image_10.png" /><figcaption>Fig. 9. Top: Skewness and kurtosis of the vertical velocity distribution and internal energy distribution, as a function of height. Bottom: Same, but in the $(K, S^2)$ plane. Each cross corresponds to one horizontal layer in the simulation.</figcaption></figure>

remains close to the Gaussian value of 3 in the bulk of simulation, the velocity distribution is distinctly platikurtic throughout – meaning that its tails are less prominent than those of a normal distribution. These results – and in particular the $(K, S^2)$ diagram shown in the bottom panel of Fig. 9 – are in good agreement with results from Large-Eddy Simulations of solar convection, both 2D (Kupka & Robinson 2007) and 3D (Cai 2018).

The fact that the kernel estimation scheme provides time- and space-dependent ensemble averages also allows to compare vertical profiles of various mean quantities in up-flows and down-flows. We show in Fig. 10 how the energy flux budget presented in Sect. 5.3 breaks down into upwards and downwards energy flow. The most striking feature is the strong positive (resp. negative) convective flux in the upflows (resp. downflows), which cancel out when the total convective flux is computed (see Fig. 7). It can also be seen that the downwards negative kinetic energy flux more than cancels out the upwards positive kinetic energy flux, thus yielding a negative total kinetic energy flux.
<figure><img src="image_11.png" /><figcaption>Fig. 10. Breakdown of the energy flux budget into upflow (solid lines) and downflow (dashed lines) energy fluxes. Blue, green, orange and black lines respectively show the convective flux $F^{(e)}$, the radiative flux $F^{(rad)}$, the kinetic energy flux $F^{(kin)}$ and the total energy flux $F^{(tot)}$, as defined by Eq. 28.</figcaption></figure>

## 5.5. Turbulent diffusivity

The Lagrangian framework adopted in our simulations also allows us to measure particle dispersion, and therefore turbulent viscosity, much more easily and accurately than in Large-Eddy Simulations. Subdividing the simulation into time series of 100 seconds each, we extracted the travel distance of each particle as a function of time increment

$$r(\tau)^2 = (x^*(t_0 + \tau) - x^*(t_0))^2 + (y^*(t_0 + \tau) - y^*(t_0))^2, \quad (33)$$

where $t_0$ refers to the initial time of the time series. We cut each individual particle trajectory when the particle gets too close to one of the domain boundaries (to avoid the effect of boundary conditions on particle dispersion). We then averaged $r(\tau)$ over all particles, for every available value of the time increment $\tau$, regardless of initial particle position. The result is shown in Fig. 11, where one can clearly make out the two expected regimes (Taylor 1922): for short travel times, particle dispersion is dominated by the initial velocity, which translates to an advective regime where $r(\tau) \propto \tau$, and the proportionality constant is the initial particle velocity $u_{drift}$ (because we average over different particle trajectories, the value $u_{drift}$ inferred from particle dispersion in the simulation must be taken to represent the RMS value of initial particle velocity). On the other hand, for long travel times, it is dominated by turbulence, which translates to a diffusive regime where $r^2(\tau) \propto \tau$, and the proportionality constant is the turbulent diffusivity $v_{turb}$.

We fit the following model to the data

$$r(\tau)^2 = \begin{cases} u_{\text{drift}}^2 \tau^2 & \text{if } \tau < \tau_0, \\ v_{\text{turb}} \left( \tau - \frac{\tau_0}{2} \right) & \text{if } \tau > \tau_0, \end{cases}$$

where $\tau_0 \equiv v_{turb}/(2u_{drift}^2)$ marks the transition between the two regimes. We find a typical initial velocity of $u_{drift} \sim 2$ km/s, and a turbulent diffusivity of $v_{turb} \sim 500$ km$^2$/s. Those are in agreement with the typical values of the convective velocity and turbulent viscosity at the top of the solar convective zone (e.g. Rincon et al. 2025). Alternatively, we can estimate the turbulent

<figure><img src="image_12.png" /><figcaption>Fig. 11. Squared particle dispersion distance $r^2$ as a function of time increment $\tau$, averaged over all particles trajectories. The black dots show the simulation data, and the solid black line shows a fit to the data, using the model described in the text. The RMS particle velocity $u_{drift}$ and turbulent diffusion coefficient $v_{turb}$ estimated from the fit are indicated on top.</figcaption></figure>

Reynolds number $Re_{turb}$, defined by

$$Re_{turb} \equiv \frac{u_{\text{conv}} L}{v_{\text{turb}}},$$

where $L$ is the typical length scale of the flow, which we take to be the vertical extent of the simulation ($L = 0.45$ Mm). With these values, we find $Re_{turb} \sim 2$. We also find that the transition between the advective and diffusive regimes occurs for $\tau_0 \sim 1$ minute, corresponding to the observed convective turnover time at the solar surface.

We can estimate the turbulent thermal diffusive time scale in a similar manner. Instead of extracting the travel distance of each particle, we extract the autocorrelation function of the internal energy fluctuations of each particle, as a function of time increment

$$C_e(t, \tau) \equiv \frac{\langle\langle(e^*(t) - \bar{e}(\mathbf{x}^*(t)))\rangle\rangle \langle\langle(e^*(t + \tau) - \bar{e}(\mathbf{x}^*(t + \tau)))\rangle\rangle}{\langle\langle(e^*(t) - \bar{e}(\mathbf{x}^*(t)))^2\rangle\rangle} . \quad (34)$$

We then average it over all particles. Figure 12 clearly shows that the autocorrelation function takes the expected form of a decaying exponential function. The decay time represents the integral Lagrangian time of the energy fluctuations (i.e. the typical time over which the thermal memory of the particle is lost). It represents the turbulent thermal diffusive time scale $\tau_{th}$: the fit shown in Fig. 12 provides $\tau_{th} = 65 \pm 5$ s, depending on which sub-timeseries of the simulation we considered. This is remarkably close to the momentum diffusivity time scale $\tau_0$ determined above, thus yielding a value $Pr \sim 1$ for the turbulent Prandtl number. This is in accordance with what is expected in the regime of strong turbulence, in line with the assumptions made in the Reynolds analogy for instance.

# 6. Conclusion

In this paper, we developed a Lagrangian stochastic model for the one-point statistics of turbulent convection in stars. This
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Autocorrelation function of particle energy fluctuations, as a function of time increment (blue), averaged over all particles trajectories. The dashed black line shows an exponential fit to the curve.</figcaption></figure>

model uses a particle representation of the flow, and is centered
around evolution equations for the position, velocity, internal en-
ergy and turbulent frequency of each particle, as a function of
time only, in the form of Eqs. 1 to 4. This set of particle-based
stochastic equations is statistically equivalent to a Fokker-Planck
equation for the Probability Density Function (PDF) of all rele-
vant turbulent variables, which itself implies transport equations
for the means and cross-variances of the velocity, internal en-
ergy, and turbulent frequency of the flow. We established this
Lagrangian stochastic model in such a way that these transport
equations are identical to the exact transport equations stemming
from first principles, using only physically motivated assump-
tions to close the equations. The main advantage of this approach
is that, thanks to the adopted Lagrangian point of view, it mod-
els the advective terms exactly, so that closure relations are not
needed there.

In a second part of the paper, we then developed a new nu-
merical code to implement this Lagrangian stochastic model in
a 2D setting. The code solves the set of stochastic differential
equations (SDE) determined in the first part of the paper, for
a large number of notional particles (several hundreds of thou-
sands to several millions). These notional particles act as real-
isations of the same flow, and interact with each other through
the mean flow fields, which are determined from the set of par-
ticle realisations themselves, through a filtering kernel estima-
tion procedure, akin to Smoothed Particle Hydrodynamics. In
classical hydrodynamic simulations, statistical averages of tur-
bulent properties can only be studied through the usual proxies
of horizontal or temporal averages. By contrast, this code mod-
els actual ensemble averages of turbulent flow properties, that is
to say both means (density, velocity, internal energy, turbulent
frequency) and second-order moments (turbulent kinetic energy,
Reynolds stress tensor, internal energy variance, internal energy
flux) as a function of time and all spatial coordinates. It is there-
fore much more suited to study how the statistical properties of
stellar turbulent convection depends on time, and especially how
they affect the propagation of waves. In addition, this code also
includes an equation for the turbulent dissipation, which is mod-
elled as a bona fide turbulent variable, which contrasts with Eu-
lerian grid-based hydrodynamic simulations where turbulent dis-

sipation is dominated by artificial viscosity, and cannot be mod-
elled physically.

The example presented in the second part of this paper aims
at showcasing the possibilities offered by this new code. In sub-
sequent studies, we aim to use this code for several purposes.
First, because statistical averages can be obtained without having
to resort to proxies such as horizontal averages or time averages,
this code is much more suited to study the interaction between
waves and convection. It will be possible to extract oscillating
modes directly from the simulation, and to study its imprint on
the various second-order moments of turbulence. In particular,
we will be able to extract directly the wave-induced fluctuation
of turbulent pressure, convective flux, internal energy variance,
etc. This will be of paramount importance to quantify surface ef-
fects (and particularly the modal part of surface effects) without
having to rely on Mixing Length Theories.

Second, as showcased in Sect. 5.5, the purely Lagrangian na-
ture of the code makes it perfect to study particle dispersion. In
classical hydrodynamic simulations, Lagrangian trackers must
be injected in order to measure turbulent diffusion, and the re-
sults usually suffer from the finite grid resolution. Here, La-
grangian trajectories of flow realisations are already followed
by construction, and the determination of diffusion coefficients
from these trajectories does not rely on any Eulerian grid. In par-
ticular, this code allows, as will be investigated in a future study,
to assess the efficiency of turbulent transport through convec-
tive/radiative interfaces, and therefore of convective penetration
in the stellar context.

Finally, it will be possible to extend the code to 3D rather
than 2D. While this is not conceptually complicated, it will be
considerably more demanding in computational power and time.
More specifically, while the appropriate number of notional par-
ticles needed to represent the flow in 2D is of the order of several
hundreds of thousands, this number will need to be increased to
tens or hundreds of millions in 3D.

*Acknowledgements.* The authors wish to thank the anonymous referee for use-
ful comments that led to improvement of the manuscript. The authors acknowl-
edge support from the "Action Thématique de Physique Stellaire" (ATPS) of
CNRS/INSU PN Astro co-funded by CEA and CNES.

References

Ahlborn, F., Kupka, F., Weiss, A., & Flaskamp, M. 2022, A&A, 667, A97

Almeida, Y. & Navarro-Martinez, S. 2021, Physics of Fluids, 33, 035155

Baglin, A., Auvergne, M., Barge, P., et al. 2006a, in ESA Special Publication, Vol. 1306, The CoRoT Mission Pre-Launch Status - Stellar Seismology and Planet Finding, ed. M. Fridlund, A. Baglin, J. Lochard, & L. Conroy, 33

Baglin, A., Auvergne, M., Boisnard, L., et al. 2006b, in 36th COSPAR Scientific Assembly, Vol. 36, 3749

Bakosi, J., Franzese, P., & Boybeyi, Z. 2008, Journal of Computational Physics, 227, 5896

Balmforth, N. J. 1992a, MNRAS, 255, 603

Balmforth, N. J. 1992b, MNRAS, 255, 632

Belkacem, K., Dupret, M. A., Baudin, F., et al. 2012, A&A, 540, L7

Belkacem, K., Kupka, F., Philidet, J., & Samadi, R. 2021, A&A, 646, L5

Belkacem, K., Kupka, F., Samadi, R., & Grimm-Strele, H. 2019, A&A, 625, A20

Belkacem, K. & Samadi, R. 2013, in Lecture Notes in Physics, Berlin Springer Verlag, ed. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green, Vol. 865 (Springer), 179

Böhm-Vitense, E. 1958, ZAp, 46, 108

Böhm-Vitense, E. 1992, Introduction to stellar astrophysics. Volume 3. Stellar structure and evolution., Vol. 3 (Cambridge University Press)

Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977

Borne, V. & Orszag, S. A. 1995, Phys. Rev. E, 51, R856

Cai, T. 2018, ApJ, 868, 12

Canuto, V. M. 1997, ApJ, 482, 827

Christensen-Dalsgaard, J. & Thompson, M. J. 1997, MNRAS, 284, 527

Das, S. K. & Durbin, P. A. 2005, Physics of Fluids, 17, 025109
Delarue, B. J. & Pope, S. B. 1997, Physics of Fluids, 9, 2704

Delarue, B. J. & Pope, S. B. 1998, Physics of Fluids, 10, 487

Dopazo, C. 1975, Physics of Fluids, 18, 397

Dreeben, T. D. & Pope, S. B. 1998, Journal of Fluid Mechanics, 357, 141

Gough, D. O. 1977, ApJ, 214, 196

Grigahcène, A., Dupret, M. A., Gabriel, M., Garrido, R., & Scuflaire, R. 2005, A&A, 434, 1055

Grigahcène, A., Dupret, M. A., Sousa, S. G., et al. 2012, MNRAS, 422, L43

Haworth, D. C. 2010, Progress in Energy and combustion Science, 36, 168

Haworth, D. C. & Pope, S. B. 1987, Physics of Fluids, 30, 1026

Heinz, S. 2003, Statistical mechanics of turbulent flows (Springer)

Houdek, G. & Dupret, M.-A. 2015, Living Reviews in Solar Physics, 12, 8

Houdek, G., Lund, M. N., Trampedach, R., et al. 2019, MNRAS, 487, 595

Houdek, G., Trampedach, R., Aarslev, M. J., & Christensen-Dalsgaard, J. 2017, MNRAS, 464, L124

Janicka, J., Kolbe, W., & Kollmann, W. 1979, Journal of Non Equilibrium Thermodynamics, 4, 47

Kippenhahn, R., Weigert, A., & Hofmeister, E. 1967, Methods in Computational Physics, 7, 129

Korzilius, S., Schilders, W. H., & Anthonissen, M. J. 2016, Journal of Applied Mathematics and Physics, 5, 168

Kupka, F., Ahlborn, F., & Weiss, A. 2022, A&A, 667, A96

Kupka, F. & Muthsam, H. J. 2017, Living Reviews in Computational Astrophysics, 3, 1

Kupka, F. & Robinson, F. J. 2007, MNRAS, 374, 305

Launder, B. & Spalding, D. 1974, Computer Methods in Applied Mechanics and Engineering

Launder, B. E. & Sharma, B. I. 1974, Letters Heat Mass Transfer, 1, 131

Lebreton, Y., Goupil, M. J., & Montalbán, J. 2014a, in EAS Publications Series, Vol. 65, EAS Publications Series, ed. Y. Lebreton, D. Valls-Gabaud, & C. Charbonnel, 99–176

Lebreton, Y., Goupil, M. J., & Montalbán, J. 2014b, in EAS Publications Series, Vol. 65, EAS Publications Series, ed. Y. Lebreton, D. Valls-Gabaud, & C. Charbonnel, 177–223

Liu, G. R., Liu, M. B., & Li, S. 2004, Computational Mechanics, 33, 491

Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10

Otí Floranes, H., Christensen-Dalsgaard, J., & Thompson, M. J. 2005, MNRAS, 356, 671

Philidet, J., Belkacem, K., & Goupil, M. J. 2021, A&A, 656, A95

Philidet, J., Belkacem, K., & Goupil, M. J. 2022, A&A, 664, A164

Platen, E. 1995, Mathematics and computers in simulation, 38, 69

Pope, S. 1979, Combustion and Flame, 35, 41

Pope, S. B. 1981, Physics of Fluids, 24, 588

Pope, S. B. 1991, Physics of Fluids A, 3, 1947

Pope, S. B. 1994a, Annual Review of Fluid Mechanics, 26, 23

Pope, S. B. 1994b, Physics of Fluids, 6, 973

Pope, S. B. 2000, Turbulent Flows (Springer)

Pope, S. B. & Chen, Y. L. 1990, Physics of Fluids A, 2, 1437

Rauer, H., Aerts, C., Cabrera, J., et al. 2025, Experimental Astronomy, 59, 26

Rincon, F., Barrère, P., & Roudier, T. 2025, A&A, 696, A143

Rodean, H. C. 1996, Stochastic Lagrangian models of turbulent diffusion, Vol. 45 (Springer)

Rodi, W. 1980, State of the Arts Paper, IAHR

Rosenthal, C. S., Christensen-Dalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689

Roxburgh, I. W. & Vorontsov, S. V. 2003, A&A, 411, 215

Samadi, R., Belkacem, K., & Sonoi, T. 2015, in EAS Publications Series, Vol. 73-74, EAS Publications Series, 111–191

Samadi, R. & Goupil, M. J. 2001, A&A, 370, 136

Samadi, R., Goupil, M. J., & Lebreton, Y. 2001, A&A, 370, 147

Schou, J. & Birch, A. C. 2020, A&A, 638, A51

Sonoi, T., Belkacem, K., Dupret, M. A., et al. 2017, A&A, 600, A31

Taylor, G. I. 1922, Proceedings of the London mathematical society, 2, 196

Unno, W. 1967, PASJ, 19, 140

Van Slooten, P. R. & Jayesh, Pope, S. B. 1998, Physics of Fluids, 10, 246

Welton, W. C. 1998, Journal of Computational Physics, 139, 410

Welton, W. C. & Pope, S. B. 1997, Journal of Computational Physics, 134, 150

Xiong, D. R., Cheng, Q. L., & Deng, L. 1997, ApJS, 108, 529

Xiong, D. R., Cheng, Q. L., & Deng, L. 2000, MNRAS, 319, 1079

Yeung, P. K. & Pope, S. B. 1989, Journal of Fluid Mechanics, 207, 531

Zahn, J. P. 1991, A&A, 252, 179

Zhou, Y., Asplund, M., & Collet, R. 2019, ApJ, 880, 13
# Appendix A: Derivation of the Lagrangian stochastic model

We proceed in two steps: first, in App. A.1, we assume a fixed turbulent frequency $\omega$ for all the fluid particles, and we derive the coefficients in the velocity and energy equations Eqs. 2 and 3. Then, in App. A.2, we introduce fluctuations of the turbulent frequency, derive the coefficients in its equation Eq. 4, and determine how the coefficients in the velocity and energy equations should be modified as a result.

## A.1. Lagrangian model with fixed turbulent frequency

In order to determine the coefficients in Eqs. 2 and 3, we first derive the equivalent Fokker-Planck equation (i.e. the transport equation for the joint velocity-energy PDF of the flow) in App. A.1.1. Then, in App. A.1.2, we deduce the equivalent transport equations for all the means (mean density, velocity and energy) and all the second-order moments (Reynolds stress tensor, internal energy variance and internal energy flux) of the turbulent flow. Finally, in App. A.1.3, we compare them to the exact equations derived from first principles, which constrains the form of the coefficients that should be adopted in our Lagrangian stochastic model.

### A.1.1. Equivalent Fokker-Planck equation

Formally, the system of stochastic differential equations can be written

$$dX_i = a_i(\mathbf{X}, t) dt + b_{ij}(\mathbf{X}, t) dW_j, \quad (A.1)$$

where $\mathbf{X}$ is the multivariate stochastic process composed of the position $\mathbf{x}^*$ of the particle, its velocity $\mathbf{u}^*$, and its specific internal energy $e^*$. The equivalent Fokker-Planck equation is

$$\frac{\partial f_L^*}{\partial t} = -\frac{\partial A_i f_L^*}{\partial X_i} + \frac{1}{2} \frac{\partial^2 B_{ij} f_L^*}{\partial X_i \partial X_j}, \quad (A.2)$$

where $A_i = a_i$, $B_{ij} = b_{ik}b_{kj}$. Using Eqs. 1, 2 and 3, and assuming that the drift coefficients in the latter two can be put in the general form of Eqs. 5 and 6, the Lagrangian Fokker-Planck equation reads

$$\frac{\partial f_L^*}{\partial t} = -\frac{\partial}{\partial x_i^*} (u_i^* f_L^*) - \frac{\partial}{\partial u_i^*} \left( \left[ G_{0i} + G_{ij} (u_j^* - \tilde{u}_j) + G_{ei} (e^* - \tilde{e}) \right] f_L^* \right) \\ - \frac{\partial}{\partial e^*} \left( \left[ K_0 + K_j (u_j^* - \tilde{u}_j) + K_e (e^* - \tilde{e}) \right] f_L^* \right) \\ + \frac{1}{2} \frac{\partial^2}{\partial u_i^* \partial u_i^*} (b_u^2 f_L^*) + \frac{1}{2} \frac{\partial^2}{\partial e^* \partial e^*} (b_e^2 f_L^*) . \quad (A.3)$$

The function $f_L^*(\mathbf{x}^*, \mathbf{u}^*, e^*, t|\mathbf{x}_0, t_0)$ is defined as the PDF of the position, velocity and specific energy at time $t$ of the fluid particle which was located at $\mathbf{x}_0$ at the initial time $t_0$. Since the marginal PDF of the initial position of the particle is simply equal to the initial gas density $\rho(\mathbf{x}_0, t_0)$, the non-conditional Lagrangian PDF $f_L(\mathbf{x}^*, \mathbf{u}^*, e^*, t)$ is given by

$$f_L(\mathbf{x}^*, \mathbf{u}^*, e^*, t) = \int d^3\mathbf{x}_0 \rho(\mathbf{x}_0, t_0) f_L^*(\mathbf{x}^*, \mathbf{u}^*, e^*, t|\mathbf{x}_0, t_0), \quad (A.4)$$

which means that the Lagrangian Fokker-Planck equation on $f_L$ is obtained by multiplying Eq. A.3 by $\rho(\mathbf{x}_0, t_0)$ and integrating over $\mathbf{x}_0$.

It can be shown (e.g. Pope 2000) that the Lagrangian PDF is related to the Eulerian PDF $f_E(\rho, \mathbf{u}, e, t; (\mathbf{x}, t))$ through

$$\int d\rho' \rho' f_E(\rho', \mathbf{u}, e; (\mathbf{x}, t)) = \int d^3\mathbf{x}_0 \rho(\mathbf{x}_0, t_0) f_L^*(\mathbf{x}, \mathbf{u}, e, t|\mathbf{x}_0, t_0) . \quad (A.5)$$

We note that the density variable is described by the Eulerian PDF but not by the Lagrangian PDF, where the information concerning the local density of the flow is contained in the dependence of $f_L^*$ on $\mathbf{x}^*$. The Fokker-Planck equation for the Eulerian PDF of the flow then reads

$$\frac{\partial}{\partial t} \left( \int d\rho' \rho' f_E \right) = -\frac{\partial}{\partial x_i} \left( \int d\rho' \rho' u_i f_E \right) \\ -\frac{\partial}{\partial u_i} \left( \int d\rho' \rho' \left[ G_{0i} + G_{ij} (u_j^* - \tilde{u}_j) + G_{ei} (e^* - \tilde{e}) \right] f_E \right) \\ -\frac{\partial}{\partial e} \left( \int d\rho' \rho' \left[ K_0 + K_j (u_j^* - \tilde{u}_j) + K_e (e^* - \tilde{e}) \right] f_E \right) \\ + \frac{1}{2} \frac{\partial^2}{\partial u_i \partial u_i} \left( \int d\rho' \rho' b_u^2 f_E \right) + \frac{1}{2} \frac{\partial^2}{\partial e^2} \left( \int d\rho' \rho' b_e^2 f_E \right) . \quad (A.6)$$

### A.1.2. Equivalent mean equations

By definition, the mean density, mean density-weighted velocity, mean density-weighted specific internal energy, Reynolds-stress tensor, energy variance, and internal energy flux are given by

$$\bar{\rho} = \int d\rho \, d^3\mathbf{u} \, de \, \rho f_E, \quad (A.7)$$

$$\tilde{\mathbf{u}} = \frac{1}{\bar{\rho}} \int d\rho \, d^3\mathbf{u} \, de \, \rho \, \mathbf{u} f_E, \quad (A.8)$$

$$R_{ij} = \frac{1}{\bar{\rho}} \int d\rho \, d^3\mathbf{u} \, de \, \rho (u_i - \tilde{u}_i) (u_j - \tilde{u}_j) f_E, \quad (A.9)$$

$$k_e = \frac{1}{\bar{\rho}} \int d\rho \, d^3\mathbf{u} \, de \, \rho (e - \tilde{e})^2 f_E, \quad (A.10)$$

$$\mathbf{F}^{(e)} = \frac{1}{\bar{\rho}} \int d\rho \, d^3\mathbf{u} \, de \, \rho (u_i - \tilde{u}_i) (e - \tilde{e}) f_E . \quad (A.11)$$

Transport equations for each of these are obtained by multiplying Eq. A.6 by the appropriate function of $\mathbf{u}$ and $e$, and integrating over all values thereof. We find

$$\frac{\partial \bar{\rho}}{\partial t} + \frac{\partial \bar{\rho} \tilde{u}_i}{\partial x_i} = 0, \quad (A.12)$$

$$\frac{\partial \bar{\rho} \tilde{u}_i}{\partial t} + \frac{\partial \bar{\rho} \tilde{u}_i \tilde{u}_j}{\partial x_j} + \frac{\partial \bar{\rho} R_{ij}}{\partial x_j} = \bar{\rho} G_{0i}, \quad (A.13)$$

$$\frac{\partial \bar{\rho} e}{\partial t} + \frac{\partial \bar{\rho} e \tilde{u}_i}{\partial x_i} + \frac{\partial \bar{\rho} F^{(e)}}{\partial x_i} = \bar{\rho} K_0, \quad (A.14)$$

$$\frac{\partial \bar{\rho} R_{ij}}{\partial t} + \frac{\partial \bar{\rho} \tilde{u}_k R_{ij}}{\partial x_k} + \frac{\partial \bar{\rho} u_j'' u_k''}{\partial x_k} + \bar{\rho} R_{ik} \frac{\partial \tilde{u}_j}{\partial x_k} + \bar{\rho} R_{jk} \frac{\partial \tilde{u}_i}{\partial x_k} = \bar{\rho} G_{ik} R_{jk} \\ + \bar{\rho} G_{jk} R_{ik} + \bar{\rho} G_{ei} F_j^{(e)} + \bar{\rho} G_{ej} F_i^{(e)} + \bar{\rho} b_u^2 \delta_{ij}, \quad (A.15)$$
$$ \frac{\partial \bar{p} k_e}{\partial t} + \frac{\partial \bar{p} \tilde{u}_j k_e}{\partial x_j} + \frac{\partial \bar{p} e''^2 \tilde{u}_j''}{\partial x_j} + 2 \bar{p} F_j^{(e)} \frac{\partial \tilde{e}}{\partial x_j} = 2 \bar{p} K_e k_e + 2 \bar{p} K_j F_j^{(e)} + \bar{p} b_e^2, \quad (A.16) $$

and

$$ \frac{\partial \bar{p} F_i^{(e)}}{\partial t} + \frac{\partial \bar{p} \tilde{u}_j F_i^{(e)}}{\partial x_j} + \frac{\partial \bar{p} \tilde{u}_i' e'' \tilde{u}_j''}{\partial x_j} + \bar{p} R_{ij} \frac{\partial \tilde{e}}{\partial x_j} + \bar{p} F_j^{(e)} \frac{\partial \tilde{u}_i}{\partial x_j} = \bar{p} G_{ij} F_j^{(e)} + \bar{p} G_{ei} k_e + \bar{p} K_j R_{ij} + \bar{p} K_e F_i^{(e)}, \quad (A.17) $$

where $\tilde{u}'' \equiv \tilde{u} - \tilde{u}'$ and $e'' \equiv e - e'$.

The diffusion coefficient in the velocity SDE Eq. 2 is related to the Lagrangian structure function

$$ D_L(\Delta t) \equiv \langle |\mathbf{u}^*(t + \Delta t) - \mathbf{u}^*(t)|^2 \rangle \quad (A.18) $$

through $D_L(\Delta t) = b_u^2 \Delta t$. But under the Kolmogorov hypotheses, the structure function is predicted to be

$$ D_L(\Delta t) = C_0 \epsilon \Delta t, \quad (A.19) $$

where $C_0$ is the Kolmogorov constant, and $\epsilon \equiv \omega k$ is the turbulent dissipation rate ($k \equiv R_{ii}/2$ is the turbulent kinetic energy). For the moment, $\omega$ is assumed constant (the effect of a fluctuating turbulent frequency is introduced in App. A.2). The correct structure function is recovered by the Lagrangian stochastic model if the velocity diffusion coefficient is

$$ b_u = \sqrt{C_0 \omega k}. \quad (A.20) $$

By analogy, we also adopt a similar expression for the energy diffusion coefficient

$$ b_e = \sqrt{C_1 \omega k_e}, \quad (A.21) $$

where $C_1$ is another dimensionless constant.

### A.1.3. Constraints from first-principles mean equations

The exact mean equations, derived from the equations of hydrodynamics, are

$$ \frac{\partial \bar{p}}{\partial t} + \frac{\partial \bar{p} \tilde{u}_i}{\partial x_i} = 0, \quad (A.22) $$

$$ \frac{\partial \bar{p} \tilde{u}_i}{\partial t} + \frac{\partial \bar{p} \tilde{u}_i \tilde{u}_j}{\partial x_j} + \frac{\partial \bar{p} R_{ij}}{\partial x_j} = - \frac{\partial \bar{p}}{\partial x_i} + g_i, \quad (A.23) $$

$$ \frac{\partial \bar{p} \tilde{e}}{\partial t} + \frac{\partial \bar{p} \tilde{e} \tilde{u}_i}{\partial x_i} + \frac{\partial \bar{p} F_i^{(e)}}{\partial x_i} = - \bar{p} \frac{\partial \tilde{u}_i}{\partial x_i} - \bar{p} \frac{\partial \tilde{u}_i''}{\partial x_i} - \overline{p' \frac{\partial \tilde{u}_i''}{\partial x_i}} + \bar{p} \kappa \frac{\partial^2 \tilde{e}}{\partial x_i \partial x_i}, \quad (A.24) $$

$$ \frac{\partial \bar{p} R_{ij}}{\partial t} + \frac{\partial \bar{p} \tilde{u}_k R_{ij}}{\partial x_k} + \frac{\partial \bar{p} \tilde{u}_i' \tilde{u}_j'' \tilde{u}_k''}{\partial x_k} + \bar{p} R_{ik} \frac{\partial \tilde{u}_j}{\partial x_k} + \bar{p} R_{jk} \frac{\partial \tilde{u}_i}{\partial x_k} \\ = - \left( \overline{u_i'' \frac{\partial \bar{p}}{\partial x_j}} + \overline{u_j'' \frac{\partial \bar{p}}{\partial x_i}} + \overline{u_j'' \frac{\partial \bar{p}'}{\partial x_i}} + \overline{u_i'' \frac{\partial \bar{p}'}{\partial x_j}} \right) \\ + \overline{u_j'' \frac{\partial \sigma_{ik}}{\partial x_k}} + \overline{u_i'' \frac{\partial \sigma_{jk}}{\partial x_k}} + \overline{u_j'' \frac{\partial \sigma_{ik}'}{\partial x_k}} + \overline{u_i'' \frac{\partial \sigma_{jk}'}{\partial x_k}}, \quad (A.25) $$

$$ \frac{\partial \bar{p} k_e}{\partial t} + \frac{\partial \bar{p} \tilde{u}_i k_e}{\partial x_i} + \frac{\partial \bar{p} e''^2 \tilde{u}_i''}{\partial x_i} + 2 \bar{p} F_i^{(e)} \frac{\partial \tilde{e}}{\partial x_i} = -2 \overline{e'' \bar{p}} \frac{\partial \tilde{u}_i}{\partial x_i} \\ - 2 \overline{p e'' \tilde{u}_i''}{\partial x_i} - 2 \overline{p' e'' \tilde{u}_i''}{\partial x_i} - 2 \overline{e'' p' \tilde{u}_j''}{\partial x_j} + 2 \overline{e'' \sigma_{ij} \tilde{u}_i''}{\partial x_j} \\ + 2 \overline{e'' \sigma_{ij} \tilde{u}_i''}{\partial x_j} + 2 \overline{\sigma_{ij} e'' \tilde{u}_i''}{\partial x_j} + 2 \overline{e'' \sigma_{ij} \tilde{u}_i''}{\partial x_j}, \quad (A.26) $$

and

$$ \frac{\partial \bar{p} F_i^{(e)}}{\partial t} + \frac{\partial \bar{p} \tilde{u}_j F_i^{(e)}}{\partial x_j} + \frac{\partial \bar{p} e'' \tilde{u}_i'' \tilde{u}_j''}{\partial x_j} + \bar{p} R_{ij} \frac{\partial \tilde{e}}{\partial x_j} + \bar{p} F_j^{(e)} \frac{\partial \tilde{u}_i}{\partial x_j} \\ = - \overline{e'' \bar{p}} \frac{\partial \tilde{p}}{\partial x_i} - \overline{e'' \frac{\partial \tilde{p}'}{\partial x_i}} + \overline{e'' \frac{\sigma_{ij}}{\partial x_j}} + \overline{e'' \frac{\partial \sigma_{ij}}{\partial x_j}} \\ - \overline{p \tilde{u}_i'' \frac{\partial \tilde{u}_j}{\partial x_j}} - \overline{p' \tilde{u}_i'' \frac{\partial \tilde{u}_j}{\partial x_j}} - \overline{p \tilde{u}_i'' \frac{\partial \tilde{u}_j''}{\partial x_j}} - \overline{p' \tilde{u}_i'' \frac{\partial \tilde{u}_j''}{\partial x_j}} \\ + \overline{u_i'' \sigma_{jk} \frac{\partial \tilde{u}_j}{\partial x_k}} + \overline{\sigma_{jk} \tilde{u}_i'' \frac{\partial \tilde{u}_j''}{\partial x_k}} + \overline{u_i'' \sigma_{jk}' \frac{\partial \tilde{u}_j}{\partial x_k}} + \overline{u_i'' \sigma_{jk}' \frac{\partial \tilde{u}_j''}{\partial x_k}}, \quad (A.27) $$

where $\kappa$ is the radiative diffusivity, and $\sigma_{ij}$ is the viscous tensor. While the mean radiative flux appears in Eq. A.24, we assumed its correlation with the velocity and energy fluctuations to be negligible, so that it does not appear in Eqs. A.26 and A.27.

The continuity equation is exactly recovered by the Lagrangian stochastic model (i.e. Eqs. A.12 and A.22 are identical). This is because in a Lagrangian description, mass is naturally conserved, and there is no need for a dedicated equation to impose that constraint. By the same token, all the advection terms are recovered exactly (i.e. the left-hand sides of Eqs. A.13, A.14, A.15, A.16 and A.17 are identical to those of Eqs. A.23, A.24, A.25, A.26 and A.27).

In the following, we make three important assumptions. First, we will assume that the turbulence is only weakly compressible, so that we may neglect the correlation of $\partial \tilde{u}_i''/\partial x_i$ with other turbulent fluctuations like $p'$ or $e''$. Second, we will assume that the viscous tensor can be neglected everywhere except when correlated with the shear tensor. This defines the turbulent dissipation rate $\epsilon$, through

$$ \sigma'_{ik} \frac{\partial \tilde{u}_j''}{\partial x_k} = - \tilde{u}_j'' \frac{\partial \sigma'_{ik}}{\partial x_k} = \frac{1}{3} \bar{p} \epsilon \delta_{ij}. \quad (A.28) $$

Finally, we assume that the fluctuations of density $\rho'$, gas pressure $p'$ and specific internal energy $e''$ are related to each other through a polytropic relation

$$ \frac{p'}{\bar{p}} = n \frac{\rho'}{\bar{\rho}} = \frac{n}{n-1} \frac{\rho e''}{\bar{\rho} e}, \quad (A.29) $$

where $n$ is the polytropic index. This polytropic relation directly entails

$$ \overline{u_i''} = - \frac{\overline{\rho' u_i''}}{\overline{\rho}} = - \frac{1}{n-1} \frac{F_i^{(e)}}{\overline{e}}, \quad (A.30) $$

$$ \overline{e''} = - \frac{\overline{\rho' e''}}{\overline{\rho}} = - \frac{1}{n-1} \frac{k_e}{\overline{e}}, \quad (A.31) $$

$$ \overline{p' u_i''} = \frac{n}{n-1} \frac{\bar{p}}{\overline{e}} F_i^{(e)} = \frac{n(\Gamma_1 - 1)}{n-1} \bar{p} F_i^{(e)}, \quad (A.32) $$

$$ \overline{p' e''} = \frac{n}{n-1} \frac{\bar{p}}{\overline{e}} k_e = \frac{n(\Gamma_1 - 1)}{n-1} \bar{p} k_e, \quad (A.33) $$
where we used the Reynolds average of the ideal gas law

$$ \overline{p} = (\Gamma_1 - 1) \overline{\rho e}, \qquad (A.34) $$

and $\Gamma_1$ is the first adiabatic exponent.

Under these assumptions, comparing the right-hand sides of Eqs. A.13, A.14, A.15, A.16 and A.17 with those of Eqs. A.23, A.24, A.25, A.26 and A.27 provides the following constraints for the coefficients of the Lagrangian stochastic model

$$ \overline{\rho} G_{0i} = -\frac{\partial \overline{p}}{\partial x_i} + \overline{\rho} g_i, \qquad (A.35) $$

$$ \overline{\rho} K_0 = -\overline{p} \frac{\partial \overline{u}_i}{\partial x_i} + \frac{\Gamma_1 - 1}{n-1} \left( \frac{\partial \overline{\rho} F_i^{(e)}}{\partial x_i} - \frac{\overline{\rho} F_i^{(e)}}{\overline{p}} \frac{\partial \overline{p}}{\partial x_i} \right) + \overline{\rho} \kappa \frac{\partial^2 \overline{e}}{\partial x_i \partial x_i} + \overline{\rho} \omega k, \qquad (A.36) $$

$$ \begin{aligned} & \overline{\rho} G_{ik} R_{jk} + \overline{\rho} G_{jk} R_{ik} + \overline{\rho} G_{ei} F_j^{(e)} + \overline{\rho} G_{ej} F_i^{(e)} + \overline{\rho} C_{0\omega} k \delta_{ij} \\ & = \frac{1}{\overline{e}(n-1)} \left( F_i^{(e)} \frac{\partial \overline{p}}{\partial x_j} + F_j^{(e)} \frac{\partial \overline{p}}{\partial x_i} \right) - \overline{u}_i'' \frac{\partial p'}{\partial x_j} - \overline{u}_j'' \frac{\partial p'}{\partial x_i} - \frac{2}{3} \overline{\rho} \omega k \delta_{ij}, \end{aligned} \qquad (A.37) $$

$$ \begin{aligned} & 2\overline{\rho} K_e k_e + 2\overline{\rho} K_i F_i^{(e)} + \overline{\rho} C_{1\omega} k_e \\ & = -2(\Gamma_1 - 1)\overline{\rho} k_e \frac{\partial \overline{u}_i}{\partial x_i} - 2\overline{\rho} \frac{2\overline{\rho} \omega k}{\overline{e}(n-1)} k_e, \end{aligned} \qquad (A.38) $$

and

$$ \begin{aligned} & \overline{\rho} G_{ij} F_j^{(e)} + \overline{\rho} G_{ei} k_e + \overline{\rho} K_j R_{ij} + \overline{\rho} K_e F_i^{(e)} = \frac{k_e}{\overline{e}(n-1)} \frac{\partial \overline{p}}{\partial x_i} \\ & \qquad - e'' \frac{\partial p'}{\partial x_i} - (\Gamma_1 - 1) \overline{\rho} F_i^{(e)} \frac{\partial \overline{u}_j}{\partial x_j} - \frac{2}{3} \overline{\rho} \omega k \delta_{ij}, \end{aligned} \qquad (A.39) $$

### A.1.4. Choice of velocity and energy coefficients

The coefficients $G_{0i}$ and $K_0$ are uniquely determined by Eqs. A.35 and A.36. By contrast, there are multiple choices for $G_{ij}$, $G_{ei}$, $K_e$ and $K_j$ that satisfy Eqs. A.37, A.38 and A.39.

First, let us treat Eqs. A.37: one possible solution is to split it into the classical constraint (e.g. Pope 1994b)

$$ \begin{aligned} & \overline{\rho} G_{ik} R_{jk} + \overline{\rho} G_{jk} R_{ik} + \overline{\rho} C_{0\omega} k \delta_{ij} \\ & = -\overline{u}_i'' \frac{\partial p'}{\partial x_j} - \overline{u}_j'' \frac{\partial p'}{\partial x_i} - \frac{2}{3} \overline{\rho} \omega k \delta_{ij}, \end{aligned} \qquad (A.40) $$

and the following constraint for $G_{ei}$

$$ \overline{\rho} G_{ei} F_j^{(e)} + \overline{\rho} G_{ej} F_i^{(e)} = \frac{1}{\overline{e}(n-1)} \left( F_i^{(e)} \frac{\partial \overline{p}}{\partial x_j} + F_j^{(e)} \frac{\partial \overline{p}}{\partial x_i} \right). \qquad (A.41) $$

Taking the trace of Eq. A.40, we get

$$ 2\overline{\rho} G_{ij} R_{ij} + 3\overline{\rho} C_{0\omega} k = -2\overline{u}_i'' \frac{\partial p'}{\partial x_i} - 2\overline{\rho} \omega k. \qquad (A.42) $$

Under our previous assumptions, the velocity-pressure-gradient tensor can be neglected, because in weakly compressible turbulence, the divergence of the acoustic flux $\partial \overline{p} \overline{u}_i'' / \partial x_i$ is a minor

source of energy transport, and we already neglected $\overline{p}' \overline{\partial u}_i'' / \partial x_i$. The simplest choice of $G_{ij}$ is to assume that it is isotropic, in which case the only possibility is

$$ G_{ij} = -\left(\frac{1}{2} + \frac{3}{4}C_0\right) \omega \delta_{ij}, \qquad (A.43) $$

which corresponds to the Simplified Langevin Model (e.g. Pope 1994b). It can be seen from Eq. A.40 that this is equivalent to the Rotta model for the pressure-rate-of-strain tensor

$$ p' \left( \frac{\partial u_i''}{\partial x_j} + \frac{\partial u_j''}{\partial x_i} \right) = -C_R \left( R_{ij} - \frac{2}{3} k \delta_{ij} \right), \qquad (A.44) $$

with $C_R = 1 + 3C_0/2$. As for Eq. A.41, it straightforwardly gives

$$ G_{ei} = \frac{1}{\overline{\rho e}(n-1)} \frac{\partial \overline{p}}{\partial x_i} = \frac{\Gamma_1 - 1}{n-1} \frac{1}{\overline{p}} \frac{\partial \overline{p}}{\partial x_i}. \qquad (A.45) $$

Let us now turn to Eq. A.38. All the terms are proportional to $k_e$, except for the second term on the left-hand side. A natural choice, then, is

$$ K_i = 0, \qquad (A.46) $$

$$ K_e = -\frac{1}{2} C_1 \omega - (\Gamma_1 - 1) \frac{\partial \overline{u}_i}{\partial x_i} - \frac{\omega k}{\overline{e}(n-1)}. \qquad (A.47) $$

Finally, using these expressions for $G_{ij}$, $G_{ei}$, $K_e$ and $K_i$, we see that Eq. A.39 reduces to the following closure relation for the energy-pressure-gradient tensor

$$ e'' \frac{\partial p'}{\partial x_i} = -C \overline{\rho} \omega F_i^{(e)}, \qquad (A.48) $$

with $C = 3C_0/4 + (1 - C_1)/2$.

### A.2. Lagrangian model including turbulent frequency

We now lift the constraint that the value of the turbulent frequency $\omega$ be fixed, and we add a stochastic equation for $\omega^*$. We follow Van Slooten & Jayesh (1998) and adopt the following stochastic differential equation

$$ d\omega^* = -\Omega(\omega^* - \bar{\omega}) dt - \omega^* \bar{\omega} S_\omega dt + \sqrt{2\sigma^2\omega^*\bar{\omega}\Omega} dW, \qquad (A.49) $$

where the conditional mean turbulent frequency $\Omega$ and the source term $S_\omega$ are given by Eq. 12 and Eq. 13 respectively.

The refined Kolmogorov hypotheses dictate (e.g. Pope & Chen 1990) that the structure function $D_L(\Delta t)$ defined by Eq. A.18 be proportional to the instantaneous turbulent dissipation rate $\epsilon^* = k\omega^*$ instead of the mean dissipation rate. As a result, the diffusion coefficients $b_u$ and $b_e$ in Eq. 2 and Eq. 3 should be given not by Eq. A.20 and Eq. A.21, but by

$$ b_u^2 = C_0 \omega^* k, \qquad (A.50) $$

$$ b_e^2 = C_1 \omega^* k_e. \qquad (A.51) $$

In order for the equivalent Fokker-Planck equation to remain identical, extra terms $A_i$ and $A_e$ should be added respectively to the coefficients $a_i$ and $a_e$ in Eqs. 2 and 3. The corresponding extra terms in the Fokker-Planck Eq. A.3 read

$$ \begin{aligned} \Delta \text{FP} = & -\frac{\partial A_i f_L^*}{\partial u_i^*} + \frac{1}{2} C_0 (\omega^* - \bar{\omega}) k \frac{\partial^2 f_L^*}{\partial u_i^* \partial u_i^*} \\ & - \frac{\partial A_e f_L^*}{\partial e^*} + \frac{1}{2} C_1 (\omega^* - \bar{\omega}) k_e \frac{\partial^2 f_L^*}{\partial e^* \partial e^*}. \end{aligned} \qquad (A.52) $$
Since we want the velocity-energy part of the Fokker-Planck equation to remain unaffected by the modification of $b_u$ and $b_e$, this must be zero. In order to determine the expressions of $A_i$ and $A_e$ leading to $\Delta$FP = 0, we will make the assumption that the local, one-time, marginal ($e''$, $u''$) PDF is joint-normal, with zero mean and a covariance matrix

$$C = \begin{bmatrix} k_e & F_x^{(e)} & F_y^{(e)} & F_z^{(e)} \\ F_x^{(e)} & R_{xx} & R_{xy} & R_{xz} \\ F_y^{(e)} & R_{xy} & R_{yy} & R_{yz} \\ F_z^{(e)} & R_{xz} & R_{yz} & R_{zz} \end{bmatrix}. \qquad (A.53)$$

Then we have

$$\frac{\partial f_L^*}{\partial X_i} = -(C^{-1})_{ij} X_j f_L^*, \qquad (A.54)$$

where $\mathbf{X} \equiv (e'', \mathbf{u}'')$, and Eq. A.52 becomes

$$\Delta \text{FP} = -\frac{\partial}{\partial u_i} \left( A_i f_L^* + \frac{1}{2} C_0 (\omega^* - \bar{\omega}) k (C^{-1})_{ij} X_j f_L^* \right) \\ -\frac{\partial}{\partial e} \left( A_e f_L^* + \frac{1}{2} C_1 (\omega^* - \bar{\omega}) k_e (C^{-1})_{0j} X_j f_L^* \right). \quad (A.55)$$

The condition that this must be zero yields the expression of the extra terms in the velocity and energy SDEs

$$A_i = -\frac{1}{2} C_0 (\omega^* - \bar{\omega}) k \left( (C^{-1})_{i0} (e^* - \bar{e}) + (C^{-1})_{ij} (u_j^* - \bar{u}_j) \right), \quad (A.56)$$

$$A_e = -\frac{1}{2} C_1 (\omega^* - \bar{\omega}) k_e \left( (C^{-1})_{00} (e^* - \bar{e}) + (C^{-1})_{0j} (u_j^* - \bar{u}_j) \right). \quad (A.57)$$

These four extra terms are added respectively to the expressions of $G_{ie}$, $G_{ij}$, $K_e$ and $K_i$ given in App. A.1.4.

## Appendix B: Corrective SPH scheme

The corrective Smoothed Particle Method (CSPM) adopted in our code works in three steps. For a given location $\mathbf{x}_i$ within the domain, we first compute the following sums through the standard SPH scheme Eq. 19

$$W = \bar{1}, \quad W_{\alpha} = \bar{1}^{\alpha}, \quad W_{\alpha\beta} = \bar{1}^{\alpha\beta}, \qquad (B.1)$$

$$X = \bar{x}, \quad X_{\alpha} = \bar{x}^{\alpha}, \quad X_{\alpha\beta} = \bar{x}^{\alpha\beta}, \qquad (B.2)$$

$$Y = \bar{y}, \quad Y_{\alpha} = \bar{y}^{\alpha}, \quad Y_{\alpha\beta} = \bar{y}^{\alpha\beta}, \qquad (B.3)$$

$$XX = \bar{x}^2, \quad XX_{\alpha} = \bar{x}^2, \quad XX_{\alpha\beta} = \bar{x}^2, \qquad (B.4)$$

$$XY = \bar{xy}, \quad XY_{\alpha} = \bar{xy}^{\alpha}, \quad XY_{\alpha\beta} = \bar{xy}^{\alpha\beta}, \qquad (B.5)$$

$$YY = \bar{y}^2, \quad YY_{\alpha} = \bar{y}^2, \quad YY_{\alpha\beta} = \bar{y}^2, \qquad (B.6)$$

where

$$\bar{Q} \equiv \sum_{j=1}^{N} \mathcal{V}^{* (j)} \mathcal{Q}^{* (j)} K(\mathbf{x}_j - \mathbf{x}_i), \qquad (B.7)$$

$$\bar{Q}^{\alpha} \equiv \sum_{j=1}^{N} \mathcal{V}^{* (j)} \mathcal{Q}^{* (j)} \left. \frac{\partial K}{\partial x_{\alpha}} \right|_{\mathbf{x}_j - \mathbf{x}_i}, \qquad (B.8)$$

$$\bar{Q}^{\alpha\beta} \equiv \sum_{j=1}^{N} \mathcal{V}^{* (j)} \mathcal{Q}^{* (j)} \left. \frac{\partial^2 K}{\partial x_{\alpha} \partial x_{\beta}} \right|_{\mathbf{x}_j - \mathbf{x}_i}. \qquad (B.9)$$

They only need to be computed once per time step. In parallel, we compute the mean of the quantity that we want to average, as well as its gradient and second order derivatives, through the standard SPH scheme Eqs. B.7, B.8 and B.9. They need to be computed separately for each mean flow quantity.

This first step requires the kernel estimation of means that are not density-weighted, so that we need to estimate the lumped volume $\mathcal{V}^{*(k)}$ occupied by each notional particle. To do this, we do a first estimation of the mean density $\bar{\rho}^{(k)}$ at each particle location, and we assign them the volume $\mathcal{V}^{*(k)} = \Delta m / \bar{\rho}^{(k)}$.

In a second step, we correct the mean and gradient through

$$\left[ \frac{\bar{Q}_{\text{CSPM}}}{\bar{Q}_{\text{CSPM}}^{y}} \quad \frac{\bar{Q}}{\bar{Q}^{x}} \quad \frac{\bar{Q}}{\bar{Q}^{y}} \right] = M_{\text{grad}}^{-1} \left[ \frac{\bar{Q}}{\bar{Q}^{y}} \right], \qquad (B.10)$$

where

$$M_{\text{grad}} = \begin{bmatrix} W & X - x_i W & Y - y_i W \\ W_x & X_x - x_i W_x & Y_x - y_i W_x \\ W_y & X_y - x_i W_y & Y_y - y_i W_y \end{bmatrix}. \qquad (B.11)$$

In a third step, we correct the hessian matrix of $\bar{Q}$ through

$$\begin{bmatrix} \bar{Q}_{\text{CSPM}}^{xx} \\ \bar{Q}_{\text{CSPM}}^{xy} \\ \bar{Q}_{\text{CSPM}}^{yy} \end{bmatrix} = M_{\text{hess}}^{-1} \begin{bmatrix} \bar{Q}^{xx} - \bar{Q}_{\text{CSPM}} W_{xx} - \bar{Q}_{\text{CSPM}}^{x} dX_{xx} - \bar{Q}_{\text{CSPM}}^{y} dY_{xx} \\ \bar{Q}^{xy} - \bar{Q}_{\text{CSPM}} W_{xy} - \bar{Q}_{\text{CSPM}}^{x} dX_{xy} - \bar{Q}_{\text{CSPM}}^{y} dY_{xy} \\ \bar{Q}^{yy} - \bar{Q}_{\text{CSPM}} W_{yy} - \bar{Q}_{\text{CSPM}}^{x} dX_{yy} - \bar{Q}_{\text{CSPM}}^{y} dY_{yy} \end{bmatrix}, \qquad (B.12)$$

where

$$M_{\text{hess}} = \frac{1}{2} \begin{bmatrix} dXX_{xx} & dXY_{xx} & dYY_{xx} \\ dXX_{xy} & dXY_{xy} & dYY_{xy} \\ dXX_{yy} & dXY_{yy} & dYY_{yy} \end{bmatrix} - \frac{1}{2} \begin{bmatrix} W_{xx} & dX_{xx} & dY_{xx} \\ W_{xy} & dX_{xy} & dY_{xy} \\ W_{yy} & dX_{yy} & dY_{yy} \end{bmatrix} M_{\text{grad}}^{-1} \begin{bmatrix} dXX & dXY & dYY \\ dXX_x & dXY_x & dYY_x \\ dXX_y & dXY_y & dYY_y \end{bmatrix}, \qquad (B.13)$$

and we have defined

$$dX_{\alpha} = X_{\alpha} - x_i W_{\alpha}, \qquad (B.14)$$

$$dY_{\alpha} = Y_{\alpha} - y_i W_{\alpha}, \qquad (B.15)$$

$$dXX_{\alpha} = XX_{\alpha} - 2x_i X_{\alpha} + x_i^2 W_{\alpha}, \qquad (B.16)$$

$$dXY_{\alpha} = XY_{\alpha} - x_i Y_{\alpha} - y_i X_{\alpha} + x_i y_i W_{\alpha}, \qquad (B.17)$$

$$dYY_{\alpha} = YY_{\alpha} - 2y_i Y_{\alpha} + y_i^2 W_{\alpha}. \qquad (B.18)$$