© K. Belkacem et al. 2021

1. Introduction

Systematic differences between observed and modelled eigenfrequencies is a long-standing problem in stellar seismology. For the Sun and solar-like stars, this has been identified and found to be associated with an improper modelling 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 the surface effects must be understood.

To this end, Rosenthal et al. (1999) used 3D hydrodynamical simulations, which allowed them 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), Magic & Weiss (2016), Houdek et al. (2017), Jørgensen & Weiss (2019), Schou & Birch (2020) for the Sun and Sonoi et al. (2015, 2017), Ball et al. (2016), Trampedach et al. (2017), Jørgensen et al. (2017, 2018, 2019, 2021), Manchon et al. (2018), Jørgensen & Angelou (2019), Houdek et al. (2019), Mosumgaard et al. (2020) 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.

For this purpose, Rosenthal et al. (1999) introduced two empirical models. The reduced-Γ₁ model (RGM), which is based on the argument that the perturbation of turbulent pressure is negligible compared to the perturbation of gas pressure, and the gas-Γ₁ model (GGM), which is based on a consideration that the perturbations of gas pressure and turbulent pressure are in phase. The authors favoured the GGM as it permits them to better reproduce the observed solar frequencies. Most of the above-mentioned works then applied the GGM, except in a few cases (e.g. Jørgensen & Weiss 2019), where the RGM was guided by the non-adiabatic calculation from 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 it demands extensive theoretical investigation to gain insights into the physics of the problem.

In this work, we aim at deriving and then assessing the theoretical validity of the GGM and RGM empirical models.

2. 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 described in Appendix A. To derive them, a first step is to express the equation for turbulent pressure.

2.1. Averaged equation for turbulent pressure

We start with the equation for the averaged Reynolds stresses (e.g. Canuto 1997; Gatski & Bonnet 2013):

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{r}_{\mathit{ij}}}{\mathrm{d}t}& =-r_{\mathit{ij}}{\partial }_k{\stackrel{\sim }{u}}_k-r_{\mathit{jk}}{\partial }_k{\stackrel{\sim }{u}}_i-r_{\mathit{ik}}{\partial }_k{\stackrel{\sim }{u}}_j-{\partial }_k\overline{\rho u_i^{″}{u}_j^{″}{u}_k^{″}}-\overline{u_i^{″}{\partial }_j{P}_{\mathrm{g}}}\hfill \\ \hfill & -\overline{u_j^{″}{\partial }_i{P}_{\mathrm{g}}}+\overline{u_j{\partial }_k{\tau }_{\mathit{ik}}}+\overline{u_i{\partial }_k{\tau }_{\mathit{jk}}},\hfill \end{array}$$(1)

where the overbar denotes the Reynolds average and the tilde denotes the Favre average (see Appendix A for a definition), u_i is 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_{\mathit{ij}}\equiv \overline{\rho u_i″u_j″}$ are the Reynolds stresses, where $u_i^{″}$ is the i-th component of the velocity fluctuation around its Favre average (${\stackrel{\sim }{u}}_i$). We also employ the notations, ∂_t ≡ ∂/∂t and ∂_i ≡ ∂/∂x_i, Einstein’s notation for repeated indices, and the pseudo-Lagrangian derivative is defined by $\mathrm{d}/\mathrm{d}t\equiv \partial /\partial t+{\stackrel{\sim }{u}}_j\partial /\partial 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. Thus, we obtain:

$$\begin{array}{cc}\hfill \frac{1}{2}\frac{\mathrm{d}{P}_{\mathrm{t}}}{\mathrm{d}t}& =-\frac{3P_{\mathrm{t}}}{2}{\partial }_{\mathrm{r}}{\stackrel{\sim }{u}}_{\mathrm{r}}-\overline{u_{\mathrm{r}}^{″}}{\partial }_{\mathrm{r}}\overline{P_{\mathrm{g}}}-{\partial }_{\mathrm{r}}(F_{\mathrm{rrr}}^{\mathrm{visc}}+F_{\mathrm{rrr}}^{\mathrm{turb}}+F_{\mathrm{r}}^p)\hfill \\ \hfill & +\overline{P_{\mathrm{g}}^{\prime }{\partial }_{\mathrm{r}}{u}_{\mathrm{r}}^{″}}-\overline{{\tau }_{\mathrm{r}k}{\partial }_k{u}_{\mathrm{r}}},\hfill \end{array}$$(2)

where $P_{\mathrm{t}}\equiv \overline{\rho u_{\mathrm{r}}″u_{\mathrm{r}}″}$ is the turbulent pressure, $F_{\mathrm{rrr}}^{\mathrm{turb}}\equiv \overline{\rho u_{\mathrm{r}}″u_{\mathrm{r}}″u_{\mathrm{r}}″}/2$, $F_{\mathrm{rrr}}^{\mathrm{visc}}\equiv -\overline{u_{\mathrm{r}}{\tau }_{\mathrm{rr}}}$, and $F_{\mathrm{r}}^p=\overline{u_{\mathrm{r}}^{″}{P}_{\mathrm{g}}\prime }$. 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 the 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_{\mathrm{rrr}}^{\mathrm{visc}}$) 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, that is, $\overline{{\tau }_{\mathrm{r}k}{\partial }_k{u}_{\mathrm{r}}}\propto \overline{{\tau }_{\mathit{ij}}{\partial }_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 ($\overline{P_{\mathrm{g}}\prime {\partial }_{\mathrm{r}}{u}_{\mathrm{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) is reduced to:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{P}_{\mathrm{t}}}{\mathrm{d}t}& =-3P_{\mathrm{t}}{\partial }_{\mathrm{r}}{\stackrel{\sim }{u}}_{\mathrm{r}}-2\overline{u_{\mathrm{r}}^{″}}{\partial }_{\mathrm{r}}\overline{P_{\mathrm{g}}}-{\partial }_{\mathrm{r}}{F}_{\mathrm{r}}^p-{\partial }_{\mathrm{r}}{F}_{\mathrm{rrr}}^{\mathrm{turb}}-\frac{2}{3}\overline{\rho }ϵ,\hfill \end{array}$$(3)

where $\overline{\rho }ϵ\equiv \overline{{\tau }_{\mathit{ij}}{\partial }_i{u}_j}$ is the dissipation rate of turbulent kinetic energy into heat.

2.2. Modelling the transport of turbulent pressure

The flux of turbulent pressure, $F_{\mathrm{rrr}}^{\mathrm{turb}}$ 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:

$$\begin{array}{c}\hfill F_{ii\mathrm{r}}^{\mathrm{turb}}=c^{-1}\frac{k}{ϵ}{F}_{\mathrm{r}}^{ϵ},\end{array}$$(4)

where $F_{ii\mathrm{r}}^{\mathrm{turb}}$ is the flux of turbulent kinetic energy, $F_{\mathrm{r}}^{ϵ}\equiv \overline{\rho ϵu_{\mathrm{r}}^{″}}$ is the flux of turbulent dissipation, $k\equiv r_{\mathit{ii}}/(2\overline{\rho })$ is the specific turbulent kinetic energy, and c is a parameter to be specified.

From Eq. (4), we then obtain:

$$\begin{array}{c}\hfill {\partial }_{\mathrm{r}}{F}_{\mathrm{rrr}}^{\mathrm{turb}}=\frac{\alpha }{c}{\omega }^{-1}{\partial }_{\mathrm{r}}{F}_{\mathrm{r}}^{ϵ}-\alpha F_{ii\mathrm{r}}^{\mathrm{turb}}{\partial }_{\mathrm{r}}ln\omega ,\end{array}$$(5)

where ω ≡ ϵ/k is the turbulence frequency and

$$\begin{array}{c}\hfill \alpha \equiv \frac{{\partial }_{\mathrm{r}}{F}_{\mathrm{rrr}}^{\mathrm{turb}}}{{\partial }_{\mathrm{r}}{F}_{ii\mathrm{r}}^{\mathrm{turb}}}.\end{array}$$(6)

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), we can rewrite Eq. (3) as:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{P}_{\mathrm{t}}}{\mathrm{d}t}& =-3P_{\mathrm{t}}{\partial }_{\mathrm{r}}{\stackrel{\sim }{u}}_{\mathrm{r}}-2\overline{u_{\mathrm{r}}^{″}}{\partial }_{\mathrm{r}}\overline{P_{\mathrm{g}}}-\frac{2}{3}\overline{\rho }ϵ-\frac{\alpha }{c}\frac{k}{ϵ}{\partial }_{\mathrm{r}}{F}_{\mathrm{r}}^{ϵ}\hfill \\ \hfill & -{\partial }_{\mathrm{r}}{F}_{\mathrm{r}}^p+\alpha F_{ii\mathrm{r}}^{\mathrm{turb}}{\partial }_{\mathrm{r}}ln\omega ,\hfill \end{array}$$(7)

where one still needs a prescription for the flux of turbulent dissipation. To do so, we use the equation governing ϵ. A standard procedure consists in constructing a parametric ϵ-equation as is done in two-equation models (e.g. k − ϵ models). Following this approach, we can write (e.g. Pope 2000; Wilcox 2006; Gatski & Bonnet 2013):

$$\begin{array}{cc}\hfill \overline{\rho }\frac{\mathrm{d}ϵ}{\mathrm{d}t}+{\partial }_{\mathrm{r}}{F}_{\mathrm{r}}^{ϵ}& =-c_1^{ϵ}\omega P_{\mathrm{t}}{\partial }_{\mathit{rmr}}{\stackrel{\sim }{u}}_{\mathrm{r}}-c_2^{ϵ}\overline{\rho }ϵ{\partial }_{\mathrm{r}}{\stackrel{\sim }{u}}_{\mathrm{r}}-c_3^{ϵ}\omega \overline{u_k^{″}}{\partial }_k\overline{P_{\mathrm{g}}}\hfill \\ \hfill & -c_4^{ϵ}\overline{\rho }\omega ϵ.\hfill \end{array}$$(8)

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^{ϵ}$ are discussed in the following section.

3. Recovering the gas-Γ₁ and reduced-Γ₁ models

Here, we derive an expression for the Lagrangian perturbation of turbulent pressure and make a number of assumptions and discuss their validity.

3.1. Perturbation of turbulent pressure

Using Eq. (7) together with Eq. (8), we obtain the desired expression for the equation governing turbulent pressure:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{P}_{\mathrm{t}}}{\mathrm{d}t}& =(3-\frac{\alpha c_2^{ϵ}}{2c}\mathrm{\Phi }-\frac{\alpha c_1^{ϵ}}{c})\frac{P_{\mathrm{t}}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}t}+(\frac{\alpha c_3^{ϵ}}{c}-2)\overline{u_{\mathrm{r}}^{″}}{\partial }_{\mathrm{r}}\overline{P_{\mathrm{g}}}\hfill \\ \hfill & +(\frac{\alpha c_4^{ϵ}}{c}-\frac{2}{3})\overline{\rho }ϵ+\frac{\alpha }{c}\frac{\overline{\rho }}{\omega }\frac{\mathrm{d}ϵ}{\mathrm{d}t}-{\partial }_{\mathrm{r}}{F}_{\mathrm{r}}^p+\alpha F_{ii\mathrm{r}}^{\mathrm{turb}}{\partial }_{\mathrm{r}}ln\omega ,\hfill \end{array}$$(9)

where we used the averaged continuity equation (Eq. (A.5)) and Φ ≡ r_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, that is, $\delta F_{\text{r}}^p=0$. (H2) The turbulence frequency (ω) is supposed to vary on a length scale that is 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, the last term of Eq. (9) vanishes. (H3) The Lagrangian perturbation of α is neglected, that is, δα = 0. (H4) The Lagrangian perturbation of the dissipation of turbulent kinetic energy into heat (δϵ) is neglected. (H5) The perturbation of the buoyancy work, $\delta (\overline{u_{\mathrm{r}}^{″}}{\partial }_{\mathrm{r}}\overline{P_{\mathrm{g}}})$, is neglected. (H6) The Lagrangian perturbation of density is assumed to be real.

Perturbing Eq. (9) and applying H1 to H6, we get:

$$\begin{array}{cc}& \frac{\delta P_{\mathrm{t}}}{P_{\mathrm{t}}}=(3-\frac{\alpha c_2^{ϵ}}{2c}\mathrm{\Phi }-\frac{\alpha c_1^{ϵ}}{c})\frac{\delta \rho }{\rho },\hfill \end{array}$$(10)

where δP_t and δρ are the pseudo-Lagrangian perturbations of turbulent pressure and density, respectively. As an introduction to the next section, we can already mention that Eq. (10) allows us 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 points 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, we find that 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 on its own is not sufficient to derive Eqs. (10) and (A.9).

The first assumption, H1, is essentially equivalent to neglecting the perturbation of the convective flux ($\delta F_{\mathrm{r}}^{\mathrm{conv}}$) because $F_{\text{r}}^p$ 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 we need to determine ϵ and the simplest possible way to do this 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), we may readily find 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 of 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 Γ₃ − 1 = 2/3, where Γ₃ ≡ (∂lnT/∂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 ought to be assessed definitively. 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 questions regarding the validity of both the RGM and GGM models because, even in the adiabatic limit, they introduce oversimplifying hypotheses regarding the properties of turbulent convection and mode physics.

3.2. Application to solar p-mode frequencies

Now we go on to investigate how Eq. (10) permits us to recover the GGM and RGM models. One prerequisite is the specification of the coefficients, $c_1^{ϵ}$ and $c_2^{ϵ}$. A commonly accepted value for the first coefficient is $c_1^{ϵ}\simeq 1.44$ (e.g. Pope 2000; Wilcox 2006; Gatski & Bonnet 2013). For $c_2^{ϵ}$, we adopt the model of Coleman & Mansour (1991) which gives:

$$\begin{array}{c}\hfill c_2^{ϵ}=\frac{1}{3}[1+3n({\mathrm{\Gamma }}_1-1)-2c_1^{ϵ}],\end{array}$$(11)

where n ≃ 0.75 is the exponent of the viscosity law on temperature. Indeed, Eq. (11) is aimed 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) was 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.

In this way, we 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 Eqs. (A.7)–(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 ± 18 K and the chemical composition is (X, Y, Z) = (0.7373, 0.2427, 0.0200) with Grevesse et al. (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 the standard solar frequency ${(G{M}_{\odot }/R_{\odot }^3)}^{1/2}/(2\pi )=99.85537\mu$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). First, we considered what guidance for c we might 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 used for such work because it does not introduce, directly or indirectly, subgrid scale viscosities: dissipation takes place by way of explicit, physical viscosity and a very small contribution from time integration only. The data were collected for a total of 12 different model configurations (five of which have been discussed in detail in Muthsam et al. 1995, 1999). They covered a range in Prandtl number Pr from 0.1 to 1 at a ‘zone Rayleigh number’ of 10⁵ to 10⁶. For a horizontal domain width eight times the size of this zone depth and with three to four granules found along each horizontal direction, this implies a ‘granulation diameter-based Rayleigh number’ that is ∼16 to 50 times larger or an Ra in the range of 10⁶ to 5 × 10⁷. This yields a product of Ra and Pr in the range of 10⁵ to 5 × 10⁷. 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 fact that the convection zones are much more shallow than in the solar case, where additionally Pr ≪ 10⁻⁶. 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 a 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 for computing $F_{\text{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 as 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.

Fig. 1. Frequency difference for radial solar modes between the observed frequencies, as inferred by Broomhall et al. (2009) and Davies et al. (2014), and the modelled frequencies computed using the GGM model, the RGM model, and the model of the present work (see text for details).

4. Conclusions

By using the Reynolds and Favre-averaged, fully compressible Navier-Stokes equations, we demonstrate 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 show the RGM and GGM models are compatible with the adiabatic approximation but this also implies more drastic and unrealistic physical assumptions regarding turbulence and mode physics.

A 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 direct 3D numerical simulations prove useful. Unfortunately, this is still out of reach for the Sun due to the limits on our current numerical capacities. Hence, extrapolations from more accessible parameter ranges remain necessary.

It is thus difficult to draw a conclusion on which of the two models is more appropriate in application. What complicates the matter further, given the above-mentioned assumptions which are needed to recover these empirical models, we can safely conclude none of them are firmly physically grounded, even in the adiabatic limit. However, we may 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 (in order 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, including issues related to free parameters and questionable physical assumptions (see Houdek & Dupret 2015, for details), which prevent them assuring that we properly grasp the physics of the problem.

Acknowledgments

F. Kupka is grateful to the Austrian Science Fund FWF for support through projects P29172-N and P33140-N and support from European Research Council (ERC) Synergy Grant WHOLESUN #810218.

References

Appendix A: The gas-Γ₁ (GGM) and reduced-Γ₁ (RGM) empirical models

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 this scenario, let us begin by considering the mass and momentum conservation equations, which read:

$$\begin{array}{cc}& {\partial }_{\mathrm{t}}\rho +{\partial }_i\left(\rho u_i\right)=0,\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& {\partial }_{\mathrm{t}}\left(\rho u_i\right)+{\partial }_k\left(\rho u_i{u}_k\right)=-{\partial }_i{P}_{\mathrm{g}}-\rho g_i+{\partial }_k{\tau }_{\mathit{jk}},\hfill \end{array}$$(A.2)

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:

$$\begin{array}{c}\hfill \overline{X}\equiv X-X^{\prime },\mathrm{with}\overline{X^{\prime }}=0,\end{array}$$(A.3)

and the Favre average is defined, for a quantity Y, by

$$\begin{array}{c}\hfill \stackrel{\sim }{Y}=\frac{\overline{\rho Y}}{\overline{Y}}\mathrm{so}\mathrm{that}Y=\stackrel{\sim }{Y}+Y^{″}\mathrm{and}\overline{\rho Y^{″}}=0.\end{array}$$(A.4)

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

$$\begin{array}{cc}& \frac{\mathrm{d}\overline{\rho }}{\mathrm{d}t}+\overline{\rho }{\partial }_i{\stackrel{\sim }{u}}_i=0,\hfill \end{array}$$(A.5)

$$\begin{array}{cc}& \overline{\rho }\frac{\mathrm{d}{\stackrel{\sim }{u}}_i}{\mathrm{d}t}+{\partial }_k{r}_{\mathit{ik}}=-{\partial }_i{\overline{P}}_{\mathrm{g}}-\overline{\rho }\overline{g_i},\hfill \end{array}$$(A.6)

where $r_{\mathit{ik}}\equiv \overline{\rho 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 $\mathrm{d}/\mathrm{d}t\equiv \partial /\partial t+{\stackrel{\sim }{u}}_j\partial /\partial 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 (${\stackrel{\sim }{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_{\mathrm{t}}\equiv \overline{\rho u_{\mathrm{r}}^{″}{u}_{\mathrm{r}}^{″}}$. In addition, we split the mean quantities such that $\overline{X}=X_0+\delta X$ (where $X_0\equiv <\overline{X}{>}_{\mathrm{t}}$ is the time-average of $\overline{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:

$$\begin{array}{cc}& \frac{\delta \rho }{{\rho }_0}+\frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2{\xi }_{\mathrm{r}}\right)=0,\hfill \end{array}$$(A.7)

$$\begin{array}{cc}& {\rho }_0{\sigma }^2{\xi }_{\mathrm{r}}=\frac{\mathrm{d}\delta P_{\mathrm{tot}}}{\mathrm{d}r}+\delta \rho g_0,\hfill \end{array}$$(A.8)

where σ is the angular frequency, ξ_r is the radial component of the eigen-displacement (${\stackrel{\sim }{u}}_{\mathrm{r}}=i\omega {\xi }_{\mathrm{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:

$$\begin{array}{c}\hfill \frac{\delta P_{\mathrm{g}}}{P_{\mathrm{g},0}}={\mathrm{\Gamma }}_1\frac{\delta \rho }{{\rho }_0},\end{array}$$(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. We note 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), we need 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:

$$\begin{array}{c}\hfill \frac{\delta P_{\mathrm{t}}}{P_{\mathrm{t}}}=\mathcal{A}\frac{\delta \rho }{\rho },\end{array}$$(A.10)

where 𝒜 is to be determined. Therefore, using Eq. (A.10), we formally write:

$$\begin{array}{cc}\hfill \frac{\delta P_{\mathrm{tot}}}{P_{\mathrm{tot}}}& ={\mathrm{\Gamma }}_1^{\mathrm{eff}}\frac{\delta \rho }{\rho },\mathrm{with}{\mathrm{\Gamma }}_1^{\mathrm{eff}}\equiv [{\mathrm{\Gamma }}_1\frac{P_{\mathrm{g}}}{P_{\mathrm{tot}}}+\mathcal{A}\frac{P_{\mathrm{t}}}{P_{\mathrm{tot}}}].\hfill \end{array}$$(A.11)

Two empirical models have previously been introduced by Rosenthal et al. (1999). The first is the RGM, which consists of the argument that δP_t = 0 or 𝒜 = 0, so that Eq. (A.11) becomes:

$$\begin{array}{c}\hfill \frac{\delta P_{\mathrm{tot}}}{P_{\mathrm{tot}}}=\frac{\delta P_{\mathrm{g}}}{P_{\mathrm{tot}}}={\mathrm{\Gamma }}_1^{\mathrm{eff}}\frac{\delta \rho }{\rho },\mathrm{where}{\mathrm{\Gamma }}_1^{\mathrm{eff}}=\frac{{\mathrm{\Gamma }}_1{P}_{\mathrm{g}}}{P_{\mathrm{tot}}}\end{array}$$(A.12)

is called the reduced Γ₁. 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, considered that the real part of δP_t can be neglected. The second model is the GGM and has been introduced by Rosenthal et al. (1999). It consists of arguing the opposite, namely, that the perturbation of gas pressure and turbulent pressure are in phase. Hence, using 𝒜 = Γ₁ together with Eqs. (A.9), (A.10), and (A.11), we obtain:

$$\begin{array}{c}\hfill \frac{\delta P_{\mathrm{tot}}}{P_{\mathrm{tot}}}=\frac{\delta P_{\mathrm{t}}}{P_{\mathrm{t}}}=\frac{\delta P_{\mathrm{g}}}{P_{\mathrm{g}}}={\mathrm{\Gamma }}_1\frac{\delta \rho }{\rho }.\end{array}$$(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, it is also necessary to replace Γ₁ by the reduced ${\mathrm{\Gamma }}_1^{\text{eff}}$ given by Eq. (A.12).

Appendix B: Relation between perturbation of gas pressure and density

We follow the derivation proposed by Rosenthal et al. (1999), which starts with the equation governing specific entropy. It reads:

$$\begin{array}{c}\hfill \rho T\frac{\mathrm{D}s}{\mathrm{D}t}=-{\partial }_k{F}_k^{\mathrm{rad}}+{\tau }_{\mathit{ik}}{\partial }_k{u}_i,\end{array}$$(B.1)

where D/Dt ≡ ∂/∂t + u_i∂_i, ρ is the density, T is the temperature, s is the specific entropy, $F_k^{\text{rad}}$ 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:

$$\begin{array}{c}\hfill \rho T\frac{\mathrm{D}s}{\mathrm{D}t}=\frac{1}{({\mathrm{\Gamma }}_3-1)}[\frac{\mathrm{D}{P}_{\mathrm{g}}}{\mathrm{D}t}-\frac{{\mathrm{\Gamma }}_1{P}_{\mathrm{g}}}{\rho }\frac{\mathrm{D}\rho }{\mathrm{D}t}],\end{array}$$(B.2)

where P_g is the gas pressure, (Γ₃ − 1)≡(∂lnT/∂lnρ)_s. Then, using Eq. (B.2) together with the conservation of mass, after averaging, Eq. (B.1) finally becomes:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\overline{P}}_{\mathrm{g}}}{\mathrm{d}t}& =-\overline{{\mathrm{\Gamma }}_1{P}_{\mathrm{g}}}{\partial }_i{\stackrel{\sim }{u}}_i-\overline{{\mathrm{\Gamma }}_1{P}_{\mathrm{g}}{\partial }_i{u}_i^{″}}-\overline{u_i^{″}{\partial }_i{P}_{\mathrm{g}}}-\overline{({\mathrm{\Gamma }}_3-1){\partial }_k{F}_k^{\mathrm{rad}}}\hfill \\ \hfill & +\overline{({\mathrm{\Gamma }}_3-1){\tau }_{\mathit{ik}}{\partial }_k{u}_i},\hfill \end{array}$$(B.3)

which is strictly equivalent to Eq. (11) of Rosenthal et al. (1999). Following the same authors, it is assumed that time-varying parts (i.e. the Lagrangian perturbation) of the last four terms of Eq. (B.3) are zero. Consequently, the perturbation of Eq. (B.3) permits us to recover Eq. (A.9), that is

$$\begin{array}{c}\hfill \frac{\delta P_{\mathrm{g}}}{P_{\mathrm{g}}}={\mathrm{\Gamma }}_1\frac{\delta \rho }{\rho },\end{array}$$(B.4)

which is the classical relation for adiabatic oscillations. We note that we assumed ${\langle \overline{{\mathrm{\Gamma }}_1{P}_g}\rangle }_t\simeq {\langle \overline{{\mathrm{\Gamma }}_1}\rangle }_t{\langle \overline{P_{\mathrm{g}}}\rangle }_{\mathrm{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 consists of assuming that the Lagrangian perturbation of the last four terms of Eq. (B.3) are zero, require 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 the underlying physical assumptions. To that end, let us recast Eq. (B.3) in the following form:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\overline{P}}_{\mathrm{g}}}{\mathrm{d}t}& =-{\mathrm{\Gamma }}_1\overline{P_{\mathrm{g}}}{\partial }_i{\stackrel{\sim }{u}}_i+({\mathrm{\Gamma }}_1-1)\overline{u_i^{″}}{\partial }_i\overline{P_{\mathrm{g}}}+({\mathrm{\Gamma }}_1-1){\partial }_i\left(\overline{u_i^{″}{P}_{\mathrm{g}}^{\prime }}\right)\hfill \\ \hfill & -{\mathrm{\Gamma }}_1{\partial }_i\left(\overline{P_{\mathrm{g}}{u}_i^{″}}\right)-({\mathrm{\Gamma }}_1-1)\overline{P_{\mathrm{g}}^{\prime }{\partial }_i{u}_i^{″}}-({\mathrm{\Gamma }}_3-1){\partial }_k\overline{F_k^{\mathrm{rad}}}\hfill \\ \hfill & +({\mathrm{\Gamma }}_3-1)\overline{{\tau }_{\mathit{ik}}{\partial }_k{u}_i},\hfill \end{array}$$(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 ($\overline{u_i^{″}}{\partial }_i\overline{P_{\mathrm{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 $\overline{P_{\mathrm{g}}{u}_i^{″}}$ and $\overline{u_i^{″}{P}_{\mathrm{g}}\prime }$ are set to zero. For the former, considering a perfect gas, it is proportional to the convective (enthalpy) flux because of the relation:

  $$\begin{array}{c}\hfill \overline{P_{\mathrm{g}}{u}_i^{″}}=\mathcal{R}\overline{\rho }\stackrel{\sim }{T^{″}{u}_i^{″}},\end{array}$$(B.6)

  where ℛ is the ideal gas constant. Concerning $\overline{u_i^{″}{P}_{\mathrm{g}}\prime }$, it is also essentially proportional to the convective flux as shown by Canuto (1997) using a polytropic relation. Therefore, we can conclude that neglecting the perturbation of those two terms is equivalent to neglecting the perturbation of the convective heat flux.

• The perturbation of the pressure-strain rate, $\overline{P_{\mathrm{g}}\prime {\partial }_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. $\overline{{\tau }_{\mathit{ik}}{\partial }_k{u}_i}$) are also considered to be negligible.

Appendix C: Averaged equation for specific entropy

Furthermore, we now look towards inferring the averaged entropy equation. To do so, rather than starting directly with the entropy equation, we consider the equation governing the enthalpy. It reads:

$$\begin{array}{c}\hfill \rho \frac{\mathrm{D}h}{\mathrm{D}t}=\frac{\mathrm{D}{P}_{\mathrm{g}}}{\mathrm{D}t}-{\partial }_k{F}_k^{\mathrm{rad}}+{\tau }_{\mathit{ik}}{\partial }_k{u}_i,\end{array}$$(C.1)

where D/Dt ≡ ∂/∂t + u_i∂_i and $F_k^{\text{rad}}$ is the k-component of the radiative flux. Now, averaging Eq. (C.1) leads to:

$$\begin{array}{cc}\hfill \overline{\rho }\frac{\mathrm{d}\stackrel{\sim }{h}}{\mathrm{d}t}& =\frac{\mathrm{d}{\overline{P}}_{\mathrm{g}}}{\mathrm{d}t}+\overline{u_k^{″}}{\partial }_k\overline{P_{\mathrm{g}}}-{\partial }_k({\overline{F}}_k^{\mathrm{rad}}+F_k^{\mathrm{conv}}-\overline{u_k^{″}{P}_{\mathrm{g}}^{\prime }})\hfill \\ \hfill & -\overline{P_{\mathrm{g}}^{\prime }{\partial }_i{u}_i^{″}}+\overline{{\tau }_{\mathit{ik}}{\partial }_k{u}_i},\hfill \end{array}$$(C.2)

where $F_k^{\mathrm{conv}}\equiv \overline{\rho h^{″}{u}_k^{″}}$ is the k-component of the convective flux. Then using the thermodynamic relation, $\overline{\rho }\stackrel{\sim }{T}\mathrm{d}\stackrel{\sim }{s}=\overline{\rho }\mathrm{d}\stackrel{\sim }{h}-\mathrm{d}{\overline{P}}_{\mathrm{g}}$, valid at leading order, and identifying the Reynolds average to the horizontal average gives:

$$\begin{array}{cc}\hfill \overline{\rho }\stackrel{\sim }{T}\frac{\mathrm{d}\stackrel{\sim }{s}}{\mathrm{d}t}& =\overline{u_{\mathrm{r}}^{″}}{\partial }_{\mathrm{r}}\overline{P_{\mathrm{g}}}-{\partial }_r({\overline{F}}_{\mathrm{r}}^{\mathrm{rad}}+F_{\mathrm{r}}^{\mathrm{conv}}-\overline{u_{\mathrm{r}}^{″}{P}_{\mathrm{g}}^{\prime }})-\overline{P_{\mathrm{g}}^{\prime }{\partial }_i{u}_i^{″}}\hfill \\ \hfill & +\overline{{\tau }_{\mathit{ik}}{\partial }_k{u}_i}.\hfill \end{array}$$(C.3)

Neglecting the pressure-strain rate as for the equation governing turbulent pressure and adopting the same notations, we finally obtain:

$$\begin{array}{cc}\hfill \overline{\rho }\stackrel{\sim }{T}\frac{\mathrm{d}\stackrel{\sim }{s}}{\mathrm{d}t}& =\overline{u_{\mathrm{r}}^{″}}{\partial }_{\mathrm{r}}\overline{P_{\mathrm{g}}}-{\partial }_{\mathrm{r}}({\overline{F}}_{\mathrm{r}}^{\mathrm{rad}}+F_{\mathrm{r}}^{\mathrm{conv}}-F_{\mathrm{r}}^p)+\overline{\rho }ϵ.\hfill \end{array}$$(C.4)

Perturbing Eq. (C.4) thus permits us to derive an expression for the perturbation of the buoyancy work. It reads:

$$\begin{array}{c}\hfill \delta \left(\overline{u_{\mathrm{r}}^{″}}{\partial }_{\mathrm{r}}\overline{P_{\mathrm{g}}}\right)=i\sigma \overline{\rho }\stackrel{\sim }{T}\delta s+\delta \left[{\partial }_{\mathrm{r}}({\overline{F}}_{\mathrm{r}}^{\mathrm{rad}}+F_{\mathrm{r}}^{\mathrm{conv}}-F_{\mathrm{r}}^p)\right]-\delta \left(\overline{\rho }ϵ\right).\end{array}$$(C.5)

Using the same approximations as described in Sect. 3.1 to derive Eq. (10), immediately leads to δs = 0.

All Figures

	Fig. 1. Frequency difference for radial solar modes between the observed frequencies, as inferred by Broomhall et al. (2009) and Davies et al. (2014), and the modelled frequencies computed using the GGM model, the RGM model, and the model of the present work (see text for details).
In the text