Solar mode damping rates: insight from a 3D hydrodynamical simulation
Key Words.:
Waves  Convection  Sun: oscillationsSpaceborne missions such as CoRoT and Kepler have provided a rich harvest of highquality photometric data for solarlike pulsators. It is now possible to measure damping rates for hundreds of mainsequence and thousands of redgiant stars with an unprecedented precision. However, among the seismic parameters, mode damping rates remain poorly understood and thus barely used for inferring the physical properties of stars. Previous approaches to model mode damping rates were based on mixinglength theory or a Reynoldsstress approach to model turbulent convection. While able to grasp the main physics of the problem, those approaches are of little help to provide quantitative estimates as well as a definitive answer on the relative contribution of each physical mechanism. Indeed, due to the high complexity of the turbulent flow and its interplay with the oscillations, those theories rely on many free parameters which inhibits an indepth understanding of the problem. Our aim is thus to assess the ability of 3D hydrodynamical simulations to infer the physical mechanisms responsible for damping of solarlike oscillations. To this end, a solar highspatial resolution and longduration hydrodynamical 3D simulation computed with the ANTARES code allows probing the coupling between turbulent convection and the normal modes of the simulated box. Indeed, normal modes of the simulation experience realistic driving and damping in the superadiabatic layers of the simulation. Therefore, investigating the properties of the normal modes in the simulation provides a unique insight into the mode physics. We demonstrate that such an approach provides constraints on the solar damping rates and is able to disentangle the relative contribution related to the perturbation (by the oscillation) of the turbulent pressure, the gas pressure, the radiative flux, and the convective flux contributions. Finally, we conclude that using the normal modes of a 3D numerical simulation is possible and is potentially able to unveil the respective role of the different physical mechanisms responsible for mode damping provided the timeduration of the simulation is long enough.
1 Introduction
The Kepler (Borucki et al. 2010) and CoRoT (Baglin et al. 2006a, b) spaceborne missions provided a wealth of highquality and longduration photometric data which allowed us to detect thousands of solarlike oscillating stars from the mainsequence up to the red clump (e.g. Chaplin & Miglio 2013). A leap forward has then been achieved concerning our understanding of the internal structure of stars and their evolution (e.g. Mosser et al. 2012, 2014; Hekker & ChristensenDalsgaard 2017). This was made possible thanks to an indepth understanding of the physics of the oscillations (e.g. Belkacem et al. 2006a, b, 2010; Dupret et al. 2009; Belkacem et al. 2012; Grosjean et al. 2014; Samadi et al. 2015; Houdek & Dupret 2015) as well as of the evolution of the properties of the power spectra along with the stellar evolution (e.g. Belkacem et al. 2011; Belkacem & Samadi 2013). However, since the discovery of solar fiveminute oscillations, our ability to understand the physical mechanisms underlying mode damping is still challenged. Indeed, mode damping occurs in the uppermost layers of solarlike stars in which convection is highly turbulent. It also corresponds to the location of the transition between convective and radiative energy transport. This makes the problem highly intricate because there is a concordance of several characteristic timescales. The modal period is found to be of the same order as both the convective and thermal timescales. Consequently, while it is notoriously difficult to model turbulent convective boundary layers (see the review by Kupka & Muthsam 2017), it is even more complex to address its coupling with the oscillations.
Several attempts to model mode damping have nevertheless been proposed. The first to investigate this issue were Ando & Osaki (1975) but no stable modes were found in the whole frequency range. Goldreich & Kumar (1991) subsequently proposed that the shear due to Reynolds stresses, modelled by an eddyviscosity, is of the same order of magnitude as the nonadiabatic component of the perturbation of gas pressure. Gough (1980) and Balmforth (1992) (see also Houdek et al. 1999) found that the damping is dominated by the modulation of turbulent pressure, while Grigahcène et al. (2005), Dupret et al. (2006), and Belkacem et al. (2012) also included the perturbation of the dissipation rate of kinetic energy into heat that acts to compensate the perturbation of turbulent pressure. Those formalisms were based on the mixinglength theory of convection (see Houdek & Dupret 2015, for a detailed discussion). It thus reduces the whole of the turbulent cascade to a single lengthscale. Therefore, the perturbation of the mixinglength cannot account for the relation between oscillations and the turbulent cascade. Xiong et al. (2000) proposed an alternative approach using a Reynolds stress formalism (Canuto 1992) to model convection and, using a perturbation method, computed mode damping rates. However, they found a number of unstable modes, which is at odds with the observational evidence. Indeed, such an approach, while being more physically grounded than the use of the mixinglength, is highly sensitive to the adopted closure models (see Kupka & Robinson 2007; Kupka & Muthsam 2007a, b, for details).
Additional constraints were thus welcomed to gain more insight into the problem of mode damping. Based on the CoRoT and Kepler observations, which have provided accurate observations of solarlike oscillations of hundreds of mainsequence stars and thousands of redgiant stars, it has been shown that damping rates follow some scaling relations. Indeed, even if it is possible to reproduce the solar damping rates by tuning some parameters, it does not ensure to reproduce the damping in other stars. The observed scaling relation of mode linewidths thus provides an important additional constraint on the modeling. First, using groundbased observations, Chaplin et al. (2009) proposed that mode linewidths follow a powerlaw dependence on effective temperature. This has been refined by Baudin et al. (2011a, b) using a homogeneous sample of stars observed by CoRoT and later by Appourchaux et al. (2012); Appourchaux (2014) for mainsequence and sub giants as well as Vrard et al. (2018) for red giants using Kepler observations. From a theoretical point of view, Chaplin et al. (2009) predicted a powerlaw of , which disagrees with CoRoT and Kepler observations. Houdek (2012) attributed the failure of this theory to a missing physical mechanism and proposed mode scattering as a possible solution. However, based on the formalism developed by Grigahcène et al. (2005), Belkacem et al. (2012) managed to reproduce the scaling relations obtained by both CoRoT and Kepler without invoking mode scattering.
Finally, despite some relative successes, much work is still needed for a proper understanding of the physics behind mode linewidths. A novel approach is needed. The use of hydrodynamical 3D numerical simulation is potentially able to offer us such an opportunity. A possible way to proceed consists in constraining the free parameters of the modelling using a 3D solar simulation as proposed by Houdek et al. (2017); Houdek (2017); Aarslev et al. (2018). However, even though the observed damping rates are reproduced, some parameters have been tuned. Therefore, it is difficult to get insight into the physics of mode damping. A more promising approach thus consists in investigating directly the normal modes of the 3D numerical simulations. Indeed, turbulent convection generates acoustic noise and, thanks to the boundary conditions, normal modes exist in the simulation. Even if the spatial structure of the modes is not realistic compared to the observed modes, they experience realistic driving and damping in the superadiabatic layers of the simulation. Therefore, investigating the properties of the normal modes in the simulation provides a unique insight into the mode physics. Indeed, contrary to the observed modes, it is possible to access the perturbations associated with the oscillations in all the physical quantities and as a function of height in the simulation. For mode driving, such an approach has already been successfully used (e.g. Stein & Nordlund 2001; Samadi et al. 2003). In this article, we aim at assessing our ability to investigate mode damping using 3D numerical simulations.
This paper is organised as follows: in Sect. 2 we present the solar 3D numerical simulation computed with the ANTARES code. The properties of the normal modes of the simulation are then described in Sect. 3. In Sect. 4, the normal mode work integrals are computed and the contributions associated with the modulation of turbulent pressure, gas pressure, radiative and convective fluxes are investigated in Sect. 5. Finally, conclusions are provided in Sect. 6.
2 The solar simulation
For our subsequent analyses we use data from a numerical, hydrodynamical simulation of the solar surface. The simulation was computed as part of a more extended research project devoted to solar physics and the adequacy of numerical tools used in modelling the solar surface (see GrimmStrele et al. 2015a). One of the simulations computed for this purpose and called cosc13 was used for the present work. Its setup and some details of the simulation code are described in the following.
2.1 The ANTARES simulation code
The solar 3D hydrodynamical simulation cosc13 has been computed with the ANTARES code (Muthsam et al. 2010) which numerically solves the NavierStokes equations (NSE) for a fully compressible fluid and accounts for radiative heat transfer, viscous processes, and realistic microphysics. ANTARES uses a conservative, 5 order weighted essentially nonoscillatory (WENO5) finite difference scheme (Shu 2003; Jiang & Shu 1996) to discretize pressure gradients and terms representing advection in the NSE and a fourth order conservative finite difference scheme (Happenhofer et al. 2013) to discretize dissipative terms. Time integration is done by an explicit threestage, second order RungeKutta scheme known as SSP RK(3,2) originally proposed by Kraaijevanger (1991) and analysed by Kupka et al. (2012). GrimmStrele et al. (2015b) demonstrated that this scheme is more efficient than traditional, total variation diminishing schemes used for this purpose (the TVD2 and TVD3 methods of Shu & Osher (1988), originally proposed by Heun (1900) and Fehlberg (1970), respectively).
For cosc13 radiative transfer was treated in the nongrey approximation as described in Muthsam et al. (2010) with a multigroup technique that assigns frequencies to one of four bins, according to the optical depth at which radiation is emitted, in order to account for the full, lineblanketed spectrum in the radiative transfer. The angular integration of the radiative transfer equation required to compute the radiative flux was performed using a threepoint GaussRadau integration along polar coordinates per hemisphere and a fourpoint equispaced integration along the azimuthal direction. Using a shortcharacteristics method (cf. Muthsam et al. 2010) the scheme hence features 18 ray directions including two vertical rays into each hemisphere. The latter was not the case for the scheme used by Muthsam et al. (2010). As described in the latter the radiative source and sink term, , is treated in a stationary approximation, since changes in the simulation occur on timescales that are orders of magnitude longer than the light crossing time of a unity optical depth (see, e.g., Mihalas & Mihalas 1984, Chapts. 6.5 and 7.2). The transition between the surface layers, for which the radiative transfer equation is solved, and the lower lying layers, where the diffusion approximation holds, is obtained from smooth interpolation between the two solutions for grid cells always located in the optically thick region: the average temperature there is typically around 12500 K and thus the mean optical Rosseland depth is in a range from about 1000 to several 1000, so even for extreme events optical thickness is ensured (cf. also Fig. 15 of Stein & Nordlund 1998).
2.2 Setup of the simulation: microphysics and simulation grid
The sum of gas and radiative pressure, the thermodynamical derivatives, and related thermodynamical quantities are obtained from interpolating in the LLNL equation of state tables (Rogers et al. 1996). Realistic (radiative) conductivities are obtained from interpolation in the opacity data of Iglesias & Rogers (1996) for interior layers. Nongrey opacities are obtained from the tables of Kurucz (1993b, a) with the binning procedure described in Muthsam et al. (2010) where also further details on the construction of opacities and equation of state from these tables are given (see GrimmStrele et al. 2015a, too). The old standard solar composition of Grevesse & Noels (1993) was assumed.
The numerical simulation has been conducted for a box with a Cartesian grid containing a volume of 3.88 Mm as measured vertically and 6 Mm as measured horizontally. The spatial location of this “box within the Sun” was chosen to contain the solar surface such that layers with an optical Rosseland depth larger than always remained inside the simulation box. Open boundaries as in GrimmStrele et al. (2015a) (type BC3b, from their Table 1) have been used in vertical directions as well as periodical boundary conditions for both horizontal directions. A uniform resolution has been chosen with a grid spacing of 11.1 km from top to bottom for cells which are 35.3 km wide. The computations for this model have been done on the VSC2 (Vienna Scientific Cluster) on 144 CPU cores in parallel. MPI parallelization was used and due to extra grid cells required to implement the boundary conditions the total grid actually contained 359 grid points in the vertical and 179 points along the two horizontal directions. From this grid 350 vertical layers can be extracted which consist of 170 by 170 cells that are used for computing horizontal averages during postprocessing of the simulation.
2.3 Initial conditions, model relaxation, and statistical evaluation
The simulation was relaxed from its initial condition for about s or sound crossing times. The latter are measured from the top of the simulation box to its bottom for the initial model. In the present case this has been the helioseismologically calibrated standard solar model S (ChristensenDalsgaard et al. 1996). Since that model is actually a bit too shallow at the top we extended the simulation domain upwards by about 200 km. In GrimmStrele et al. (2015a) the procedures to do so have been explained. It has also been shown there that the specific 1D starting model is not crucial, as simulations started from different solar structure models which agree in effective temperature and input entropy in the quasiadiabatic layers near the bottom all yield the same mean stratification after relaxing the simulation with respect to kinetic energy, i.e., after about one solar hour.
The relaxed simulation was then evaluated every 1/20 of a sound crossing time, i.e. at an interval of 15.84 sec. This way 2527 samples have been produced covering a time of about 40012 sec or slightly more than 11 hours. The samples are not strictly equidistant in time, since the samples are picked from the dynamically varying timestepping of the simulation (variations introduced this way are in any case less than about 0.2 sec and most of the time below 0.1 sec). Statistical quantities have been computed for this output in a post processing step to calculate both horizontally averaged variables for each output step, i.e. for each of the 2527 samples, as well as ensembles averages from the horizontally averages quantities over the entire integration time of more than 126.3 sound crossing times. The data generated this way was used in the analysis presented in the following and the difference between the temperature gradient and the adiabatic gradient, as well as the Mach number, are displayed in Fig. 1 as a function of depth in the simulation. From averaging the vertical component of the radiative flux at the surface over the entire simulation time an effective temperature can be determined for the numerical simulation. Over the entire simulation time is found to be on average , with a negligibly small drift of per hour over that time and rare extrema not exceeding . The drift has been computed from a least square fit of a linear function to as determined at 0.1 Mm below the top of the simulation, where within less than 0.1%. Visual inspection and a comparison with a straight mean, which is more than lower than the mean of the least square fit line, confirm that the drift becomes smaller with time and it is in any case well within the range of fluctuations of determined for our simulation. We hence conclude the simulation to show an agreement with solar effective temperature. A value for the effective temperature of the Sun which is commonly used in recent literature is K. It can be derived from the results discussed in ChristensenDalsgaard (2009) and is also provided in Table 18.2 of Cox & Giuli (1968). However, values differing from this result by a few K depending on the specific measurements for radius and luminosity used are also common (see, e.g., Lebreton et al. 2008). This is certainly sufficiently accurate to investigate solar p modes. The surface gravity naturally remains fixed at its initial value of as does the chemical composition, , with its Grevesse & Noels (1993) metallicity distribution.
3 Normal radial modes of the simulation
The first step is to identify the normal modes of the simulation and thus to determine their characteristics throughout the simulation. To this end, we consider radial modes only so that, as described in Sect. 2, we consider horizontal averages of the simulation for the physical quantities.
3.1 Mode profiles and fitting procedure
We Fourier transform the timeseries described in Sect. 2 using a Fast Fourier Transform algorithm. As shown by Fig. 2, one can clearly distinguish three normal modes with a Lorentzian profile that is characteristic for solarlike oscillations. Indeed, in the time series, the vertical velocity of a radial solarlike mode can be written as
(1) 
where is the is the amplitude at . The observed signal is a sum of many such terms, each with their own amplitude and zeropoint of time, and also phase. is the radial displacement eigenfunction, is the pulsational eigenfrequency, is the time, is the damping rate, and the phase.
In the power spectrum, for , the Fourier transform of Eq. (1) can be approximated by a Lorentzian function such as (e.g. Baudin et al. 2005; Appourchaux 2014)
(2) 
and stands for the mode height, is the mode linewidth that is related to the mode damping rate through . Mode height is subsequently related to the mode amplitude and mode linewidth by (e.g. Samadi 2011)
(3) 
Therefore, a normal mode in the Fourier spectrum can be characterized by several quantities. First, through the global quantities (that do not vary with depth), i.e. frequency and mode linewidth. Second, through the mode height and phase which depend on the location in the simulation.
Mode  [Hz]  [Hz]  [Hz]  [Hz]  [Hz] 

1  2397.95  0.42       
2  3540.41  0.98  50.00  2.17  2.08 
3  4955.14  4.39  319.28  14.08  13.48 
We then fitted the power spectrum of the vertical component of the velocity by means of the maximumlikelihood estimator (see e.g. Toutain & Appourchaux 1994; Appourchaux et al. 1998). Each mode is fitted separately using a constant background. The global parameters, i.e. mode frequencies and linewidths, are obtained by performing a simultaneous fitting in several layers. For the fundamental mode of the simulation (hereafter mode 1), we consider all layers except near the upper and bottom boundaries.
The results for the frequencies and linewidths are summarized in Table 1.^{1}^{1}1Internal errors are computed from the Hessian matrix as described in Press et al. (2002). We note that the linewidth of mode 1 is not provided since it is lower than the actual resolution (which is about Hz). We also emphasize that although the values of the global parameters in Table 1 are precise, due to the relatively short duration of the simulation, it is difficult to obtain accurate results. Indeed, as shown for instance by Appourchaux (2014), the determination of mode linewidth is subject to many biases. For example, an overestimation of the mode height leads to an underestimation of the linewidth. Therefore, the values provided in Table 1 should be considered with care and are only intended to provide order of magnitudes.
To obtain the mode height and phase we consider the frequency bin with the largest power, near the eigenfrequency. We note that inferring those quantities from fits at each layer within the simulation is possible, but it provides much more noisy results. Figure 3 displays the mode velocity power densities (top panel) and mode phases (bottom panel) as a function of depth. The first mode corresponds to the fundamental mode with no node, while the second and third mode exhibit one and two nodes, respectively. We note that both for the mode velocity power densities and phases, there is a rapid variation near the peak of the superadiabatic gradient which is typical for nonadiabatic effects. Indeed, this feature is the result of a rapid variation of the entropy perturbations (e.g Belkacem et al. 2011; Samadi et al. 2012).
3.2 Comparison with adiabatic oscillations of a solar 1D model
We can even go a step further and identify the normal modes of the simulation with the observed modes. Indeed, as already shown by Stein & Nordlund (2001), it is possible to compare and identify the 3D modes with modes computed with a standard 1D model. To this end, one has to identify modes that exhibit a node at the bottom of the simulation and compare the mode velocity profiles. To do so, the velocity of the normal modes of the simulation is computed in the Fourier domain using Eq. (3) for the amplitude. Note that for the first mode, we used the spectral resolution instead of the linewidth since it is not resolved.
For consistency, we thus compute a 1D solar model that matches both the solar gravity and effective temperature but also the mean temperature at the bottom of the 3D numerical simulation. This model has been obtained using the CESTAM evolutionary code (Marques et al. 2013) assuming standard physics: Convection was included according to Canuto et al. (1996), with a mixinglength parameter , and turbulent pressure is ignored. Microscopic diffusion was included. The OPAL equation of state is assumed. The chemical mixture of the heavy elements is similar to that of Grevesse & Noels (1993)’s mixture. Subsequently, we constructed a 1D model following Trampedach (1997) as detailed in Samadi et al. (2008, 2010) in such a way that their outer layers are replaced by the averaged 3D simulations. Finally, 1D adiabatic oscillations are computed using the ADIPLS code (ChristensenDalsgaard 2011) and the ”gas ” hypothesis to account for the turbulent pressure (see Rosenthal et al. 1999; Sonoi et al. 2015, 2017).
Figure 4 (top panel) shows the comparison of normalized velocities and the resulting mode identification for modes 1 and 2. The eigenvelocities from the 1D computation have been normalized so that their kinetic energy equals the kinetic energy of the corresponding mode in the 3D simulation. It turns out that the three modes of the 3D simulation correspond to the observed modes with respective radial orders , , and . There is a very good match between the velocity profiles of both modes 1 and 2 as a function of radius compared to the adiabatic 1D computation. The main differences occur in the atmospheric layers since the upper boundary of the 1D model and the 3D simulation are different as well as near the peak of the superadiabatic gradient (i.e. at the bottom of the photosphere). In the latter region, the sharp variations exhibited by the 3D modes are the result of purely nonadiabatic effects so that adiabatic computations are unable to reproduce those features. As it will be shown in the following sections, these patterns are important for investigating the physics of mode damping. In the uppermost region, the difference between the 1D and 3D models can be attributed to the boundary condition of the 3D simulation, which forces a node for the normal modes.
In Fig. 4 (bottom panel), we illustrate the mode identification by showing the solar damping rates as a function of frequency. The modes that correspond to the normal modes of the simulation are overplotted in red and it appears that modes 1 and 2 of the 3D simulation bracket the frequency of the maximum height (Hz) while mode 3 corresponds to a mode near the cutoff frequency. Since the relative importance of the various physical contributions to the damping rates is expected to vary with frequency (e.g. Balmforth 1992; Belkacem et al. 2012), this will allow us to probe different physical regimes with respect to the mode damping. We also note that the frequencies of the modes in the simulation and the frequencies of the observed modes are comparable but not exactly the same. This is a consequence of the location of the bottom of the simulation and, to a lesser extent, to nonadiabatic effects (e.g. Sonoi et al. 2017).
3.3 Scaling for mode amplitudes
Going beyond the shape of the eigenfunctions, one can expect that the mode physics in the simulation is realistic enough to gain some insight into the physical mechanisms responsible for mode damping. Indeed, mode damping occurs in the uppermost layers of the stars and more precisely when the thermal timescale becomes equal or higher than the modal period (see for instance Belkacem et al. 2011, 2012). This occurs in the superadiabatic layers and in the atmosphere, but in the quasiadiabatic layers, mode damping and driving are negligible. Consequently, the extent of the simulation appears to be sufficient to investigate the physical mechanisms responsible for mode damping. However, at first sight, the normal modes of the 3D numerical simulation cannot be directly compared to the observed solar oscillations or to the computed solar oscillations from a 1D model. The fundamental difference is the size of the resonant cavity, which modifies the mode masses (see Eq. 5) as well as the large separation. Indeed, while for the Sun the large separation is Hz, we found Hz for the 3D simulation by computing the firstorder asymptotic expression of the large separation (e.g. Tassoul 1980).
Using order of magnitude estimates, it is nevertheless possible to obtain a relation between the mode amplitude in the 3D simulation and in the Sun. To do so, let us first define the energy of a mode as (e.g. Samadi 2011; Belkacem & Samadi 2013; Samadi et al. 2015)
(4) 
where is the mode velocity observed at the surface of the Sun and the mode mass defined as
(5) 
where is the eigenvelocity.
To go further, we note that the eigenvelocity amplitude becomes important near the surface layers so that its integration over the vertical coordinate is similar within the 3D box and in the 1D model. More precisely, using the 1D model, the relative contribution of the mode mass in the upper most layers corresponding to the 3D simulation domain is about , , and for modes 1, 2, and 3 respectively. Consequently, in the 3D and 1D models, they mainly differ due to the horizontal integral, so that the mode masses per unit surface area are approximately similar in the Sun and in the 3D simulation. Therefore
(6) 
with
(7) 
and
(8) 
where is the solar radius, is the horizontal surface of the simulation, is the horizontally and timeaveraged density, the depth of the 3D simulation, and the subscripts and stand for the modes of the 3D simulation and the solar modes, respectively. By considering the second mode of the 3D simulation (which is resolved and has a sufficient signal to noise ratio) and the corresponding mode in our solar 1D model (see Sect. 3.2), we found that Eq. (6) is verified to a few percent.
Now using Eq. (4) together with Eq. (6), one obtains
(9) 
We also note that the mode energy can be written as (see Samadi et al. 2015, for details)
(10) 
where is the excitation rate. Consequently, mode energy is independent of the mode mass (or mode inertia) since both the excitation and damping rates can be written to be inversely proportional to the mode mass. Now, if the 3D simulation is realistic enough, one can hence expect that the energies of a mode in the 3D simulation and its equivalent in the Sun are similar (i.e. ). We checked this relation by computing directly the mode energies in the 3D simulation and in the solar model using Eq. (4). For the 3D simulation, the mode velocity at the surface is computed directly using the Fourier transform of the velocity and for the solar surface velocity we consider from Baudin et al. (2005). It gives, using Eq. (3) and Eq. (4), J and J (for ), which is in a reasonable agreement given the uncertainties of our fit^{2}^{2}2We also note that the observation made with the GOLF instrument is obtained at an altitude much higher than the photosphere. Therefore, the mode energy of the observed mode at the photosphere is expected to be lower (see Kjeldsen et al. 2008, at the frequency of the maximum power). (see Sect. 3.1).
Finally, using the aforementioned argument and using Eq. (9), we get
(11) 
Solar mode amplitudes are typically between m.s and m.s near the frequency of the maximum amplitude and using Eq. (11) one recovers the order of magnitude of the amplitudes of the normal modes in the 3D simulation, which are about m.s. Obviously, this estimate yields a very rough order of magnitude and differences in mode physics between the 3D model and the observations are still to be expected, for instance, due to the boundary conditions of the simulation. Nevertheless, this estimate favors the idea that mode damping and mode driving occurring within the 3D simulation is worth being considered for constraining the underlying physical processes.
4 Computation of the work integral
Mode damping can be directly inferred by fitting the modes in the power spectrum of the vertical velocity. However, if one wants to go further and decipher the contributions to this damping, it is necessary to compute the work integral. Indeed, such an approach permits us to explicitly split the contributions but also to infer information on the location of both driving and damping regions in the simulations.
4.1 Averaged equations
The first step consists in averaging the primitive equations. To do so, as we are considering a compressible flow, we will consider both Reynolds and mass weighted averages (e.g. Canuto 1997a; Nordlund & Stein 2001). In the following, we will approximate the Reynolds average by the straight horizontal average in the 3D simulation. Therefore, any quantity , can be decomposed such that
(12) 
This average will be applied to the density, pressure, radiative and convective fluxes.
The second average is named as the mass weighted or Favre average (Favre 1969), so that for a quantity , it is defined as
(13) 
The quantity can thus be decomposed such that
(14) 
It immediately follows
(15) 
A more detailed description of the properties of the Favre average is provided in Appendix. A.1. Notice that in our case, the mass average will be applied to the velocity, and entropy fields (see Canuto 1997a; Nordlund & Stein 2001, for details). Because we consider a compressible flow, this choice permits us to simplify the equations. Indeed, many correlation terms involving density fluctuations are incorporated into the Favre mean quantities and consequently are no longer present in the governing equations by using the Favre average for the abovementioned quantities.
Subsequently, we average the mass conservation equation (see Appendix. A.2.1 for a detailed derivation) to obtain
(16) 
where is the density, is the vertical component of the velocity. We also introduced the pseudoLagrangian derivative such that
(17) 
Applying the same procedure, and omitting the terms that cancel from hydrostatic equilibrium, we get from the momentum conservation equation (see Appendix. A.2.1 for details)
(18) 
where is the turbulent pressure, is the gas pressure, is the gravitational field that is considered constant as in the 3D simulation. In addition, for any quantity , we have defined with the time average of . is therefore the pseudoLagrangian perturbation of . Finally, we introduced the notation .
4.2 Integral expression of the damping rates
To get insight into the physics of mode damping, it is necessary to determine the phase lag between the Lagrangian perturbations of pressure and density for a given mode within the simulation and the integral of this phase lag provides the total damping rate for that mode (e.g. Samadi et al. 2015). This is what we call the integral approach. Therefore, we will work in the time Fourier space. To start, we thus consider the time Fourier transformation of Eq. (18) and we multiply it by to obtain
(19) 
where is the total fluctuation of pressure, the symbol stands for the temporal Fourier transform, the symbol stands for the complex conjugate, and is the cyclic frequency.^{3}^{3}3Notice that (the cyclic frequency in the time Fourier domain) is not to be confused with (the modal frequency). Note that the tilde has been omitted for ease of notation. This equation can be further simplified if one keeps only the dominant order in the LHS of Eq. (18). This is equivalent to assuming in the RHS of Eq. (19) such that
(20) 
To go further, we integrate over the mean mass column density (defined as ), perform integration by part and multiply by to finally obtain
(21) 
where
(22) 
Note that at the modal frequencies, can be identified with the mode inertia.
It is now useful to take advantage of the Fourier transform of the mass conservation equation that reads
(23) 
where again , so that Eq. (21) leads to
(24) 
where
(25)  
(26)  
(27) 
We assumed that near the modal frequency (i.e. ) the coherent contributions (associated with the oscillation) of the pressure and density fluctuations dominate over the random contributions (associated with the turbulence). The eigenfrequency is complex , with the damping rate, and in Eq. (24) we assumed since we are in the situation for which .
At this step, several comments are necessary. The second term of Eq. (24), i.e. , vanishes because is null at the bottom of the simulation box and both and tend to zero at the upper boundary. We checked numerically that this contribution is negligible as well as the contribution of the third term (). Consequently, Eq. (24) reduces to
(28) 
Thus, the damping rate is determined by the phase lag between mode compressibility and pressure fluctuations. This is a classical result in 1D nonadiabatic calculations (e.g. Ledoux & Walraven 1958; Unno et al. 1989).
An alternative approach, as proposed by Nordlund & Stein (2001), is to split the pressure fluctuations in an adiabatic and nonadiabatic component, i.e.
(29) 
where
(30) 
with the horizontal and time averaged squared sound speed (see Appendix in Samadi et al. 2007). The adiabatic pressure fluctuations are dominant over the nonadiabatic ones but do not contribute to the damping. This can be easily seen by inserting Eq. (29) into Eq. (28). Consequently, splitting explicitly the adiabatic and nonadiabatic contribution can help to improve the accuracy of the computation. However, we checked that using either the total or the nonadiabatic pressure fluctuations makes almost no difference when computing Eq. (28).
4.3 Computation of the cumulative damping
We then computed Eq. (28). To do so, the pseudoLagrangian quantities () are computed by interpolating the physical variables () onto the timeaveraged mean column density () before subtracting their timeaveraged values. Then, the amplitudes and phases of each quantity () are taken from the Fourier transforms. Eventually, the damping rates are obtained at the frequency of the modes. More precisely;

For the first mode, since it is not resolved, we consider the bin at the maximum height of the mode. The result is Hz, which is consistent with our upper limit (i.e. the resolution of Hz) because of the relatively short duration of the simulation.

For the second mode, we adopt the same approach by selecting the bin with the highest amplitude in the Fourier spectrum because this bin is the least affected by the noise. We find Hz, which is quite close to the value found from fitting the modepeak (see Table 1).

For the third mode, the situation is more complex because of its large width. In this case, and due to the low signal to noise ratio, is highly varying from bin to bin preventing us from making conclusions.
Finally, given the inherent limitations of the simulation, one can conclude that for the first two modes the damping rates computed from Eq. (28) and the measured damping rates are roughly consistent (see Table 2). However, for the third one, it is not possible to draw conclusions. Indeed, one can easily expect large uncertainties on the fitted mode linewidths because of the relatively short duration of the simulation compared to the fifteen years of continuous observation of the Sun by the GOLF instrument. The same is true for the damping rate computed using the integral expression (Eq. 28) because near the bottom of the simulated box the imaginary part of mode perturbations become smaller than stochastic convective fluctuations so that the estimate of the phase lag between the nonadiabatic pressure and the density is highly affected by the noise.
Mode  [Hz]  (fit) [Hz]  (work integral) [Hz] 

1  2397.95  25  5.13 
2  3540.41  50.00  59.75 
3  4955.14  319.28   
Figure 5 displays the cumulative damping for modes 1 and 2. We stress that when cumulative damping increases (decreases) outward there is a stabilizing effect (destabilizing effect). This convention will be used in the following. For both modes one can distinguish three regions, namely;

The inner quasiadiabatic region, for which , stabilizes the modes. In this region, however, the cumulative damping is noisy since it exhibits some oscillations. This is due to the Fourier transform (at the modal frequencies) of the pseudoLagrangian density which is affected by the noncoherent turbulent fluctuations in this region. This effect is small but as we are considering phase differences it has an important impact on the results.

The atmosphere stabilizes the mode 2 while the cumulative damping is almost constant for mode 1. We also note that near the upper boundary () the cumulative damping suddenly increases. This is related to the rapid decrease of the eigenfunctions (see Fig. 4, top panel). It is likely an effect of the boundary condition of the simulation.

Finally, the superadiabatic region, which corresponds to the region of hydrogen ionisation (), destabilizes. Indeed, this region corresponds to the region where the opacity mechanism related to the ionisation of hydrogen is effective. In the eigenfunctions, this can be seen to be responsible for the rapid variations of the eigenvelocities (see Fig. 4, top panel)
Despite of the effect of the noncoherent turbulent fluctuations in the quasiadiabatic region, the behavioural patterns of the cumulative damping (as described above) are in qualitative agreement with 1D nonadiabatic calculation (e.g. Balmforth 1992; Belkacem et al. 2011).
5 Contributions to the mode damping
As shown in the previous section, the measured and computed damping rates are found consistent even if the inherent limitations of the simulation prevents us from getting a perfect match. This encourages us to go further by disentangling the different physical contributions to this damping rate.
5.1 The role of gas and turbulent pressure fluctuations
Let us start by splitting the perturbation of total pressure appearing in Eq. (28) as
(31) 
where is the perturbation of gas pressure and the perturbation of turbulent pressure. Therefore, Eq. (28) can be rewritten such that
(32) 
where
(33)  
(34) 
For these individual contributions to , shown in Fig. 6, we stress again that for the results must be considered with care because of the fluctuations of the density perturbation but also of the perturbation of the turbulent pressure. The Fourier spectra of turbulent pressure are rather noisy, even at modal frequencies, making it difficult to conclude that the signal is dominated by the coherent oscillating signal. Notwithstanding these words of caution, we suggest to distinguish between mainly two regions (for ).
First, in the atmospheric layers both modes and both contributions behave the same way, i.e. the cumulative damping is almost neutral or only slightly damping. This is due to the fact that, for the contribution associated with the gas pressure, the modal period is much longer than the local thermal timescale. Consequently, the medium adapts almost instantaneously to any perturbation so that the total energy flux is frozen (see Samadi et al. 2015, for a detailed explanation). For the contribution associated with the turbulent pressure, its contribution is very small in the uppermost layers and turns destabilizing in the innerlayers. This behavior is in qualitative agreement with previous findings by Balmforth (1992) and Sonoi et al. (2017). The latter had shown that the turbulent pressure contribution is mainly controlled by both the Mach number and the ratio of the local convective timescale to the modal frequency. In the atmospheric layers, both factors are small but increase towards the inner layers to become nonnegligible for .
Secondly, in the superadiabatic layers (and more precisely for ), the modal period is of the same order of magnitude as the local thermal timescale. Consequently, it corresponds to the region where the destabilization effect due to the mechanism associated with the ionisation of hydrogen is at work. It also corresponds to the location of the rapid variation of the eigenvelocity shown by Fig. 4. This explains the destabilizing effect of the gas pressure contribution. For the work associated with the turbulent pressure modulation, the situation is equivalent in this region, i.e. it destabilizes the modes. In contrast, in the inner layers (i.e. for mode 1 and for mode 2), excitation by competes against damping by , resulting in a net stabilizing effect for both modes. We note that those behavioural patterns support previous findings using a timedependent, mixinglength formulation of convection (e.g. Fig. 14 in Balmforth 1992).
Figures 7 show the contributions to the mode damping from the perturbation of the gas pressure (Eq. 33) and the turbulent pressure (Eq. 34). It is worth mentioning that, for the two considered normal modes, both the gas and turbulent pressure contributions have an overall stabilizing effect. For mode 1, the contribution of the turbulent pressure is slightly dominant while the two contributions are of the same order of magnitude for mode 2. This is in line with previous findings (see Houdek & Dupret 2015, for a detailed discussion) that the role of the perturbation of turbulent pressure is an essential ingredient for stabilizing the modes. However, our result further suggests that the contribution related to the perturbation of the gas pressure also stabilizes the modes and that its contribution becomes dominant, compared to the contribution of turbulent pressure, for highfrequency modes.
5.2 The contributions to the gas pressure
To disentangle the various terms that contribute to the perturbation of gas pressure, it is first necessary to use the perturbed equation of state
(35) 
with
(36) 
Since the adiabatic part of the gas pressure perturbation (second term of Eq. 35) cancels in the integrand of Eq. (28), the gas component of the damping rate becomes
(37) 
As a consistency check, we computed Eq. (37) and recover the values obtained for the gas contribution of the damping as provided by Eq. (28).
To go further, it is necessary to write down the equation governing the perturbation of entropy. This is obtained from the equation of conservation of the internal energy, which after some manipulations permits us to get to the lowest order (see Appendix. A.2.3 for details)
(38) 
where is the specific entropy, is the vertical component of the radiative flux, is the vertical component of the convective flux defined by
(39) 
where is the specific internal energy and the perturbation of dissipation rate of turbulent kinetic energy into heat () is defined by
(40) 
with standing for viscous dissipation. Equation (38) is the perturbation of the energy equation and is similar to what is obtained from the linear perturbation theory (e.g. Ledoux & Walraven 1958; Ledoux 1958; Grigahcène et al. 2005). It emphasizes that the perturbation of gas pressure is the combination of three main physical ingredients, namely the perturbation of the radiative flux, the perturbation of the convective flux and the perturbation of dissipation rate of turbulent kinetic energy into heat.
Using Eq. (37) together with Eq. (38), finally gives
(41) 
where is given by Eq. (34) and
(42)  
(43)  
(44) 
The last term () balances the contribution of the turbulent pressure in the quasiadiabatic region as shown by Ledoux & Walraven (1958) and Grigahcène et al. (2005). It is thus a nonnegligible contribution to the total damping rate but is difficult to estimate directly from the simulation. There are two main reasons, one which is method dependent while the other one is related to the physics of dissipation of turbulent kinetic energy. To compute the former in the case of ANTARES would require to evaluate the dissipation through the (SmagorinskyLilly type) subgrid scale model (see Muthsam et al. 2010, and in particular Sect. 2.6 of Mundprecht et al. 2013). But this were only a lower limit, since shocks and steep gradients in general would be dealt with by a local, nonlinear viscosity which is built into the weighted essentially nonoscillatory scheme used by ANTARES (Muthsam et al. 2010). One of the advantages of the latter is that it aims at minimizing the amount of viscosity added by the scheme to the numerical solution. On the other hand, it is difficult to accurately quantify the exact amount of numerical viscosity introduced this way, at least from a single simulation. A more fundamental, physical problem is that the dissipation of turbulent kinetic energy is dominated by variations of the numerical solution close to the grid scale. Computed values are thus very sensitive to both the numerical method used and to the resolution of the simulation. Since kinetic energy is dissipated down to scales orders of magnitudes smaller than can reasonably be resolved in a simulation of solar convection (a classical result, see also Kupka & Muthsam 2017, for references and estimates), any direct computation of this quantity is necessarily inaccurate. To get a better insight would require an expensive series of simulations down to very high resolutions to see if some simple scaling estimates were adequate.
An alternative which provides a first approximation and is available from the present calculations is to estimate from its definition, i.e.,
(45) 
We use this relation in the following, although we will focus on the radiative and convective flux contributions to the damping rates, i.e. Eq. (42) and Eq. (43), respectively.
Figures 8 display the contributions of the divergence of the radiative and convective fluxes to the damping for modes 1 and 2. For both modes, the behavioral patterns are the same for both the perturbation of radiative and convective fluxes. In the superadiabatic region, the perturbation of the divergence of the radiative flux stabilizes the modes while the perturbation of the divergence of the convective flux destabilizes them. In the atmospheric layers, the situation is reversed, except for the upper layers, where contributions from both fluxes are close to zero, and the layers at the very top, which are influenced by the boundary conditions. It is worth noting that both contributions are quite large in absolute values but they compensate each other so that the resulting contribution remains small. In addition, if we do not consider the uppermost layers in which the effect of the boundary conditions of the numerical simulation are important, the convective contribution dominates over the radiative one. The total contribution to the damping rate is thus a residual that results from a balance between the two flux divergences (see Houdek & Dupret 2015, for a detailed discussion on this issue). Finally, we note that since the contribution of the two fluxes nearly cancels each other, the role of the perturbation of the dissipation rate of turbulent kinetic energy into heat becomes important. This is shown in Figs. 8. Indeed, these contributions seem to dominate the effect of the perturbation of gas pressure on the mode damping as they essentially destabilize the modes. But as already mentioned, it is difficult to estimate it directly from the numerical simulation, therefore Eq. (45) is used to circumvent the problem and as such must be considered with care. Nevertheless, the fact that the contributions of the perturbation of the radiative and convective fluxes nearly cancel each other still emphasizes the essential role of the perturbation of the dissipation rate into heat.
6 Concluding remarks
Using a 3D hydrodynamical simulation of the Sun computed with the ANTARES code and with a timeduration of about 11 hrs, we have shown that it is possible to identify at least three radial normal modes in the simulation, two of which were usable for our purposes. Those modes have been shown to have properties similar to the normal modes of 1D solar models that happen to have a node at the bottom of the simulation domain. In contrast, the amplitudes of simulation modes are found to be much higher than in the Sun due to the relatively very small horizontal area of the simulation. Assuming that the physical background experienced by those modes is realistic enough, we have demonstrated that it is possible to gain some insight into the physics governing the mode damping rates.
For the two first normal modes of the simulation, it has been possible to investigate the work related to their damping. Except for the quasiadiabatic region of the simulation, for which the ratio of the mode amplitude compared to the turbulent noise is not large enough, we have been able to exhibit the different regions in which the modes are destabilized and stabilized. Going further, we disentangled the respective role of both the perturbation of the gas and turbulent pressure. From a qualitative point of view our results are in good agreement with previous findings (e.g. Balmforth 1992; Belkacem et al. 2012). However, from a quantitative point of view, it appears that both contributions have an overall stabilizing effect. Indeed, in contrast to previous results based on timedependent extensions of the 1D, mixinglength formulation, the relative contribution associated with the perturbation of the turbulent pressure is not found to be always dominant over the contribution associated with the perturbation of the gas pressure. In addition, it has also been possible to gain insight into the contributions of the perturbation of the divergence of both the radiative and convective fluxes. It appears that those two contributions nearly cancel each other both in the atmospheric and in the superadiabatic layers. Indeed, while each of them has an important absolute contribution, the sum of the two remains small since they tend to cancel each other. However, it appears that while the radiative contribution destabilizes the modes, the convective one stabilizes the mode. The latter is found to be dominant, even if only slightly.
Consequently, we have shown that investigating the properties of the normal modes of a 3D simulation is of great help to gain some physical insight into the physics of mode damping. From this work, such an approach is found to be feasible. However, the main limitation has been found to be the duration of the simulation. Indeed, it is a crucial point which affects our ability to provide reliable quantitative estimates. A much longer simulation is thus highly desirable. First, it would ensure that the mode of the lowest frequency is resolved, i.e. that the linewidth of the mode is higher than one bin in the Fourier domain. It would also ensure that high frequency modes with a very large linewidth could be fitted properly. Second, a longer duration of the simulation is important to improve the ratio between the amplitudes of the normal modes and the turbulent noise in the Fourier spectrum. This is particularly important in the innermost layers of the simulation. Indeed, mode 2 can already be successfully analysed in detail with the present simulations though of course a longer simulation would be desirable to do the same for mode 1 at and possibly also for mode 3.
We conclude that using 3D hydrodynamical simulations to investigate the physics of mode damping reveals to be a promising approach. A drawback is that such an approach is highly demanding in terms of computational efforts because, to get precise and accurate constraints, one would need simulations with a timeduration of several days. This is certainly an objective to attain for our future works on the issue of mode damping. Another important issue that will certainly deserve further work is related to the top boundary conditions. Indeed, in the very top layers of the simulation the normal modes must be considered with care. This emphasises the need of investigating systematically the influence of boundary conditions on the normal modes of the simulation.
Finally, we emphasize that the estimate of the influence of the dissipation on the oscillation remains to be consolidated. Indeed, due to the ENO scheme employed by ANTARES, it is not possible to directly and properly estimate the dissipation within the simulation. A future work dedicated to this issue is thus highly desirable. A possible approach would be to include artificial viscosities in ANTARES and compute the advective fluxes from standard centred stencil. Such a procedure, while being less accurate for the modeling of the dissipation, presents the advantage of permitting to quantify directly its impact on the oscillations and thus to verify that our indirect estimate is valid.
Acknowledgements.
F. Kupka is grateful for support through Austrian Science Fund (FWF) projects P21742N16, P25229N27, and P29172N27. K.B. and R.S. acknowledge support from the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, France.References
 Aarslev et al. (2018) Aarslev, M. J., Houdek, G., Handberg, R., & ChristensenDalsgaard, J. 2018, MNRAS, 478, 69
 Ando & Osaki (1975) Ando, H. & Osaki, Y. 1975, PASJ, 27, 581
 Appourchaux (2014) Appourchaux, T. 2014, A crash course on data analysis in asteroseismology, ed. P. L. Pallé & C. Esteban, 123
 Appourchaux et al. (2012) Appourchaux, T., Benomar, O., Gruberbauer, M., et al. 2012, A&A, 537, A134
 Appourchaux et al. (1998) Appourchaux, T., Gizon, L., & RabelloSoares, M.C. 1998, A&AS, 132, 107
 Baglin et al. (2006a) Baglin, A., Auvergne, M., Barge, P., et al. 2006a, in ESA Special Publication, Vol. 1306, ESA Special Publication, ed. M. Fridlund, A. Baglin, J. Lochard, & L. Conroy, 33
 Baglin et al. (2006b) Baglin, A., Auvergne, M., Boisnard, L., et al. 2006b, in COSPAR Meeting, Vol. 36, 36th COSPAR Scientific Assembly, 3749
 Balmforth (1992) Balmforth, N. J. 1992, MNRAS, 255, 603
 Baudin et al. (2011a) Baudin, F., Barban, C., Belkacem, K., et al. 2011a, A&A, 529, A84
 Baudin et al. (2011b) Baudin, F., Barban, C., Belkacem, K., et al. 2011b, A&A, 535, C1
 Baudin et al. (2005) Baudin, F., Samadi, R., Goupil, M.J., et al. 2005, A&A, 433, 349
 Belkacem et al. (2012) Belkacem, K., Dupret, M. A., Baudin, F., et al. 2012, A&A, 540, L7
 Belkacem et al. (2011) Belkacem, K., Goupil, M. J., Dupret, M. A., et al. 2011, A&A, 530, A142
 Belkacem & Samadi (2013) Belkacem, K. & Samadi, R. 2013, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 865, Lecture Notes in Physics, Berlin Springer Verlag, ed. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green, 179
 Belkacem et al. (2010) Belkacem, K., Samadi, R., Goupil, M. J., et al. 2010, A&A, 522, L2
 Belkacem et al. (2006a) Belkacem, K., Samadi, R., Goupil, M. J., & Kupka, F. 2006a, A&A, 460, 173
 Belkacem et al. (2006b) Belkacem, K., Samadi, R., Goupil, M. J., Kupka, F., & Baudin, F. 2006b, A&A, 460, 183
 Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
 Canuto (1992) Canuto, V. M. 1992, ApJ, 392, 218
 Canuto (1997a) Canuto, V. M. 1997a, ApJ, 482, 827
 Canuto (1997b) Canuto, V. M. 1997b, ApJ, 489, L71
 Canuto et al. (1996) Canuto, V. M., Goldman, I., & Mazzitelli, I. 1996, ApJ, 473, 550
 Chaplin et al. (2009) Chaplin, W. J., Houdek, G., Karoff, C., Elsworth, Y., & New, R. 2009, A&A, 500, L21
 Chaplin & Miglio (2013) Chaplin, W. J. & Miglio, A. 2013, ARA&A, 51, 353
 ChristensenDalsgaard (2009) ChristensenDalsgaard, J. 2009, in IAU Symposium, Vol. 258, The Ages of Stars, ed. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse, 431–442
 ChristensenDalsgaard (2011) ChristensenDalsgaard, J. 2011, ADIPLS: Aarhus Adiabatic Oscillation Package (ADIPACK), astrophysics Source Code Library
 ChristensenDalsgaard et al. (1996) ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286
 Cox & Giuli (1968) Cox, J. P. & Giuli, R. T. 1968, Principles of stellar structure, ed. Cox, J. P. & Giuli, R. T.
 Dupret et al. (2009) Dupret, M., Belkacem, K., Samadi, R., et al. 2009, A&A, 506, 57
 Dupret et al. (2006) Dupret, M. A., Barban, C., Goupil, M.J., et al. 2006, in ESA Special Publication, Vol. 624, Proceedings of SOHO 18/GONG 2006/HELAS I, Beyond the spherical Sun
 Favre (1969) Favre, A. 1969, SIAM, Philadelphia, 231
 Fehlberg (1970) Fehlberg, E. 1970, Computing, 6, 61
 Goldreich & Kumar (1991) Goldreich, P. & Kumar, P. 1991, ApJ, 374, 366
 Gough (1980) Gough, D. 1980, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 125, Nonradial and Nonlinear Stellar Pulsation, ed. H. A. Hill & W. A. Dziembowski, 273–299
 Grevesse & Noels (1993) Grevesse, N. & Noels, A. 1993, in Origin and Evolution of the Elements, ed. N. Prantzos, E. VangioniFlam, & M. Casse, 15–25
 Grigahcène et al. (2005) Grigahcène, A., Dupret, M.A., Gabriel, M., Garrido, R., & Scuflaire, R. 2005, A&A, 434, 1055
 GrimmStrele et al. (2015a) GrimmStrele, H., Kupka, F., LöwBaselli, B., et al. 2015a, New A, 34, 278
 GrimmStrele et al. (2015b) GrimmStrele, H., Kupka, F., & Muthsam, H. J. 2015b, Comp. Phys. Comm., 188, 7
 Grosjean et al. (2014) Grosjean, M., Dupret, M.A., Belkacem, K., et al. 2014, ArXiv eprints
 Happenhofer et al. (2013) Happenhofer, N., GrimmStrele, H., Kupka, F., LöwBaselli, B., & Muthsam, H. 2013, Journal of Computational Physics, 236, 96
 Hekker & ChristensenDalsgaard (2017) Hekker, S. & ChristensenDalsgaard, J. 2017, A&A Rev., 25, 1
 Heun (1900) Heun, K. 1900, Z. Math. Phys., 45, 23
 Houdek (2012) Houdek, G. 2012, ArXiv eprints (arXiv:1201.0194), to be published in: Proc. 61st Fujihara Seminar ”Progress in solar/stellar physics with helio and asteroseismology”, Shibahashi H., ed.
 Houdek (2017) Houdek, G. 2017, in European Physical Journal Web of Conferences, Vol. 160, European Physical Journal Web of Conferences, 02003
 Houdek et al. (1999) Houdek, G., Balmforth, N. J., ChristensenDalsgaard, J., & Gough, D. O. 1999, A&A, 351, 582
 Houdek & Dupret (2015) Houdek, G. & Dupret, M.A. 2015, Living Reviews in Solar Physics, 12, 8
 Houdek et al. (2017) Houdek, G., Trampedach, R., Aarslev, M. J., & ChristensenDalsgaard, J. 2017, MNRAS, 464, L124
 Iglesias & Rogers (1996) Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943
 Jiang & Shu (1996) Jiang, G.S. & Shu, C.W. 1996, Journal of Computational Physics, 126, 202
 Kjeldsen et al. (2008) Kjeldsen, H., Bedding, T. R., Arentoft, T., et al. 2008, ApJ, 682, 1370
 Kraaijevanger (1991) Kraaijevanger, J. 1991, BIT, 31, 482
 Kupka et al. (2012) Kupka, F., Happenhofer, N., Higueras, I., & Koch, O. 2012, Journal of Computational Physics, 231, 3561
 Kupka & Muthsam (2007a) Kupka, F. & Muthsam, H. J. 2007a, in IAU Symposium, Vol. 239, Convection in Astrophysics, ed. F. Kupka, I. Roxburgh, & K. L. Chan, 80–82
 Kupka & Muthsam (2007b) Kupka, F. & Muthsam, H. J. 2007b, in IAU Symposium, Vol. 239, Convection in Astrophysics, ed. F. Kupka, I. Roxburgh, & K. L. Chan, 83–85
 Kupka & Muthsam (2017) Kupka, F. & Muthsam, H. J. 2017, Living Reviews in Computational Astrophysics, 3, 1
 Kupka & Robinson (2007) Kupka, F. & Robinson, F. J. 2007, MNRAS, 374, 305
 Kurucz (1993a) Kurucz, R. 1993a, ATLAS9 Stellar Atmosphere Programs and 2 km/s grid. Kurucz CDROM No. 13. Cambridge, Mass.: Smithsonian Astrophysical Observatory, 1993., 13
 Kurucz (1993b) Kurucz, R. 1993b, Opacities for Stellar Atmospheres: [+0.0],[+0.5],[+1.0]. Kurucz CDROM No. 2. Cambridge, Mass.: Smithsonian Astrophysical Observatory, 1993., 2
 Lebreton et al. (2008) Lebreton, Y., Monteiro, M. J. P. F. G., Montalbán, J., et al. 2008, Ap&SS, 316, 1
 Ledoux (1958) Ledoux, P. 1958, Handbuch der Physik, 51, 605
 Ledoux & Walraven (1958) Ledoux, P. & Walraven, T. 1958, Handbuch der Physik, 51, 353
 Marques et al. (2013) Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74
 Mihalas & Mihalas (1984) Mihalas, D. & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics (Oxford University Press)
 Mosser et al. (2014) Mosser, B., Benomar, O., Belkacem, K., et al. 2014, A&A, 572, L5
 Mosser et al. (2012) Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10
 Mundprecht et al. (2013) Mundprecht, E., Muthsam, H. J., & Kupka, F. 2013, MNRAS, 435, 3191
 Muthsam et al. (2010) Muthsam, H. J., Kupka, F., LöwBaselli, B., et al. 2010, New Astronomy, 15, 460
 Nordlund & Stein (2001) Nordlund, Å. & Stein, R. F. 2001, ApJ, 546, 576
 Press et al. (2002) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical recipes in C++ : the art of scientific computing
 Rogers et al. (1996) Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902
 Rosenthal et al. (1999) Rosenthal, C. S., ChristensenDalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689
 Samadi (2011) Samadi, R. 2011, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 832, Lecture Notes in Physics, Berlin Springer Verlag, ed. J.P. Rozelot & C. Neiner, 305
 Samadi et al. (2012) Samadi, R., Belkacem, K., Dupret, M.A., et al. 2012, A&A, 543, A120
 Samadi et al. (2008) Samadi, R., Belkacem, K., Goupil, M. J., Dupret, M.A., & Kupka, F. 2008, A&A, 489, 291
 Samadi et al. (2015) Samadi, R., Belkacem, K., & Sonoi, T. 2015, in EAS Publications Series, Vol. 73, EAS Publications Series, 111–191
 Samadi et al. (2007) Samadi, R., Georgobiani, D., Trampedach, R., et al. 2007, A&A, 463, 297
 Samadi et al. (2010) Samadi, R., Ludwig, H.G., Belkacem, K., Goupil, M. J., & Dupret, M.A. 2010, A&A, 509, A15
 Samadi et al. (2003) Samadi, R., Nordlund, Å., Stein, R. F., Goupil, M. J., & Roxburgh, I. 2003, A&A, 404, 1129
 Shu (2003) Shu, C.W. 2003, International Journal of Computational Fluid Dynamics, 17, 107
 Shu & Osher (1988) Shu, C.W. & Osher, S. 1988, Journal of Computational Physics, 77, 439
 Sonoi et al. (2017) Sonoi, T., Belkacem, K., Dupret, M.A., et al. 2017, A&A, 600, A31
 Sonoi et al. (2015) Sonoi, T., Samadi, R., Belkacem, K., et al. 2015, A&A, 583, A112
 Stein & Nordlund (1998) Stein, R. F. & Nordlund, A. 1998, ApJ, 499, 914
 Stein & Nordlund (2001) Stein, R. F. & Nordlund, Å. 2001, ApJ, 546, 585
 Tassoul (1980) Tassoul, M. 1980, ApJS, 43, 469
 Toutain & Appourchaux (1994) Toutain, T. & Appourchaux, T. 1994, A&A, 289, 649
 Trampedach (1997) Trampedach, R. 1997, PhD thesis, Aarhus University (1997)
 Unno et al. (1989) Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars (University of Tokyo Press, 1989, 2nd ed.)
 Vrard et al. (2018) Vrard, M., Kallinger, T., Mosser, B., et al. 2018, ArXiv eprints
 Xiong et al. (2000) Xiong, D. R., Cheng, Q. L., & Deng, L. 2000, MNRAS, 319, 1079
Appendix A Averaged equations
a.1 Reynolds and mass (or Favre) averages
Two types of horizontal averages are defined, the first of which is commonly known as the Reynolds average. In the following we will approximate the Reynolds average by the straight horizontal average in the 3D simulation. Therefore, any quantity , is decomposed such that
(46) 
so that we obviously have
(47) 
The second average is called the mass or Favre average, so that for a quantity , it is defined as
(48) 
with
(49) 
It immediately follows
(50)  
(51) 
From the above equations, it is quite straightforward to derive the following relations
(52)  
(53)  
(54) 
which will be used in the following.
a.2 Mean equations
The Reynolds averages will be applied to the density and pressure while mass average will be applied to the velocity, temperature and entropy fields.
a.2.1 Mass conservation
(55) 
where stands for the density and the component of the velocity field. After averaging it gives
(56) 
To go further, we assume that there is no large scale horizontal flow, like meridional circulation, so that there is no mean horizontal momentum flux. Therefore, Eq. (56) reduces to
(57) 
If we now introduce the pseudoLagrangian derivative as
(58) 
we finally get
(59) 
This equation is the same as used by Nordlund & Stein (2001).
a.2.2 Momentum conservation
(60) 
with the gas pressure and the gravitational acceleration. After averaging, one gets
(61) 
Using the same approximation as for Eq. (57), it reduces to
(62) 
Because we are interested in radial modes, we consider to obtain
(63) 
In the pseudoLagrangian frame this simplifies to
(64) 
with
(65) 
From Eq. (64), the stationary solution in the pseudoLagragian frame reads
(66) 
where denotes the temporal average. Thus, one can use Eqs. 64 and 66 to finally get
(67) 
where .
a.2.3 Energy equation
To go further, and disentangle the various contributions of the perturbation of gas pressure, it is necessary to write the equation governing the perturbation of entropy. To this end,
(68) 
where is the specific internal energy, is the radiative flux, and stands for viscous dissipation. After averaging, Eq. (68) becomes
(69) 
Finally, after some rearrangements the pseudoLagrangian energy equation is
(70) 
with
(71) 
Now, to the lowest order, we use the relation
(72) 
where is the temperature and the specific entropy. Consequently, Eq. (70) becomes to the lowest order
(73) 
where we have defined
(74) 
Finally, we note that while contributions by viscosity can be neglected in Eq. (60) and thus also in Eq. (67), the same does not hold true for Eq. (68) and hence Eq. (73)Eq. (74), as has been critically discussed by Canuto (1997b). Hence, in Canuto (1997a) the viscous flux was dropped in the dynamical equation for the large scale velocity field, whereas the dissipation rate of turbulent kinetic energy into heat, , was kept in the dynamical equations for the Reynolds stresses and for the temperature field, i.e., for the energy equation.