On the computation of eigenfrequencies for equilibrium models including turbulent pressure

TL;DR

This paper investigates how to accurately include turbulent pressure effects in the computation of stellar oscillations, improving agreement with observations and proposing a formalism for broader application.

Contribution

It evaluates the impact of turbulent pressure perturbations on stellar oscillation frequencies and compares different modeling approaches, proposing a new formalism for better accuracy.

Findings

TDC formalism improves frequency agreement with observations.

Perturbation of density and advection terms significantly affect frequency shifts.

Proposed formalism allows extrapolation of results to other stars.

Abstract

The space-borne missions have provided a wealth of highly accurate data. However, our inability to properly model the upper-most region of solar-like stars prevents us from making the best of these observations. This problem is called "surface effect" and a key ingredient to solve it is turbulent pressure for the computation of both the equilibrium models and the oscillations. While 3D hydrodynamic simulations help to include properly the turbulent pressure in the equilibrium models, the way this surface effect is included in the computation of stellar oscillations is still subject to uncertainties. We aim at determining how to properly include the effect of turbulent pressure and its Lagrangian perturbation in the adiabatic computation of the oscillations. We also discuss the validity of the gas-gamma model (GGM) and reduced gamma model (RGM) approximations, which have been used to compute adiabatic oscillations of equilibrium models including turbulent pressure. We use a patched model of the Sun with an inner part constructed by a 1D stellar evolution code (CESTAM) and an outer part by the 3D hydrodynamical code (CO$^5$BOLD). Then, the adiabatic oscillations are computed using the ADIPLS code for the GGM and RGM and with the MAD code imposing the adiabatic condition on an existing time-dependent convection (TDC) formalism. We show that the computation of the oscillations using the TDC formalism in the adiabatic limit improves significantly the agreement with the observed frequencies compared to the GGM and RGM. Of the components of the turbulent pressure perturbation, the perturbation of the density and advection term is found to contribute most to the frequency shift. We propose a formalism to evaluate the frequency shift due to the inclusion of the term with the turbulent pressure perturbation in the variational principle in order to extrapolate our result to other stars.

Tables

Table 1. Table 1: Characteristics of the patched model (PM)

$T_{\mathrm{eff}}$ [K]	$\log g$ [g cm-2]	$T_b$ [K]	$p_b$ [g cm-1 s-2]	$M$ [$M_{\odot }$]	Age [Gyr]	$\alpha$
5775	4.44	$1.53\times {10}^4$	$3.66\times {10}^6$	1.01	4.61	1.65

Equations

$$\frac{\delta p_{\mathrm{turb}}}{p_{\mathrm{turb}}}\simeq \frac{\delta p_{\mathrm{tot}}}{p_{\mathrm{tot}}}\simeq \frac{\delta p_{\mathrm{th}}}{p_{\mathrm{th}}}={\Gamma }_1\frac{\delta \rho }{\rho },$$

$$\frac{\delta p_{\mathrm{turb}}}{p_{\mathrm{turb}}}=0.$$

$$\frac{\delta p_{\mathrm{tot}}}{p_{\mathrm{tot}}}=\frac{p_{\mathrm{th}}}{p_{\mathrm{tot}}}\frac{\delta p_{\mathrm{th}}}{p_{\mathrm{th}}}={\Gamma }_1^r\frac{\delta \rho }{\rho }.$$

$$\frac{\delta p_{\mathrm{turb},\mathrm{l}}}{p_{\mathrm{turb},\mathrm{l}}}=\frac{\delta \rho }{\rho }+2\frac{\overline{V_r\delta V_r}}{\overline{V_r^2}},$$

$$\frac{\overline{V_r\delta V_r}}{\overline{V_r^2}}$$

$$B$$

$$C$$

$$D$$

$$A=\frac{\overline{\rho V_r^2}}{\overline{\rho (V_{\theta }^2+V_{\varphi }^2)}},$$

$$p_{\mathrm{turb},\mathrm{nl}}({\zeta }_0)={\int }_{-\infty }^{+\infty }{p}_{\mathrm{turb},\mathrm{l}}{\mathrm{e}}^{-b|\zeta -{\zeta }_0|}\mathrm{d}\zeta ,$$

$$F_{\mathrm{c},\mathrm{nl}}({\zeta }_0)={\int }_{-\infty }^{+\infty }{F}_{\mathrm{c},\mathrm{l}}{\mathrm{e}}^{-a|\zeta -{\zeta }_0|}\mathrm{d}\zeta ,$$

$${\mathrm{d}}^2{p}_{\mathrm{turb},\mathrm{nl}}/\mathrm{d}{\zeta }^2=b^2(p_{\mathrm{turb},\mathrm{nl}}-p_{\mathrm{turb},\mathrm{l}}),$$

$${\mathrm{d}}^2{F}_{\mathrm{c},\mathrm{nl}}/\mathrm{d}{\zeta }^2=a^2(F_{\mathrm{c},\mathrm{nl}}-F_{\mathrm{c},\mathrm{l}}).$$

$${\nu }^2$$

$$N_{\mathrm{turb}}(m)=\frac{1}{8{\pi }^2\nu }{\left({\int }_0^M|\xi {|}^2\mathrm{d}{m}^{\prime }\right)}^{-1}{\int }_0^m\mathrm{Re}\left[\frac{\delta {\rho }^{\ast }}{\rho }\frac{\delta p_{\mathrm{turb}}}{\rho }\right]\mathrm{d}{m}^{\prime },$$

$$\frac{\delta p_{\mathrm{turb},\mathrm{l}}}{p_{\mathrm{tot}}}=\Pi \frac{\delta p_{\mathrm{th}}}{p_{\mathrm{tot}}}+\Xi \frac{\xi }{R}$$

$$\Pi =\mathcal{G}\left[\frac{1}{{\Gamma }_1}+\mathcal{F}({\Pi }_1+{\Pi }_2+{\Pi }_3+{\Pi }_4)\right],$$

$$\Xi =\mathcal{FH}[{\Xi }_1+{\Xi }_2+{\Xi }_3+{\Xi }_4],$$

$$\mathcal{F}=\frac{2}{B+[(i\Omega +\beta )\sigma {\tau }_c+1]D},$$

$$\mathcal{G}=\frac{p_{\mathrm{turb},\mathrm{l}}}{p_{\mathrm{th}}}{\left[1-\frac{\mathcal{F}}{1+(\sigma {\tau }_c{)}^2}\frac{p_{\mathrm{turb},\mathrm{l}}}{p_{\mathrm{tot}}}(2{\omega }_R{\tau }_{c}D+1)(D+1)\right]}^{-1},$$

$$\mathcal{H}$$

$$N_{\mathrm{turb}}(m)$$

$$\mathcal{F}\to \frac{\Omega \Lambda }{i\sigma {\tau }_c},\mathcal{G}\to \frac{p_{\mathrm{turb},\mathrm{l}}}{p_{\mathrm{th}}},$$

$$\frac{\delta p_{\mathrm{turb},\mathrm{l}}}{p_{\mathrm{tot}}}\to \frac{p_{\mathrm{turb},\mathrm{l}}}{p_{\mathrm{th}}}\frac{1}{{\Gamma }_1}\left(1+\frac{2A}{A+1}\right)\frac{\delta p_{\mathrm{th}}}{p_{\mathrm{tot}}}\sim \frac{p_{\mathrm{turb},\mathrm{l}}}{p_{\mathrm{th}}}\frac{\delta p_{\mathrm{th}}}{p_{\mathrm{tot}}},$$

$$N_0(m)=\frac{1}{8{\pi }^2\nu }{\left({\int }_0^M|\xi {|}^2\mathrm{d}{m}^{\prime }\right)}^{-1}{\int }_0^m\frac{p_{\mathrm{th}}}{{\Gamma }_1\rho }{\left|\frac{\delta p_{\mathrm{th}}}{p_{\mathrm{th}}}\right|}^2\mathrm{Re}\left(\frac{\mathcal{G}}{{\mathrm{\Gamma }}_1}\right)\mathrm{d}{m}^{\prime }$$

$$N_i(m)$$

$$\Delta {\nu }_{\mathrm{turb}}$$

$$\nabla \cdot \overrightarrow{V}=0,$$

$$\overline{\rho }\frac{\mathrm{d}\overrightarrow{V}}{\mathrm{d}t}=\frac{\Delta \rho }{\overline{\rho }}\nabla \overline{p_{\mathrm{tot}}}-\nabla \Delta p_{\mathrm{tot}}-\rho \overrightarrow{V}\cdot \nabla \overrightarrow{u}-\Lambda \frac{\overline{\rho }\overrightarrow{V}}{{\tau }_c},$$

$$\frac{\Delta (\rho T)}{\overline{\rho T}}\frac{\mathrm{d}\overline{s}}{\mathrm{d}t}+\overrightarrow{V}\cdot \nabla \overline{s}=-\frac{{\omega }_R{\tau }_c+1}{{\tau }_c}\Delta s,$$

$$\Lambda \frac{\overline{\rho }\overrightarrow{V}}{{\tau }_c}$$

$$\overline{\rho }\overline{T}\frac{\Delta s}{{\tau }_c}$$

$$\nabla \cdot \Delta {\overrightarrow{F}}_{\mathrm{R}}$$

$$l$$

$$\Gamma (\Gamma +1)=\mathcal{A}(\nabla -{\nabla }_{\mathrm{ad}}),$$

Full text

*Aims. *We aim at determining how to properly include the effect of turbulent pressure and its Lagrangian perturbation in the adiabatic computation of the oscillations. We also discuss the validity of the gas-gamma model and reduced gamma model approximations, which have been used to compute adiabatic oscillations of equilibrium models including turbulent pressure.

*Methods. *We use a patched model of the Sun with an inner part constructed by a 1D stellar evolution code (CESTAM) and an outer part by the 3D hydrodynamical code (CO5BOLD). Then, the adiabatic oscillations are computed using the ADIPLS code for the gas-gamma and reduced gamma model approximations and with the MAD code imposing the adiabatic condition on an existing time-dependent convection formalism. Finally, all those results are compared to the observed solar frequencies.

*Results. *We show that the computation of the oscillations using the time-dependent convection formalism in the adiabatic limit improves significantly the agreement with the observed frequencies compared to the gas-gamma and reduced gamma model approximations. Of the components of the perturbation of the turbulent pressure, the perturbation of the density and advection term is found to contribute most to the frequency shift.

*Conclusions. *The turbulent pressure is certainly the dominant factor responsible for the surface effects. Its inclusion into the equilibrium models is thus necessary but not sufficient. Indeed, the perturbation of the turbulent pressure must be properly taken into account for computing adiabatic oscillation frequencies. We propose a formalism to evaluate the frequency shift due to the inclusion of the term with the turbulent pressure perturbation in the variational principle in order to extrapolate our result to other stars at various evolutionary stages. Although this work is limited to adiabatic oscillations and the inclusion of the turbulent pressure, future works will have to account for the nonadiabatic effect and convective backwarming.

††offprints: takafumi.sonoi@obspm.fr

1 Introduction

As shown by the space missions CoRoT (Baglin et al. 2006b, a; Michel et al. 2008) and Kepler (Borucki et al. 2010; Bedding et al. 2010; Chaplin et al. 2011), solar-like oscillations are ubiquitous to low-mass stars from the main-sequence to the red giant branch. They have been widely used to infer the internal structure of those stars and have permitted us to dramatically improve our knowledge as well as to put stringent constraints on stellar structure and evolution (e.g. Chaplin & Miglio 2013).

However, there are still some fundamental difficulties to overcome so as to exploit the full potential of the asteroseismic observations. Surface effects are likely to be the most important. This generic term names the systematic differences between the observed and model frequencies due to our deficient physical description of the upper-most layers of solar-like stars (e.g. Christensen-Dalsgaard 2016). One of the key ingredients of those surface effects is the turbulent pressure. In standard models of stellar equilibrium structure and oscillations, it is generally neglected because its modelling is difficult. Nevertheless, it is a key factor to obtain accurate frequencies of stellar models and particularly for $p$ modes that are very sensitive to the surface layers. The crucial role of turbulent pressure in computing stellar oscillations has been emphasized in many studies (Brown 1984; Zhugzhda & Stix 1994; Schlattl et al. 1997; Petrovay et al. 2007; Houdek 2010). More recently, analyses of surface effects have been carried out using 3D hydrodynamical models (Stein & Nordlund 1991; Rosenthal et al. 1995, 1999; Yang & Li 2007; Piau et al. 2014; Bhattacharya et al. 2015; Sonoi et al. 2015; Magic & Weiss 2016; Ball et al. 2016; Houdek et al. 2017; Trampedach et al. 2017), because these models provide a realistic description of the equilibrium structure including the turbulent pressure.

However, when turbulent pressure is included in the equilibrium model, the computation of the related stellar oscillations becomes tricky and we have to care about the possible inconsistency between the oscillation formalism and the equilibrium models. To consider this problem, Rosenthal et al. (1995, 1999) proposed two approximations, the gas-gamma model (GGM), for which the Lagrangian perturbation of the turbulent pressure equals to the perturbation of the gas pressure, and the reduced gamma model (RGM), for which the Lagrangian perturbation of the turbulent pressure vanishes. They have shown that the GGM frequencies better reproduce the observed frequencies compared to ones obtained with the RGM assumption. Their result implies that it is important to take the perturbation of the turbulent pressure into account in order to obtain accurate frequencies. However, the GGM assumption does not rely on a convincing principle and deserves more investigation.

In this work, we consider the computation of adiabatic oscillations for an equilibrium model including turbulent pressure. To do so, we use a time-dependent convection (TDC) formalism that enables us to account for the perturbation of turbulent pressure. We also discuss the validity of the GGM approximation. We use the TDC formalism developed by Grigahcène et al. (2005), which originates from the work of Unno (1967) and was generalized for nonradial oscillations by Gabriel et al. (1975). This formalism has been so far adopted for the computation of the full nonadiabatic oscillations in order to explain the excitation of the classical pulsators (e.g. Dupret et al. 2005, 2008), or to fit to the damping rates of the solar-like oscillations (e.g. Dupret et al. 2006a; Belkacem et al. 2012; Grosjean et al. 2014). Dupret et al. (2006b) developed it for treating the nonlocal convection. For our purpose, we will impose the adiabatic condition on this formalism to see the validity of the GGM approximation. Moreover, such an approach allows us to consider the effect of the turbulent pressure separately from the nonadiabatic effect, which is also expected to affect eigenfrequencies (Houdek 2010).

The paper is organized as follows: Section 2 introduces how to compute eigenfrequencies with turbulent pressure. Section 3 discusses the dominant causes of the frequency shift due to the perturbation of the turbulent pressure. Section 4 gives the conclusion.

2 Modelling eigenfrequencies with turbulent pressure

2.1 Equilibrium model

We use the solar “patched” model (PM) described in Samadi et al. (2007) and Sonoi et al. (2015). The inner part of this model was constructed using the 1D stellar evolution code CESTAM (Marques et al. 2013) while the near-surface layers have been obtained using temporal and horizontal averages of the 3D hydrodynamical simulation by the CO5BOLD (Freytag et al. 2012) code with the CIFIST grid (Ludwig et al. 2009). The turbulent pressure is thus included only in the 3D upper layers. The matching between the inner and outer layers have been computed through an optimization of the 1D model with a Lenvenberg-Marquardt algorithm. The constraints for the optimization are the effective temperature ($T_{\rm eff}$) of the 3D model, the gravity acceleration at the photosphere ($g$), and the temperature at the bottom of the 3D model ($T_{\rm b}$). As for the last one, the temperature at the level having the same total pressure with the bottom of the 3D model ($p_{\rm b}$) is matched with $T_{\rm b}$. The free parameters are the stellar age, total mass ($M$), and mixing length parameter ($\alpha$). The resulting values are provided in Table 1. Our matching point is deep enough since, at the bottom of our 3D model, the fraction of the turbulent pressure to the total pressure is small enough ($\simeq 0.014$) that it does not affect frequencies of acoustic modes, of which amplitude is concentrated in the upper layers.

PM is constructed by replacing the outer part of the optimized 1D model, which we call “unpatched” model (UPM), with the averaged 3D model. The additional support by turbulent pressure modifies the hydrostatic equilibrium so that, at the photosphere, the radius of PM is larger than UPM by about 0.02% (see also Table 2 and Fig. 2 in Sonoi et al. 2015).

2.2 Computation of adiabatic oscillations: the gas-gamma and reduced gamma approximations

Following the work of Rosenthal et al. (1995, 1999), two approximations can be adopted to account for the turbulent pressure in the equilibrium model, namely the gas-gamma model (hereafter GGM) and the reduced gamma model (hereafter RGM).

The GGM assumes that the relative Lagrangian perturbation of turbulent pressure equals the relative Lagrangian perturbation of thermal pressure, which is the sum of gas and radiation pressures, and hence is equal to that of the total pressure,

[TABLE]

where $\delta$ denotes the Lagrangian perturbation, $p_{\rm turb}$ is the turbulent pressure, $p_{\rm th}$ is the thermal pressure, and $p_{\rm tot}(=p_{\rm th}+p_{\rm turb})$ is the total pressure.

The RGM approximation introduces the reduced $\Gamma_{1}$, defined as $\Gamma^{r}_{1}\equiv(p_{\rm th}/p_{\rm tot})\Gamma_{1}$. In this approximation, the Lagrangian perturbation of turbulent pressure is neglected:

[TABLE]

We have thus

[TABLE]

Figure 1 shows the adiabatic exponent in UPM and PM. We note that there is no difference between GGM and RGM for UPM, which does not include any turbulent pressure. In the 3D layer of PM, $\Gamma^{r}_{1}$ has lower values than $\Gamma_{1}$ as a consequence of the presence of turbulent pressure.

Figure 2 shows the difference between the observed frequencies as given by Broomhall et al. (2009) and the computed frequencies obtained using the PM and UPM described in Sect. 2.1. The larger radius of PM makes frequencies lower than for UPM. Then, the value of $(\nu_{\rm obs}-\nu_{\rm model})$ is higher for PM. On the other hand, the GGM treatment, namely including the perturbation of the turbulent pressure, oppositely increases the frequencies, and reduces the deviation of the RGM frequencies from the observed frequencies. As a result, the GGM frequencies are in better agreement with the observation than the RGM ones. This result confirms the result of Rosenthal et al. (1999).

The deviation of the GGM frequencies from the observation is at most $\sim 6\mu$Hz in our analysis. This is of the same order but a little larger than those of the other analyses using the other 3D hydrodynamical models ($\sim 4\mu$Hz in Rosenthal et al. 1999 and Magic & Weiss 2016, and $\sim 3\mu$Hz in Ball et al. 2016). The deviation of the RGM frequencies is at most $\sim 10\mu$Hz, similarly to Houdek et al. (2017).

Finally, we note that the GGM and RGM approximations are easily implemented in an adiabatic oscillation code. However, the underlying assumptions are rather crude and deserve more attention. Particularly, the perturbation of the turbulent pressure should be out of phase with that of the gas pressure and density (Houdek 2000; Houdek et al. 2017). Then, computations only with the real part of the eigenfrequency such as those performed by the ADIPLS code are not valid. Therefore, one has to provide a modelling of the perturbation of the turbulent pressure and this is permitted by using a time-dependent modelling of convection (TDC) as provided in the following section. Indeed, the phase lag between the turbulent pressure and the other variables takes place in computation with TDC.

2.3 Computation of adiabatic oscillations: the TDC treatment for nonlocal convection

In the following, we adopt the TDC formalism developed by Grigahcène et al. (2005) and Dupret et al. (2006b) to compute frequencies of PM. While it is usually used to compute non-adiabatic oscillations, we consider the limit of adiabatic oscillations by setting $\delta s=0$, where $s$ is the specific entropy. Such an approach allows us to properly consider both the effects of turbulent pressure on the equilibrium structure and of its perturbation in the adiabatic limit. Moreover, this clarifies the individual contribution of turbulent pressure to the surface effects separately from the nonadiabatic effects.

As mentioned at the end of Sect. 2.2, the phase lag occurs between the perturbation of the turbulent pressure and the other variables when we adopt a TDC formalism. The phase lag leads to excitation or damping of oscillation amplitude. Namely, the eigenfrequencies of the oscillation become complex. In this work, however, we only pay attention to the real part of the eigenfrequencies, since we need nonadiabatic treatment to exactly investigate the damping rates. We are aware that the nonadiabatic effects would be important not only for the damping rates, but also for the oscillation frequencies. This is however out of the scope of the present article and will be considered in a following work.

Therefore, we start by considering the expression of the perturbation of the turbulent pressure,

[TABLE]

where $p_{\rm turb,l}$ is the turbulent pressure as obtained in the framework of a local theory of convection and $V_{r}$ is the radial component of the convective velocity. The overbar indicates averaging in the coarse grain, which is much larger than most convective eddies but much smaller than the scale of the oscillation wavelength. To go further, we consider the perturbation of the radial convective velocity (Eq. 44 or Eq. 21 in Dupret et al. 2006b) in the adiabatic limit, but for the sake of simplicity, we limit ourselves to the case of radial oscillations ($\ell=0$). This gives

[TABLE]

with

[TABLE]

where $\sigma$ $(\equiv 2\pi\nu)$ is the oscillation frequency in unit of rad s*-1*, $\tau_{c}$ is the convective timescale, $c_{p}$ is the specific heat capacity at constant pressure, $\rho$ is the density, $\xi$ is the displacement, $\omega_{R}$ is the inverse of the radiative cooling timescale of convection eddies, $l$ is the mixing length defined by Eq. (34), and $Q[\equiv-(\partial\ln\rho/\partial\ln T)_{p_{\rm th}}]$ is the volume expanding rate.

The free parameters $\beta$ and $\Omega$ are related to the closure of the TDC theory. The parameter $\beta$ is a complex value and is introduced in Eq. (39). The parameter $\Omega$ is an adjusting function introduced in the closure terms of the momentum and energy equations for the convective fluctuations (Eqs. 40 and 41). For stationary convection, it has the same meaning as in the formalism of Canuto & Mazzitelli (1991). This quantity is determined by matching with the results given by the 3D simulation using Eqs. (38b), (42) and (43).

The parameter $A$ stands for the anisotropy of the turbulence and is defined as

[TABLE]

where $V_{\theta}$ and $V_{\phi}$ are the horizontal components of the convective velocity. In this work, this parameter is obtained directly from the 3D simulation. For the layers extracted from the 1D model, we fix the value as given at the bottom of the 3D simulation. This quantity is displayed in Fig. 3.

For taking the non-locality into account, we adopt the approach of Spiegel (1963). It consists in using an analogy with radiative transfer. The local values, as given by the mixing length theory (MLT), are considered as source terms and then the nonlocal values are obtained by performing an average, that is,

[TABLE]

where ${\rm d}\zeta={\rm d}r/H_{p}$, and $a$ and $b$ are free parameters as introduced by Balmforth (1992). The temporally and horizontally averaged values of turbulent pressure and convective flux in the 3D model are substituted into $p_{\rm turb,nl}$ and $F_{\rm c,nl}$, respectively. The quantities $p_{\rm turb,l}$, $F_{\rm c,l}$ stand for their local counterparts. These equations can be recast by taking the second order derivative

[TABLE]

Equations (12) and (13) are then used to infer the values of $a$ and $b$ as well as the local values of the turbulent pressure and convective flux from the 3D numerical simulation. In the overshooting region, the two local quantities, $p_{\rm turb,l}$ and $F_{\rm c,l}$, vanish, so that $a$ and $b$ are obtained by fitting an exponential function to the turbulent pressure and convective flux as given by the 3D simulation. From our model, we get $a=6.975$ and $b=1.697$. Subsequently, the local counterparts ($p_{\rm turb,l}$ and $F_{\rm c,l}$) are easily obtained by solving Eqs. (12) and (13) in the convective region (Fig. 4 and top panel of Fig. 6). With the equations for stationary convection (Eqs. 35, 38b and 43), we can evaluate $\Omega$, $\alpha$ and $\Gamma$ as functions of the depth, where $\Gamma\,[=(\omega_{R}\tau_{c})^{-1}]$ is the convective efficiency.

The perturbations of the nonlocal turbulent pressure and convective flux ($\delta p_{\rm turb,nl}$ and $\delta F_{\rm c,nl}$) are obtained by solving the eigenvalue problem of the perturbed hydrodynamical equations of mean flow combined with the perturbed equations of (12) and (13). On the other hand, their local counterparts ($\delta p_{\rm turb,l}$ and $\delta F_{\rm c,l}$) can be evaluated with the linear combination of the eigenfunctions, which is given by Eqs. (4) and (5) for $\delta p_{\rm turb,l}$, and Eq. (45) for $\delta F_{\rm c,l}$.

2.4 Comparison of frequencies among GGM, RGM, and TDC

For the computation with the TDC formalism, we need to give the value of the free parameter $\beta$. The calibrated $\beta$ values have been of the order of unity in the previous studies (Dupret et al. 2005, 2006a, 2008; Belkacem et al. 2011; Grosjean et al. 2014). In this work, the real part were ranged from 0.2 to 2.0, while the imaginary part from $-$2.0 to 2.0 at 0.2 intervals. The top panel of Fig. 5 shows the results with the different values of $\beta$ (the black lines, at 0.4 intervals in both the real and imaginary parts for visibility). Evaluating $\chi^{2}=\sum_{n}(\nu^{n}_{\rm model}-\nu^{n}_{\rm obs})^{2}$ for each $\beta$, we found that the value of $\beta=0.2-1.2i$ gives the smallest deviation from the observed frequencies (the red line).

The bottom panel compares the case of $\beta=0.2-1.2i$ with the GGM and RGM. First, compared to the GGM, the TDC treatment improves the agreement with the observations, particularly for the intermediate radial order modes. The deviation from the observed frequencies is at most $\sim$ 4$\mu$Hz. Although our analysis is adiabatic, it provides results of the same order as the nonadiabatic analysis of Houdek et al. (2017), who used another TDC formalism (Gough 1977b, a) and PM with a 3D model of Trampedach et al. (2013) and reported $\sim 3\mu$Hz deviation from the observed frequencies.

Secondly, the GGM frequencies are closer to the TDC ones than the RGM ones are. Although Rosenthal et al. (1995, 1999) and our results in Sect. 2.2 (Fig. 2) show that the GGM reproduced the observations better than the RGM, this result implies that the GGM is superior to the RGM also from the theoretical point of view. Namely, it would be worth taking the perturbation of the turbulent pressure into account even for the adiabatic computations from both observational and theoretical viewpoints.

We note some difference between the TDC and GGM for $n\gtrsim 10$, while their frequencies are almost identical for the lower radial orders. This difference implies that the turbulent pressure perturbation is not as simple as provided by Eq. (1), and that the GGM cannot reproduce the influence of the turbulent pressure with enough precision. We discuss such effects in the following section.

3 Contribution to the frequency shift introduced by turbulent pressure perturbation

In Sect. 2.2 and Fig. 2, we have shown the frequency shift due to the elevation of the upper layer due to the turbulent pressure in the equilibrium model, comparing the PM and UPM. Here, we discuss the contribution to the frequency shift due to the perturbation of the turbulent pressure. First, we determine which region in the star contributes to the frequency shift (Sect. 3.1). As shown in Eqs. (4) and (5), the perturbation of the turbulent pressure consists in different perturbative processes. Secondly, we determine which perturbative process in convection is dominant (Sect. 3.2).

3.1 Contributing region to the frequency shift

To see the contribution of the turbulent pressure perturbation, we adopt the variational principle. Multiplying $\xi_{r}^{*}$ in both sides of the equation of movement (48), using Eqs. (49) and (50), integrating over the mass of the star and taking the real part, we obtain

[TABLE]

Except for low-order modes, the terms in the second line of Eq. (14) hardly contribute since $|\xi/r|\ll|d\xi/dr|$ and $p_{\rm turb}\ll p_{\rm th}$. Here, we discuss the term with the turbulent pressure perturbation, $\delta p_{\rm turb}$. We introduce

[TABLE]

so that the integral to the surface, $N_{\rm turb}(m=M)$, represents the frequency shift to which the turbulent pressure perturbation contributes. To be exact, this term includes some part of the effect of the upper layer elevation, which appears in the equilibrium variables, $\rho$ and ${\rm d}m$. Nevertheless, it is useful to see the contribution of the turbulent pressure perturbation. Indeed, since $p_{\rm th}\gg p_{\rm turb}$ and hence —$\delta p_{\rm th}|\gg|\delta p_{\rm turb}|$, most part of the elevation effect is included in the term with $\delta p_{\rm th}$.

The bottom panel of Fig. 6 shows the profiles of $N_{\rm turb}$ for four radial modes. Here, the perturbation of the nonlocal turbulent pressure given by the MAD code is substituted into Eq. (15). The integral $N_{\rm turb}$ increases mainly at $\log\,T\simeq$ 4.0–4.4, just below the peak of the $p_{\rm turb}/p_{\rm tot}$ ratio, shown in the top panel. By the way, it slightly increases even in the overshooting region above the boundary determined by the Schwarzschild criterion, since the nonlocal turbulent pressure contributes there. As the radial order $n$ increases, $N_{\rm turb}$ increases more substantially. The inertia gives the major contribution to this tendency. With the increasing radial order, the amplitude becomes confined in the near-surface region. Because of the low density of this region, the frequency becomes easier to shift. We discuss the dominant causes of the frequency shift in the following section.

3.2 Dominant perturbative process

In the previous section, we have confirmed that the zone just below the peak of the $p_{\rm turb}/p_{\rm tot}$ ratio dominantly contributes to the frequency shift. Here we identify the respective contribution of the different processes to the total perturbation of turbulent pressure. For this purpose, we recast the expression of the perturbation of the local turbulent pressure (Eqs. 4 and 5). More precisely, we express it as the linear combination of the thermal pressure perturbation $\delta p_{\rm th}$ and the displacement $\xi$. The detailed procedure is described in Appendix B. Then, the perturbation of the local turbulent pressure (Eq. 4) is expressed as

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where the definitions of $\Pi_{1,2,3,4}$ and $\Xi_{1,2,3,4}$ are given by Eqs. (55) to (62). The coefficients $\Pi_{1}$ and $\Xi_{1}$ correspond to the advection term in the equation of movement, $\Pi_{2}$ and $\Xi_{2}$ to the perturbation of the mixing length, $\Pi_{3}$ and $\Xi_{3}$ to the perturbation of the radiative cooling timescale of convection eddies, and $\Pi_{4}$ and $\Xi_{4}$ to the remaining parts.

Panel a) of Fig. 7 shows the absolute values of the eigenfunctions obtained with the MAD code. The perturbation of the local turbulent pressure reduces to zero toward the boundary determined by the Schwarzschild criterion (red line). However, the perturbation of the nonlocal turbulent pressure has amplitude even in the overshooting zone due to the nonlocal effects expressed as Eq. (12) (cyan line).

As shown in panel a), the second term of Eq. (16) is negligible (green line). Then, Eq. (16) simplifies to $\delta p_{\rm turb,l}/p_{\rm tot}\simeq\Pi\delta p_{\rm th}/p_{\rm tot}$. Using the adiabatic relation $\delta\rho/\rho=\delta p_{\rm th}/p_{\rm th}/\Gamma_{1}$, Eq. (15) becomes

[TABLE]

Although $\Pi$ is a complex number, we should pay attention to only its real part to discuss the frequency shift. Panel b) shows that the real part of $\Pi$ has a peak ($\log\,T\simeq 3.96$) located deeper than $p_{\rm turb,l}/p_{\rm th}$, which corresponds to the GGM approximation. However, their values are of the same order. It implies that the GGM treatment gives a good prediction to some extent.

Indeed, we can analytically understand that the GGM is valid in the bottom part of the convection zone. Since $\sigma\tau_{c}\gg 1\gg\omega_{R}\tau_{c}$ in such a part, we have

[TABLE]

and ${\cal G}/\Gamma_{1}$, which corresponds to the density perturbation in Eq. (4), and ${\cal GF}\Pi_{1}$ are much larger than the other terms in Eq. (17). Therefore, we can derive

[TABLE]

which implies that the situation is close to the GGM (Eq. 1) in the bottom part of the convection zone.

To see the contribution of each component of $\Pi$ to the frequency shift, we introduce the variational principle like Eq. (15):

[TABLE]

and

[TABLE]

for $i=1,2,3,4$. Panel c) shows that the terms with ${\cal G}/\Gamma_{1}$ and $\Pi_{1}$ dominantly contribute to the frequency shift. Although Fig. 7 shows the case of $n=28$, the contributions of $\Pi_{2}$, $\Pi_{3}$ , and $\Pi_{4}$ are even more negligible for the other lower-order modes since the mode amplitude is distributed in the inner region. Panel d) shows that all the terms except $\Pi_{2}$, related to the perturbation of the mixing length, certainly contribute to the perturbation of the turbulent pressure in the top part of the convection zone. Particularly, the low convective efficiency, namely the low value of $\Gamma\,[=(\omega_{R}\tau_{c})^{-1}]$, makes $\Pi_{3}$ contributive near to the upper boundary of the convection zone. However, the physical processes in the top part of the convection zone hardly contribute to the frequency shift. Then, Eq. (15) would be written as

[TABLE]

We note that we restricted our analysis to radial oscillations in this work. For nonradial oscillations, we should adopt Eq. (44) instead of Eq. (5). Besides, Eqs. (49) and (50) are no longer valid in the derivation. For high $\ell$ modes, the quantity $\ell(\ell+1)\xi_{h}$ may become important.

4 Conclusion

Previous studies (Rosenthal et al. 1995, 1999) have found that the frequencies obtained with the gas-gamma model (GGM) approximation better agree with the observations than those obtained with the reduced gamma model (RGM) approximation. This treatment is easy to adopt for computing the adiabatic oscillations of models including the turbulent pressure. However this crude approximation has no clear physical background. In this study, we computed the frequencies with a TDC formalism imposing the adiabatic condition. We found that the GGM provides closer frequencies to the TDC ones compared to the RGM. It implies that the GGM is superior to the RGM from not only observational but also theoretical viewpoints. Besides, the TDC computations reproduced the frequencies closer to the observation than did the GGM, regardless of the values of the free parameter $\beta$. Although our work is limited to the Sun, it is worth extrapolating our results obtained by the TDC to other stars. Using the variational principle, we found that the perturbation of the density and advection term mainly contribute to the frequency shift due to the perturbation of the turbulent pressure. Equation (27) can be then used to evaluate the frequency shift for adiabatic radial oscillations.

As discussed in previous studies (e.g. Brown 1984; Rosenthal et al. 1999; Sonoi et al. 2015), the turbulent pressure in the equilibrium model affects the frequencies because of the elevation of the outer layers. However its perturbation is also important for the frequencies, as discussed in this paper. Although this subject has been already shown by Houdek (2010) using the equilibrium convection models and the TDC formalism based on the theory of Gough (1977a, b), our study used the convection profiles obtained with the 3D simulations. As a first step, we limited ourselves to adiabatic oscillations and the effect of turbulent pressure. However, future works should consider nonadiabatic effects as well as the effect of convective backwarming. As for the latter, Trampedach et al. (2013, 2017) reported that the high temperature sensitivity of the opacity in the top of the convection zone causes warming by upflows of convection surpassing cooling by the downflows coupled with the non-linear nature of radiative transfer. The resultant net warming leads to the elevation of the outer layers as well as turbulent pressure. They also reported that the contribution of the backwarming has a similar magnitude as that of the turbulent pressure.

For the asteroseismology of solar-like stars, we need correct model frequencies. Since stellar ages are substantially affected by the surface effect, many studies have adopted the empirical relation based on the solar frequencies proposed by Kjeldsen et al. (2008). However Sonoi et al. (2015) found that this solar-calibrated relation has difficulty in correcting the frequencies in different stellar models and at different evolutionary stages. Then, it may be preferable to find a method of the correction based on a strong physical approach. Especially, the convective effects both in the equilibrium state and perturbation and also the nonadiabatic effect may be important for this problem. Therefore, we will extend the work of Sonoi et al. (2015) including these effects.

Acknowledgements.

T.S. has been supported by the ANR (Agence Nationale de la Recherche) program IDEE (Interaction Des Étoiles et des Exoplanètes) and CNES (Centre National d’Études Spatiales). H.G.L. acknowledges financial support by the Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject A4) of the German Research Foundation (DFG).

Appendix A Time-dependent convection formalism for nonlocal convection

Following Dupret et al. (2006b), we introduce the way to adopt the results given by the 3D simulations of nonlocal convection to the time-dependent convection (TDC) formalism of Grigahcène et al. (2005).

A.1 Hydrodynamical equations for local convection

The TDC formalism of Grigahcène et al. (2005) originates from the one proposed by Unno (1967). Later, Unno’s formalism was developed for nonradial oscillations by Gabriel et al. (1975). The classical mixing length theory (MLT) of Böhm-Vitense (1958) is the description for convection in the hydrostatic equilibrium state. On the other hand, the TDC formalism includes variation of convection on the dynamical timescale. However, if we impose the stationary condition on the formalism, we can obtain consistent results with the MLT.

First, we derive the equation of convection in the equilibrium state. We thus decompose the physical variables in the hydrodynamical equations into the mean flow and convective fluctuations as $y=\overline{y}+\Delta y$ for the scalars and $\overrightarrow{v}=\overrightarrow{u}+\overrightarrow{V}$ for the velocity. In Unno’s formalism, the convective fluctuation parts of the hydrodynamical equations of the continuity, movement, and energy conservation are given by

[TABLE]

where the notations follow the definitions introduced in Sect. 2.3 of this paper. To obtain the above equations, the following approximations have been made for the closure of Eqs. (29) and (30):

[TABLE]

We adopt the Boussinesq approximation, in which the pressure fluctuations are neglected except in the equation of movement (Eq. 29) and the density fluctuations are neglected in the equation of continuity (Eq. 28). Besides, the spatial variation in the density is assumed to be much smaller than that in the convective velocity in Eq. (28). The closure approximations (31), (32), and (33) follow the assumption that turbulent viscosity and thermal conductivity due to smaller eddies are expressed with the typical scale given by a representative convective element including them. Equation (34) is the usual closure equation of the MLT. Assuming constant coefficients and $\Lambda=8/3$, the above equations give the stationary solution consistent with the MLT:

[TABLE]

where ${\cal A}=P_{T}p_{\rm tot}/(2P_{\rho}\rho)[\kappa c_{p}\rho^{2}gl^{2}/(12acT^{3}p_{\rm tot})]^{2}$ and $\Gamma=(\omega_{\rm R}\tau_{c})^{-1}$. As mentioned above, the Boussinesq approximation includes the neglect of the density fluctuations in Eq. (28). However, this assumption is not valid in near-surface layers of solar-like stars since convective velocity can be comparable with sound speed of surrounding materials. Besides, the assumption that the spatial variation in the density is much smaller than that in the convective velocity is invalid in the deep part of a convection zone, where the surrounding structure is no longer homogeneous in the representative scale of convective eddies. However, this is a standard hypothesis made in most TDC approaches. Without such assumption, it is difficult to build a TDC formalism. Besides it is a consequence of the adoption of the MLT.

A.2 Perturbative theory for local convection

To consider the behaviour of convection with the oscillations, we perturb the above formalism, Eqs. (28)–(30), which allows us to evaluate the perturbation of correlated quantities of the convective fluctuations. However the closure described above is crude, and many complex physical processes are neglected including the whole cascade of energy. Then, uncertainty cannot be avoided when perturbing the closure terms. Because of such uncertainties, the unphysical, short wavelength oscillations appear in the eigenfunctions of the differential equations for the oscillations. To deal with this problem, Grigahcène et al. (2005) proposed to introduce a free complex parameter $\beta$ in the perturbation of the thermal closure equations:

[TABLE]

Introducing this parameter leads to phase lags between the oscillations and the way the turbulence cascade adapts to them.

Then, we search for the solutions of the perturbed convective fluctuation equations of the form $\delta(\Delta X)=\delta(\Delta X)_{\overrightarrow{k}}e^{i\overrightarrow{k}\cdot\overrightarrow{r}}e^{i\sigma t}$, assuming constant coefficients within the coarse grain, which is much larger than most of convective eddies but much smaller than the scale of the perturbation wavelength. Next, we integrate these particular solutions over all values of $k_{\theta}$ and $k_{\phi}$ so that $k^{2}_{\theta}+k^{2}_{\phi}=Ak^{2}_{r}$, keeping $A$ constant and that every direction of the horizontal component of $\overrightarrow{k}$ has the same probability. The value of $A$ is the free parameter, given by Eq. (9) based on the 3D simulation in this study. We have to introduce this distribution of $\overrightarrow{k}$ values to obtain an expression for the perturbation of the Reynolds tensor which allows the proper separation of the variables in terms of spherical harmonics in the equation of motion. Finally, the obtained values of the perturbation of the correlated values are implemented into the differential equations of the oscillations.

A.3 Procedure for taking equilibrium values given by 3D simulation into account

3D hydrodynamic simulations (e.g. Stein & Nordlund 1991, 1998; Rosenthal et al. 1999; Yang & Li 2007; Piau et al. 2014) provide much more realistic profiles of the convection zones than with the MLT. Here we discuss how to extend the above formalism to the nonlocal case following Dupret et al. (2006b).

As discussed in Sect. A.1, most of the uncertainties are included in the closure terms (Eqs. 31–34). Then, we modify these terms introducing a free function varying with depth, $\Omega$, which has the same meaning as in the formalism of Canuto & Mazzitelli (1991). It can be assumed to be a function of the convective efficiency $\Gamma$ following Canuto & Mazzitelli. We also set the usual mixing length $\alpha$ as an additional free function varying with the depth or $\Gamma$. More precisely, we multiply the left hand side of Eq. (31) by $\Omega(\Gamma)$ and the left hand side of Eqs. (32) and (33) by $1/\Omega(\Gamma)$. Then, Eqs. (29) and (30) become

[TABLE]

In the stationary case, these new equations have a form similar to the old ones (Eqs. 29 and 30). Equation (35) remains unchanged, giving the same meaning to $\Gamma$ as in the previous case. Equation (38) is still verified (with varying $\alpha$), but Eqs. (36) and (37) are slightly modified:

[TABLE]

[TABLE]

[TABLE]

By adjusting $\Omega$ and $\alpha$, we can fit these equations to the results given by 3D simulations in combination with Eqs. (12) and (13). The quantities $F_{\rm c,nl}$, $p_{\rm turb,nl}$, $(\nabla-\nabla_{\rm ad})$ and other thermodynamic quantities are deduced from the 3D simulations, and we take appropriate horizontal and time averages. Using Eqs. (12) and (13), the local counterparts of turbulent pressure and convective flux, $p_{\rm turb,l}$ and $F_{\rm c,l}$, are obtained based on $p_{\rm turb,nl}$ and $F_{\rm c,nl}$. Using these local counterparts, we obtain appropriate values of $\Omega$ and $\alpha$ at each location with Eqs. (43) and (38b), respectively.

To generalize the perturbative theory presented in Sect. A.2, we replace the equations of movement and energy conservation for the local case, (29) and (30), with the ones for the 3D case, (40) and (41). We follow the same procedure as in Sect. A.2. Assuming again constant coefficient and searching for solutions in the form of plane waves, we obtain the new expressions for the perturbed local convective quantities such as the convective flux and turbulent pressure.

The main uncertainties in this approach appear in the way to perturb $\Omega$ and $\alpha$. The free parameter $\beta$ introduced in Eq. (39) is also somehow related to these uncertainties. At present, we have no theoretical prescription how to perturb $\Omega$ and $\alpha$, and then we neglect their perturbations. However there is no reason to expect them to be small, and we should not be too optimistic when using this new perturbative treatment.

Here, we do not discuss the derivation which is very similar to those of Grigahcène et al. (2005). The final results of the perturbation of the radial components of the local convective velocities and convective flux are given as Eqs. (44) and (45). They are not so different from the former expressions (Eqs. 12 and 18 in Grigahcène et al. 2005):

[TABLE]

[TABLE]

We note that Eq. (44) becomes Eq. (5) for adiabatic radial oscillations.

On the other hand, the perturbation of the nonlocal turbulent pressure and convective flux is obtained by solving the eigenvalue problem of the differential equations of oscillations combining the perturbed equations of (12) and (13).

Appendix B Recasting the expression of turbulent pressure perturbation

Here, we describe the procedure for recasting the expression of turbulent pressure perturbation, which is required for the discussion in Sect. 3.2. We begin with Eqs. (4) and (5) and aim to express them as the linear combination of the thermal pressure perturbation, $\delta p_{\rm th}$, and the displacement, $\xi$.

For the perturbation of the mixing length, we adopt the expression,

[TABLE]

where the perturbation of the pressure scale height is described by

[TABLE]

To cancel the term ${\rm d}\delta p_{\rm tot}/{\rm d}p_{\rm tot}$, we adopt the perturbed equation of movement (Eq. D.3 in Grigahcène et al. 2005), neglecting the perturbation of the divergence of the Reynolds tensor:

[TABLE]

and the perturbed equation of continuity (Eq. D.1 in Grigahcène et al. 2005),

[TABLE]

For radial oscillation, the Poisson equation becomes

[TABLE]

With the adiabatic condition $\delta\rho/\rho=\delta p_{\rm th}/p_{\rm th}/\Gamma_{1}$, we can express Eq. (5) as a linear combination of $\delta p_{\rm th}$, $\xi$ and $\delta p_{\rm turb,l}$. We categorize the terms in Eq. (5) into four parts as follows. The first part is the term which stems from the advection term in the equation of movement:

[TABLE]

The second part is the perturbation of the mixing length:

[TABLE]

The third part is the perturbation of the inverse of the radiative cooling timescale of convection eddies:

[TABLE]

The last part corresponds to the remaining terms:

[TABLE]

Finally, we obtain the expression of $\delta p_{\rm turb,l}/p_{\rm tot}$ as Eq. (16). The coefficient of $\delta p_{\rm th}/p_{\rm tot}$ ($\Pi$) consists in $\Pi_{i}$’s ($i=1,2,3,4$) as shown in Eq. (17). The expression of $\Pi_{i}$’s is

[TABLE]

On the other hand, the coefficients of $\xi/R$, $\Xi$, are

[TABLE]