© The Authors 2024

1. Introduction

Comparisons between the first asteroseismic models and solar surveys revealed a deviation between theoretical and observational eigenfrequencies, which affected modes of various degrees. These deviations were attributed to the presence of a highly turbulent region in the superficial layer of the Sun. The interaction or coupling between the stellar oscillations and the turbulent convective motion has several consequences. A major consequence are the discrepancies between theory and observations regarding the eigenfrequencies of stellar oscillations, referred to as “surface effects,” which are directly linked to the presence of this superficial convective zone (Dziembowski et al. 1988; Christensen-Dalsgaard et al. 1996; Rosenthal et al. 1999). They are not solely related to the Sun, but affect cold stars with a convective region in their superficial layer (e.g., low-mass main sequence stars and red giant stars). These discrepancies constitute an important limitation in accurate seismic probing of those stars.

Multiple attempts have been made to correct the surfaces effects. On one hand, the structural impact has been modeled using patched models followed by empirical laws to account for the surface effects (Kjeldsen et al. 2008; Ball & Gizon 2014; Sonoi et al. 2015; Trampedach et al. 2017). After computing the adiabatic frequencies of the patched model, free parameters in the empirical laws had to be calibrated or optimized to fit the observations.

On the other hand, the modal impact is processed via the time-dependent convection models based on perturbation of the mixing-length approach for turbulence (Unno 1967; Gough 1977). These latter models (e.g., Gabriel et al. 1975; Houdek 1996; Grigahcène et al. 2005; Sonoi et al. 2017; Houdek et al. 2017, 2019) suffer from the underlying hypotheses behind the mixing-length theory (MLT) and their one-dimensional approach of the problematic (along the radial direction of the star). In the MLT, the turbulent cascade is reduced to a single scale, the mixing length, leading to misrepresentation of turbulent characteristics such as the turbulent pressure. Furthermore, the actual coupling occurs over the entire turbulent spectrum, which possesses features that vary from one region to another (downflows or upflows). Due to horizontal averages, those MLT models are dependent on a closure equation for the Reynolds stress and rely on free parameters – the calibration of which is not always physical – to account for all the features erased by the averaging process. More complex analytical approaches beyond the mixing length have been proposed; these include the use of Reynolds stress models (Xiong et al. 2000), but these suffer from similar drawbacks. Another family of models is based on the use of stellar atmosphere simulation for a direct measure of the mode line widths and damping rates (Nordlund & Stein 2001; Stein & Nordlund 2001; Belkacem et al. 2019; Zhou et al. 2020). As a consequence of the boundary conditions used in the hydrodynamic simulations, normal modes appear to be present in the simulation box in addition to the turbulent acoustic noise. The properties of the normal modes are then examined in order to study the mode physics. However, there are several downsides to studying modes present in a simulation box. Distinguishing those normal modes from convective fluctuations in the simulation box requires a very long simulation of several tens of hours or days and the normal modes identified actually correspond to modes of the box and not to global modes of the star. Therefore, due to the relative difference between the size of the simulation box and the star, the inertia of those modes is much smaller than global stellar modes, which leads to the obtention of modes of a much larger amplitude than those observed. The physics behind those modes is modified compared to reality; for example, nonlinearities resulting from their excessive amplitude could be responsible for an increased importance of certain physical mechanisms, such as the role of dissipation of turbulent kinetic energy in this frame.

In solar-like stars, the mechanism behind their oscillations is the stochastic excitation and damping induced by the turbulence around the photosphere of the star. The coupling between convection and oscillation impacts the energetics of the oscillations, not only in the stars affected by surface effects but also in others such as white dwarfs (Van Grootel et al. 2012), γ Dor, and δ Scuti stars (Dupret et al. 2005). Therefore, the damping of the modes is also directly related to turbulence. Several attempts have been made to model the driving (e.g., Goldreich & Keeley 1977a,b; Samadi & Goupil 2001; Chaplin et al. 2005; Belkacem et al. 2006, 2008, 2010) and damping (e.g., Goldreich & Kumar 1991; Grigahcène et al. 2005; Dupret et al. 2006; Belkacem et al. 2012) of the modes but these suffer the same limitations as the attempts to model surface effects. Moreover, damping models are often computed and adjusted via different methods from those used to deal with surface effects models, which constitutes an anomaly in itself, as both are related in the same manner to the convective motion in the superficial layer of the star; more details regarding the modal effect of turbulence are found in reviews by Houdek & Dupret (2015) and Samadi et al. (2015).

Looking solely at the structure of turbulence is an extremely challenging topic. This is especially stressed in the stellar astrophysical context, with turbulent flows with Reynolds numbers of around 10¹⁰ and a very large spectrum of space- and timescales. For those reasons, the first numerical simulations of turbulent flows relied on oversimplification, with models such as the Reynolds-averaged Navier-Stokes (RANS), in which the well-known Prandtl mixing length can be used as a closure equation. However, in recent decades, tremendous progress has been made in computational fluid dynamics with the rise of the Large Eddy Simulation (LES). These simulations are able to solve the Navier-Stokes equations down to a spatial limit, the smallest resolved length scale. These improvements in turbulence modeling also benefit stellar flows (Schmidt 2015; Kupka & Muthsam 2017). For instance, a state-of-the-art LES for stellar surface convection, called “COnservative COde for the COmputation of COmpressible COnvection in a BOx of L Dimensions with l = 2,3” (CO⁵BOLD), can be used to study convection in solar-type stars, red supergiants, substellar objects, and white dwarfs (Freytag et al. 2012).

As mentioned above, the interaction between convection and stellar oscillations, the modal effect of turbulence, is still an open question in asteroseismology. Despite improvements brought to stellar oscillations models, the current approaches have reached limits regarding questions as to the impact of convection, whether that be in terms of surface effects or damping of the mode. Even with the most recent models based on hydrodynamic simulations, the corrections of the surface effects are unsatisfactory and the damping remains poorly understood. In the current time-dependent convection models, the excitation of the modes, for instance in the γ Dor (Dupret et al. 2005) and white dwarfs (Van Grootel et al. 2012), is not in agreement with the observations and therefore requires further investigation. These points still constitute the main limitations to accurate seismic probing of solar-like and red giants stars. Alongside the progress made in the structure of turbulence, this motivated us to change the current paradigm by moving to a full 3D quasi-radial global stellar oscillation model. In our approach, we propose to remain as close as possible to the complex 3D nature of turbulence which cannot be accurately represented in one dimension.

In this paper, we define a new formalism for solving the quasi-radial global nonadiabatic oscillation equations in a full 3D framework. In Sect. 2, we define the specific simulation used to develop our formalism and establish the relevance of a three-dimensional approach. The derivation of the main equations characterizing the stellar oscillation in a 3D environment is presented in Sect. 3. With our formalism, we introduce a spectral decomposition that allows us to take into account the full spectrum of turbulence and to assess the impact of all the different turbulent scales on the stellar oscillations. Based on this Fourier decomposition, in Sect. 4, we separate the impact of the average flow and its fluctuations. This allows us to derive wave equations describing the average oscillations as well as corrections – called fluctuation oscillations – due to the turbulent motion; consequently, the eigenvalue problem – to be solved numerically – is presented in Sect. 5. We also discuss the extension of our formalism beyond the simulation box. Due to the peculiar nature of our formalism, a specific derivation is required for the work integral and the damping of the mode, the development of which is presented in Sect. 6. We also discuss possibilities offered by a 3D stellar oscillation model to precisely locate regions of driving and damping of the mode and to characterize their physical origin. Finally, we draw our conclusions regarding our formalism in Sect. 7. We note that this first paper of the series is dedicated to the mathematical developments of our formalism. The numerical implementation of this formalism and a presentation of the results will follow in a subsequent paper of this series.

2. Three-dimensional hydrodynamic simulation

To develop a 3D formalism depicting the interaction between convection and oscillations in low-mass stars, a 3D structure model has to be defined. We set the 3D structure of the convective superficial layer of the star via the CO⁵BOLD code. Given the high accuracy of the observational data alongside the well-known driving mechanism behind its oscillations, a star of great interest for our formalism is the Sun. Hence, a custom simulation of the Sun was produced and serves as reference for the derivation of our formalism in Sect. 3. This simulation was run in a 3D cartesian box where the z-coordinate is aligned with the radial direction of the Sun and the x − y plane is defined normal to this radial or z direction. With its dimensions (18.5 × 18.5 × 8.4 × 10⁸ cm) and location, the bottom of the box reaches deep enough into the convective zone, where turbulent fluctuations are almost nonexistent (see Fig. 4), while the top of the box is located above the photosphere in the atmosphere of the star. The heart of the convection zone is therefore fully represented.

Several hypotheses and specifications are also accounted for and adhered to in simulation runs. As the box is relatively small in comparison to the size of the Sun, a plan parallel approach is used and therefore the curvature of the envelope is neglected. Periodic conditions are implemented for the four sides of the box. An entropy flux is imposed at the bottom of the simulation (“inoutflow” condition) in order to obtain an effective temperature in the simulation corresponding to that of the Sun. Finally, a transmitting boundary condition is imposed at the top of the box, which implies an exponential decrease in density. Regarding the radiative transport, the optically thick approximation (also known as diffusion approximation) is used up to the photosphere, while the remaining part uses gray opacities. Every other parameter of the simulation is set according to the classical prescription for a solar hydrodynamic simulation. We note that, based on these specifications, two simulations were run, that is, at high and low resolution. The only distinction between those simulations is the number of nodes present in the box and therefore the spatial resolution of the simulation. A grid of 534 × 534 × 424 nodes is used in the high-resolution simulation, while the low-resolution simulation uses a grid of 189 × 189 × 150 nodes.

Using the high-resolution simulation, we can properly illustrate our statement regarding the difficulties faced by 1D stellar oscillation models to fully capture the nature of the turbulent flow and the usefulness of a 3D formalism. Figure 1 shows an internal-energy snapshot of a solar CO⁵BOLD simulation where we can clearly see the intricacy of the solar convection zone just below the photosphere. We also see clear asymmetries in the downflows (blue shades) and upflows (orange shades) in terms of effective area but also intensity. We also see the presence of a complex velocity field whose interaction with the stellar oscillations can only be fully assessed in a 3D framework. We set the foundation of this 3D formalism – which fully accounts for the 3D nature of turbulence – via the use of a Large Eddy Simulation. We therefore start by establishing a new set of equations that describe the stellar oscillations in this nonadiabatic context.

Fig. 1. Horizontal cut of the CO5BOLD 3D simulation box below the photosphere of the normalized temperature. Blue shades correspond to downflows while orange shades identify upflows. The snapshot shown here corresponds to a high-resolution simulation of the Sun performed with CO5BOLD (Grid nodes: 534 × 534 × 424).

3. Mathematical model

3.1. Starting point

As the formalism is developed in a 3D environment associated with a hydrodynamic simulation, the starting point of our model is the well-known Navier-Stokes equations; given the high Reynolds number present in astrophysical flows, viscosity is left out in this instance. Here we reiterate the continuity, momentum, and conservation of energy equations in a Cartesian reference frame:

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \underline{v})& =0,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial \underline{v}}{\partial t}+\underline{v}·\mathrm{\nabla }\underline{v}& =-\frac{\mathrm{\nabla }P}{\rho }+\underline{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \rho T\frac{\mathrm{D}s}{\mathrm{D}t}& =-\mathrm{\nabla }·\underline{F},\hfill \end{array}$$(3)

where ρ is the density, t the time, $\underline{v}$ the velocity vector, P the pressure, $\underline{g}$ the gravitational acceleration vector, T the temperature, s the entropy and finally $\underline{F}$ is the radiative heat flux vector.

In the conservation of energy Eq. (3), we note the presence of the radiative heat flux $\underline{F}$. Therefore, coupled with the Navier-Stokes equations (Eqs. (1)–(3)), we define an additional expression to express this physical quantity. We decided to use the classical diffusion equation for the radiative transport, which is valid where the mean free path of photons is negligible, as follows:

$$\begin{array}{c}\hfill \underline{F}=-\mathcal{X}\mathrm{\nabla }T,\end{array}$$(4)

where 𝒳 stands for the radiative conductivity.

We apply the perturbative approach to the 3D medium, here applying a small Eulerian perturbation directly to the time- and space-dependent Navier-Stokes equations (Eqs. (1)–(3)). We obtain the following system of coupled equations:

$$\begin{array}{cc}\hfill \frac{\partial {\rho }^{\prime }}{\partial t}+\mathrm{\nabla }·({\rho }^{\prime }\underline{v}+\rho {\underline{v}}^{\prime })& =0,\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill \frac{\partial {\underline{v}}^{\prime }}{\partial t}+{\underline{v}}^{\prime }·\mathrm{\nabla }\underline{v}+\underline{v}·\mathrm{\nabla }{\underline{v}}^{\prime }& =-\frac{\mathrm{\nabla }{P}^{\prime }}{\rho }-\mathrm{\nabla }P{\left(\frac{1}{\rho }\right)}^{\prime }+{\underline{g}}^{\prime },\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill (\frac{{\rho }^{\prime }}{\rho }+\frac{T^{\prime }}{T})\left(\rho T\frac{\mathrm{D}s}{\mathrm{D}t}\right)+\left(\rho T\right){\left(\frac{\mathrm{D}s}{\mathrm{D}t}\right)}^{\prime }& =-\mathrm{\nabla }·{\underline{F}}^{\prime },\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill {\underline{F}}^{\prime }& =-{\mathcal{X}}^{\prime }\mathrm{\nabla }T-\mathcal{X}\mathrm{\nabla }{T}^{\prime },\hfill \end{array}$$(8)

where X′ corresponds the Eulerian perturbation of the physical quantity X.

The physical quantities appearing in these equations (Eqs. (5)–(8)) depend on the three spatial coordinates and time. In particular, the coefficients are directly obtained from 3D hydrodynamic simulations. We aim to solve these for the entire star considered and in this 3D environment. It is worth insisting that a major change of paradigm with the current stellar oscillation models can already be observed. We do not restrict ourselves to the horizontal average, which allows us to avoid introducing challenging parameters or unknowns such as the convective flux in the conservation of energy equation (Eq. (7)) or the Reynolds stress in the conservation of momentum equation (Eq. (6)). This automatically means that the most difficult part of turbulence modeling is avoided, namely the closure equation.

A fundamental distinction can also be observed between our formalism and the various existing approaches based on hydrodynamic simulation (Nordlund & Stein 2001; Belkacem et al. 2019; Zhou et al. 2020). First, we aim to solve the system of equations (Eqs. (5)–(8)) for the entire star and we do not limit our approach to the simulation box. This leads to global oscillation modes and frequencies corresponding to the studied star. Second, in our formalism, as mentioned earlier, turbulence and convection are considered as known data. More specifically, the medium in which the waves propagate is an input provided by a previously computed 3D hydrodynamics simulation and we compute the 3D eigenfunctions corresponding to the oscillation modes (the primed quantities in Eqs. (5)–(8)). Our approach enables us to avoid significant difficulties and limitations associated with other existing 3D approaches. As we do not have to search the oscillation modes present in a 3D simulation, we completely avoid the difficult and problematic task of having to distinguish convection from oscillations in a precomputed 3D simulation. Moreover, our formalism allows us to compute global oscillation modes. Therefore, those global modes are not subject to major inaccuracies associated with box modes present in the simulation: excessively large amplitudes of oscillations due to their overly low inertia artificially enhances nonlinear effects and leads to inaccurate frequencies, and so on. Therefore, our formalism enables us to accurately quantify the 3D modal surface effects on frequencies, which is not possible with box modes. We can illustrate this difference between our formalism and the current 3D approaches based on hydrodynamic simulation modes. Without loss of generality, let us take a physical quantity X as follows:

$$\begin{array}{c}\hfill X(\underline{r},t)=X_{\mathrm{conv}}(\underline{r},t)+X^{\prime }(\underline{r},t).\end{array}$$(9)

In our model, we want to solve a small perturbation X′ (unknown stellar oscillation) around an equilibrium position, which is given by the hydrodynamic simulation X_conv (considered data, a coefficient of our system of equations). It immediately follows that convection and stellar oscillations are naturally separated in our formalism.

Nevertheless, a solution cannot be found numerically and further development is required to obtain a solvable system. As we are making use of the perturbative approach and wish to isolate an oscillation mode, we assume a small Eulerian perturbation; we can decompose it into a space-related component $a^{\prime }(\underline{r})$ ($\underline{r}\to r(x,\mathit{y},z)$) and a time-related component e^i σt, where σ is the oscillation frequency:

$$\begin{array}{c}\hfill A^{\prime }(\underline{r},t)=a^{\prime }(\underline{r}){\mathrm{e}}^{i\sigma t}.\end{array}$$(10)

While one might be tempted to use spherical coordinates, we ought to comply to the reference frame used in the 3D hydrodynamic simulation, which is 3D Cartesian (see Sect. 2). After the introduction of this decomposition into our system, we take the Fourier transform of this set of equations at the (unknown) oscillation frequency and over a shorter time interval than the total duration of the 3D simulation:

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }t}{\displaystyle {\int }_t}{A}^{\prime }(\underline{r},t)B(\underline{r},t){\mathrm{e}}^{-i\sigma t}\mathrm{d}t=\frac{a^{\prime }(\underline{r})}{\mathrm{\Delta }t}{\displaystyle {\int }_t}B(\underline{r},t){\mathrm{e}}^{i\sigma t}{\mathrm{e}}^{-i\sigma t}\mathrm{d}t=a^{\prime }(\underline{r})\overline{B}(\underline{r}).\end{array}$$(11)

Equation (11) shows an example of this time integration, where we highlight the case where a coefficient B, that is, any physical quantity of a 3D hydrodynamic simulation, is multiplied by a perturbation A′, any unknown to be solved in our model. In combination with our decomposition of the perturbation (Eq. (10)), this time integration turns the coefficient into a time-averaged coefficient $\overline{B},$ which then allows us, without loss of generality, to use the local time averages of the physical quantities available in our Large Eddy Simulation as coefficients in our stellar oscillation model. This step is essential in our new formalism as it allows us to remove the time dependency of all the variables, unknowns, and coefficients alike, as illustrated in Eq. (11). After some algebra, we obtain the following set of equations for the continuity, momentum, energy, and radiative transport equations. We note that in the remainder of the paper, the time average symbol “$\overline{A}$” is omitted for the sake of clarity and therefore any coefficients present in the equations are time averaged, unless clearly specified otherwise.

$$\begin{array}{cc}\hfill i\sigma {\rho }^{\prime }+\mathrm{\nabla }·({\rho }^{\prime }\overline{\underline{v}}+\overline{\rho }{\underline{v}}^{\prime })& =0,\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill i\sigma {\underline{v}}^{\prime }+{\underline{v}}^{\prime }·\overline{\mathrm{\nabla }\underline{v}}+\overline{\underline{v}}·\mathrm{\nabla }{\underline{v}}^{\prime }& =-\frac{\mathrm{\nabla }{P}^{\prime }}{\overline{\rho }}-\overline{\mathrm{\nabla }P}{\left(\frac{1}{\rho }\right)}^{\prime }+{\underline{g}}^{\prime },\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill (\frac{{\rho }^{\prime }}{\rho }+\frac{T^{\prime }}{T})\overline{\left(\rho T\frac{\mathrm{D}s}{\mathrm{D}t}\right)}+\overline{\left(\rho T\right)}{\left(\frac{\mathrm{D}s}{\mathrm{D}t}\right)}^{\prime }& =-\mathrm{\nabla }·{\underline{F}}^{\prime },\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill {\underline{F}}^{\prime }& =-{\mathcal{X}}^{\prime }\overline{\mathrm{\nabla }T}-\overline{\mathcal{X}}\mathrm{\nabla }{T}^{\prime }.\hfill \end{array}$$(15)

A new stellar oscillation model must be closer to reality and should avoid relying on free parameters by staying in a 3D environment and capturing the nature of turbulence. This could potentially be undermined if the length of the integration time chosen removes those turbulent characteristics. Indeed, as the length of the integration time is increased, we can lose some features of turbulence such as the amplitude of the motion as well as the smallest timescale of the simulation. Due to the nature of the solar flow in its superficial layer (Rayleigh-Bénard convection), with an ill-chosen time of integration, one might even be able to remove convective cells. Figure 2 shows a comparison between a snapshot and several time-averaged intervals performed over 3, 5 and 10 min with CO⁵BOLD for the case of the Sun. Comparing the normalized internal energy between a snapshot and a three-minute time average, we can clearly see that the differences in terms of turbulence pattern and intensity of the downflows and upflows are almost negligible. A quick look at the scales and the shape of the convective cells allows us to state that, over three minutes, the features of turbulence have been conserved. With a five-minute time average, the turbulent structure is still similar to the snapshot. While the upflows remain at the same amplitude, we note a small decrease in the maximum amplitude of the downflows. Nonetheless, the features of turbulence are also conserved. Theses similarities between the snapshot and time average strengthen the time processing step of our model, given that the turbulence does not drastically evolve over the typical oscillation period of the Sun. We note that taking a time average that is longer than the oscillation period and convective turn-over time does not lead to the complete removal of turbulence, as the largest scales are stable and preserved, meaning that the average medium in which the modes propagate does not tend to zero. However, with longer time averaging intervals, as shown by the ten-minutes average in Fig. 2, the turbulent features are reduced. Indeed, only the largest scales remain and the amplitudes of the downflows have been significantly diminished. Therefore, in our formalism, shorter time averaging intervals of approximately one oscillation period are recommended in order to avoid the loss of some features of turbulence.

Fig. 2. Comparison between a snapshot of the normalized internal energy (top left panel) and time averages over 3 min (top right), 5 min (bottom left), and 10 min (bottom right). Blue shades correspond to downflows while orange shades signal upflows. The snapshots shown here corresponds to low resolution simulations of the Sun performed with CO5BOLD (Grid nodes: 189 × 189 × 150).

3.2. Fourier approach

Solving the system in its current state (Eqs. (12)–(15)) is still out of reach for most numerical schemes and algorithms. Consequently, we use the Fourier series and transforms to derive the final state of our new formalism. We note that we work with a spatial approach in the Fourier domain. Without loss of generality, one can express the variables of our system of equations, unknown perturbations (e.g., ρ′, T′), and 3D simulation quantities and coefficients (e.g., ρ, 𝒳) alike as a 2D Fourier series in the horizontal plane $\underline{x}=(x,\mathit{y})$:

$$\begin{array}{cc}& A(\underline{r})={\displaystyle \sum _{m=-\frac{N_x}{2}}^{\frac{N_x}{2}-1}\sum _{n=-\frac{N_y}{2}}^{\frac{N_y}{2}-1}}{a}_{m,n}(z){\mathrm{e}}^{i\underline{k}·\underline{x}}={\displaystyle \sum _k}{a}_k(z){\mathrm{e}}^{i\underline{k}·\underline{x}},\hfill \end{array}$$(16)

$$\begin{array}{cc}& A^{\prime }(\underline{r})={\displaystyle \sum _k}{a}_k^{\prime }(z){\mathrm{e}}^{i\underline{k}·\underline{x}},\hfill \end{array}$$(17)

where m,n are summation indexes, N_x is the number of nodes in the x direction, N_y is the number of nodes in the y direction, $\underline{k}=\frac{2\pi m}{N_x\mathrm{\Delta }x}{\underline{e}}_x+\frac{2\pi n}{N_{\mathit{y}}\mathrm{\Delta }\mathit{y}}{\underline{e}}_{\mathit{y}}$, Δx,Δy are the mesh spatial resolutions in the horizontal plane, $A(\underline{r})$ stands for a 3D simulation quantity (coefficient), and $A^{\prime }(\underline{r})$ represents a perturbation (unknown).

Before further presentation of our formalism, it might be useful to briefly discuss the reasons for and relevance of using the spatial Fourier transforms as a tool to solve our system of coupled partial differential equations established in real space (Eqs. (12)–(15)). First, we insist on the fact that the Fourier transforms are only a tool used to resolve the oscillation equations in 3D. Once a solution is found, an inverse Fourier transform is used to return to real space, enabling us to verify that Eqs. (12)–(15) have been precisely solved and to interpret the results. Second, as detailed in Sect. 4, introducing the Fourier series (Eqs. (16) and (17)) allows us to obtain a natural distinction between the horizontal average value and the fluctuations around this horizontal average value. Indeed, looking at a Fourier series of a coefficient of our system of equations, $A(\underline{r})$ (Eq. (16)), we can immediately see that the term k = 0 in the sum corresponds to the horizontal average value of A at a given z-coordinate, while the remaining terms (k ≠ 0) represent fluctuations around this horizontal average value. This natural distinction between horizontal average and fluctuations is crucial in order to obtain the desired state of our system of equations and to solve them numerically. This distinction can just as easily be illustrated in real space via Eq. (9). Without loss of generality, X_conv and X′ alike can be decomposed into their horizontal average and fluctuations as follows:

$$\begin{array}{c}\hfill X(\underline{r},t)=X_{\mathrm{conv},\mathrm{avg}}(r)+X_{\mathrm{conv},\mathrm{fluc}}(\underline{r})+X_{\mathrm{avg}}^{\prime }(r){\mathrm{e}}^{i\sigma t}+X_{\mathrm{fluc}}^{\prime }(\underline{r}){\mathrm{e}}^{i\sigma t}.\end{array}$$(18)

We reexpressed X_conv as X_conv, avg and X_conv, fluc where X_conv, fluc is defined as X_conv − X_conv, avg. After a similar decomposition of the oscillations X′, we obtained Eq. (18). The parallel with the Fourier space is then straightforward. We recall here that this decomposition into horizontal average and fluctuations is performed on both 3D hydrodynamic quantities and perturbations. Whether horizontal average or fluctuation, we note that any quantity coming from the hydrodynamic simulation is considered known and is used as a coefficient in the equations.

Using Fourier series (Eqs. (16) and (17)) to express physical quantities from the simulation is possible due to the periodic boundary conditions used for the four sides of the CO⁵BOLD box. Regarding the perturbation, the periodicity condition in the horizontal plane is also naturally implied and a Fourier series expansion is then possible. Indeed, assuming a nonzero resulting horizontal perturbation, an aperiodic perturbation in the horizontal plane would result in a global stellar motion in one of the horizontal directions, which is incompatible with the observations. Applying this decomposition in 2D Fourier series (Eqs. (16) and (17)) to our system of equations (Eqs. (12)–(15)) leads to the following in terms of continuity, momentum, energy, and radiative transport:

$$\begin{array}{cc}& {\displaystyle \sum _{k_1}}\{i\sigma {\rho }_{k_1}^{\prime }{\mathrm{e}}^{i{\underline{k}}_1·\underline{x}}\}+{\displaystyle \sum _{k_1}}{\displaystyle \sum _{k_2}}\{[i({\underline{k}}_1+{\underline{k}}_2)]·{\rho }_{k_1}^{\prime }{\underline{v}}_{{h,k}_2}{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}+[i({\underline{k}}_1+{\underline{k}}_2)]·{\rho }_{k_1}{\underline{v}}_{{h,k}_2}^{\prime }{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}\}\hfill \\ \hfill & +{\displaystyle \sum _{k_1}}{\displaystyle \sum _{k_2}}\{\frac{\partial }{\partial z}\left({\rho }_{k_1}^{\prime }{v}_{{z,k}_2}\right){\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}+\frac{\partial }{\partial z}\left({\rho }_{k_1}{v}_{{z,k}_2}^{\prime }\right){\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}\}=0,\hfill \end{array}$$(19)

$$\begin{array}{cc}& {\displaystyle \sum _{k_1}}\{i\sigma {\underline{v}}_{k_1}^{\prime }{\mathrm{e}}^{i{\underline{k}}_1·\underline{x}}\}+{\displaystyle \sum _{k_1}}{\displaystyle \sum _{k_2}}\{{\underline{v}}_{{h,k}_1}^{\prime }·(i\underline{k_2})\otimes {\underline{v}}_{k_2}+{\underline{v}}_{{h,k}_1}·(i\underline{k_2})\otimes {\underline{v}}_{k_2}^{\prime }{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}\}\hfill \\ \hfill & +{\displaystyle \sum _{k_1}}{\displaystyle \sum _{k_2}}\{v_{z,k_1}\prime \frac{\partial }{\partial z}\left({\underline{v}}_{k_2}\right){\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}+v_{z,k_1}\frac{\partial }{\partial z}\left({\underline{v}}_{k_2}^{\prime }\right)\}=-{\displaystyle \sum _{k_1}}{\displaystyle \sum _{k_2}}\{i\underline{k_1}{P}_{k_1}^{\prime }{\left(\frac{1}{\rho }\right)}_{k_2}{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}+i\underline{k_1}{P}_{k_1}{\left(\frac{1}{\rho }\right)}_{k_2}^{\prime }{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}\hfill \\ \hfill & +{\underline{e}}_z[\frac{\partial P_{k_1}^{\prime }}{\partial z}{\left(\frac{1}{\rho }\right)}_{k_2}+\frac{\partial P_{k_1}}{\partial z}{\left(\frac{1}{\rho }\right)}_{k_2}^{\prime }]{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}\}+{\displaystyle \sum _{k_1}}{\underline{g}}_{k_1}{\mathrm{e}}^{i{\underline{k}}_1·\underline{x}},\hfill \end{array}$$(20)

$$\begin{array}{cc}& {\displaystyle \sum _{k_1}}{\displaystyle \sum _{k_2}}\{{(\frac{{\rho }^{\prime }}{\rho }+\frac{T^{\prime }}{T})}_{k_1}{\left(\rho T\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_{k_2}{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}+{\left(\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_{k_1}^{\prime }{\left(\rho T\right)}_{k_2}{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}\}=-{\displaystyle \sum _{k_1}}\{i{\underline{k}}_1·{\underline{F}}_{{{h,k}_1}^{\prime }}{\mathrm{e}}^{i{\underline{k}}_1·\underline{x}}+\frac{\partial }{\partial z}(F_{z,k_1}\prime ){\mathrm{e}}^{i{\underline{k}}_1·\underline{x}}\},\hfill \end{array}$$(21)

$$\begin{array}{cc}& {\displaystyle \sum _{k_1}}\{{\underline{F}}_{k_1}^{\prime }{\mathrm{e}}^{i{\underline{k}}_1·\underline{x}}\}=-{\displaystyle \sum _{k_1}}{\displaystyle \sum _{k_2}}\{{\mathcal{X}}_{k_1}^{\prime }\left(i\underline{k_2}\right)T_{k_2}{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}+{\mathcal{X}}_{k_2}\left(i\underline{k_2}\right)T_{k_2}^{\prime }{\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}\}\hfill \\ \hfill & -{\underline{e}}_z{\displaystyle \sum _{k_1}}{\displaystyle \sum _{k_2}}\{{\mathcal{X}}_{k_1}^{\prime }\frac{\partial }{\partial z}\left(T_{k_2}\right){\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}+{\mathcal{X}}_{k_2}\frac{\partial }{\partial z}\left(T_{k_1}^{\prime }\right){\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}\},\hfill \end{array}$$(22)

where ⊗ denotes the dyadic product, ${\underline{e}}_z$ is the unit vector in the z/radial direction, and k₁ and k₂ are arbitrary summation indexes.

In Eqs. (19)–(22), we note several major consequences from the introduction of our 2D Fourier decomposition (Eqs. (16) and (17)). First, any product turns into double summations. Second, any horizontal gradient or divergence is replaced by a product following the well-known properties of the Fourier series or transforms. Third and finally, the only remaining derivatives are with respect to the radial coordinate z, which is hereafter the only independent variable of our system of equations.

After introducing the 2D Fourier decomposition, double summations appeared. To remove this inconvenience, we can project our system of equations (Eqs. (19)–(22)) in a specific direction of the Fourier space (horizontal plane). Let us use ${\underline{k}}_0$ to denote an arbitrary direction in the 2D Fourier space. Projecting any quantity A, an equation in our case, is simply done in the Fourier space by applying the following integral, where Σ stands for the horizontal surface:

$$\begin{array}{c}\hfill A\to \frac{1}{\mathrm{\Sigma }}{\displaystyle \int \int A}{\mathrm{e}}^{-i{\underline{k}}_0·\underline{x}}\mathrm{d}\mathrm{\Sigma }.\end{array}$$(23)

Applying this projection (Eq. (23)) to our system of equations (Eqs. (19)–(22)) leads to specific relations between the three parameters ${\underline{k}}_0$ (direction of the projection), ${\underline{k}}_1,$ and ${\underline{k}}_2$ (both parameters of the decomposition in 2D Fourier series). Those relations can be established based on the two possibilities that can already be observed in the continuity equation (Eq. (19)). The whole purpose of this projection is to provide a link between our previously independent parameters ${\underline{k}}_1$ and ${\underline{k}}_2$ from the Fourier series decomposition.

The first relation can be derived based on the projection of the first term of the continuity equation. As we apply a surface integral, only the imaginary exponential remains as an integrand, which leads to the following relation in order to avoid a trivial result:

$$\begin{array}{c}\hfill {\displaystyle \int \int {\mathrm{e}}^{i{\underline{k}}_1·\underline{x}}}{\mathrm{e}}^{-i{\underline{k}}_0·\underline{x}}\mathrm{d}\mathrm{\Sigma }\ne 0\iff {\underline{k}}_0={\underline{k}}_1.\end{array}$$(24)

The second relation can be derived similarly based on the projection of the second term of the continuity equation this time. Once again applying the surface integral, only the imaginary exponential remains as an integrand, which leads us to the following relation in order to avoid a trivial result:

$$\begin{array}{c}\hfill {\displaystyle \int \int {\mathrm{e}}^{i({\underline{k}}_1+{\underline{k}}_2)·\underline{x}}}{\mathrm{e}}^{-i{\underline{k}}_0·\underline{x}}\mathrm{d}\mathrm{\Sigma }\ne 0\iff {\underline{k}}_0={\underline{k}}_1+{\underline{k}}_2.\end{array}$$(25)

Using these two established relations (Eqs. (24) and (25)) derived from the projection in the specific direction ${\underline{k}}_0$, after some algebra, we can rewrite our system of equations as a convolution product in the Fourier space and obtain the following for the continuity, momentum, energy, and radiative transport (we note that we replaced the summation index k₁ with k to simplify the notations):

$$\begin{array}{cc}& i\sigma {\rho }_{k_0}^{\prime }+{\displaystyle \sum _k}\{i{\underline{k}}_0·{\rho }_k^{\prime }{\underline{v}}_{{h,k}_0-k}+i{\underline{k}}_0·{\rho }_k{\underline{v}}_{{h,k}_0-k}\prime +\frac{\partial }{\partial z}({\rho }_k^{\prime }{v}_{{z,k}_0-k}+{\rho }_k{v}_{{z,k}_0-k}\prime )\}=0,\hfill \end{array}$$(26)

$$\begin{array}{cc}& i\sigma {\underline{v}}_{k_0}^{\prime }+{\displaystyle \sum _k}\{i{\underline{v}}_{h,k}^{\prime }·(\underline{k_0}-\underline{k})\otimes {\underline{v}}_{k_0-k}+i{\underline{v}}_{h,k}·(\underline{k_0}-\underline{k})\otimes {\underline{v}}_{k_0-k}\prime +v_{z,k}\prime \frac{\partial }{\partial z}\left({\underline{v}}_{k_0-k}\right)+v_{z,k}\frac{\partial }{\partial z}({\underline{v}}_{k_0-k}\prime )\}\hfill \\ \hfill & =-{\displaystyle \sum _k}\{i\underline{k}{P}_k^{\prime }{\left(\frac{1}{\rho }\right)}_{k_0-k}+i\underline{k}{P}_k{\left(\frac{1}{\rho }\right)}_{k_0-k}\prime +{\underline{e}}_z[\frac{\partial P_k^{\prime }}{\partial z}{\left(\frac{1}{\rho }\right)}_{k_0-k}+\frac{\partial P_k}{\partial z}{\left(\frac{1}{\rho }\right)}_{k_0-k}\prime ]\}+{\underline{g}}_{k_0}^{\prime },\hfill \end{array}$$(27)

$$\begin{array}{cc}& {\displaystyle \sum _k}\{{(\frac{{\rho }^{\prime }}{\rho }+\frac{T^{\prime }}{T})}_k{\left(\rho T\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_{k_0-k}+{\left(\rho T\right)}_{k_0-k}{\left(\frac{\mathit{Ds}}{\mathit{Dt}}\right)}_k^{\prime }\}=-i{\underline{k}}_0·{\underline{F}}_{{h,k}_0}^{\prime }-\frac{\partial }{\partial z}\left(F_{{z,k}_0}^{\prime }\right),\hfill \end{array}$$(28)

$$\begin{array}{cc}& {\underline{F}}_{k_0}^{\prime }=-i{\displaystyle \sum _k}\{{\mathcal{X}}_k^{\prime }(\underline{k_0}-\underline{k})T_{k_0-k}+{\mathcal{X}}_{k_0-k}\left(\underline{k}\right)T_k^{\prime }\}-{\underline{e}}_z{\displaystyle \sum _k}\{{\mathcal{X}}_k^{\prime }\frac{\partial }{\partial z}\left(T_{k_0-k}\right)+{\mathcal{X}}_{k_0-k}\frac{\partial }{\partial z}\left(T_k^{\prime }\right)\}.\hfill \end{array}$$(29)

The advantages of this 2D Fourier decomposition and then projection are as follows. First and foremost, the 2D Fourier space turns a system of coupled partial differential equations into a system of differential equations with only one independent variable, the radial coordinate z, for each value of k₀ (we note that k and k₀ are evaluated at the nodes of the simulation box). Once its value is set, the algebraic form of our system of equations is now in an ideal state to be solved by many different numerical schemes. Presently, a large quantity of unknowns or data are no longer considered an inconvenience with the capabilities of modern computers, which are able to store massive amounts of data and carry out efficient parallelization. However, reaching a simpler version of an equation can be a key factor in the stability and efficiency of the numerical solver. Second, we obtained a natural distinction between the horizontal average equations (comparable to the 1D model currently in use) when the value of k₀ is set to zero in our system of equations, and the different scales of turbulence, any value of k₀ different from 0. Hence, our model allows us to assess the impact of all the different scales of the turbulent spectrum on the stellar oscillations and ultimately their individual contributions to the eigenfrequencies and damping rates of the modes. Finally, being in Fourier space does not induce any drawbacks to our model, because once a solution is found, we can easily access the real space by a simple inverse Fourier transform.

4. Scales of our formalism

We now discern the horizontal average (k₀ = 0) from fluctuation (k₀ ≠ 0) equations and derive the final form of our new global nonadiabatic oscillation equations in a full 3D framework. While our formalism sets the mathematical foundation of a new 3D stellar oscillation model, it has to be borne in mind that a solution to this system of equations is only achievable via the use of a numerical scheme and an iterative method. Although a specific numerical scheme to solve the problem is out of the scope of this paper and will be presented in a subsequent paper of this series, it is nevertheless implied that, during the iterative process, a decoupling will be performed between the average unknowns and the many fluctuations (see Sect. 5). In the hydrodynamic simulation, the average (k = 0) of physical quantities is larger than their fluctuations (k ≠ 0), with the exception of the velocities. In Figs. 1 (temperature) and 2 (internal energy), we can clearly observe that the fluctuations are indeed small in comparison to the horizontal average values, being around 10%–20% of the horizontal average values for the most part and very locally up to 40%. Hence, we assume that this trend will be followed by the 3D eigenfunctions. This hypothesis will have to be verified a posteriori and, in the first instance, leads us to first limit our approach to quasi-radial modes and, later, to nonradial modes of low l degrees. Consequently, a small change of any fluctuation is not expected to lead to a significant modification of the Fourier series, which is dominated by the average term. Therefore, in the following sections when dealing with the average equations, it is assumed that the fluctuations (unknowns for each k ≠ 0) were determined in a previous step of the algorithm and are then fixed when average equations are to be solved. The same assumption is applied when fluctuation equations are considered, although in that case the average unknowns (k = 0) are fixed and assumed solved in a previous step. We note that, although the fluctuation oscillations are small, they impact the average oscillations, and their contribution as a whole has a dominant effect on the damping and therefore on the modal contribution to the surface effects. Indeed, these fluctuation oscillations constitute, amongst other things, the perturbation of the turbulent pressure, which is known to be a major contribution to the damping of the modes in solar-like oscillations (Houdek & Dupret 2015; Samadi et al. 2015).

4.1. Average equations k₀ = 0

The equations for the horizontal average perturbations correspond to k₀ = 0 in the general form of our system (Eqs. (26)–(29)). After some algebra, we get the following for the horizontal average perturbations:

$$\begin{array}{cc}& i\sigma {\rho }_0\prime +{\displaystyle \sum _k}\{\frac{\partial }{\partial z}({\rho }_k^{\prime }{v}_{z,k}^{\ast }+{\rho }_k^{\ast }{v}_{z,k}\prime )\}=0,\hfill \end{array}$$(30)

$$\begin{array}{cc}& i\sigma v_{z,0}\prime +{\displaystyle \sum _k}\{-i{\underline{v}}_{h,k}\prime ·\underline{k}{v}_{z,k}^{\ast }+i{\underline{v}}_{h,k}^{\ast }·\underline{k}{v}_{z,k}\prime +\frac{\partial }{\partial z}({\underline{v}}_{z,k}\prime {\underline{v}}_{z,k}^{\ast })\}=-{\displaystyle \sum _k}\{\frac{\partial P_k^{\prime }}{\partial z}{\left(\frac{1}{\rho }\right)}_k^{\ast }+\frac{\partial P_k^{\ast }}{\partial z}{\left(\frac{1}{\rho }\right)}_k^{\prime }\}-g_0^{\prime },\hfill \end{array}$$(31)

$$\begin{array}{cc}& i\sigma {\underline{v}}_{h,0}^{\prime }+{\displaystyle \sum _k}\{-i{\underline{v}}_{h,k}\prime ·\underline{k}\otimes {\underline{v}}_{h,k}^{\ast }+i{\underline{v}}_{h,k}^{\ast }·\underline{k}\otimes {\underline{v}}_{h,k}\prime +v_{z,k}\prime \frac{\partial }{\partial z}\left({\underline{v}}_{h,k}^{\ast }\right)+v_{z,k}^{\ast }\frac{\partial }{\partial z}({\underline{v}}_{h,k}\prime )\}={\displaystyle \sum _k}\{-i\underline{k}{P}_k^{\prime }{\left(\frac{1}{\rho }\right)}_k^{\ast }+i\underline{k}{P}_k^{\ast }{\left(\frac{1}{\rho }\right)}_k^{\prime }\},\hfill \end{array}$$(32)

$$\begin{array}{cc}& {\displaystyle \sum _k}\{{(\frac{{\rho }^{\prime }}{\rho }+\frac{T^{\prime }}{T})}_k{\left(\rho T\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_k^{\ast }+{\left(\rho T\right)}_k^{\ast }{\left(\frac{\mathit{Ds}}{\mathit{Dt}}\right)}_k^{\prime }\}=-\frac{\partial }{\partial z}(F_{z,0}\prime ),\hfill \end{array}$$(33)

$$\begin{array}{cc}& F_{z,0}\prime =-{\displaystyle \sum _k}\{{\mathcal{X}}_k^{\prime }\frac{\partial }{\partial z}\left(T_k^{\ast }\right)+{\mathcal{X}}_k^{\ast }\frac{\partial }{\partial z}\left(T_k^{\prime }\right)\},\hfill \end{array}$$(34)

$$\begin{array}{cc}& {\underline{F}}_{h,0}^{\prime }=-i{\displaystyle \sum _k}\{-{\mathcal{X}}_k^{\prime }\left(\underline{k}\right)T_k^{\ast }+{\mathcal{X}}_k^{\ast }\left(\underline{k}\right)T_k^{\prime }\}.\hfill \end{array}$$(35)

As we are working with real data from the 3D hydrodynamic simulation, A_−k is simply – according to the properties of Fourier transform – the complex conjugate of A_k (noted $A_k^{*}$). Here, we note that we have separated the equations of momentum conservation (Eqs. (31) and (32)) as well as those of radiative transport (Eqs. (34) and (35)) into a radial and a horizontal component. Following a direct investigation of the horizontal radiative flux ${\underline{F}}_{h,0\prime }$ in this system of equations, we can immediately note that this latter only appears in Eq. (35) and is fully determined by the fluctuations, because the first term of the sum (k = 0) is trivial. Regarding the horizontal velocity ${\underline{v}}_{h,0}^{\prime }$, a closer examination of the system leads to the conclusion that this unknown is decoupled from the remaining horizontal average equations. Indeed, the horizontal average velocity only appears in its own equation (Eq. (32)) and is fully determined by the average radial velocity and the fluctuations, as shown in the following, where fluctuations are gathered in one variable ${\underline{M}}_{h0k}$:

$$\begin{array}{c}\hfill i\sigma {\underline{v}}_{h,0}^{\prime }+v_{z,0}\prime \frac{\partial }{\partial z}\left({\underline{v}}_{h,0}^{\ast }\right)+v_{z,0}^{\ast }\frac{\partial }{\partial z}\left({\underline{v}}_{h,0}^{\prime }\right)+{\underline{M}}_{h0k}=0.\end{array}$$(36)

Therefore, this equation can be solved in a later step once the radial velocity has been solved or even ignored if the perturbation of the horizontal velocity is not important for a specific use of this formalism.

When dealing with the average equations of our Fourier expansion of the perturbed Navier-Stokes equations, as shown above, two of the six equations can be decoupled from the four remaining ones. However, in those four remaining equations, we highlight the presence of more unknowns than equations even when fixing the values of all the fluctuations. However, those unknowns are not totally independent and can be reexpressed via an equation of state. For this work, we use the pressure $P_0^{\prime }$, entropy $s_0^{\prime }$, radial velocity ${\underline{v}}_{z,0\prime }$, and the radial radiative flux $F_{z,0}^{\prime }$ as independent variables. Every other unknown perturbation is then expressed via an equation of state in terms of pressure and entropy. Appendix A provides an explanation of how each equation of state was obtained. Here, we only reiterate the general form of the equation of state used in the case k₀ = 0 (average perturbations):

$$\begin{array}{c}\hfill X_0\prime =X_{P_0}{P}_0\prime +X_{s_0}{s}_0\prime +X_{\mathit{eosk}},\end{array}$$(37)

where X_eosk regroups all the fluctuations as detailed in Appendix A.

We also restate here that, in the following developments, we consider any order k (k ≠ 0) unknown fluctuation of our problem to be fixed and determined in the previous step of the numerical scheme. Therefore, we regroup all those terms under one variable as demonstrated in the horizontal velocity equation (Eq. (36)) or in the following continuity equation (see Sect. 4.1.1). To obtain an eigenvalue problem, we introduce the following decomposition of the eigenfrequency σ:

$$\begin{array}{c}\hfill \sigma ={\sigma }_0+\delta \sigma ,{\sigma }^2={\sigma }_0^2+2{\sigma }_0\delta \sigma ,{\sigma }^{-1}={\sigma }_0^{-1}-\frac{\delta \sigma }{{\sigma }_0^2},\end{array}$$(38)

where σ₀ is the initial eigenfrequency guess and δσ corresponds to the correction.

With this decomposition, the system of equations takes the canonical form of an eigenvalue problem:

$$\begin{array}{c}\hfill \underline{\underline{A}}\underline{x\prime }=\delta \sigma \underline{\underline{B}}\underline{x^{\prime }}.\end{array}$$(39)

Under this form (Eq. (39)), the global eigenvalue problem can be solved by various numerical algorithms, such as the inverse iteration algorithm, to compute the eigenvectors. This algorithm and the decomposition of the eigenvalue are commonly used in nonadiabatic codes for stellar oscillations, where the adiabatic eigenvalue is typically taken as the initial guess. As the surface effects are small (Christensen-Dalsgaard et al. 1996; Houdek & Dupret 2015; Samadi et al. 2015), δσ is much smaller than σ₀ and the numerical convergence is fast. Let us now take a closer look at each equation describing the average perturbations in order to derive their final form in our formalism.

4.1.1. Average continuity equation

From the general form of the average continuity equation (Eq. (30)), we can split the summation and separate the first term (k = 0) from all the remaining terms (k ≠ 0). By doing so, we can regroup all the fluctuations under one variable C_0bk and simplify the form of the equation as:

$$\begin{array}{c}\hfill i\sigma {\rho }_0^{\prime }+\frac{\partial }{\partial z}({\rho }_0^{\prime }{v}_{z,0}^{\ast }+{\rho }_0^{\ast }{v}_{z,0}\prime )+\underset{C_{0k}}{\underbrace{{\sum }_{k,k\ne 0}\{\frac{\partial }{\partial z}({\rho }_k^{\prime }{v}_{z,k}^{\ast }+{\rho }_k^{\ast }{v}_{z,k}\prime )\}}}=0.\end{array}$$(40)

Using the equations of state for reexpressing the density from Appendix A and our decomposition of σ (Eqs. (38)), we obtain the final form of the average continuity equation:

$$\begin{array}{cc}& C_{0\mathrm{d}P}\frac{\partial P_0\prime }{\partial z}+C_{0P}{P}_0\prime +C_{0\mathrm{d}s}\frac{\partial s_0\prime }{\partial z}+C_{0s}{s}_0\prime +C_{0\mathrm{d}vz}\frac{\partial v_{z,0}\prime }{\partial z}+C_{0vz}{v}_{z,0}\prime +C_0=\delta \sigma [-i{\rho }_{P_0}{P}_0^{\prime }-i{\rho }_{s_0}{s}_0^{\prime }-i{\rho }_{\mathit{eosk}}].\hfill \end{array}$$(41)

4.1.2. Average radial momentum equation

From the radial component of the momentum equation (Eq. (31)), we immediately have

$$\begin{array}{c}\hfill i\sigma v_{z,0}^{\prime }+\frac{\partial }{\partial z}({\underline{v}}_{z,0}\prime {\underline{v}}_{z,0}^{\ast })+\frac{\partial P_0^{\prime }}{\partial z}{\left(\frac{1}{\rho }\right)}_0^{\ast }+\frac{\partial P_0^{\ast }}{\partial z}{\left(\frac{1}{\rho }\right)}_0\prime +g_0^{\prime }+M_{z0k}=0.\end{array}$$(42)

In this expression (Eq. (42)), the Eulerian perturbation of a body force, the gravity $g_0^{\prime }$, appears. Unlike state variables, where we have a set of equations with which to substitute them in terms of our independent variables, a body force cannot be expressed in such a manner. Nonetheless, an expression for the gravity is available:

$$\begin{array}{c}\hfill g=\frac{Gm}{r^2}.\end{array}$$(43)

For this specific case, a Lagrangian perturbation will be implied instead of the Eulerian one used so far. Indeed the Lagrangian perturbation is more convenient when dealing with the mass, because we are in the superficial layer of the star, and the mass in the gravity expression can be approximated by the total mass of the star, whose envelope has a negligible mass. Consequently, the Lagrangian perturbation of the mass is trivial and the Eulerian perturbation of the gravity can be recovered via the well-known Euler–Lagrange relation. After some algebra, an expression for the perturbation of gravity is derived as follows:

$$\begin{array}{c}\hfill g^{\prime }=-\frac{G}{r^2}\frac{\partial m}{\partial z}{\xi }_r=i\frac{4\pi G\rho }{\sigma }{v}_z^{\prime }\Rightarrow g_0^{\prime }=\frac{1}{\sigma }{M}_{0g}{v}_{z,0}^{\prime }.\end{array}$$(44)

Using Eq. (44) in combination with the equation of state for the specific volume (see Appendix A) and our decomposition of σ (Eq. (38)) in Eq. (42), the radial component of the momentum reads as

$$\begin{array}{cc}& M_{z0\mathrm{d}vz}\frac{\partial v_{z,0}^{\prime }}{\partial z}+M_{z0vz}{v}_{z,0}^{\prime }+M_{z0\mathrm{d}P}\frac{\partial P_0^{\prime }}{\partial z}+M_{z0P}{P}_0^{\prime }+M_{z0s}{s}_0^{\prime }+M_{z0}=\delta \sigma \left[(-i+{\sigma }_0^{-2}{M}_{0g})v_{z,0}^{\prime }\right].\hfill \end{array}$$(45)

4.1.3. Average horizontal momentum equation

As already demonstrated in the beginning of this section, the horizontal component of the momentum equation (Eq. (36)) is decoupled from our system of equations and can be solved once the radial velocity is computed. Nevertheless, we reexpress it in a similar form to the remaining equations. Our decomposition of σ, as in Eq. (38), leads, after some algebra, to

$$\begin{array}{c}\hfill M_{h0\mathrm{d}vh}\frac{\partial }{\partial z}\left({\underline{v}}_{h,0}^{\prime }\right)+M_{h0vh}{\underline{v}}_{h,0}^{\prime }+{\underline{M}}_{h0vz}{v}_{z,0}^{\prime }+{\underline{M}}_{h0k}=\delta \sigma [-i{\underline{v}}_{h,0}^{\prime }].\end{array}$$(46)

4.1.4. Average radial radiative flux equation

Given that the horizontal component of the radiative flux is entirely determined by the fluctuations (Eq. (35)), no further developments are required. However, this is not the case for the radial component of the radiative flux (Eq. (34)), which can be written under the following form:

$$\begin{array}{c}\hfill F_{z,0}\prime =-{\mathcal{X}}_0\prime \frac{\partial }{\partial z}\left(T_0^{\ast }\right)+{\mathcal{X}}_0^{\ast }\frac{\partial }{\partial z}(T_0\prime )-R_{z0k}.\end{array}$$(47)

In this expression, we still need to express ${\mathcal{X}}_0^{\prime }$ and $T_0^{\prime }$ in terms of our chosen set of independent variables via the use of equations of state (see Appendix A). After some algebra, Eq. (34) becomes

$$\begin{array}{c}\hfill F_{z,0}\prime =R_{z0\mathrm{d}P}\frac{\partial P_0^{\prime }}{\partial z}+R_{z0P}{P}_0^{\prime }+R_{z0\mathrm{d}s}\frac{\partial s_0^{\prime }}{\partial z}+R_{z0s}{s}_0^{\prime }+R_{z0}.\end{array}$$(48)

4.1.5. Average energy equation

Based on the energy equation (Eq. (33)) and grouping all the fluctuation (k ≠ 0) terms, we have

$$\begin{array}{c}\hfill {(\frac{{\rho }^{\prime }}{\rho }+\frac{T^{\prime }}{T})}_0{\left(\rho T\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_0^{\ast }+{\left(\rho T\right)}_0^{\ast }{\left(\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_0\prime +\frac{\partial }{\partial z}(F_{z,0}\prime )+E_{0k}=0.\end{array}$$(49)

We resort to equations of state (see Appendix A) to obtain an equation expressed solely in terms of our independent variables. So far, we have kept the material derivative of the entropy in our system of equations. However, a material derivative is not compatible with a numerical scheme, because all other derivatives are partial ones. To substitute this term in the energy equation (Eq. (49)), we use the expression (Eq. (B.3)) developed in Appendix B. Finally, applying our decomposition of σ (Eq. (38)) leads to

$$\begin{array}{cc}& E_{0\mathrm{d}s}\frac{\partial s_0^{\prime }}{\partial z}+E_{0s}{s}_0^{\prime }+E_{0P}{P}_0^{\prime }+E_{0\mathrm{d}vz}\frac{\partial v_{z,0}^{\prime }}{\partial z}+\frac{\partial }{\partial z}(F_{z,0}\prime )+E_0=\delta \sigma [-i{\left(\rho T\right)}_0^{\ast }{s}_0^{\prime }].\hfill \end{array}$$(50)

4.2. Fluctuation equations k₀ ≠ 0

We have now established the final form of the average equations of our Fourier expansion of the perturbed Navier-Stokes equations. The second part of the problem consists in describing the fluctuation equations following our formalism in the 2D Fourier space (Eqs. (26)–(29)). For each fluctuation, we have an associated value of k₀ leading to a system of equations to solve for this peculiar fluctuation. Moreover, due to the convolution products, each fluctuation is coupled to the average term and all the other fluctuations. Trying to solve such a system of equations would simply be too costly in terms of computation time in the current form of the fluctuation equations. With a full coupling between the scales, this corresponds to the perturbed full 3D Navier-Stokes equations expressed in the Fourier space. Therefore, in order to obtain a solvable system, we have to make some simplifying assumptions.

In Sect. 4, we stated an important hypothesis according to which the average term of the Fourier series dominates the others. Based on this assumption, we split our system into average and fluctuation equations. In the fluctuation equations, this assumption, $\Vert x_0\Vert ,\Vert x_0^{\prime }\Vert \gg \Vert x_k\Vert ,\Vert x_k^{\prime }\Vert$, leads to a simplification of the convolution product in our system of equations (Eqs. (26)–(29)). This is also related to the real space where the average quantities are much larger than the fluctuations. In consequence, for each fluctuation equation (k₀ fixed) in each convolution product, all terms are neglected, except for the dominating one corresponding to the coupling between a fluctuation and an average value. With this working hypothesis, two important observations can be made. On the one hand, we decouple each fluctuation (except through the average values) and, on the other hand, we have reached a state for our system of equations where parallel algorithms can be used to obtain a numerical solution. Following our working hypothesis, we can derive a system of equations for each fluctuation where, as in the average equations, the momentum equation and radiative transport have been split into a radial and a horizontal component.

$$\begin{array}{cc}& i\sigma {\rho }_k^{\prime }+i\underline{k}·({\rho }_0^{\prime }{\underline{v}}_{h,k}+{\rho }_0{\underline{v}}_{h,k}\prime +{\rho }_k^{\prime }{\underline{v}}_{h,0}+{\rho }_k{\underline{v}}_{h,0}^{\prime })+\frac{\partial }{\partial z}({\rho }_0^{\prime }{v}_{z,k}+{\rho }_0{v}_{z,k}\prime +{\rho }_k^{\prime }{v}_{z,0}+{\rho }_k{v}_{z,0}^{\prime })=0,\hfill \end{array}$$(51)

$$\begin{array}{cc}& i\sigma v_{z,k}\prime +\frac{\partial }{\partial z}(v_{z,0}^{\prime }{v}_{z,k}+v_{z,0}{v}_{z,k}\prime )+i{\underline{v}}_{h,0}^{\prime }·\underline{k}{v}_{z,k}+i{\underline{v}}_{h,0}·\underline{k}{v}_{z,k}\prime =-[\frac{\partial P_0^{\prime }}{\partial z}{\left(\frac{1}{\rho }\right)}_k+\frac{\partial P_0}{\partial z}{\left(\frac{1}{\rho }\right)}_k^{\prime }+\frac{\partial P_k^{\prime }}{\partial z}{\left(\frac{1}{\rho }\right)}_0+\frac{\partial P_k}{\partial z}{\left(\frac{1}{\rho }\right)}_0\prime ],\hfill \end{array}$$(52)

$$\begin{array}{cc}& i\sigma {\underline{v}}_{h,k}\prime +v_{z,0}^{\prime }\frac{\partial }{\partial z}\left({\underline{v}}_{h,k}\right)+v_{z,0}\frac{\partial }{\partial z}({\underline{v}}_{h,k}\prime )+v_{z,k}\prime \frac{\partial }{\partial z}\left({\underline{v}}_{h,0}\right)+v_{z,k}\frac{\partial }{\partial z}\left({\underline{v}}_{h,0}^{\prime }\right)+i{\underline{v}}_{h,0}^{\prime }·\underline{k}\otimes {\underline{v}}_{h,k}+{\underline{v}}_{h,0}·\underline{k}\otimes {\underline{v}}_{h,k}\prime =-i\underline{k}[P_k^{\prime }{\left(\frac{1}{\rho }\right)}_0+P_k{\left(\frac{1}{\rho }\right)}_0\prime ],\hfill \end{array}$$(53)

$$\begin{array}{cc}& {(\frac{{\rho }^{\prime }}{\rho }+\frac{T^{\prime }}{T})}_0{\left(\rho T\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_k+{(\frac{{\rho }^{\prime }}{\rho }+\frac{T^{\prime }}{T})}_k{\left(\rho T\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_0+{\left(\rho T\right)}_k{\left(\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_0^{\prime }+{\left(\rho T\right)}_0{\left(\frac{\mathrm{D}s}{\mathrm{D}t}\right)}_k^{\prime }=-i\underline{k}·{\underline{F}}_{h,k}\prime -\frac{\partial }{\partial z}(F_{z,k}\prime )\hfill \end{array}$$(54)

$$\begin{array}{cc}& F_{z,k}\prime =-[{\mathcal{X}}_0\prime \frac{\partial }{\partial z}\left(T_k\right)+{\mathcal{X}}_0\frac{\partial }{\partial z}\left(T_k^{\prime }\right)+{\mathcal{X}}_k^{\prime }\frac{\partial }{\partial z}\left(T_0\right)+{\mathcal{X}}_k\frac{\partial }{\partial z}(T_0\prime )],\hfill \end{array}$$(55)

$$\begin{array}{cc}& {\underline{F}}_{h,k}\prime =-i[{\mathcal{X}}_0\prime \left(\underline{k}\right)T_k+{\mathcal{X}}_0\left(\underline{k}\right)T_k^{\prime }].\hfill \end{array}$$(56)

We note that the index of the fluctuation is now denoted by k and no longer k₀ in order to avoid confusion and to simplify the notation, because a given fluctuation is only coupled with the average (index 0) from now on.

While there is only one set of equations for the average perturbations in every scenario possible, depending on the grid resolution of the hydrodynamic simulation, the number of turbulent scales can increase significantly. If we assume a 3D simulation box with N by M nodes in the horizontal plane, applying the 2D Fourier transform would lead to N × M fluctuations and consequently the same number of equations in the system to be solved. A robust example is provided by the case of the low-resolution CO⁵BOLD simulation shown in Fig. 2, which involves 35 721 different fluctuations, compared to the high-resolution simulation in Fig. 1, which implies 264 196. With such an important and rapidly increasing number of fluctuations to be solved, reducing the number of equations describing a fluctuation can lead to a drastic reduction in computation time and improvement of the parallelization of our formalism.

We propose a substitution of the equations of momentum (Eqs. (52) and (53)) in the continuity equation (Eq. (51)), as well as a substitution of the radiative flux (Eqs. (55) and (56)) in the conservation of energy equation (Eq. (54)), giving rise to, respectively, the continuity-momentum and energy–radiation equations (see Sects. 4.2.1 and 4.2.2). Similarly to the average equations, at this point, we still have more unknowns than equations and some unknowns will be reexpressed via an equation of state, whose general form, when fluctuations are considered, is

$$\begin{array}{c}\hfill X_k^{\prime }=X_P{P}_k^{\prime }+X_s{s}_k^{\prime }+X_{eos0}.\end{array}$$(57)

For more details regarding the equations of state and there developments, please refer to Appendix A. When dealing with the equations for the fluctuations of our Fourier expansion of the perturbed Navier-Stokes equations, we keep the pressure $P_k^{\prime }$ and entropy $s_k^{\prime }$ as independent variables. In the following developments, as mentioned in Sect. 3.2, we consider that any average unknown or perturbation of our problem is fixed and was determined in the previous step when we solved for the average equations. Therefore, we regroup all those under one variable.

4.2.1. Continuity–momentum equation

To obtain the continuity–momentum equation, we start from the continuity equation (Eq. (51)), which is multiplied by iσ, and where all known terms, average and simulation data, are regrouped in C_k0. After some algebra, we have

$$\begin{array}{c}\hfill -({\sigma }^2+\sigma \underline{k}·{\underline{v}}_{h,0}){\rho }_k^{\prime }+i\sigma \frac{\partial }{\partial z}\left({\rho }_k^{\prime }{v}_{z,0}\right)+\sigma C_{k0}+i\underline{k}·({\rho }_{0}i\sigma {\underline{v}}_{h,k}\prime )+\frac{\partial }{\partial z}({\rho }_{0}i\sigma v_{z,k}\prime )=0.\end{array}$$(58)

Regarding the radial component of the momentum equation (Eq. (52)), we assume the average radial velocity v_z, 0 to be small due to a compensation between the upflows and downflows, which automatically leads to $v_{z,0}^{\prime }{v}_{z,k}\gg v_{z,0}{v}_{z,k}^{\prime }$. Moreover, in the simulation box, with periodic boundary conditions and the type of turbulent flows in the superficial layers of stars (Rayleigh-Bénard convection), there is no physical reasons for a resulting average horizontal motion (${\underline{v}}_{h,0}$). From this observation, it is safe to assume that $\sigma \gg {\underline{v}}_{h,0}·\underline{k}$. Based on those observations and hypotheses, alongside injecting equations of state (see Appendix A), the following expression for the radial velocity fluctuation can be derived:

$$\begin{array}{cc}& i\sigma {\underline{v}}_{z,k}\prime +\frac{\partial }{\partial z}\left(v_{z,0}^{\prime }{v}_{z,k}\right)+i{v}_{z,k}({\underline{v}}_{h,0}^{\prime }·\underline{k})=-[\frac{\partial P_0^{\prime }}{\partial z}{\left(\frac{1}{\rho }\right)}_k+\frac{\partial P_0}{\partial z}{\left(\frac{1}{\rho }\right)}_k^{\prime }+\frac{\partial P_k^{\prime }}{\partial z}{\left(\frac{1}{\rho }\right)}_0+\frac{\partial P_k}{\partial z}{\left(\frac{1}{\rho }\right)}_0\prime ],\hfill \end{array}$$(59)

$$\begin{array}{cc}& ⟹i\sigma v_{z,k}\prime =M_{zkb0}+M_{\mathit{zkP}}{P}_k^{\prime }+M_{\mathit{zks}}{s}_k^{\prime }+M_{zk\mathrm{d}P}\frac{\partial P_k^{\prime }}{\partial z}.\hfill \end{array}$$(60)

As for the horizontal momentum equation (Eq. (53)), similar hypotheses can be made regarding v_z, 0 and ${\underline{v}}_{h,0}$, which means that the following terms can be neglected: $v_{z,0}\frac{\partial }{\partial z}\left({\underline{v}}_{h,k\prime }\right)$,$v_{z,k\prime }\frac{\partial }{\partial z}\left({\underline{v}}_{h,0}\right)$ and ${\underline{v}}_{h,0}·\underline{k}\otimes {\underline{v}}_{h,k\prime }$. Finally, for the horizontal component of the momentum equation, we obtain

$$\begin{array}{cc}& i\sigma {\underline{v}}_{h,k}\prime +v_{z,0}^{\prime }\frac{\partial }{\partial z}\left({\underline{v}}_{h,k}\right)+v_{z,k}\frac{\partial }{\partial z}\left({\underline{v}}_{h,0}^{\prime }\right)+i{\underline{v}}_{h,0}^{\prime }·\underline{k}\otimes {\underline{v}}_{h,k}=-i\underline{k}[P_k^{\prime }{\left(\frac{1}{\rho }\right)}_0+P_k{\left(\frac{1}{\rho }\right)}_0\prime ],\hfill \end{array}$$(61)