Issue 
A&A
Volume 625, May 2019



Article Number  A20  
Number of page(s)  15  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201834223  
Published online  01 May 2019 
Solar pmode damping rates: Insight from a 3D hydrodynamical simulation
^{1}
LESIA, Observatoire de Paris, CNRS, PSL Research University, Université Pierre et Marie Curie, Université Denis Diderot, 92195 Meudon, France
email: kevin.belkacem@obspm.fr
^{2}
Institute for Astrophysics, Faculty of Physics, Univ. Göttingen, FriedrichHundPlatz 1, 37077 Göttingen, Germany
^{3}
Department of Flow and Material Simulation, Fraunhofer Institute for Industrial Mathematics ITWM, FraunhoferPlatz 1, 67663 Kaiserslautern, Germany
Received:
11
September
2018
Accepted:
12
March
2019
Spaceborne 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 they can be used to grasp the main physics of the problem, such 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.
Key words: waves / convection / Sun: oscillations
© K. Belkacem et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
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, 2012; Dupret et al. 2009; 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 include 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 modelling. 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 subgiants 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 to properly understand 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, 5th 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 Runge–Kutta 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 Gauss–Radau 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, ∇ ⋅ F_{rad}, 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, Chapters 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 12 500 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 (1993a,b) 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 10^{−4} 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 5884 s or 18.58 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, that is, after about one solar hour.
The relaxed simulation was then evaluated every 1/20 of a sound crossing time, that is at an interval of 15.84 s. This way 2527 samples have been produced covering a time of about 40 012 s or slightly more than 11 h. 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 s and most of the time below 0.1 s). Statistical quantities have been computed for this output in a post processing step to calculate both horizontally averaged variables for each output step, that is for each of the 2527 samples, as well as ensembles averaged 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 T_{eff} is found to be on average 5750 K ± 18 K, with a negligibly small drift of −1.1 K per hour over that time and rare extrema not exceeding 52 K. The drift has been computed from a least square fit of a linear function to T_{eff} as determined at 0.1 Mm below the top of the simulation, where F_{tot} = F_{rad} within less than 0.1%. Visual inspection and a comparison with a straight mean, which is more than 1.6 K 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 T_{eff} 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 5779 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 log(g) = 4.4377 as does the chemical composition, (X, Y, Z) = (0.7373, 0.2427, 0.0200), with its Grevesse & Noels (1993) metallicity distribution.
Fig. 1. Superadiabatic gradient (∇ − ∇_{ad}) and Mach number (defined as the ratio between the temporally and horizontally averaged convective vertical velocity and sound speed) versus the simulation depth. The zeropoint depth is chosen where the temporally and horizontally averaged temperature equals the effective temperature. 
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
where A is the amplitude at t = 0. The observed signal is a sum of many such terms, each with their own amplitude and zeropoint of time, and also phase. ξ_{r} is the radial displacement eigenfunction, ω_{0} = 2πν_{0} is the pulsational eigenfrequency, t is the time, η is the damping rate, and ϕ the phase.
Fig. 2. Power spectrum of the vertical velocity at the photosphere, normalized by its maximum value, as a function of the frequency. 
In the power spectrum, for ν ≈ ν_{0}, the Fourier transform of Eq. (1) can be approximated by a Lorentzian function such as (e.g. Baudin et al. 2005; Appourchaux 2014)
and H stands for the mode height, Γ is the mode linewidth that is related to the mode damping rate through Γ = η/π. The mode height is subsequently related to the mode amplitude and mode linewidth by (e.g. Samadi 2011)
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), that is frequency and mode linewidth. Second, through the mode height and phase which depend on the location in the simulation.
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, that is 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}. We note that the linewidth of mode 1 is not provided since it is lower than the actual resolution (which is about 25 μ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).
Fig. 3. Top panel: mode velocity power density as a function of the depth for the three modes identified in Table 1. As for Fig. 1, the zeropoint depth is chosen where the temporally and horizontally averaged temperature equals the effective temperature. Bottom panel: same as for top panel except that the mode phases, computed directly in the Fourier space, are displayed. The zeropoint phases are chosen at a depth equal to 3 Mm. 
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. 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 α = 0.67, 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 Γ_{1}” 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 n = 16, n = 24, and n = 33. 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.
Fig. 4. Top panel: mode velocity as a function of the radius normalized by the total radius of the Sun obtained using the 1D solar model (i.e., where the temperature equals the effective temperature). For the normal modes of the 3D simulation, the velocity profiles have been obtained as described in Sect. 3.1 and using Eq. (3). For the normal modes computed using the 1D solar model, the velocities are obtained as described in Sect. 3.3 and they are normalized so that their kinetic energy equals the kinetic energy of the corresponding mode in the simulation. Bottom panel: mode damping rates observed by the GOLF instrument as a function of frequency. The observations are taken from Baudin et al. (2005). The overplotted red vertical dashed lines correspond to the frequencies of the modes in the 3D simulation (see Table 1) and are used to emphasize the correspondence with the observed solar 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 (ν_{max} = 3100 μ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 Δν_{sun} ≃ 135 μHz, we found Δν_{3D} ≃ 1930 μ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)
where v_{osc, s} is the mode velocity observed at the surface of the Sun and ℳ the mode mass defined as
where v_{osc} 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 11%, 8%, and 7% 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
with
and
where R is the solar radius, S is the horizontal surface of the simulation, ρ_{0} is the horizontally and timeaveraged density, z the depth of the 3D simulation, and the subscripts 3D and Sun 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 per cent.
Now using Eq. (4) together with Eq. (6), one obtains
We also note that the mode energy can be written as (see Samadi et al. 2015, for details)
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. E_{osc, 3D} ≃ E_{osc, Sun}). 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 v_{osc, s, Sun} from Baudin et al. (2005). It gives, using Eq. (3) and Eq. (4), E_{osc, 3D} ≃ 2.5 × 10^{19}J and E_{osc, Sun} ≃ 1.9 × 10^{20}J (for ℓ = 1), which is in a reasonable agreement given the uncertainties of our fit^{2} (see Sect. 3.1).
Finally, using the aforementioned argument and using Eq. (9), we get
Solar mode amplitudes are typically between 0.1 m s^{−1} and 0.3 m s^{−1} 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 100 m s^{−1}. 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 favours 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 X, can be decomposed such that
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 X, it is defined as
The quantity X can thus be decomposed such that
It immediately follows
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
where ρ is the density, u_{z} velocity. We also introduced the pseudoLagrangian derivative such that
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)
where P_{t} is the turbulent pressure, P_{g} is the gas pressure, g is the gravitational field that is considered constant as in the 3D simulation. In addition, for any quantity X, we have defined with the time average of . δX is therefore the pseudoLagrangian perturbation of X. 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
where δP = δP_{g} + δP_{t} 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}. We 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
To go further, we integrate over the mean mass column density τ_{0} ≡ 〈τ〉_{t} (defined as dτ_{0} = ρ_{0}dz), perform integration by parts and multiply by −iω to finally obtain
where
We note that at the modal frequencies, E 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
where again , so that Eq. (21) leads to
where
We assumed that near the modal frequency (i.e. ω ≃ ω_{0}) 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 σ_{0} = σ_{R} + iη, with η the damping rate, and in Eq. (24) we assumed since we are in the situation for which σ_{R} ≫ η.
At this step, several comments are necessary. The second term of Eq. (24), that is 𝒯_{2}, 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 (𝒯_{3}). Consequently, Eq. (24) reduces to
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, that is
where
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 (δX) are computed by interpolating the physical variables (X) onto the timeaveraged mean column density (τ_{0}) before subtracting their timeaveraged values. Then, the amplitudes and phases of each quantity (δX) 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 η_{1} = 5.13 μHz, which is consistent with our upper limit (i.e. the resolution of 25 μ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 η_{2} = 59.75 μ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 becomes 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.
Global characteristics of the three first normal modes.
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;
Fig. 5. Cumulative damping computed using Eq. (28) (normalized to the total mode damping rate) starting from the bottom of the simulation, as a function of the logarithm of temperature. Note that due to boundary condition effects, the uppermost layers (log T < 3.64) must be considered with care (see text for details). 
– The inner quasiadiabatic region, for which log T_{0} ≳ 4.1, 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 (log T_{0} ≃ 3.65) 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 (3.9 ≲ log T_{0} ≲ 4.0), 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
where δP_{g} is the perturbation of gas pressure and δP_{t} the perturbation of turbulent pressure. Therefore, Eq. (28) can be rewritten such that
where
For these individual contributions to η, shown in Fig. 6, we stress again that for log T ≳ 4.1 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 log T ≲ 4.1).
Fig. 6. Cumulative work integrals contributions (see Eq. (32)) associated with the gas pressure (top panel) and the turbulent pressure (bottom panel), integrated from the bottom of the simulation, as a function of the logarithm of temperature. The total mode damping of each mode is used to normalize the work integrals. We note that due to boundary condition effects, the uppermost layers (log T < 3.64) must be considered with care (see text for details). 
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 behaviour 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 log T ≳ 3.9.
Secondly, in the superadiabatic layers (and more precisely for 3.9 ≲ log T ≲ 4.0), 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. log T ≳ 4.0 for mode 1 and log T ≳ 4.05 for mode 2), excitation by δP_{g} competes against damping by δP_{t}, 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).
Figure 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.
Fig. 7. Cumulative work integrals contributions (see Eq. (32)) associated with mode 1 (top panel) and mode 2 (bottom panel), starting from the bottom of the simulation, as a function of the logarithm of temperature. The total mode damping of each mode is used to normalize the work integrals. We note that due to boundary condition effects, the uppermost layers (log T < 3.64) must be considered with care (see text for details). 
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
with
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
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)
where s is the specific entropy, F_{rad, z} is the vertical component of the radiative flux, F_{conv, z} is the vertical component of the convective flux defined by
where e is the specific internal energy and the perturbation of dissipation rate of turbulent kinetic energy into heat (ϵ) is defined by
with Q_{diss} 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 the dissipation rate of turbulent kinetic energy into heat.
Using Eq. (37) together with Eq. (38), finally gives
where η_{turb} is given by Eq. (34) and
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, that is,
We use this relation in the following, although we will focus on the radiative and convective flux contributions to the damping rates, i.e. Eqs. (42) and (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 behavioural 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 Fig. 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.
Fig. 8. Cumulative work integrals contributions of the gas pressure associated with mode 1 (top panel) and mode 2 (bottom panel), starting from the bottom of the simulation, as a function of the logarithm of temperature. The total mode damping of each mode is used to normalize the work integrals. We note that due to boundary condition effects, the uppermost layers (log T < 3.64) must be considered with care (see text for details). 
6. Concluding remarks
Using a 3D hydrodynamical simulation of the Sun computed with the ANTARES code and with a timeduration of about 11 h, 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, that is 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 log T > 4.1 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 modelling 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.
Internal errors are computed from the Hessian matrix as described in Press et al. (2002).
We 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).
Acknowledgments
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, M. J., Houdek, G., Handberg, R., & ChristensenDalsgaard, J. 2018, MNRAS, 478, 69 [CrossRef] [Google Scholar]
 Ando, H., & Osaki, Y. 1975, PASJ, 27, 581 [NASA ADS] [Google Scholar]
 Appourchaux, T. 2014, in A Crash Course on Data Analysis in Asteroseismology, eds. P. L. Pallé, & C. Esteban, 123 [Google Scholar]
 Appourchaux, T., Benomar, O., Gruberbauer, M., et al. 2012, A&A, 537, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Appourchaux, T., Gizon, L., & RabelloSoares, M.C. 1998, A&AS, 132, 107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baglin, A., Auvergne, M., & Barge, P. 2006a, in ESA Spec. Publ., eds. M. Fridlund, A. Baglin, J. Lochard, & L. Conroy, 1306, 33 [Google Scholar]
 Baglin, A., Auvergne, M., & Boisnard, L. 2006b, in 36th COSPAR Scientific Assembly, COSPAR Meeting, 36, 3749 [Google Scholar]
 Balmforth, N. J. 1992, MNRAS, 255, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Baudin, F., Samadi, R., Goupil, M.J., et al. 2005, A&A, 433, 349 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baudin, F., Barban, C., Belkacem, K., et al. 2011a, A&A, 529, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baudin, F., Barban, C., Belkacem, K., et al. 2011b, A&A, 535, C1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., & Samadi, R. 2013, in Lecture Notes in Physics, eds. M. Goupil, K. Belkacem, C. Neiner, F. Ligniéres, & J. J. Esteban, 865, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., & Kupka, F. 2006a, A&A, 460, 173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., Kupka, F., & Baudin, F. 2006b, A&A, 460, 183 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., et al. 2010, A&A, 522, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Goupil, M. J., Dupret, M. A., et al. 2011, A&A, 530, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Dupret, M. A., Baudin, F., et al. 2012, A&A, 540, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Canuto, V. M. 1992, ApJ, 392, 218 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M. 1997a, ApJ, 482, 827 [Google Scholar]
 Canuto, V. M. 1997b, ApJ, 489, L71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Canuto, V. M., Goldman, I., & Mazzitelli, I. 1996, ApJ, 473, 550 [NASA ADS] [CrossRef] [Google Scholar]
 Chaplin, W. J., & Miglio, A. 2013, ARA&A, 51, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Chaplin, W. J., Houdek, G., Karoff, C., Elsworth, Y., & New, R. 2009, A&A, 500, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ChristensenDalsgaard, J. 2009, in The Ages of Stars, eds. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse, IAU Symp., 258, 431 [NASA ADS] [Google Scholar]
 ChristensenDalsgaard, J. 2011, Astrophysics Source Code Library [record ascl:1109.002] [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [CrossRef] [Google Scholar]
 Cox, J. P., & Giuli, R. T. 1968, in Principles of Stellar Structure, eds. J. P. Cox, & R. T. Giuli (New York: Gordon and Breach) [Google Scholar]
 Dupret, M., Belkacem, K., Samadi, R., et al. 2009, A&A, 506, 57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dupret, M. A., Barban, C., & Goupil, M. J. 2006, in Proceedings of SOHO 18/GONG 2006/HELAS I, Beyond the Spherical Sun, ESA Spec. Publ., 624, 78 [Google Scholar]
 Favre, A. 1969 (Philadelphia: SIAM), 231 [Google Scholar]
 Fehlberg, E. 1970, Computing, 6, 61 [CrossRef] [Google Scholar]
 Goldreich, P., & Kumar, P. 1991, ApJ, 374, 366 [NASA ADS] [CrossRef] [Google Scholar]
 Gough, D. 1980, in Nonradial and Nonlinear Stellar Pulsation, eds. H. A. Hill, C. Barban, & W. A. Dziembowski (Berlin: Springer Verlag), Lecture Notes in Physics, 125, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, eds. N. Prantzos, E. VangioniFlam, & M. Casse, 15 [Google Scholar]
 Grigahcène, A., Dupret, M.A., Gabriel, M., Garrido, R., & Scuflaire, R. 2005, A&A, 434, 1055 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 GrimmStrele, H., Kupka, F., LöwBaselli, B., et al. 2015a, New A, 34, 278 [NASA ADS] [CrossRef] [Google Scholar]
 GrimmStrele, H., Kupka, F., & Muthsam, H. J. 2015b, Comp. Phys. Comm., 188, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Grosjean, M., Dupret, M. A., Belkacem, K., et al. 2014, A&A, 572, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Happenhofer, N., GrimmStrele, H., Kupka, F., LöwBaselli, B., & Muthsam, H. 2013, J. Comput. Phys., 236, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Hekker, S., & ChristensenDalsgaard, J. 2017, A&ARv, 25, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heun, K. 1900, Z. Math. Phys., 45, 23 [Google Scholar]
 Houdek, G. 2012, ASP Conf. Proc., 462, 7 [NASA ADS] [Google Scholar]
 Houdek, G. 2017, Eur. Phys. J. Web Conf., 160, 02003 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Houdek, G., & Dupret, M.A. 2015, Liv. Rev. Sol. Phys., 12, 8 [CrossRef] [Google Scholar]
 Houdek, G., Balmforth, N. J., ChristensenDalsgaard, J., & Gough, D. O. 1999, A&A, 351, 582 [NASA ADS] [Google Scholar]
 Houdek, G., Trampedach, R., Aarslev, M. J., & ChristensenDalsgaard, J. 2017, MNRAS, 464, L124 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, G.S., & Shu, C.W. 1996, J. Comput. Phys., 126, 202 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kjeldsen, H., Bedding, T. R., Arentoft, T., et al. 2008, ApJ, 682, 1370 [NASA ADS] [CrossRef] [Google Scholar]
 Kraaijevanger, J. 1991, BIT, 31, 482 [CrossRef] [Google Scholar]
 Kupka, F., Happenhofer, N., Higueras, I., & Koch, O. 2012, J. Comput. Phys., 231, 3561 [NASA ADS] [CrossRef] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2007a, in Convection in Astrophysics, eds. F. Kupka, I. Roxburgh, & K. L. Chan, IAU Symp., 239, 80 [NASA ADS] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2007b, in Convection in Astrophysics, eds. F. Kupka, I. Roxburgh, & K. L. Chan, IAU Symp., 239, 83 [NASA ADS] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2017, Liv. Rev. Comput. Astrophys., 3, 1 [Google Scholar]
 Kupka, F., & Robinson, F. J. 2007, MNRAS, 374, 305 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. 1993a, in Opacities for Stellar Atmospheres: [+0.0], [+0.5], [+1.0] (Cambridge, MA: Smithsonian Astrophysical Observatory), Kurucz CDROM, 2 [Google Scholar]
 Kurucz, R. 1993b, in ATLAS9 Stellar Atmosphere Programs and 2 km/s grid (Cambridge, MA: Smithsonian Astrophysical Observatory), Kurucz CDROM, 13 [Google Scholar]
 Lebreton, Y., Monteiro, M. J. P. F. G., Montalbán, J., et al. 2008, ApS&S, 316, 1 [Google Scholar]
 Ledoux, P. 1958, Handbuch der Physik, 51, 605 [NASA ADS] [CrossRef] [Google Scholar]
 Ledoux, P., & Walraven, T. 1958, Handbuch der Physik, 51, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (Oxford University Press) [Google Scholar]
 Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mosser, B., Benomar, O., Belkacem, K., et al. 2014, A&A, 572, L5 [NASA ADS] [CrossRef] [EDP Sciences] [PubMed] [Google Scholar]
 Mundprecht, E., Muthsam, H. J., & Kupka, F. 2013, MNRAS, 435, 3191 [NASA ADS] [CrossRef] [Google Scholar]
 Muthsam, H. J., Kupka, F., LöwBaselli, B., et al. 2010, New Astron., 15, 460 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, Å., & Stein, R. F. 2001, ApJ, 546, 576 [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical Recipes in C++: the Art of Scientific Computing [Google Scholar]
 Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenthal, C. S., ChristensenDalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689 [NASA ADS] [Google Scholar]
 Samadi, R. 2011, in Lecture Notes in Physics, eds. J.P. Rozelot, & C. Neiner (Berlin: Springer Verlag), 832, 305 [NASA ADS] [CrossRef] [Google Scholar]
 Samadi, R., Nordlund, Å., Stein, R. F., Goupil, M. J., & Roxburgh, I. 2003, A&A, 404, 1129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Georgobiani, D., Trampedach, R., et al. 2007, A&A, 463, 297 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Belkacem, K., Goupil, M. J., Dupret, M.A., & Kupka, F. 2008, A&A, 489, 291 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Ludwig, H.G., Belkacem, K., Goupil, M. J., & Dupret, M.A. 2010, A&A, 509, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Belkacem, K., Dupret, M.A., et al. 2012, A&A, 543, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Belkacem, K., & Sonoi, T. 2015, EAS Publ. Ser., 73, 111 [CrossRef] [Google Scholar]
 Shu, C.W. 2003, Int. J. Comput. Fluid Dyn., 17, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, C.W., & Osher, S. 1988, J. Comput. Phys., 77, 439 [Google Scholar]
 Sonoi, T., Samadi, R., Belkacem, K., et al. 2015, A&A, 583, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sonoi, T., Belkacem, K., Dupret, M.A., et al. 2017, A&A, 600, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stein, R. F., & Nordlund, A. 1998, ApJ, 499, 914 [Google Scholar]
 Stein, R. F., & Nordlund, Å. 2001, ApJ, 546, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Tassoul, M. 1980, ApJS, 43, 469 [NASA ADS] [CrossRef] [Google Scholar]
 Toutain, T., & Appourchaux, T. 1994, A&A, 289, 649 [NASA ADS] [Google Scholar]
 Trampedach, R. 1997, PhD Thesis, Aarhus University [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., Shibahashi, H. 1989, Nonradial Oscillations of Stars, 2nd edn. (University of Tokyo Press) [Google Scholar]
 Vrard, M., Kallinger, T., Mosser, B., et al. 2018, A&A, 616, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Xiong, D. R., Cheng, Q. L., & Deng, L. 2000, MNRAS, 319, 1079 [CrossRef] [Google Scholar]
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 X, is decomposed such that
so that we obviously have
The second average is called the mass or Favre average, so that for a quantity X, it is defined as
with
It immediately follows
From the above equations, it is quite straightforward to derive the following relations
which will be used in the following.
A.2. Mean equations
The Reynolds averages will be applied to the density and pressure while mass averages will be applied to the velocity, temperature and entropy fields.
A.2.1. Mass conservation
where ρ stands for the density and u_{j} is the j component of the velocity field. After averaging it gives
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. (A.11) reduces to
If we now introduce the pseudoLagrangian derivative as
we finally get
This equation is the same as used by Nordlund & Stein (2001).
A.2.2. Momentum conservation
with P_{g} the gas pressure and g the gravitational acceleration. After averaging, one gets
Using the same approximation as for Eq. (A.12), it reduces to
Because we are interested in radial modes, we consider j = z to obtain
In the pseudoLagrangian frame this simplifies to
with
From Eq. (A.19), the stationary solution in the pseudoLagragian frame reads
where 〈 〉_{t} denotes the temporal average. Thus, one can use Eqs. (A.19) and (A.21) to finally get
where .
A.2.3. Energy equation
To go further, and disentangle the various contributions of the perturbation of gas pressure, it is necessary to derive the equation governing the perturbation of entropy. To this end,
where e is the specific internal energy, F_{rad} is the radiative flux, and Q_{diss} stands for viscous dissipation. After averaging, Eq. (A.23) becomes
Finally, after some rearrangements the pseudoLagrangian energy equation is
with
Now, to the lowest order, we use the relation
where T is the temperature and s the specific entropy. Consequently, Eq. (A.25) becomes to the lowest order
where we have defined
Finally, we note that while contributions by viscosity can be neglected in Eq. (A.15) and thus also in Eq. (A.22), the same does not hold true for Eq. (A.23) and hence Eqs. (A.28)–(A.29), 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, that is, for the energy equation.
All Tables
All Figures
Fig. 1. Superadiabatic gradient (∇ − ∇_{ad}) and Mach number (defined as the ratio between the temporally and horizontally averaged convective vertical velocity and sound speed) versus the simulation depth. The zeropoint depth is chosen where the temporally and horizontally averaged temperature equals the effective temperature. 

In the text 
Fig. 2. Power spectrum of the vertical velocity at the photosphere, normalized by its maximum value, as a function of the frequency. 

In the text 
Fig. 3. Top panel: mode velocity power density as a function of the depth for the three modes identified in Table 1. As for Fig. 1, the zeropoint depth is chosen where the temporally and horizontally averaged temperature equals the effective temperature. Bottom panel: same as for top panel except that the mode phases, computed directly in the Fourier space, are displayed. The zeropoint phases are chosen at a depth equal to 3 Mm. 

In the text 
Fig. 4. Top panel: mode velocity as a function of the radius normalized by the total radius of the Sun obtained using the 1D solar model (i.e., where the temperature equals the effective temperature). For the normal modes of the 3D simulation, the velocity profiles have been obtained as described in Sect. 3.1 and using Eq. (3). For the normal modes computed using the 1D solar model, the velocities are obtained as described in Sect. 3.3 and they are normalized so that their kinetic energy equals the kinetic energy of the corresponding mode in the simulation. Bottom panel: mode damping rates observed by the GOLF instrument as a function of frequency. The observations are taken from Baudin et al. (2005). The overplotted red vertical dashed lines correspond to the frequencies of the modes in the 3D simulation (see Table 1) and are used to emphasize the correspondence with the observed solar modes. 

In the text 
Fig. 5. Cumulative damping computed using Eq. (28) (normalized to the total mode damping rate) starting from the bottom of the simulation, as a function of the logarithm of temperature. Note that due to boundary condition effects, the uppermost layers (log T < 3.64) must be considered with care (see text for details). 

In the text 
Fig. 6. Cumulative work integrals contributions (see Eq. (32)) associated with the gas pressure (top panel) and the turbulent pressure (bottom panel), integrated from the bottom of the simulation, as a function of the logarithm of temperature. The total mode damping of each mode is used to normalize the work integrals. We note that due to boundary condition effects, the uppermost layers (log T < 3.64) must be considered with care (see text for details). 

In the text 
Fig. 7. Cumulative work integrals contributions (see Eq. (32)) associated with mode 1 (top panel) and mode 2 (bottom panel), starting from the bottom of the simulation, as a function of the logarithm of temperature. The total mode damping of each mode is used to normalize the work integrals. We note that due to boundary condition effects, the uppermost layers (log T < 3.64) must be considered with care (see text for details). 

In the text 
Fig. 8. Cumulative work integrals contributions of the gas pressure associated with mode 1 (top panel) and mode 2 (bottom panel), starting from the bottom of the simulation, as a function of the logarithm of temperature. The total mode damping of each mode is used to normalize the work integrals. We note that due to boundary condition effects, the uppermost layers (log T < 3.64) must be considered with care (see text for details). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.