Surface effects and turbulent pressure. Assessing the Gas-$\Gamma_1$ and Reduced-$\Gamma_1$ empirical models

The use of the full potential of stellar seismology is made difficult by the improper modeling of the upper-most layers of solar-like stars and their influence on the modeled frequencies. Our knowledge on these \emph{surface effects} has improved thanks to the use of 3D hydrodynamical simulations but the calculation of eigenfrequencies relies on empirical models for the description of the Lagrangian perturbation of turbulent pressure: the reduced-$\Gamma_1$ model (RGM) and the gas-$\Gamma_1$ model (GGM). Starting from the fully compressible turbulence equations, we derive both the GGM and RGM models using a closure to model the flux of turbulent kinetic energy. It is found that both models originate from two terms: the source of turbulent pressure due to compression produced by the oscillations and the divergence of the flux of turbulent pressure. It is also demonstrated that they are both compatible with the adiabatic approximation but also imply a number of questionable assumptions mainly regarding mode physics. Among others hypothesis, one has to neglect the Lagrangian perturbation of the dissipation of turbulent kinetic energy into heat and the Lagrangian perturbation of buoyancy work.


Introduction
Systematic differences between observed and modeled eigenfrequencies is a long-standing problem in stellar seismology. For the Sun and solar-like stars, this has been identified and recognized to be associated to the improper modeling of their uppermost layers (e.g. Brown 1984; Dziembowski et al. 1988;Christensen-Dalsgaard & Thompson 1997). To circumvent this problem, well-chosen combinations of frequencies (e.g. Roxburgh & Vorontsov 2003) or ad-hoc corrections (Kjeldsen et al. 2008;Ball & Gizon 2014;Sonoi et al. 2015) are commonly used. Nevertheless, to be able to exploit all the information contained in the observed frequencies, the physics underlying surface effects must be understood.
To do so, Rosenthal et al. (1999) used 3D hydrodynamical simulations. This allowed to account for the mean turbulent pressure (as well as convective backwarming, see Trampedach et al. 2017) in the equilibrium structure. Their work was followed by Piau et al. (2014) Mosumgaard et al. (2020); Jørgensen et al. (2021) for solar-like stars. A drawback of this approach is the need to compute adiabatic eigenfrequencies using empirical models to describe the Lagrangian perturbation of turbulent pressure.
To this end, Rosenthal et al. (1999) introduced two empirical models. The reduced-Γ 1 model (RGM), which consists in arguing that the perturbation of turbulent pressure is negligible com-pared to the perturbation of gas pressure, and the gas-Γ 1 model (GGM), which consists in considering that perturbations of gas pressure and turbulent pressure are in phase. The authors favored the GGM as it permits them to better reproduce the observed solar frequencies. Most of the above-mentioned works then used the GGM, except for a few (e.g. Jørgensen & Weiss 2019) who considered the RGM guided by the non-adiabatic calculation by Houdek et al. (2017). However, the latter was recently challenged by Schou & Birch (2020) using the eigenfunctions as inferred directly from 3D numerical simulations. Therefore, the issue of surface effects remains and demands extensive theoretical investigation to gain insight into the physics of the problem.
Hence, the recent progress on surface effects relies on the use of empirical models for which the physics is not understood well enough. In this article, we aim at deriving and then assessing the theoretical validity of the GGM and RGM empirical models.

Equation governing the turbulent pressure
When including the mean turbulent pressure in the adiabatic oscillation equations, the Lagrangian perturbation of turbulent pressure is to be prescribed. This is the role of the GGM and RGM models, which are extensively introduced in Appendix A. To derive them, a first step is to express the equation for turbulent pressure.
where the overbar denotes the Reynolds average and the tilde denotes the Favre average (see Appendix A for a definition), u i the i component of the velocity field, P g is the gas pressure (it includes contributions from the radiative field as well as body forces), ρ is the density, τ jk is the viscous stress tensor, and r i j ≡ ρu i u j are the Reynolds stresses where u i is the i-th component of the velocity fluctuation around its Favre average (ũ i ). We also employ the notation ∂ t ≡ ∂/∂t and ∂ i ≡ ∂/∂x i , Einstein's notation for repeated indices, and the pseudo-Lagrangian derivative is defined by d/dt ≡ ∂/∂t + u j ∂/∂x j . We then consider the rr component of Eq.
(1) (where r is the radial coordinate) and identify the Reynolds average with the horizontal average. We obtain 1 2 where P t ≡ ρu r u r is the turbulent pressure, F turb rrr ≡ ρu r u r u r /2, F visc rrr ≡ −u r τ rr , and F p r = u r P g . Except for notational differences, Eq. (2) is strictly equivalent to Eq. (16b) in Canuto (1997). The first term of the right-hand side of Eq. (2) is a source of turbulent pressure due to compression produced by radial oscillations. The second term corresponds to the pressure work (or the buoyancy work), which is a source of turbulent pressure in convective regions. The three following terms in brackets correspond to transport terms. Finally, the last two terms are dissipative terms.
Several assumptions are now needed. The molecular diffusion flux (F visc rrr ) is neglected compared to the other fluxes appearing in Eq. (2). This is justified as we are considering fully developed turbulence with high Reynolds numbers (e.g. Nordlund et al. 2009). The last term of Eq. (2) is assumed to be proportional to the rate of dissipation of turbulent kinetic energy into heat, i.e. τ rk ∂ k u r ∝ τ i j ∂ i u j . This standard assumption is made possible because dissipation by molecular forces occurs at almost isotropic small scales (e.g. Pope 2000). Finally, the pressure-dilatation term (P g ∂ r u r ) is neglected because it scales as the square of the turbulent Mach number (Sarkar 1992) (the maximum is about 0.3 in the Sun). By adopting those approximations, Eq. (2) reduces to where ρ ≡ τ i j ∂ i u j is the dissipation rate of turbulent kinetic energy into heat.

Modelling the transport of turbulent pressure
The flux of turbulent pressure, F turb rrr in Eq. (3), is difficult to model in the uppermost layers of solar-like stars for which the down-gradient approximation fails (see Canuto 2009;Kupka & Muthsam 2017). To overcome this issue, we adopt the closure initially proposed by Canuto (1992) (see also Canuto 2011) and supported by 3D direct numerical simulations (Kupka & Muthsam 2007). It reads where F turb iir is the flux of turbulent kinetic energy, F r ≡ ρ u r is the flux of turbulent dissipation, k ≡ r ii /(2ρ) is the specific turbulent kinetic energy, and c is a parameter to be specified.
From Eq. (4), we then obtain where ω ≡ /k is the turbulence frequency and The parameter α, which can be understood as the degree of anisotropy of the flux of turbulent kinetic energy, will be supposed to be known and obtained directly from the solar 3D numerical simulation (see Sect. 3.2). Now, using Eq. (5), Eq.
(3) can be rewritten where one still needs a prescription for the flux of turbulent dissipation.
The first term in the right-hand side of Eq. (8) is related to the production of turbulent kinetic energy by compression. The second term represents for the effect of bulk compressions and expansions onto (Coleman & Mansour 1991). The third term is related to production of turbulent kinetic energy by the buoyancy work, and the last term is related to both the viscous destruction and production due to vortex stretching (Gatski & Bonnet 2013). The coefficients c 1 , c 2 , c 3 , and c 4 will be discussed in the following section.

Recovering the gas-Γ 1 and reduced-Γ 1 models
We will now derive an expression for the Lagrangian perturbation of turbulent pressure and make a number of assumptions whose validity will be discussed.

Perturbation of turbulent pressure
Using Eq. (7) together with Eq. (8), we obtain the desired expression for the equation governing turbulent pressure where we used the averaged continuity equation (Eq. A.5), and Φ ≡ 2r ii /P t is the anisotropy factor.
To go further, we assume the following hypothesis: (H1) the Lagrangian perturbation of the pressure-velocity fluctuations is neglected, i.e. δF p r = 0. (H2) the turbulence frequency (ω) is supposed to vary on a length-scale much larger than the length-scale of F iir . Accordingly, the second term of the right-hand-side of Eq. (5) can be neglected and thus last term of Eq. (9) vanishes. (H3) the Lagrangian perturbation of α is neglected, i.e. δα = 0. (H4) the Lagrangian perturbation of the dissipation of turbulent kinetic energy into heat (δ ) is neglected. (H5) the perturbation of the buoyancy work, δ(u r ∂ r P g ), is neglected. (H6) the Lagrangian perturbation of density is assumed to be real.
Perturbing Eq. (9) and applying H1 to H6, one gets where δP t and δρ are the pseudo-Lagrangian perturbations of turbulent pressure and density, respectively. Anticipating on the next section, one can already mention that Eq. (10) allows one to recover the RGM and GGM models.
Tracing back the origin of Eq. (10) shows that it results from two terms in Eq. (2): the source of turbulent pressure due to compression produced by the oscillations and the divergence of the flux of turbulent pressure. It is also important to mention that H1 to H6 are fully consistent with the assumptions made to obtain the equation for the perturbation of gas pressure (δP g , see A.9 and Appendix B). Moreover, given H1, H3 and H6 and further neglecting the perturbation of the radiative flux, H5 is equivalent to adopting the adiabatic limit (see Appendix C). Hence, both equations for δP t (Eq. 10) and δP g (Eq. A.9) are compatible with the adiabatic approximation. Conversely, the adiabatic approximation only is not suffisent to derive Eq. (10) and Eq. (A.9).
The first assumption, H1, is essentially equivalent to neglecting the perturbation of the convective flux (δF conv r ) because F p r is proportional to the convective flux as shown by Canuto (1997) (see his Eq.35c). Such an assumption is however not strictly valid because the perturbation of the convective flux does not vanish even in the adiabatic limit (e.g. Sonoi et al. 2017). H2 is difficult to properly assess because one needs to determine and the simplest possible way is to consider ∝ k 3/2 /H p (e.g. Pope 2000), where the dissipation length-scale is assumed to scale as the pressure scale-height (H p ). Using, the solar 3D simulation described in Belkacem et al. (2019) one readily finds that H2 is valid near the super-adiabatic peak but not in the quasi-adiabatic regions. For H3, it is equivalent to assuming that perturbations of horizontal and vertical contributions of turbulent kinetic energy adjust instantaneously to each other and are thus in phase. A look at the third-order equation on fluxes (see Canuto 1997) demonstrates that the situation is much more complex because many terms are capable to introduce some redistribution, and in turn phase shifts, when perturbed. Concerning H4 and H5, it consists in neglecting the Lagrangian perturbation of both the turbulent dissipation and the buoyancy work. For mode damping, these two contributions exactly compensate the contribution of turbulent pressure in the limit of a vanishing flux of turbulent kinetic energy and with Γ 3 −1 = 2/3 (where Γ 3 ≡ (∂ ln T/∂ ln ρ) s ) (e.g. Ledoux & Walraven 1958;Grigahcène et al. 2005). More recently, Belkacem et al. (2019) demonstrated using the normal modes of a 3D solar hydrodynamical simulation that the perturbation of both terms plays an essential role for the mode damping rates. Hence, as modal surface effects also partly rely on the phase differences between the perturbations of density and turbulent pressure, the impacts of those contributions to surface effects are definitively to be assessed. Finally, H6 assumes that adding the mean turbulent pressure to the hydrostatic equilibrium only introduces a negligible phase shift to the perturbation of density in the adiabatic limit.
These remarks lead to question the validity of both the RGM and GGM models because, even in the adiabatic limit, they introduces oversimplifying hypotheses regarding the properties of turbulent convection and mode physics.

Application to solar p-mode frequencies
We will now investigate how Eq. (10) permits us to recover the GGM and RGM models. A prerequisite is to specify the coefficients c 1 and c 2 . A commonly accepted value for the first coefficient is c 1 1.
where n 0.75 is the exponent of the viscosity law on temperature. Indeed, Eq. (11) aims at accounting for the effect of bulk compressions and expansions onto due to the dependence of the viscosity on temperature (Coleman & Mansour 1991). We note that Eq. (11) has been obtained using the rapid distorsion theory (see Hunt & Carruthers 1990;Cambon et al. 1993, for reviews) and under the assumption of adiabatic compression. In the quasi-adiabatic regions of the solar convection zone, those approximations are relatively applicable because the modal period is much shorter than both the typical turn-over time-scale and the thermal time-scale. However, in the super-adiabatic layers, those time scales become of comparable magnitude and the validity of Eq. (11) becomes questionable. Therefore, the adopted value of c 2 used in this work must be considered a guideline rather than a firm value. We thus computed the frequency differences (for radial modes) between the observed frequencies, taken from Broomhall et al. (2009) and Davies et al. (2014), and theoretical frequencies computed with a classical shooting method. For the latter, we used the RGM, GGM, and the model developed in this work using several values of the closure coefficient c (see Eq. 4). Theoretical frequencies have been obtained by integrating equations Eqs. (A.7), (A.8), (A.9), complemented by either Eq. (A.12) for the RGM, Eq. (A.13) for the GGM, and Eq. (10) for the model proposed in this article. The equilibrium model had been obtained by patching a CESTAM model together with a solar ANTARES 3D simulation (see Belkacem et al. 2019, for details) using the same methodology as described in Sonoi et al. (2015). The 3D model is however not exactly solar because the effective temperature is 5750 ± 18K and the chemical composition is (X, Y, Z) = (0.7373, 0.2427, 0.0200) with Grevesse & Noels (1993) mixture. To allow for a comparison with observed solar frequencies, the 3D model is patched with a CESTAM model with a helium abundance of 0.2485 and the resulting frequencies are rescaled to match with the standard solar frequency (GM /R 3 ) 1/2 /(2π) = 99.85537 µHz (G is the gravitational constant, M the solar mass, and R the solar radius). Such a procedure is sufficient for our purposes, even if a frequency shift remains at low frequencies, because we are interested in differential effects. From this patched model, all the equilibrium quantities have been inferred. Figure 1 shows that both the frequency differences obtained using the GGM and the RGM models can be recovered by adjusting the value of the closure coefficient c introduced in Eq. (4). What guidance for c can we take from existing 3D numerical simulations? A value of c ≈ 0.6 was derived from four direct numerical simulations of fully compressible convection, two of which were discussed in Kupka & Muthsam (2007). Those simulations were done with the ASCIC3 code (Muthsam et al. 1995) which is distinguished for such work by not introducing, directly or indirectly, subgrid scale viscosities: dissipation is by explicit, physical viscosity and a very small contribution from time integration only. Data were collected for a total of 12 different model configurations (five of which had been discussed in detail in Muthsam et al. 1995 andin Muthsam et al. 1999). They covered a range in Prandtl number Pr from 0.1 to 1 at a "zone Rayleigh number" of 10 5 to 10 6 . For a horizontal domain width eight times the size of this zone depth and with 3 to 4 granules found along each horizontal direction, this implies a "granulation diameter based Rayleigh number" ∼ 16 to 50 times larger or an Ra in the range of 10 6 to 5 · 10 7 . This yields a product of Ra and Pr in the range of 10 5 to 5 · 10 7 . That (squared) ratio of the thermal diffusion time scale to the buoyancy time scale agrees with results for the upper part of the solar convection zone (Kupka et al. 2020) despite the convection zones are much more shallow than in the solar case, where additionally Pr 10 −6 . Kupka & Muthsam (2007) found c ≈ 0.6 in all cases where sufficient numerical resolution had been ensured, also for cases not shown therein. Confirming those results for much lower values of Pr and higher ratio of total flux to radiative flux would be useful. Recalling the discussion in Belkacem et al. (2019) on computing the dissipation rate from realistic solar simulations we have to point out here that large eddy simulations, whether they use hyperviscosity, a Smagorinsky-type subgrid scale model, or a Riemann solver, are not the best tools to compute F r or : those quantities depend on viscosity related processes which peak near small multiples of the grid scale. This could heavily bias computations of those quantities by the numerical method used. We thus consider conclusions based on low Pr direct numerical simulations for the physical range of Ra · Pr of interest the safer way to estimate c. In comparison, the ratio α from Eq. (5) can be safely estimated from solar granulation simulations, as the contributions to this quantity peak at length scales resolved in those simulations.

Conclusions
By using the Reynolds and Favre averaged fully compressible Navier-Stokes equations, we have demonstrated that it is possible to develop a model which recovers both the RGM and GGM empirical models. Interestingly enough, this is based on a relation that has been shown to originate from a compensation between the source of turbulent pressure due to compression produced by the oscillations and the divergence of the flux of turbulent pressure. We then showed the RGM and GGM models are compatible with the adiabatic approximation but also imply more drastic and unrealistic physical assumptions regarding turbulence and mode physics.
Comparison with solar frequencies shows that, while recovering the RGM and GGM, the results are sensitive to the closure coefficients (c in Eq. 4 but also c 2 appearing in Eq. 10). To consolidate the value of these parameters, only a direct 3D numerical simulations could help. Unfortunately, this is still outof-reach for the Sun due to our current numerical capacities and  hence extrapolations from more accessible parameter ranges remain necessary.
It is thus difficult to draw a conclusion on which is the more appropriate to use among the two models. Even worse, given the above-mentioned assumptions which are needed to recover these empirical models, one can safely conclude none of them are firmly physically grounded, even in the adiabatic limit. However, one shall still quantify the hypotheses (H1 to H6) on which both RGM and GGM models rely, and more precisely, their individual effect on mode frequencies. This will be made possible by either using a realistic treatment of turbulent convection based on a time-dependent and non-adiabatic treatment or by using normal modes of direct 3D numerical simulations. For the latter, dedicated long-duration simulations (to have a sufficient statistic and to resolve the normal modes) with a large spatial extension (to have a sufficient number of normal modes) need to be computed and this must be done by resolving all spatial scales to obtain an accurate estimate of turbulent dissipation. For the former, a qualitative leap forward is needed because current 1D formalisms based on the mixing-length theory all have their shortcomings among which are free parameters and questionable physical assumptions (see Houdek & Dupret 2015, for details), which prevent them from guaranteeing to properly grasp the physics of the problem.
When turbulent pressure is included in the mean stratification for computing the classical adiabatic oscillations, it is necessary to determine the perturbation of turbulent pressure and subsequently the perturbation of total pressure (e.g. Rosenthal et al. 1999). To illustrate it let us begin by considering the mass and momentum conservation equations that read where P g is the gas pressure, g i is the i-th component of the gravitational acceleration, τ jk is the viscous stress tensor, ρ is the density, and u i the i-th component of the velocity field. We also employ the notation ∂ t ≡ ∂/∂t and ∂ i ≡ ∂/∂x i as well as Einstein's notation for repeated indices. Equations (A.1) and (A.2) are averaged using both a classical Reynolds average and a density-weighted average (also commonly named to as Favre average, see Favre 1969). For a quantity X, the Reynolds average is defined as and the Favre average is defined, for a quantity Y, by As we consider a compressible flow, using a Reynolds average for density and gas pressure and a Favre average for the velocity field greatly simplifies the averaged equations (e.g. Canuto 1997;Nordlund & Stein 2001;Belkacem et al. 2019). Applying this procedure for Eqs. (A.1) and (A.2) gives where r ik ≡ ρ u i u k are the Reynolds stresses (where u k is the k-th component of the velocity fluctuation around its Favre average), the pseudo-Lagrangian derivative is defined by d/dt ≡ ∂/∂t + u j ∂/∂x j . To derive Eqs. (A.5) and (A.6), gravity fluctuations have been neglected and it has been assumed that viscosity does not affect the mean momentum equation and thus the large-scale flow ( u j ). It is important to notice that this approximation does not mean that viscous dissipation is neglected, because it appears in the equations governing the turbulent quantities (Canuto 1997). We now identify the Reynolds average with the horizontal average so that in Eq. (A.6) only the rr component of the Reynolds stress remains. Turbulent pressure is thus defined by P t ≡ ρu r u r . In addition, we split the mean quantities such that X = X 0 + δX (where X 0 ≡< X > t is the time-average of X and δX is the pseudo-Lagrangian perturbation corresponding to the radial oscillations). Therefore, from Eqs. (A.5) and (A.6) we get the desired oscillation equations in the pseudo-Lagrangian frame where σ is the angular frequency, ξ r is the radial component of the eigen-displacement ( u r = iω ξ r ) and δP tot = δP g + δP t is the Lagrangian perturbation of total pressure. In addition, in the adiabatic limit, we also have the thermodynamic relation δP g P g,0 = Γ 1 δρ ρ 0 , (A.9) which has been derived by using a number of approximations that are explicitly stated in Appendix B. Therefore, except for the boundary conditions, the system is not closed because we must specify δP t or equivalently δP tot . Notice that in the following the subscript "0" for denoting the temporally and horizontally averaged quantities will be dropped for ease of notation. To solve Eqs. (A.7) and (A.8), one needs to specify a prescription for the perturbation of total pressure. It is thus necessary to express the perturbation of turbulent pressure with the perturbation of density so that where A is to be determined. Therefore, using Eq. (A.10), one formally writes Two empirical models have then been introduced by Rosenthal et al. (1999). The first is the reduced-Γ 1 model (RGM), which consists in arguing that δP t = 0 or A = 0 so that Eq. (A.11) becomes is called the reduced Γ 1 . This approximation was initially introduced by Rosenthal et al. (1995) based on the observation that some non-local mixing-length theory shows that density and gas pressure perturbations are almost in phase quadrature with the perturbation of turbulent pressure. The authors therefore consider that the real part of δP t can be neglected. The second model is the gas-Γ 1 model (GGM) and has been introduced by Rosenthal et al. (1999). It consists in arguing the opposite, i.e. that perturbation of gas pressure and turbulent pressure are in phase. Hence, using A = Γ 1 together with Eqs. (A.9), (A.10), and (A.11), one obtains δP tot P tot = δP t P t = δP g P g = Γ 1 δρ ρ . (A.13) As can be seen with Eqs. (A.12) and (A.13), both the RGM and GGM present the major advantage of being easily implemented for computing adiabatic oscillation while including the effect of turbulent pressure, provided that the mean turbulent pressure is prescribed. For GGM it is sufficient to replace the gas pressure by the total pressure in the classical adiabatic oscillation equations while for the RGM one has to also replace Γ 1 by the reduced Γ eff 1 given by Eq. (A.12).
where D/Dt ≡ ∂/∂t + u i ∂ i , ρ is the density, T is the temperature, s is the specific entropy, F rad k is the k-component of the radiative flux, and τ ik is the viscous stress tensor. Equation (B.1) can be recast in terms of gas pressure and density by using the thermodynamic identity where P g is the gas pressure, (Γ 3 − 1) ≡ (∂ ln T/∂ ln ρ) s . Then, using Eq. (B.2) together with the conservation of mass, Eq. (B.1) finally becomes after averaging which is strictly equivalent to Eq. (11) of Rosenthal et al. (1999 which is the classical relation for adiabatic oscillations. Note that we assumed Γ 1 P g t Γ 1 t P g t , which is quite an accurate approximation as verified by our solar numerical 3D simulation. The assumptions made to derive Eq. (B.4) from Eq. (B.3), which consist in assuming that the Lagrangian perturbation of the last four terms of Eq. (B.3) are zero, need further discussion. Indeed, as recognized by Rosenthal et al. (1999), those approximations are quite radical. It is nevertheless useful to go a step further and to explain what are the underlying physical assumptions. To that end, let's recast Eq. (B.3) in the following form dP g dt = −Γ 1 P g ∂ i u i + (Γ 1 − 1) u i ∂ i P g + (Γ 1 − 1) ∂ i u i P g − Γ 1 ∂ i P g u i − (Γ 1 − 1) P g ∂ i u i − (Γ 3 − 1) ∂ k F rad k + (Γ 3 − 1) τ ik ∂ k u i , (B.5) where, for sake of simplicity and without loss of meaning, we assumed that thermodynamic quantities are time-independent. To recover Eq. (B.4) from Eq. (B.5), the Lagrangian perturbations of the last six terms must be neglected. More precisely; the perturbation of the buoyancy work (u i ∂ i P g ), which also appears as a source of turbulent kinetic energy and thus a sink of thermal energy, is considered to be null. the perturbation of the divergence of the terms P g u i and u i P g are set to zero. For the former, considering a perfect gas, it is proportional to the convective (enthalpy) flux because of the relation where R is the ideal gas constant. Concerning u i P g , it is also essentially proportional to the convective flux as shown by Canuto (1997) using a polytropic relation. Therefore, one can conclude that neglecting the perturbation of those two terms is equivalent to neglect the perturbation of the convective heat flux.
the perturbation of the pressure-strain rate, P g ∂ i u i , is neglected. This can be justified as it scales as the squared turbulent Mach number of the rate of dissipation of turbulent kinetic energy into heat (Sarkar 1992). Hence, it can be neglected because we are considering turbulent flow at relatively low turbulent Mach numbers. the perturbations of the radiative flux as well as the dissipation rate of turbulent energy into heat (i.e. τ ik ∂ k u i ) are also considered to be negligible.