R. Samadi ^1,2 - Å. Nordlund ³ - R. F. Stein ⁴ - M. J. Goupil ² - I. Roxburgh ^1,2

1 - Astronomy Unit, Queen Mary, University of London, London E14NS, UK
2 - Observatoire de Paris, LESIA, CNRS UMR 8109, 92195 Meudon, France
3 - Niels Bohr Institute for Astronomy Physics and Geophysics, Copenhagen, Denmark
4 - Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan, USA

Abstract
Analyses of a 3D simulation of the upper layers of a solar convective envelope provide constraints on the physical quantities which enter the theoretical formulation of a stochastic excitation model of solar p modes, for instance the convective velocities and the turbulent kinetic energy spectrum. These constraints are then used to compute the acoustic excitation rate for solar p modes, P. The resulting values are found $\sim$ 5 times larger than the values resulting from a computation in which convective velocities and entropy fluctuations are obtained with a 1D solar envelope model built with the time-dependent, nonlocal Gough (1977) extension of the mixing length formulation for convection (GMLT).
This difference is mainly due to the assumed mean anisotropy properties of the velocity field in the excitation region. The 3D simulation suggests much larger horizontal velocities compared to vertical ones than in the 1D GMLT solar model. The values of P obtained with the 3D simulation constraints however are still too small compared with the values inferred from solar observations.
Improvements in the description of the turbulent kinetic energy spectrum and its depth dependence yield further increased theoretical values of P which bring them closer to the observations. It is also found that the source of excitation arising from the advection of the turbulent fluctuations of entropy by the turbulent movements contributes $\sim$ $65{-}75 \%$ to the excitation and therefore remains dominant over the Reynolds stress contribution. The derived theoretical values of P obtained with the 3D simulation constraints remain smaller by a factor $\sim$ 3 compared with the solar observations. This shows that the stochastic excitation model still needs to be improved.

Key words: convection - turbulence - stars: oscillations - Sun: oscillations

1 Introduction

Solar-type oscillations are believed to be stochastically excited by turbulent convection in the near-surface layers of the star. The excitation is caused by turbulent convective motions which generate acoustic energy which in turn is injected into the p modes (e.g. Goldreich & Keeley 1977). Measurements of the acoustic energy injected into solar-like oscillations are among the goals of future space seismic missions such as the COROT (Baglin & The Corot Team 1998) and Eddington (Favata et al. 2000) missions. These seismic data will make it possible to constrain the theory of the oscillation excitation and damping, to provide valuable information about the properties of stellar convection, and hence to severely constrain stellar models.

Models for stochastic excitation of stellar p modes have been proposed by several authors, (e.g. Goldreich & Keeley 1977; Osaki 1990; Balmforth 1992a; Goldreich et al. 1994; Samadi & Goupil 2001). These theoretical approaches yield the acoustic energy injected into solar-like oscillations. This offers the advantage of testing separately several properties entering the excitation mechanism which are not well understood or modeled.

Such approaches require simplifying assumptions which need to be validated before they can be used with confidence. They require an accurate knowledge of the properties of turbulent convection and, unfortunately, current observations of the solar granulation cannot provide a determination of the turbulent spectrum precise enough in the present context (Rieutord et al. 2000; Nordlund et al. 1997). On the theoretical side, theoretical models of turbulent convection, such as the mixing-length approaches or multiple size eddies approaches (e.g. Canuto & Mazzitelli 1991; Canuto et al. 1996), provide a too limited description of the characteristic scale length of the solar turbulent spectrum.

These theoretical formulations of stochastic excitation also involve scaling parameters which are determined to by the requirement that the computed values of the oscillation amplitudes give the best fit to the solar seismic measurements (e.g. Houdek et al. 1999; Samadi et al. 2001, Paper II hereafter). When the scaling parameters are so adjusted, constraints and validation on the turbulent stellar medium can only come from seismic observations of other stars. Such accurate data on the excitation rates for other stars than the Sun are not yet available.

An alternative way is then to consider results from 3D numerical simulations. They indeed enable one to compute directly the rate at which p modes are excited (e.g. this was undertaken for the Sun by Stein & Nordlund 2001). Such methods are time consuming and do not easily allow massive computations of the excitation rate for stars with different temperatures and luminosities. They can provide quantities which can be implemented in a formulation for the excitation rate P. In any case we cannot avoid to use a 1D model for computing accurate eigenfrequencies for the whole observed frequency range.

The purpose of the present paper is to provide a better insight into the excitation model with a semianalytical approach but using a model of turbulence and values of the scaling parameters derived from a 3D simulation of the solar outer layers. We consider in this work the theoretical formulation of stochastic excitation by Samadi & Goupil (2001, hereafter Paper I, see also Samadi 2001 for a detailed summary) which includes a detailed treatment of turbulent convection. This formulation involves two scaling parameters which are related to the spatial and temporal characteristics of the turbulence model. Our final goal is to test the excitation model without adjusting these parameters and without the use of the mixing-length approach for estimating convective velocities and entropy fluctuations.

The paper is organized as follows: in Sect. 2 we briefly recall the adopted formulation for estimating the rate at which turbulent convection supplies energy to the p modes (excitation rate $P(\nu)$ ). We emphasize some assumptions and approximations entering this formulation.

In Sect. 3, a 3D numerical simulation of the upper part of the solar convection zone is used in order to determine the time averaged properties of turbulent convection: this provides constraints on the ingredients involved in the theoretical expression of the excitation rate, such as scaling parameters, velocity anisotropy factor, the values of convective velocities and entropy fluctuations and the k (wavenumber) dependence of the kinetic turbulent spectrum.

These constraints are then used in Sect. 4 to compute the excitation rate $P(\nu)$ , for radial solar p modes. The results are compared with solar seismic observations as given in Chaplin et al. (1998) and with a 1D mixing-length model built according to Gough (1977)'s non-local formulation of the mixing-length theory (GMLT hereafter). In Sect. 5 we summarize our results and discuss some possible origins of the remaining discrepancies with solar seismic observations and results by Stein & Nordlund (2001).

2 Stochastic excitation

2.1 The excitation model

The rate at which turbulent motions of the convective elements supply energy to acoustic oscillation modes is computed as in Paper I. For a given mode with eigenfrequency $\omega_0$ , the excitation rate can be written as (Eqs. (58) and (59) of Paper I):

The above expression for P is mainly based on the assumption that the medium is incompressible. In other words, we adopt the Boussinesq approximation i.e. assume a homogeneous model for the turbulence and the excitation mechanism. We therefore neglect effects of the stratification in the excitation process.

2.2 The turbulence model

Let k₀(r) be the wavenumber at which energy is injected into the turbulent cascade and the energy E(k) is maximum. k₀(r) is related to the mixing-length $\Lambda \equiv \alpha ~ H_{\rm p}$ by (Paper I):

The Gaussian function is usually assumed for modeling $\chi_k(\omega)$ (e.g. Stein 1967; Goldreich & Keeley 1977) as a consequence of the turbulent nature of the medium where the stochastic excitation occurs. The Gaussian function takes the form

Let $\tau_k$ be the characteristic time correlation length of an eddy of wavenumber k. Equation (13) corresponds in the time domain to a Gaussian function with linewidth equal to $2 / \omega_k$ . Then $\omega_k= {2 / \tau_k}$ for a Gaussian time spectrum.

The energy supply rate P crucially depends on the correlation time-scale $\tau_k$ (see Paper II). Following Balmforth (1992a) we define it as:

The parameter $\lambda$ in Eq. (14) accounts for our lack of precise knowledge of the time correlation $\tau_k$ in stellar conditions. In the present paper, we assume $\lambda=1$ while $\beta$ (Eq. (12)) and $\Phi (z)$ (Eq. (4)) are given by a simulation of the upper part of the solar convective zone in Sect. 3.4 below.

2.3 Computations of the excitation rate ${P(\omega )}$

In practice, we compute the excitation rate ${P(\omega )}$ according to Eq. (1). The calculation requires the knowledge of several quantities which can be obtained either from a 1D model (Paper II) or at least partly from a 3D simulation. Comparison of the results using both options yields insights in the excitation mechanism and its modelling. Hence in the following:

$\bullet$ The velocity, entropy fluctuations, anisotropy and turbulent spectra E's are obtained from a 3D simulation as described in the next section.

$\bullet$ The mean density, $\rho_0$ , the thermodynamic quantity $\alpha_{\rm s}$ , the oscillation properties - eigenfrequencies and eigenfunctions - are calculated from a solar envelope equilibrium model and Balmforth (1992b)'s pulsation code. The envelope model is built with a treatment of convection as prescribed by the GMLT formulation and is computed in the manner of Balmforth (1992b) and Houdek et al. (1999). This solar envelope model (hereafter GMLT solar model) is identical to the one considered in Samadi et al. (2002, hereafter Paper III). In particular, it incorporates turbulent pressure (momentum flux) in the equilibrium model envelope. The entire envelope is integrated using the equations appropriate to the nonlocal mixing-length formulation by Gough (1976) and to the Eddington approximation to radiative transfer (Unno & Spiegel 1966). The equation of state included a detailed treatment of the ionization of C, N, and O, and a treatment of the ionization of the next seven most abundant elements (Christensen-Dalsgaard 1982), as well as "pressure ionization'' by the method of Eggleton, Faulkner & Flannery (Eggleton et al. 1973). In this generalization of the mixing-length approach, two additional parameters, namely a and b, are introduced which control the spatial coherence of the ensemble of eddies contributing to the total heat and momentum fluxes (a), and the degree to which the turbulent fluxes are coupled to the local stratification (b). These convection parameters are calibrated to a solar model to obtain the helioseismically inferred depth of the solar convection zone of 0.287 of the solar radius (Christensen-Dalsgaard et al. 1991). The adopted value for the shape factor $\Phi=1.3745$ , a value which provides the best fit between computed solar damping rates and measurements by Chaplin et al. (1998) (see Houdek et al. 2001). The detailed equations describing the equilibrium and pulsation models were discussed by Balmforth (1992b) and by Houdek (1996).

For implementation in Eqs. (1)-(11), the quantities from the 3D simulation are interpolated at the GMLT model mesh points. The grid of mesh points of the simulated domain is matched with the GMLT one such that w in the 3D simulation has its maximum at the same layer as in the GMLT model. In the simulation, w peaks $\sim$ 40 km above the layer at which the mean optical depth $<\tau>$ is unity while in the GMLT model, w peaks $\sim$ 130 km below the photosphere ( $<\tau>$ = 2/3).

3 Constraints from the 3D simulation

We consider a 3D simulation of the upper part of the solar convective zone obtained with the 3D numerical code developed at the Niels Bohr Institute for Astronomy, Physics and Geophysics (Copenhagen, Denmark).

The simulated domain is 3.2 Mm deep and its surface is 6 $\times$ 6 ${\rm Mm}^2$ . The grid of mesh points is 256 $\times$ 256 $\times$ 163, the total duration 27 min and the sampling time 30 s. Physical assumptions are described in Stein & Nordlund (1998).

Output of the simulation considered here are the velocity field $\vec u(x,y,z,t)$ and the entropy s(x,y,z,t). They are used to determine the quantities $\tilde s^2$ , $\Phi (z)$ , w(z), $E_{\rm s}(k,z)$ , E(k,z)which enter the excitation rate through Eqs. (2), (7), (8), (10), (11).

3.1 Fourier transforms and averaging

We compute the 2D Fourier transform, along horizontal planes, of the velocity field $\vec u$ and the entropy s, at each layer z. This provides $\vec {\hat u}(\vec k,z,t)$ and $\hat s(\vec k,z,t)$ where $\vec k$ is the wavenumber along the horizontal plane. Next we integrate $\vec {\hat u}^2(\vec k,z,t)$ and ${\hat s}^2(\vec k,z,t)$ over circles with radius k at each given layer z. Finally take a time average of the various quantities over the time series. This yields $\vec {\hat u}(k,z)$ and ${\hat s}(k,z)$ where $k = \Vert \vec k \Vert$ is the wavenumber norm.

3.2 Convective velocities and entropy fluctuations

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{3271f1.eps}\end{figure}$

Figure 1: The root mean square of the vertical component of the velocity ( $w=\sqrt {<u_z^2>-<u_z>^2}$ ) in the upper layers of a solar model is plotted versus depth for the 3D simulation (solid line) and for the 1D GMLT model (dashed line). The abscissa is the depth $z = r -R_\odot$ where $R_\odot$ is the radius at the photosphere. The w maximum corresponds to the top of the superadiabatic region and is reached at the depth $z \simeq - 130~{\rm km}$ in the GMLT model. The grid of mesh points of the simulated domain is adjusted so that the w maxima of the 3D simulation and the GMLT coincide at the same layer.

Figures 1 and 2 present w(z), $\tilde s^2 (z)$ versus depth for the 3D simulation. For comparison purpose, the plots also show w and $\tilde s^2$ obtained with the GMLT solar model.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{3271f2.eps}\end{figure}$

Figure 2: Same as Fig. 1 for the mean square of the entropy fluctuations ( $\tilde s^2$ ). The peak is narrower than for the velocity w because $\tilde s^2$ scales approximatively as w⁴.

The vertical velocity GMLT w is larger at the top of superadiabatic region but smaller just beneath compared to values from the simulation. The GMLT $\tilde s^2$ is larger than in the simulation ( $\sim$ 20%). This explains that the relative contribution of the entropy source term to the excitation is overestimated with the GMLT model (see Sect. 4.1.1). Differences between the GMLT and the simulation are likely to be related to differences in the convective efficiency: GMLT is less efficient than the 3D simulation. Indeed as pointed out by Houdek & Gough (2002), a single eddy approach such as the GMLT results in a larger peak for the superadiabatic gradient.

3.3 Velocity anisotropy at large scale

As it will be shown in Sect. 4.1.2, the value of $\Phi$ plays a crucial role in controlling the depth of the excitation region and therefore the total amount of acoustic energy injected into the oscillation modes.

Figure 3 displays the anisotropy factor $\Phi$ versus depth z for the 3D simulation. $\Phi (z)$ sharply decreases from the value $\Phi=3$ at the top of the CZ down to $\Phi =2$ and then slowly decreases to reach the value $\Phi\simeq 1.3$ at the bottom of the simulation. The decrease of $\Phi (z)$ with depth is explained first by the onset of the convection and the formation of convective plumes at $z \sim 0$ and then by the relative increase in number of the plumes inward in the simulation. Indeed, plumes are highly anisotropic structures whereas turbulent cells are quite isotropic. The turbulent Mach number increases with z and reaches its maximum value at the top of the CZ. Therefore the fluid is more turbulent outward in the atmosphere. Consequently the number of turbulent isotropic cells increases with z up to the top of the CZ whereas the number of plumes remains roughly constant. The medium is thus more isotropic outward than inward.

In most of the excitation region, the value of $\Phi =2$ consistent with the BV-MLT is in better agreement with the values of $\Phi (z)$ inferred from the simulation compared to the value $\Phi=1.3745$ which must be imposed for the GMLT solar model in order to match the observed solar damping rates.

3.4 Turbulent kinetic energy spectrum E(k)

Variations of E and $E_{\rm s}$ with k at different depths z are depicted in Fig. 4. The spectra clearly show two regimes: at large scale (small values of k), the spectra increase approximately as k⁺¹ which can roughly be explained using dimensional analysis. At small scale (large values of k), the spectra decrease very rapidly with k. The Kolmogorov law (k ^-5/3) is observed only over a small k-range. Departures of the computed spectra from a Kolmogorov law at high values of k can be explained by the finite resolution of the simulation spatial grid.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{3271f4.eps}\par\includegraphics[width=8.8cm,clip]{3271f5.eps}\end{figure}$

Figure 4: Turbulent kinetic energy spectra E (top) and $E_{\rm s}$ (bottom) from the simulation are plotted versus k and for different depths z in the simulation. The straight solid lines delimitate the slopes k¹ and k^-5/3 of the EKS spectrum (Eq. (21)). Intersection of the slopes determines k₀^E, the scale of maximum energy at each depth z.

The main characteristics of the kinetic spectrum E(k,z) - k dependency - derived from the 3D simulation are approximatively reproduced by an analytical expression which was considered by Musielak et al. (1994), namely the "Extended Kolmogorov Spectrum'' (EKS hereafter) defined in Musielak et al. (1994) as:

At each layer, k₀^E is determined by imposing that the EKS, as defined above, matches the turbulent spectrum E(k,z) calculated from the simulation as well as possible. This then fixes the z dependency of k₀^E. A similar procedure is applied for $E_{\rm s}(k,z)$ for which we introduce $k_0^{E_{\rm s}}$ . All spectra satisfy their respective normalisation condition as given in Eq. (20).

For comparison, in Fig. 7, the "Nesis Kolmogorov Spectrum'' (NKS hereafter) determined from solar observations of Nesis et al. (1993) is also shown. The NKS scales as k^-5 in the energy injection region for k < k₀^E and down to $k_{\rm min}=0.7~k_0^E$ . This spectrum does not agree with turbulent spectrum E(k,z) calculated from the simulation. In particular, the NKS underestimates the velocity of the small size turbulent elements in the cascade (k>k₀) and overestimates the velocity of the turbulent with wavenumber $k \sim k_0^E$ . As we will show in Sect. 4.2, differences between the EKS and the NKS have an important impact on ${P(\omega )}$ .

If we assume that $\displaystyle k_0^{E} = k_0^{E_{\rm s}}$ , one can show that ${P(\omega )}$ scales as k₀^-4. ${P(\omega )}$ is therefore very dependent on the values reached by k₀(z) in the excitation region. Variation of k₀(z) with depth is thus shown in Fig. 5 for E and $E_{\rm s}$ : $\displaystyle k_0^{E_{\rm s}}$ and $\displaystyle k_0^{E}$ vary slowly within the excitation region.

For comparison, in Fig. 5 we have also plotted $\displaystyle k_0^{\rm MLT}(z)$ , the MLT value for k₀(z) according to Eq. (12). The scaling parameter $\beta$ in the definition of $k_0^{\rm MLT}$ is determined such that $k_0^{\rm MLT}$ and $k_0^{\rm E}$ take the same value at the layer $z \simeq - 130~{\rm km}$ where w reaches its maximum (and consequently the layer where the excitation is maximum). The derived value is $\beta =3.48$ .

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{3271f6.eps}\end{figure}$

Figure 5: The wavenumbers k₀^E (solid line) and $k_0^{E_{\rm s}}$ (dashed line) are plotted versus z (k₀^E and $k_0^{E_{\rm s}}$ are obtained by fitting, at each layer, the EKS (Eq. 21) to the computed spectra E and $E_{\rm s}$ of Fig. 4 resp., see text for details). The dotted line corresponds to $k_0^{\rm MLT}(z)$ obtained according to Eq. (12). In computing $k_0^{\rm MLT}(z)$ , we assume $\beta =3.48$ in order for $k_0^{\rm MLT}$ to match the value reached by $\displaystyle k_0^{E}$ (solid line) at the layer $z \simeq - 130~{\rm km}$ where w reaches its maximum ( $k_0^{E}=3.62~{\rm Mm}^{-1}$ at that layer).

$k_0^{\rm MLT}$ varies slowly with depth below the top of the superadiabatic region ( $z \simeq - 130~{\rm km}$ ) but increases very rapidly above. Such a behavior is explained by the rapid decrease with z of the pressure scale height $H_{\rm p}$ (which enters in the definition of $k_0^{\rm MLT}$ , Eq. (12)) in the atmosphere.

Comparison between $k_0^{\rm MLT}(z)$ and k₀^E(z) shows that the mixing-length approach does not model satisfactorily the behavior of k₀^E(z) in particular just above the layer at which w reaches its maximum value. Consequences in terms of mode excitation are investigated in Sect. 4.3.

4 Consequences in term of p modes excitation

The acoustic energy supply rate P injected into the solar oscillations is related to the rms value $v_{\rm s}$ of surface velocity as:

We derive the "observed'' P from Chaplin et al. (1998)'s seismic data according to Eq. (22) where the mode damping rate, $\eta$ , and the mode surface velocity, $v_{\rm s}$ , are obtained from Chaplin et al. (1998)'s data. The mode mass ${I / \xi_r^2(r_{\rm s})}$ is given by the GMLT model and we adopt $r_{\rm s} = R_\odot + 200$ km consistent with the observations.

Theoretical values of P are computed according to Eq. (1). In Eqs. (10) and (11) the integrations over k are performed from $k=k_{\rm min}$ (where $k_{\rm min}$ depends on the adopted turbulent spectra E and $E_{\rm s}$ ) to k=20 k₀. We checked numerically that contributions to the excitation rate from turbulent elements with $k \ga 20~k_0$ are negligible.

For the other quantities (w, $\tilde s^2$ , $\phi$ , k₀, E(k/k₀) and $E_{\rm s}(k/k_0)$ ) involved in the expression for P we investigate several possible assumptions.

4.1 Convective velocities and large scale anisotropy

In this section, the excitation rate P (Eq. (1)) is computed with the following assumptions:
- the k-dependency E and $E_{\rm s}$ is given by the analytical form of Eq. (21), also called the EKS.
- $k_0^{E}= k_0^{E_{\rm s}} = k_0^{\rm MLT}$ where $k_0^{\rm MLT}(z)$ is given in Eq. (12) with $\beta =3.48$ so that $k_0^{\rm MLT}$ takes the value reached by k₀^E(z) (solid line, Fig. 5) at the layer $z \simeq - 130~{\rm km}$ where w reaches its maximum.

For the quantities $\Phi$ , w and $\tilde s^2$ we investigate the effects of using either the values derived from the 3D simulation (see Sects. 3.2 and 3.3) or calculated with the GMLT solar model.

4.1.1 Convective velocities and entropy fluctuations

The values of w, $\tilde s^2$ and $\Phi (z)$ are fixed by the 3D simulation inside the simulation domain and by the 1D equilibrium model outside this domain. Either, if we impose zero values or if we assume quantities from the 1D MLT model, no sensitivity on the calculation of P is found.

Results are shown in Fig. 6 for P and for the relative contribution of the Reynolds stress to the total energy supply rate P. When the excitation rate $P(\nu)$ is computed with quantities derived from the 3D simulation as described in Sect. 2.3, the resulting excitation rate at maximum is found too small by a factor $\sim$ 4 compared with the observations.

Provided the appropriate value for $\Phi$ is given in the GMLT estimations (see Sect. 4.1.2 below), no significant difference is found in the excitation rate when computed with the values of w and $\tilde s^2$ from the simulation or their respective GMLT estimations.

The main effect is illustrated in the bottom panel of Fig. 6: the 3D simulation generates a larger relative contribution of the Reynolds stress to P than the GMLT model. This is explained as follows: within most part of the excitation region - except at the top of superadiabatic region - the values reached by w are larger whereas values reached by $\tilde s^2$ are smaller than their corresponding GMLT estimations.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{3271f7.eps}\par\includegraphics[width=8.8cm,clip]{3271f8.eps}\end{figure}$

Figure 6: Top: Rate P at which acoustic energy is injected into the solar radial modes. The filled dots represent P computed from Chaplin et al. (1998)'s solar seismic data according to Eq. (22). The curves represent theoretical values of P computed according Eq. (1) and for different computations of w, $\tilde s^2$ and $\Phi (z)$ : Solid line: values of w, $\tilde s^2$ and $\Phi (z)$ are fixed by the 3D simulation inside the simulated domain and by the 1D equilibrium model outside this domain. Dashed line (- - -): same as solid line but fixing $\Phi$ to BV's value $\Phi =2$ . Dot dashed line ( $-~\cdot~-~\cdot-$ ):values of w, $\tilde s^2$ and $\Phi~$ (=1.37) are fixed by the 1D equilibrium model (GMLT). Three dots dashed line ( $-\cdot \cdot \cdot -\cdot \cdot \cdot -$ ): same as the dot dashed line but fixing $\Phi$ to BV's value $\Phi =2$ . Bottom: Same as top panel for the relative contribution of the Reynolds stress, $P_{\rm R}$ , to the total acoustic energy P.

4.1.2 Velocity anisotropy at large scales

The main consequence (in term of p modes excitation) of the differences between the time averaged properties of the convective region inferred from the 3D simulation and from the GMLT solar model (Fig. 6) is due to differences in their respective anisotropy factor $\Phi$ values (Fig. 3).

Within most of the excitation region, $\Phi (z)$ is found close to $\sim$ 2 and thus larger than the value $\Phi=1.37$ assumed for the 1D equilibrium model (see Sect. 3.3 and Fig. 3). Smaller values of $\Phi$ decrease the rms total convective velocity which results in larger values of $\tau_k$ (see Eqs. (14), (15)) and therefore in a smaller depth of excitation for a given mode frequency (see Paper III for more details). Smaller values of the rms total convective velocity also induce smaller values of E in the integrand of Eq. (1). Consequently, as it is illustrated in Fig. 6, the total amount of acoustic energy injected into the modes is $\sim$ 5 times smaller for the constant value $\Phi=1.37$ compared to the constant value $\Phi =2$ (the relative contribution of the Reynolds stress to P is found $\sim$ 2 times larger in the simulation).

The effect of the depth dependency of $\Phi$ on the mode excitation is small except at high frequency. This is illustrated in the bottom panel of Fig. 6 (compare the solid line with the dashed line). Just above the top of the superadiabatic region ( $z \ga -130~$ km), $\Phi (z)$ increases rapidly with z until the value $\simeq$ 3 (Fig. 3). Most of the injection of acoustic energy into the high frequency modes occurs at the top of superadiabatic region. The high frequency modes are therefore more sensitive to this rapid increase of $\Phi (z)$ . As a consequence, the relative contribution of the Reynolds stress is larger for the high frequency modes than it is when assuming the constant value $\Phi =2$ .

4.2 Turbulent spectra

In this section we compare the excitation rate obtained assuming, for the turbulent spectra (E and $E_{\rm s}$ ) either - the EKS spectrum (Eq. (21)) with slopes given by the 3D simulation as in the previous section - or assuming the NKS spectrum from solar observations of Nesis et al. (1993) (see also Paper I).

As in Sect. 4.1 we compute P using w, $\tilde s^2$ and $\Phi$ derived from the 3D simulation and assuming that $k_0^{E}= k_0^{E_{\rm s}} = k_0^{\rm MLT}$ . The results are plotted in Fig. 8.

The NKS overestimates the maximum in P by a factor $\sim$ 1.5 while the EKS underestimates it by a factor $\sim$ 4. This is because most of the kinetic energy in the NKS is concentrated at $k \sim k_0^E$ whereas in the EKS a large part of the energy is concentrated both at large scales (k < k₀^E) and at small scales (k > k₀^E).

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{3271f10.eps}\end{figure}$

Figure 8: Acoustic energy supply rate P computed according to Eq. (1) and assuming for E(k) the EKS and the NKS plotted in Fig. 7. Dots represent the energy supply rate injected into the oscillations derived from the solar observations with the help of Eq. (22).

4.3 Effects due to the stratification of the turbulent spectrum at large scale

In Sect. 3.4 we showed that the variations of k₀^E and $k_0^{E_{\rm s}}$ with z deduced from the 3D simulation differ from the MLT estimation as given by Eq. (12) (see Fig. 5). Figure 9 presents the consequences of the z variations of k₀^E on the oscillation amplitudes (as the variations of k₀^Eand $k_0^{E_{\rm s}}$ with depth are quite similar we assume for the sake of simplicity that $k_0^{E_{\rm s}}(z)$ is equal to k₀^E(z)).

The z dependency of k₀^E causes the maximum of P to be larger than when assuming $k_0^E=k_0^{\rm MLT}$ ( $\sim$ 50% larger). This is due to the fact that in most part of the excitation region $k_0^{\rm MLT}$ is smaller than k₀^E and $k_0^{E_{\rm s}}$ except above the top of the superadiabatic region (see Fig. 5). A larger k₀^E results in a larger linewidth $\omega_k \equiv (k u_k) / \lambda$ for $\chi_k(\omega)$ hence in a larger amount of acoustic energy injected to the mode (see Eq. (13)).

Furthermore, at high frequency, P decreases with $\nu$ more rapidly than when assuming $k_0^{E}= k_0^{E_{\rm s}} = k_0^{\rm MLT}$ . Taking into account the actual variation k₀^E with z instead of assuming $k_0^E=k_0^{\rm MLT}$ makes then the $\nu$ -dependency of P at high frequency closer to that of the observed excitation spectrum. This is because, above the top of the superadiabatic region, k₀^E decreases with z whereas $k_0^{\rm MLT}$ increases with z. Indeed, the excitation of the high frequency modes occurs predominantly in the upper most part of the top of the superadiabatic region. As mentionned above the line width of $\chi_k(\omega)$ decreases with decreasing k₀. Therefore the contribution of the term $\chi_k(\omega)$ to the excitation of high frequency mode is smaller when assuming the actual variation of k₀^E with z than when assuming that k₀^E varies as $k_0^{\rm MLT}$ .

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{3271f11.eps}\end{figure}$

Figure 9: Same as Fig. 6. Solid line: the variation of k₀^E and $k_0^{E_{\rm s}}$ with z are obtained from the simulation (see Fig. 5 and Sect. 3.4). The dashed line is identical to the solid line of Fig. 6 where $k_0^E=k_0^{E_{\rm s}}=k_0^{\rm MLT}(z)$ .

5 Summary and discussion

An analysis of a 3D simulation of the upper part of the solar convective zone provides time averaged constraints upon several physical parameters which enter the theoretical expression for the supply rate of energy, P, injected into the solar p modes. These constraints are:

1) the depth dependency: of u², the mean square velocity - of w², the mean square vertical component of the velocity - of $\tilde s^2$ , the mean square values of entropy fluctuations.

2) the wavenumber (k) dependency of E and $E_{\rm s}$ the turbulent kinetic energy spectrum and the turbulent entropy spectrum respectively.

3) the depth dependency of the wavenumbers k₀^E and $k_0^{E_{\rm s}}$ , the wavenumbers at which convective energy is maximum and is injected into the turbulent inertial ranges of the turbulent kinetic energy spectra E, $E_{\rm s}$ respectively.

Differences between w² - and $\tilde s^2$ - and their respective GMLT estimations have only small consequences on the profile of the excitation rate ${P(\omega )}$ . However the values reached by w² and $\tilde s^2$ with the 3D simulation are responsible for an increase of the relative contribution of the Reynolds stress to ${P(\omega )}$ by a factor $\sim$ 1.5 at low frequency $\nu \lesssim 3$ mHz compared to the one obtained with the GMLT solar model. This is because the GMLT model overestimates $\tilde s^2$ by $\sim$ 20% at the top of excitation region and underestimates w² within most part of the excitation region by up to $\sim$ 15%.

The energy distributions E and $E_{\rm s}$ over eddies with wavenumber k obtained in the simulation scale approximately as k⁺¹in the domain $k \lesssim k_0^E$ . They therefore have approximately the same behavior as the "Extended Kolmogorov Spectrum'' (EKS) defined in Musielak et al. (1994). In contrast, their k-dependencies significantly differ from those assumed in the Nesis Kolmogorov Spectrum (NKS) which scales as k^-5 below k₀^E. The NKS predicts much larger maximum values for P than does the EKS. This is because the NKS concentrates kinetic energy in the vicinity of k₀^E.

The 3D simulation indicates that $k_0^{E}\simeq 3.6~ {\rm Mm}^{-1}$ at the top of excitation region. This corresponds to the horizontal size of the granulation ( $\sim$ 2 Mm). It is worth noting that taking the depth dependency of k₀(z) into account results in a increase of the maximum of P by as much as $\sim$ 50% and brings these values even closer to the observations. On the other hand, only minor differences are seen on the frequency dependence of the excitation rates P when using the depth dependency of k₀(z) from the simulation or assuming the form $k_0^{\rm MLT}=2 \pi/\beta \Lambda$ with $\Lambda= \alpha H_{\rm p}$ provided that $\beta$ is adjusted in order for $k_0^{\rm MLT}$ to match the value reached by k₀^E at the layer where w reaches its maximum.

The excitation rate ${P(\omega )}$ is very sensitive to the value of $\Phi$ . In the GMLT formulation, the quantity $\Phi$ is a parameter which is adjusted in order to obtain the best fit between computed solar damping rates and the solar measurements: the adopted value is $\Phi=1.37$ (see Houdek et al. 2001). On the other hand, the 3D simulation suggests a higher value within the excitation region ( $\Phi \simeq2$ ). Larger values of $\Phi$ result in an increase of the mode driving by the turbulent motions. We find that using the value $\Phi=1.37$ underestimates ${P(\omega )}$ by a factor $\sim$ 5 relatively to ${P(\omega )}$ computed with $\Phi (z)$ in the excitation region from the 3D simulation. On the other hand, using the GMLT formulation for the convective velocity with a value $\Phi\sim2$ , as suggested by the 3D simulation, yields a power ${P(\omega )}$ close to the one obtained by the 3D simulation. To fix ideas, the maximum amplitude is $\sim$ 4 cm/s, $\sim$ 8 cm/s, $\sim$ 10 cm/s when calculated with GMLT and $\Phi=1.37$ , with GMLT and $\Phi =2$ and when using velocities and $\Phi (z)$ derived from the 3 D simulation respectively. These figures must be compared to the observed maximum amplitude $\sim$ 23 cm/s.

This shows that the values of $\Phi$ found for the solar GMLT model when adjusted to the damping rates is not compatible with the actual properties of the turbulent medium in the excitation region. An improvement could come from a consistent calculation which would assume a depth dependent $\Phi (z)$ , as suggested by the simulations, in both damping rate and excitation rate computations. Damping rates are indeed expected to be sensitive to depths deeper than the excitation rate where smaller values of $\Phi$ are encountered and the simulation shows that the velocity anisotropy factor $\Phi$ decreases from 2 down to 1.3 from top of the superadiabatic region to bottom of the simulated solar region.

Without any adjustment of scaling parameters but using all the constraints inferred from the 3D simulation considered here, we find a maximum of P much larger ( $\sim$ 5 times larger) than the P maximum obtained using a 1D GMLT solar model when $\Phi$ is fixed by the observed damping rates. It is also found that the so-called entropy source term, which arises from the advection of the turbulent fluctuations of entropy by the turbulent motions, is still the dominant source of the excitation. However its contribution to the excitation is now $\sim$ 65-75% instead of $\sim$ 95% as found in Paper II.

Our computation still underestimates by a factor $\sim$ 3 the maximum value of P compared with the one derived from the solar seismic observations by Chaplin et al. (1998). Moreover the decrease of P with $\nu$ at high frequency ( $\nu \ga 3.5$ mHz) is found to be significantly smaller than the one inferred from the solar seismic observations, indicating a deficiency in the present modelling at high frequency.

As a final point, we discuss the model for the turbulent kinetic energy spectrum:

In Paper II the parameter $\lambda$ and k₀^E were adjusted - given a turbulent spectrum E(k) - so as to obtain the best possible agreement between computed and measured values of the maximum solar oscillation amplitude and its frequency position, as well as the frequency-dependence of the oscillation amplitudes. Adjustments of these scaling parameters led to a better agreement between computed values of P and the seismic observations when using the NKS than the EKS. However, the present results from a 3D simulation strongly suggest that the EKS is a better model for the solar turbulent kinetic energy spectrum.

The better agreement obtained with NKS than EKS when adjusting the free parameters is due to the fact that the NKS concentrates most of the kinetic energy in the vicinity of k₀^E. Indeed, the NKS predominantly excites the modes whose period are close to the characteristic lifetime of the eddies of wavenumber k₀^E, i.e. the modes with frequency close to the frequency at which P peaks ( $\nu \sim 3$ mHz). As a consequence, the amount of energy going into the high frequency modes is relatively smaller with the NKS than it is with the EKS. This explains why the NKS reproduces better the steep decrease with $\nu$ of P at high frequency and results in a value for k₀^E identical to that inferred from the simulation ( $k_0^{E}\simeq 3.6~ {\rm Mm}^{-1}$ ). In contrast, whatever the adjustment, the EKS reproduces neither the frequency dependence of P at high frequency nor the value $k_0^{E}\simeq 3.6~ {\rm Mm}^{-1}$ . Hence assuming that the 3D simulation yields the proper behavior of the solar kinetic energy spectrum, well modelled by the EKS, one is led to conclude that the excitation as given by the present stochastic excitation model is not efficient enough at large scales ( $k \sim k_0$ ) and too efficient at small scales (k > k₀).

Discrepancies between our calculations and the observed excication rates or the results by Stein & Nordlund (2001) are likely due to dynamic properties of turbulence which are not properly taken into account in the excitation model. Indeed, the dynamic properties of turbulence are modeled by the function $\chi_k$ . All current theoretical calculations of the excitation rates assume a Gaussian function for $\chi_k$ (e.g. Goldreich & Keeley 1977; Balmforth 1992a). The Gaussian model is likely to be at the origin of the current under-estimate of the rates at which solar p-modes are excited (see forthcoming paper Samadi et al. 2003). This may also explain the fact that we find that the entropy source term is dominant over the Reynolds stress contribution whereas Stein & Nordlund (2001) in their direct computations found the reverse. In a recent study based on a frequency analysis of the present simulation we investigate what model can correctly reproduce model $\chi_k$ in the frequency range where the acoustic energy injected into the solar p-modes is important (see forthcoming paper Samadi et al. 2003).

In the manner of Rosenthal et al. (1999) constraints from 3D simulation can be imposed to the 1D model. According to the authors, such constraints result in a better agreement between the observed frequencies of the solar p-modes and the eigenfrequencies of the computed adiabatic oscillations. An improvement in the calculation of the excitation rates at solar-type oscillations could then also come from a more consistent calculation of the eigenmodes which would use such constrained 1D model.

Contents

Numerical constraints on the model of stochastic excitation of solar-type oscillations