Stellar turbulence and mode physics

A review on this topic was recently given by Appourchaux et al. (2009). Therefore I shall summarize only the most important matters. According to Eq. (1) we need to model the (linear) damping rate $\eta$ and the energy supply rate $P$ for estimating the mode height $H$. In this paper I shall address only the modelling of the energy supply rate. One of the problems that one faces in deriving a stochastic excitation model is the separation of the total fluid motion into (large-scale) oscillatory motion (oscillation modes), with displacement $\boldsymbol{\xi}(\mathchoice{\mbox{\boldmath$\displaystyle r$}}{\mbox{\boldmath$\textstyle r$}}{\mbox{\boldmath$\scriptstyle r$}}{\mbox{\boldmath$\scriptscriptstyle r$}},t)$, and small-scale turbulent motion with the convective velocity field $\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}(\mathchoice{\mbox{\boldmath$\displaystyle r$}}{\mbox{\boldmath$\textstyle r$}}{\mbox{\boldmath$\scriptstyle r$}}{\mbox{\boldmath$\scriptscriptstyle r$}},t)=(u,v,w)$. This separation is superficially straightforward for radial p modes (Gough, 1977a); however, for nonradial modes this separation is a much more difficult task and can possibly not be accomplished for the very-high-degree p modes and also very-high-order g modes (see e.g. Appourchaux et al., 2009). For simplicity we shall discuss here only radial p modes.

Following the basic principle of Lighthill’s acoustic analogy (Lighthill, 1952, see also Crighton 1975) one typically arranges the linearized global mode variables on the left and the nonlinear fluctuating terms associated with the turbulent convection on the right of the equations of motion, which may then be written as

$$\rho\left(\frac{\partial^{2}\boldsymbol{\xi}}{\partial t^{2}}+2\eta\frac{\partial\boldsymbol{\xi}}{\partial t}+\mathcal{L}\boldsymbol{\xi}\right)=\boldsymbol{\mathcal{F}}(\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}),$$

(7)

in which the linear spatial differential operator $\cal L$ satisfies the (unforced) homogeneous wave equation ${\cal L}\boldsymbol{\xi}_{i}=\omega_{i}^{2}\boldsymbol{\xi}$ for the mode eigenfunctions $\boldsymbol{\xi}_{i}$ with frequencies $\omega_{i}$, and $\boldsymbol{\mathcal{F}}$ is the nonlinear stochastic forcing and damping term that depends only on the turbulent velocity field $\textstyle u$. Here we consider only the Reynolds stress driving term, which dominates over the entropy driving term in the Sun and in most solar-type stars (Balmforth, 1992b; Stein & Nordlund, 2001; Stein et al., 2004; Belkacem et al., 2006; Samadi et al., 2007). It should, however, be noted that there is partial cancellation between the fluctuating Reynolds stress and entropy source terms, which could lead to smaller oscillation amplitudes than adopting the Reynolds stress contribution alone (Osaki, 1990; Stein et al., 2004; Houdek, 2006). Because the mode amplitudes of solar-type oscillations are small it is additionally assumed that they do not interact with the nonlinear right hand side of the wave equation, and consequently one can solve Eq. (7) by a nonsingular perturbation method (Goldreich & Keeley, 1977; Balmforth, 1992b; Samadi & Goupil, 2001; Chaplin et al., 2005). For radial modes only the vertical component $\mathcal{F}_{r}$ of the fluctuating Reynolds stress source term $\boldsymbol{\mathcal{F}}$ is important,

$$\mathcal{F}_{r}(w)=\partial_{r}(\rho ww-\langle\rho ww\rangle)\,,$$

(8)

which depends on the $(r,r)$ component of a two-point correlation function, $R_{rr}:=\langle ww\rangle$, where angular brackets denote an ensemble average. For incompressible, homogeneous isotropic turbulence the Fourier transform $\hat{R}_{rr}$ of $R_{rr}$ is proportional to the turbulent energy spectrum function $E(k,\omega)$, i.e.

$$\hat{R}_{rr}\propto k^{-2}E(k,\omega)$$

(9)

(Batchelor, 1953), where $k$ is a wavenumber.

Following Stein (1967) we factorize the energy spectrum function into $E(k,\omega)=E(k\ell/\pi)\Omega(\omega,\tau_{k};r)$, where

$$\tau_{k}=\lambda/ku_{k}$$

(10)

is the correlation time-scale of eddies with vertical size $\ell=\pi/k$ and velocity $u_{k}$; the correlation factor $\lambda$ is of order unity and accounts for uncertainties in defining $\tau_{k}$. For statistically stationary turbulence the energy supply rate is then given by (see e.g., Chaplin et al. 2005 for details)

$$P=\frac{\pi}{9I}\int_{0}^{R}\ell^{3}\left(\Phi\Psi rp_{\rm t}\frac{\partial{\xi}_{ri}}{\partial r}\right)^{2}{\cal S}(\omega_{i};r)\,{\rm d}r\,,$$

(11)

with

$${\cal S}(\omega_{i};r)=\int_{0}^{\infty}\kappa^{-2}\tilde{E}^{2}(\kappa)\tilde{\Omega}(\omega_{i},\tau_{k};r)\,{\rm d}\kappa\,,$$

(12)

where $p_{\rm t}=\langle\rho ww\rangle$ is the $(r,r)$ component of the (mean) Reynolds stress tensor, $\kappa=k\ell/\pi$, $R$ is surface radius, $\xi_{ri}$ is the normalized radial part of $\boldsymbol{\xi}_{i}$, and the product of $\Phi$ and $\Psi$ is a factor of unity accounting for the anisotropy of the turbulent velocity field. The spectral function ${\cal S}$ accounts for contributions to $P$ from the small-scale turbulence. For the normalized spatial turbulence energy spectrum $\tilde{E}(\kappa)$ it has been common to adopt, for example, the Kolmogorov (1941) spectrum. For the frequency-dependent factor $\Omega(\omega,\tau_{k};r)$, however, which is used for evaluating the self-convolution $\tilde{\Omega}(\omega_{i},\tau_{k};r)=\int\Omega(\omega,\tau_{k};r)\Omega(\omega_{i}\!-\!\omega,\tau_{k};r)\,{\rm d}\omega$ in Eq. (12), no satisfactory theory exists. Various functional forms were proposed in the past, which lead, however, to rather different values for the modelled oscillation amplitudes $H$ (Samadi et al., 2003; Chaplin et al., 2005; Appourchaux et al., 2009). We shall therefore discuss the most commonly adopted frequency-dependent factors for stellar turbulence spectra in the next section.

Differences in the (properly normalized) spatial spectrum $\tilde{E}(\kappa)$ create only minor differences in the energy supply rate $P$ (e.g., Balmforth, 1992b; Samadi & Goupil, 2001; Chaplin et al., 2005). But differences in the frequency-dependent (temporal) spectrum $\Omega(\omega,\tau_{k};r)$ would create rather large differences in the p-mode and also g-mode amplitudes (Samadi et al., 2003; Chaplin et al., 2005; Belkacem et al., 2009; Appourchaux et al., 2009). Those differences are produced by the convection-oscillation interactions in the deeper convectively unstable stellar layers, which are off resonance, and whose magnitude depends crucially on the adopted frequency dependence of the turbulence spectrum in the high-frequency tail of the cascade. We consider the following factors:

Fig. 3 compares the dimensionless quantity $\Omega(\omega,\tau_{k};r)\tau^{-1}_{k}$ for the three frequency factors (13)–(15). All three factors are normalized according to $\int\Omega(\omega,\tau_{k};r)\,{\rm d}\omega=1$, the integral being from $-\infty$ to $\infty$. Furthermore, the full-width at half-maximum, i.e. $2\sqrt{2\ln 2}\tau^{-1}_{k}$, has been chosen to be the same for all three factors, in order to have a meaningful comparison. The Gaussian frequency factor decreases more rapidly with $\omega\tau_{k}$, than both the Exponential and Lorentzian factors, i.e. the Gaussian factor produces the smallest magnitude in the high-frequency tail of the spectrum. Convective eddies situated in the deeper convectively unstable stellar layers have longer characteristic time scales $\tau_{k}$, i.e. $\tau_{k}$ [see Eq. (10)] increases with depth in the star (see for example Fig. 5 of Chaplin et al., 2005). The Lorentzian factor decreases more slowly with depth at constant frequency and consequently a larger fraction of the integrands of Eqs (11) and (12) arises from large off-resonant eddies situated deep in the star, whereas the Gaussian frequency factor gives less weight to the large off-resonance eddies. As a result of this the modelled oscillation amplitudes are larger with a Lorentzian time-correlation function and therefore in better agreement with the observations than with a Gaussian (Samadi et al., 2003, 2007). However, Chaplin et al. (2005) reported that the Lorentzian time-correlation function leads to overestimated heights $H$ at low frequencies for solar p modes. Also for solar g modes Belkacem et al. (2009) concluded that the best fit for reproducing the observations was obtained with a Gaussian factor at the lowest and a Lorentzian at the highest frequencies. Moreover, Belkacem et al. (2009) reported that the correlation parameter $\lambda$ [see Eq. (10)] has to increase with stellar depth in order to reproduce the 3D numerical simulations by Miesch et al. (2008), a result that is consistent with the findings by Chaplin et al. (2005). The effect of $\lambda$ on the time-correlation functions is demonstrated in Fig. 4 for a Lorentzian factor. Increasing $\lambda$ leads to a more rapidly declining tail and consequently to a smaller contribution to the energy supply rate $P$ from the off-resonant eddies situated in the deep convection zone. Finding the correct time-correlation function for stellar turbulence remains an open issue.

From the theoretical point of view a Lorentzian factor is a result predicted for the largest, most-energetic eddies by the time-dependent mixing-length formulation of Gough (1977a). Moreover, in agreement with Kolmogorov (1941) theory Sawford (1991) demonstrated that in the limit of high-Reynolds-number turbulence the two-point velocity autocorrelation, normalized by the velocity variance $\sigma^{2}:=\langle ww\rangle-\langle w\rangle\langle w\rangle$, is

$$\tilde{R}(\tau)=\langle w(t)w(t+\tau)\rangle/\sigma^{2}={\rm e}^{-|\tau|/T_{\rm L}}\,,$$

(16)

where

$$T_{\rm L}=\int_{0}^{t}\tilde{R}(\tau)\,{\rm d}\tau$$

(17)

is the Lagrangian integral time scale (see also Legg & Raupach, 1982). Consequently the normalized velocity autocorrelation $\tilde{R}$ is proportional to an exponential decrease with separation time $\tau$ in the inertial range, and the velocity spectrum becomes with Eq. (16)

	$\displaystyle F(\omega)$	$\displaystyle=$	$\displaystyle\int_{0}^{\infty}\tilde{R}(\tau)\cos(\omega\tau)\,{\rm d}\tau$		(18)
		$\displaystyle=$	$\displaystyle\frac{T_{\rm L}}{1+(\omega T_{\rm L})^{2}}\,,$

i.e. it is of Lorentzian shape with $F(\omega)\propto\omega^{-2}$ for $\omega\rightarrow\infty$ and at infinite Reynolds number $R_{\rm e}$. Sawford extended this analysis for arbitrary values of the Reynolds number $R_{\rm e}:=(t_{\rm E}/t_{\eta})^{2}$, where $t_{\rm E}$ and $t_{\eta}$ are the time scales of the energy-containing eddies and the Kolmogorov timescale respectively (Tennekes & Lumley, 1972). For this case the velocity spectrum is (Sawford, 1991)

$$F(\omega)=\frac{(1+R_{\rm e}^{1/2})R_{\rm e}^{1/2}T_{\rm L}}{(R_{\rm e}^{1/2}-\omega^{2}T^{2}_{\rm L})^{2}+(1+R_{\rm e}^{1/2})^{2}\omega^{2}T^{2}_{\rm L}}\,.$$

(19)

Fig. 5 displays the spectrum $F(\omega)$ of the normalized velocity autocorrelation $\tilde{R}$ for two different values of the Reynolds number, $R_{\rm e}=(100,\infty)$. It is interesting to note that there is essentially no inertial range with a frequency decay of $\omega^{-2}$, i.e. a Lorentzian frequency dependence [Eq. (18)], for $R_{\rm e}=100$; the spectrum decreases as $\omega^{-4}$ in the dissipation subrange, i.e. for turbulent scales close to the Kolmogorov dissipation scale.

These findings are supported by laboratory experiments. Fig. 6 displays the measured velocity power spectrum in a rotating shear turbulence experiment (Kàrmàn swirling flow) with water. The solid curve displays the measurements and the dot-dashed curve is a fit of expression (19) to the data with $R^{1/2}_{\rm e}=T_{\rm L}/T_{2}$, where $T_{\rm L}\simeq T_{\rm E}$ and $T_{2}\simeq t_{\eta}$. There is a well defined inertial range for frequencies between $(2\pi T_{\rm L})^{-1}$ and $(2\pi T_{2})^{-1}$ with a decay of about $(\omega/2\pi)^{-2}$, i.e. a Lorentzian spectrum. The high-frequency tail, however, decays more rapidly. The measurements are superficially similar to the second-order model by Sawford (1991) depicted in Fig. 5. It should, however, be noted that the molecular Prandtl number of about 6.8 in the water experiment is several magnitudes larger than the molecular Prandtl number in stellar convection zones, which is typically of order $10^{-6}$.

Georgobiani et al. (2004) used, as did Samadi et al. (2003), the 3D Large-Eddy-Simulations by Stein & Nordlund (2001) to investigate the time-correlation function $\Omega(\omega,\tau_{k};r)$ in the Sun. Fig. 7 compares the power spectrum (squared fast Fourier transform), $|w(\nu)|^{2}$, of the vertical component of the turbulent velocity field, ${\boldsymbol{u}}=(u,v,w)$, at fixed horizontal wavenumber $k$ and depth, between the simulation results (solid curve) and various analytical time-correlation functions [Exponential (EF), Gaussian (GF) and Lorentzian (LF) frequency factors]. The results are shown for a wavenumber $k=4\,$Mm^-1 computed 250 km below the solar surface. Neither of the three analytical functions fit the simulation results, particularly in the high-frequency tail of the spectrum. Moreover, Georgobiani et al. (2004) reported that the functional form (frequency dependence) of the simulation results changes substantially with wavenumber and depth, and consequently that the turbulent energy spectrum function $E(k,\nu)$ is not separable into wavenumber $E(k)$ and frequency $\Omega(\nu,\tau_{k};r)$. The authors suggest using instead the following empirical function for the temporal part of the vertical velocity power spectrum

$$|\hat{w}(\nu)|^{2}\propto\frac{1}{\left(1+\Lambda^{2}\nu^{2}\right)^{\beta}}\,,$$

(20)

where the two coefficients $\Lambda$ and $\beta$ are determined from fitting Eq. (20) to the 3D simulation results. Note that for $\beta=1$ a Lorentzian frequency factor is recovered. Results from solar simulations are given in Fig. 8 for various horizontal wavenumbers and depths. In agreement with Chaplin et al. (2005) and Belkacem et al. (2009) the coefficient $\Lambda\propto\lambda$ increases with depth, i.e. the width of the frequency factor decreases with depth (see also Fig. 4), particularly at smaller scales (i.e. at large wavenumbers). Except near the surface the values of the exponent $\beta$ are rather larger than unity.

Numerical simulations of stellar convection resolve only the largest scales of the turbulent motion. The smaller scales still need to be approximated by a so-called sub-grid model. Such simulations are called Large-Eddy-Simulations. The total number of scales that are numerically resolved, i.e. the ratio of the largest scale $l_{\rm L}$ to the smaller scale $l_{\eta}$, is related to the Reynolds number by (e.g. Tennekes & Lumley, 1972) ${l_{\rm L}}/{l_{\eta}}\sim R^{3/4}_{\rm e}\,,$ and consequently the total number $N$ of grid points in the simulation is $N=({l_{\rm L}}/{l_{\eta}})^{3}\sim R^{9/4}_{\rm e}\,.$ A solar Reynolds number of $R_{\rm e}\simeq 10^{12}$ would require a total meshpoint number of $N\simeq 10^{27}$. With today’s super computers the maximum achievable number of meshpoints is $N_{\rm max}\simeq 10^{12}$, and is therefore some 15 magnitudes too small for what is required to resolve all the turbulent scales of solar convection. Consequently a sub-grid model is necessary for describing the dynamics of the numerically unresolved smaller scales of the turbulent cascade. Various models are available. The most commonly used models are hyperviscosity models and the Smagorinsky model. All sub-grid models assume that turbulent transport is a diffusive process. Hyperviscosity models, for example, use higher derivatives for the diffusion operator in the momentum equation, thereby extending the inertial range, which also leads to a better representation of the dynamics of the larger scales. An overview of sub-grid models was recently presented by Miesch (2005).
An obvious question to ask is what are the effects of using different sub-grid models on the mode properties in simulations of solar-type stars? This question was recently addressed, in part, by Jacoutot et al. (2008), who compared simulated solar energy supply rates $P$ using various sub-grid models. Their results are summarized in Fig. 9. Agreement with observations is generally satisfactory except for the classical Smagorinsky model.

It is also important to note that the Prandtl number in 3D simulations is currently about 0.01 $-$ 0.25 (e.g. Miesch et al., 2008). It is therefore substantially larger than the Prandtl number in the Sun and in solar-type stars. Although numerical simulations provide an important input to our understanding of stellar convection and are very important for calibrating semi-analytical convection models, we must remain aware of the shortcomings of the currently used 3D numerical simulations.