© J. Philidet et al. 2020

1. Introduction

Solar-like oscillations are known to be stochastically excited and damped by turbulence occurring close to the surface of low-mass stars (see e.g. Goldreich & Keeley 1977a,b or Samadi et al. 2015 for a review). The power spectral density of such oscillations is expected to feature a Lorentzian-shaped peak centred around their eigenfrequencies. This idealised line profile has been extensively used to fit observations (see e.g. Jefferies et al. 1991). However, as the resolution reached in helioseismic measurements (both ground-based and space-borne) has increased, it has become apparent that the observed line profiles feature a certain degree of asymmetry (see e.g. Duvall et al. 1993 for observations made at the geographic South Pole; Toutain et al. 1998 for data from the MDI and SPM instruments aboard the SOHO spacecraft).

Since the discovery of this skew symmetry in solar p-mode line profiles, several studies have been devoted to explaining this feature. In particular, it had been recognised early on that a source of excitation that is highly localised compared to the mode wavelength (which we refer in the rest of the paper as “source localisation”) could lead to a certain degree of mode asymmetry, depending on the position of the source (Gabriel 1992, 1993; Duvall et al. 1993; Abrams & Kumar 1996; Roxburgh & Vorontsov 1995, 1997). Line profile asymmetries have then been used to infer some properties of the turbulent source, especially its radial location and its multipolar nature (see e.g. Roxburgh & Vorontsov 1997; Nigam et al. 1998).

Furthermore, Duvall et al. (1993) noticed an inversion of the sense of asymmetry between spectrometric and photometric measurements, with line profiles in the velocity spectrum featuring more power in their low-frequency wing than in their high-frequency wing and vice-versa for line profiles in the intensity spectrum. Since intensity perturbations were expected to be proportional to velocity perturbations, one would have expected the asymmetries to be the same. Many hypotheses were posited to explain this puzzling result. Duvall et al. (1993) suggested that it was due to non-adiabatic effects lifting the proportionality relationship between the two kinds of perturbations (fluid displacement and temperature) but this hypothesis was later contradicted by Rast & Bogdan (1998). Non-adiabaticity was brought up again later on by Georgobiani et al. (2003) who suggested that the explanation resided in radiative transfer between the mode and the medium. Indeed, the observed radiation temperature corresponds to the gas temperature at local optical depth τ = 1. But optical depth depends on opacity, which non-linearly depends on temperature. Therefore, the temperature fluctuations due to the oscillating mode entails opacity fluctuations, which in turn impacts the “observed” radiation temperature. Given the non-linear nature of the κ − T relation, this modulation decreases the observed temperature fluctuations more significantly in the low-frequency wing of the mode than in its high-frequency wing. Since this radiative transfer does not impact the velocity measurements, this could explain the asymmetry reversal between velocity and intensity spectra. Using 3D simulations, Georgobiani et al. (2003) computed mode line profiles in both the velocity and the intensity power spectrum alternatively at mean unity optical depth and instantaneous unity optical depth. Their results indeed show that the modulation of the “observed” intensity fluctuations due to radiative transfer could be significant enough to reverse the sense of mode asymmetry. One of the hypothesis enjoying the most support for asymmetry reversal, however, is based on the effect of turbulent perturbations partially correlated with the mode, which thus impact its line profile (Nigam et al. 1998; Roxburgh & Vorontsov 1997; Rast & Bogdan 1998; Kumar & Basu 1999). Indeed, a part of these perturbations is coherent with the mode and, thus, leads to interference. This interference term may be constructive or destructive, depending on the phase difference between the mode and the coherent turbulent perturbations. For frequencies at which the interference is constructive, the power spectral density is slightly elevated, whereas it drops slightly for frequencies at which it is destructive. Typically, in the vicinity of a resonant mode, the dependence of the phase difference between mode and turbulent perturbation is such that the interference term is constructive for frequencies located in one wing of the mode and destructive in the other. Therefore, as a result of this interference behaviour, one of the wings falls off more slowly and the other more rapidly, leading to mode asymmetry. It has been suggested that the degree of correlation between the turbulent perturbations and the oscillation it excites is higher in intensity than in velocity, so that it changes the sign of mode asymmetry only in the intensity spectrum. While it is widely accepted that correlated turbulent perturbations must be taken into account to explain asymmetries in the intensity spectrum, the question of whether it has a significant impact on the velocity spectrum remains an open issue (see e.g. Jefferies et al. 2003).

The possibility that correlated turbulent fluctuations have an affect on mode asymmetry has led many authors to include them in their models for the power spectrum. Even though correlated noise was introduced to explain the particular puzzle of asymmetry reversal between velocity and intensity measurements, several models include correlated noise in the velocity spectrum as well as in the intensity spectrum. This is the case for the model developed by Severino et al. (2001) and later used, for instance, by Barban et al. (2004), which includes three types of noise (coherent-correlated, coherent-uncorrelated and incoherent) in both the velocity spectrum, the intensity spectrum, and the velocity-intensity cross-spectrum. They considered, however, that the “pure oscillation” (without the noise) has a Lorentzian shape, thus discarding the contribution of source localisation. This model was later refined by Wachter & Kosovichev (2005) to take this contribution into account.

These prior studies have one thing in common, however, and that is that they all treat the various sources of asymmetry (mainly source localisation and correlated noise) in a simplified, parametrised fashion. Indeed, the excitation is consistently modelled as a point-like source, with radial position and multipolar development left as free parameters. This prescription remains somewhat unsatisfactory in the sense that it does not take into account the finer properties of the source of excitation, such as its spatial extent or its dependence on frequency, for instance. As such, these prior models lack a realistic description of the source of excitation. Likewise, for models including the effect of noise on the power spectrum, the various relative amplitudes and phase differences between modal oscillation and correlated noise in both spectra are also left as free parameters. The approach followed by these studies is to find best-fit values for all their free parameters by fitting their model to observations in order to localise the source.

In contrast, in the present paper, we follow a different approach: we model both the source of excitation and the correlated background by constraining their properties using an analytical model of stochastic excitation, coupled with a 3D simulation of the solar atmosphere. The novelty of our approach lies in the fact that we do not fit a parametrised model to the observations but, instead, we predict the dependence of mode asymmetry on frequency, which we then compare to observations in order to validate our model. Our model of mode asymmetry is, therefore, both more realistic (in its description of the source of excitation) and more complete (in its lack of freely adjustable parameters). It can then be used to deepen our understanding of the underlying physical mechanisms behind mode asymmetry by exploring how varying a given physical constraint impacts the results yielded by our model. Finally, our model allows for a much higher predictability of mode asymmetry, which is essential when it comes to applying these results to other solar-like oscillators. We note that this paper is devoted to the modelling of the velocity power spectrum only and, as a result, we do not address the problem of asymmetry reversal, which is a separate challenge altogether.

These efforts to model the line profiles of solar-like oscillations are also necessary in order to correctly infer mode properties from observations. Indeed, it was discovered early on that using a Lorentzian shape to fit skew symmetric line profiles led to a significant bias in the eigenfrequency determination, which may be higher than the frequency resolution achieved by helioseismic measurements (Duvall et al. 1993; Abrams & Kumar 1996; Chaplin et al. 1999; Thiery et al. 2000; Toutain et al. 1998). Such eigenfrequency determination bias has also been revealed for solar-like oscillations in stars other than the Sun by Benomar et al. (2018). Inversion methods used to infer the internal structure of solar-like oscillators, whether they be analytical or numerical, require a very accurate determination of the mode eigenfrequencies. For spectra extracted from very long time series, the resolution is high enough that this bias in eigenfrequencies impacts the results obtained by inversion methods (see e.g. Toutain et al. 1998, who show that the difference between the sound speed squared inferred from symmetric and asymmetric fits can reach 0.3% in the core). When fitting these observations, mode asymmetry must, therefore, be taken into account. Since it has proven very difficult to determine accurate mode eigenfrequency without prior knowledge on its line shape, obtaining an a priori model of p-mode line profiles is of primary importance.

In this paper, we present a predictive model of solar radial p-mode line profile in the velocity spectrum. In particular, we use a realistic model for stochastic excitation, following a method similar to that of Samadi & Goupil (2001). Furthermore, we include the effect of correlated turbulent perturbations in the model in a non-parametrised way, unlike what was done in previous works (see e.g. Severino et al. 2001). The paper is structured as follows: we present the analytical model of the Sun’s velocity power spectral density in Sect. 2 and its numerical implementation in Sect. 3. We then present the results yielded by our model concerning mode asymmetry in Sect. 4. In Sect. 5, we briefly describe the development of a toy model to describe the impact of source localisation on mode asymmetry and use it to interpret our results; we also investigate the matter of the influence of correlated turbulent perturbations. We then confront our results with the related observations in Sect. 6 and discuss the issue of eigenfrequency determination bias entailed by the skewness of the mode line profiles.

2. Modelling the p-mode line profiles

To extract the asymmetries of solar radial p-modes, we first need to model their line profile in the velocity power spectrum. In this section, we present the analytical developments that led us to this model. First, we define the disc-integrated velocity power spectrum in terms of the radial fluid displacement. We then present the inhomogeneous, radial wave equation associated to the acoustic modes and detail how convolving its Green’s function with its inhomogeneous part gives us access to the velocity power spectral density.

2.1. Definition of the velocity power spectral density

Before embarking on a discussion of the actual modelling of the line profiles, the spectrum from which they are extracted needs to be defined. In this paper, we restrict ourselves to the study of radial acoustic modes in the Sun. Furthermore, as part of the definition of the spectrum, we include the effect of limb-darkening and of disk integration that affect the Sun-as-a-star measurements. We note, however, that other instrumental effects – in particular mode leakage – are not accounted for.

To derive an expression for the observed power spectral density, we separate the total surface velocity into an oscillatory part v_osc and a turbulent part u, where it is understood that the modes are described by the oscillatory part. The observations made for the Sun-as-a-star are obtained by integrating the velocities over the entire solar disk. Neither the mode velocity (for radial modes), nor the turbulent perturbations depend on the point of the disk at which it is estimated; however, the projection on the line of sight n does. This integration over the solar disk is affected by limb-darkening h(μ) (where μ refers to the cosine of the angle between the local radial direction and the line of sight). Furthermore, since it is the turbulent perturbations that excite the mode, a certain fraction of the former must be correlated with the latter, so that the contribution of turbulent perturbations to the velocity spectrum must be considered.

With these considerations, the observed velocity power spectral density can be expressed as

$$\begin{array}{c}\hfill P(\omega )={\displaystyle \frac{1}{{\displaystyle \int \mathrm{d}\mathrm{\Omega }h(\mu )}}}\langle {\left|{\displaystyle \int \mathrm{d}\mathrm{\Omega }h(\mu )(\widehat{{\mathit{v}}_{\mathrm{osc}}}(\omega )+\widehat{\mathit{u}}(\omega )).\mathit{n}}\right|}^2\rangle ,\end{array}$$(1)

where the integration is performed over the solar disk, Ω refers to the solid angle, n is the unit vector along the line of sight, v_osc is the mode velocity, u represents the fluctuations of the turbulent velocity around its mean value, ω is the angular frequency, the notation $\left(\widehat{.}\right)$ refers to temporal Fourier transform, and ⟨.⟩ refers to ensemble average. Since we are only considering radial modes, v_osc is exclusively radial. Thus, Eq. (1) becomes

$$\begin{array}{c}\hfill P(\omega )=\langle {\left|{\displaystyle \widehat{v_{\mathrm{osc}}}(\omega )\int \mathrm{d}\mathrm{\Omega }\mu \stackrel{\sim }{h}(\mu )+\int \mathrm{d}\mathrm{\Omega }\stackrel{\sim }{h}(\mu )\widehat{u_n}(\omega )}\right|}^2\rangle ,\end{array}$$(2)

where u_n is the component of the turbulent velocity along the line of sight. We introduced the reduced limb-darkening $\stackrel{\sim }{h}(\mu )$ so that its integral over the solar disk is normalised to unity.

We expand the square in the above expression and we consider that the term containing $⟨|\widehat{u_n}{|}^2⟩$ is negligible compared to the terms that contain $⟨|\widehat{{\mathit{v}}_{\mathrm{osc}}}{|}^2⟩$ and $\mathrm{Re}(⟨\widehat{u_n}{\widehat{{\mathit{v}}_{\mathrm{osc}}}}^{\star }⟩)$, respectively. Indeed, the power spectral density is several orders of magnitude higher for the mode velocity than for the turbulent velocity (typically, the former is of order 10⁵ m² s⁻² Hz⁻¹, while the latter is of order 10 m² s⁻² Hz⁻¹, e.g. Turck-Chièze et al. 2004, Fig. 2), so that

$$\begin{array}{c}\hfill ⟨|\widehat{u_n}{|}^2⟩\ll \mathrm{Re}(⟨\widehat{u_n}{\widehat{v_{\mathrm{osc}}}}^{\star }⟩)\ll ⟨|\widehat{v_{\mathrm{osc}}}{|}^2⟩,\end{array}$$(3)

where the notation Re refers to the real part of a complex quantity, and ⋆ refers to its complex conjugate. Finally, we obtain

$$\begin{array}{cc}\hfill P(\omega )=& {\left({\displaystyle \int \mathrm{d}\mathrm{\Omega }\mu \stackrel{\sim }{h}(\mu )}\right)}^2\langle {\left|\widehat{v_{\mathrm{osc}}}(\omega )\right|}^2\rangle \hfill \\ \hfill & {\displaystyle +2\int \mathrm{d}\mathrm{\Omega }\mu \stackrel{\sim }{h}(\mu )\mathrm{Re}\left({\displaystyle \int \mathrm{d}\mathrm{\Omega }\stackrel{\sim }{h}(\mu )\langle \widehat{v_{\mathrm{osc}}}(\omega ){\widehat{u_n}}^{\star }(\omega )\rangle }\right).}\hfill \end{array}$$(4)

The first term corresponds to the spectral power density of the mode velocity v_osc. In itself, the line profile generated by this term is already asymmetric; indeed, it has been known for a long time that source localisation can generate line profile asymmetry (see e.g. Abrams & Kumar 1996; Roxburgh & Vorontsov 1997; Chaplin & Appourchaux 1999). The second term corresponds to what the literature commonly refers to as correlated turbulent perturbations and which is also expected to significantly impact mode asymmetry in photometric measurements (see e.g. Nigam et al. 1998; Roxburgh & Vorontsov 1997; Kumar & Basu 1999), although its importance in velocity measurements is not as clear.

2.2. The inhomogeneous wave equation

Going further, we write the radial wave equation associated to v_osc with the same formalism as Unno et al. (1989). We detail its derivation in Appendix A. Although we included both the source terms due to Reynolds stress fluctuations and non-adiabatic pressure fluctuations in the computation detailed in Appendix A, we only consider the former in the following. Indeed, it is the dominant source of excitation for acoustic modes in the Sun (e.g. Belkacem et al. 2008). When it is temporally Fourier transformed, the inhomogeneous wave equation for radial modes reads:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\Psi }}_{\omega }}{\mathrm{d}{r}^2}}+({\displaystyle \frac{{\omega }^2}{c^2}}-V(r)){\mathrm{\Psi }}_{\omega }=S(r),\end{array}$$(5)

where c is the sound speed, the wave variable Ψ_ω(r) is related to the radial fluid displacement ξ_r(r) through

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\omega }(r)=rc(r)\sqrt{{\rho }_0(r)}{\xi }_r(r),\end{array}$$(6)

and the acoustic potential and source term are given by

$$\begin{array}{c}\hfill \begin{array}{c}V(r)={\displaystyle \frac{N^2-4\pi G{\rho }_0}{c^2}}+{\displaystyle \frac{2}{x^2}}{\left({\displaystyle \frac{\mathrm{d}x}{\mathrm{d}r}}\right)}^2-{\displaystyle \frac{1}{x}}{\displaystyle \frac{{\mathrm{d}}^{2}x}{\mathrm{d}{r}^2}},\hfill \\ x(r)={\displaystyle \frac{r\sqrt{I}}{c}},\hfill \\ S(r)={\displaystyle \frac{r}{c\sqrt{{\rho }_0(r)}}}{\displaystyle \frac{\mathrm{d}{p}_t^{\prime }}{\mathrm{d}r}},\hfill \\ I(r)=exp\left({\displaystyle {\int }_0^r{\displaystyle \frac{N^2}{g_0}}-{\displaystyle \frac{g_0}{c^2}}\mathrm{d}{r}^{\prime }}\right),\hfill \end{array}\end{array}$$(7)

where r is the radial coordinate, ρ₀ is the density, N is the Brunt-Väisälä frequency, g₀ is the gravitational acceleration, G is the gravitational constant, and $p_t^{\prime }$ refers to the fluctuations of the Reynolds stress around its mean value. Indeed, only the fluctuating part of the Reynolds stress contributes to the source term S(r) and its mean value only modifies the equilibrium structure. The subscript 0 refers to the equilibrium structure and all the above quantities are dependent on the radius at which they are estimated, even when not explicitly specified. We note that we only model radial modes in this paper, so that the wave equation (Eq. (5)) is of the second order despite the fact that we did not use the Cowling approximation.

Mode damping is not included in Eq. (5). Indeed, we did not take into account the feedback of modal oscillations on the equilibrium state through modulations in the fluid density, pressure, opacity, etc. Such feedback allows mechanical work and thermal transfer to occur from the mode to the medium in which it develops; depending on the phase-lag between these different modulations energy can be exchanged with the surrounding medium. However, the modelling of damping rates of solar-like oscillations is extremely difficult (Samadi et al. 2015). Thus, we directly introduce damping in the wave equation in the form of a mode lifetime, or, equivalently, by a (frequency-dependent) linewidth Γ_ω, so that the wave equation takes the following form

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\Psi }}_{\omega }}{\mathrm{d}{r}^2}}+({\displaystyle \frac{{\omega }^2+j\omega {\mathrm{\Gamma }}_{\omega }}{c^2}}-V(r)){\mathrm{\Psi }}_{\omega }=S(r),\end{array}$$(8)

where j denotes the imaginary unit and the linewidths Γ_ω are inferred from observations. We used the line-widths presented in Baudin et al. (2005) (see their Table 2), which were inferred from GOLF data. Note, however, that our definition of Γ_ω corresponds to their Γ multiplied by 2π, or equivalently to twice their damping rate η. We completed these data with low-frequency line-widths obtained by Davies et al. (2014) through BiSON, which go as low as ∼900 μHz (see their Table 1). We reproduce the dependence of the linewidth we used on frequency in Table 1. We note that damping can potentially be a source of mode asymmetry. However, the impact of damping on mode asymmetry is very weak compared to the other sources of asymmetry (Abrams & Kumar 1996), so that the direct introduction of observed line-widths in our model is unlikely to have an impact on our results.

2.3. Expression of the velocity power spectral density

By definition, the Green’s function G_ω(r_o, r_s) is the value taken by the function Ψ_ω at the radius r = r_o (the variable r_o refers to the height in the atmosphere at which the spectrum is observed and the variable r_s refers to the position of the point-like source term), where Ψ_ω is the solution to the inhomogeneous wave equation,

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\Psi }}_{\omega }}{\mathrm{d}{r}^2}}+({\displaystyle \frac{{\omega }^2+j\omega {\mathrm{\Gamma }}_{\omega }}{c^2}}-V(r)){\mathrm{\Psi }}_{\omega }=\delta (r-r_{\mathrm{s}}),\end{array}$$(9)

and δ refers to the Dirac function. Once the Green’s function is known, it can be used to express explicitly v_osc in Eq. (4). Indeed, on the one hand, the general solution to the inhomogeneous wave equation with a source term S(r_s) is

$$\begin{array}{c}\hfill {\displaystyle {\mathrm{\Psi }}_{\omega }(r_{\mathrm{o}})=\int \mathrm{d}{r}_{\mathrm{s}}{G}_{\omega }(r_{\mathrm{o}},r_{\mathrm{s}})S(r_{\mathrm{s}}),}\end{array}$$(10)

where the source term is given by Eq. (7). The pulsational velocity v_osc is related to the variable Ψ_ω through

$$\begin{array}{c}\hfill v_{\mathrm{osc}}(r_{\mathrm{o}})={\displaystyle \frac{j\omega }{r_{\mathrm{o}}c(r_{\mathrm{o}})\sqrt{{\rho }_0(r_{\mathrm{o}})}}}{\mathrm{\Psi }}_{\omega }(r_{\mathrm{o}}).\end{array}$$(11)

Using the source term given by Eq. (7) in Eqs. (10) and (11) and after finally performing an integration by part, we write the velocity Fourier transform at angular frequency ω as

$$\begin{array}{cc}\hfill \widehat{v_{\mathrm{osc}}}(\omega ,r_{\mathrm{o}})=& -{\displaystyle \frac{j\omega }{r_{\mathrm{o}}c(r_{\mathrm{o}})\sqrt{{\rho }_0(r_{\mathrm{o}})}}}\hfill \\ \hfill & {\displaystyle \times \int {\mathrm{d}}^3{\mathit{r}}_{\mathrm{s}}\mathbf{\nabla }\left(G_{\omega }({\mathit{r}}_{\mathrm{s}},r_{\mathrm{o}}){\displaystyle \frac{||{\mathit{r}}_{\mathrm{s}}||}{c({\mathit{r}}_{\mathrm{s}})\sqrt{{\rho }_0({\mathit{r}}_{\mathrm{s}})}}}\right).\left({\rho }_0\widehat{u_r\mathit{u}})({\mathit{r}}_{\mathrm{s}})\right).}\hfill \end{array}$$(12)

In the following, the observation height r_o will be fixed, so that we drop it for ease of notation. However, since the observation height depends on the transition line used for the observations and on whether the observations rely on spectrometric or photometric measurements, it significantly varies from instrument to instrument (see Sect. 6 for more details).

Using Eq. (12) in Eq. (4) then gives an expression for the velocity power spectral density in terms of Green’s function G_ω(r_s):

$$\begin{array}{c}\hfill P(\omega )={\left({\displaystyle \int \mathrm{d}\mathrm{\Omega }\mu \stackrel{\sim }{h}(\mu )}\right)}^2[\langle {\left|\widehat{v_{\mathrm{osc}}}(\omega )\right|}^2\rangle +C(\omega )],\end{array}$$(13)

where $\langle {\left|\widehat{{\mathit{v}}_{\mathrm{osc}}}(\omega )\right|}^2\rangle$ and C(ω) are given, respectively, by Eqs. (B.19) and (B.28). We note that the effects of limb-darkening and disk integration are now contained in a single factor and, thus, these will only have an effect on mode amplitude. Since the asymmetry of a mode does not depend on its amplitude, it is not impacted by such a factor.

The calculations leading from Eqs. (4)–(13) are detailed in Appendix B. In the following, we only provide the main steps and assumptions. We split the calculations two ways, focussing separately on the first term inside the brackets of Eq. (13), which we hereby refer to as the leading term, and on its second term, which we hereby refer to as the cross term.

2.3.1. Closure models

The calculations leading from Eqs. (4)–(13) involve the evaluation of fourth-order and third-order two-point correlation moments of the turbulent velocity. Therefore, an appropriate closure model is needed to express these high-order moments as a function of second-order moments. We devote the following subsection to presenting and developing these closure models.

Fourth-order moments. To describe the fourth-order correlation moments of the turbulent velocity, we make use of the quasi-normal approximation (QNA). This closure model consists in considering that all turbulent quantities are normally distributed, in which case fourth-order moments can be analytically expressed as a combination of second-order moments (Lesieur 2008):

$$\begin{array}{c}\hfill ⟨abcd⟩=⟨ab⟩⟨cd⟩+⟨ac⟩⟨bd⟩+⟨ad⟩⟨bc⟩,\end{array}$$(14)

where a, b, c, and d refer to any turbulent scalar quantity. Applying the QNA to isotropic, homogeneous turbulence inhibits energy transfers among modes of different wave numbers, thus leading to violations of the energy conservation principle (Kraichnan 1957). This is due to the fact the QNA entails vanishing third-order correlation moments. When it comes to estimating the fourth-order moments, however, the picture is different. Belkacem et al. (2006a) have studied the validity of the QNA for two-points, fourth-order correlation moments of the vertical turbulent velocity, in the form of $\langle u_{r,1}^2{u}_{r,2}^2\rangle$ (where the indices 1 and 2 refer to two different points in space), using 3D simulations of the solar atmosphere. They found that the dependence of this correlation moment on the distance ΔX between the two points is correctly estimated by the QNA but that its absolute value (which can be taken as the corresponding one-point moment) is not. Consequently, the amplitude of the modes are largely underestimated when the QNA is used. However, the asymmetry of the modes does not depend on their amplitude, so that mode asymmetry will be unaffected by a discrepancy in the absolute value of the two-points, fourth-order moments. As such, the decomposition given by Eq. (14) can be considered valid when it comes to studying mode asymmetry.

Third-order moments. While the QNA provides an adequate closure relation for fourth-order moments, as mentioned above, it assumes vanishing third-order moments. Therefore, in order to estimate these third-order moments, we make use of another closure model, the Plume closure model (PCM), which was developed by Belkacem et al. (2006b). The idea behind this closure model is to separate the flows directed upwards from those directed downwards (the latter being referred to as plumes) and to apply the QNA to both separately. The anisotropy between the two types of flow – in particular, turbulence is more prominent in the downwards plumes (e.g. Goode et al. 1998) – yields non-vanishing third-order correlation moments:

$$\begin{array}{cc}\hfill ⟨u_r{(\mathit{R},t)}^2{u}_r(\mathit{R}+\mathit{r},t+\tau )⟩=& [a{(1-a)}^3-a^3(1-a)]\delta u^3\hfill \\ \hfill & -a(1-a)[2⟨\stackrel{\sim }{u_d}(\mathit{R},t)\stackrel{\sim }{u_d}(\mathit{R}+\mathit{r},t+\tau )⟩\hfill \\ \hfill & +⟨\stackrel{\sim }{u_d}{(\mathit{R},t)}^2⟩]\delta u,\hfill \end{array}$$(15)

where u_r is the vertical component of the turbulent velocity, a is the relative horizontal section of the upflows, δu is the difference between the mean velocity of the two types of flows (considering their respective signs, it actually is the sum of their absolute values), and $\stackrel{\sim }{u_d}$ is the fluctuation of the vertical velocity around its mean value in the downflows.

We note that, strictly speaking, the third-order moment given by Eq. (15) and yielded by the PCM are centred. However, we consider that the mean value of the overall vertical velocity of the flow is sufficiently low (compared to its standard deviation for instance) to be neglected. Therefore, these moments may interchangeably refer here either to centred or non-centred moments.

We also note that this closure relation is written here in terms of $\stackrel{\sim }{u_d}$ (i.e. the turbulent fluctuations in the downflows only). It would be more practical to rewrite it in terms of u_r (i.e. the total turbulent fluctuations). The two are related through

$$\begin{array}{c}\hfill ⟨\stackrel{\sim }{u_d}(\mathit{R},t)\stackrel{\sim }{u_d}(\mathit{R}+\mathit{r},t+\tau )⟩={\displaystyle \frac{1}{1-a}}⟨u_r(\mathit{R},t)u_r(\mathit{R}+\mathit{r},t+\tau )⟩-a\delta u^2.\end{array}$$(16)

2.3.2. The leading term

In the following, we detail the derivation of the first term of Eq. (13). This term corresponds to the pulsational velocity itself, without correlated turbulent perturbations. As such, any asymmetry featured by this term alone represents the effect of source localisation. The first step consists in separating the scales relevant to the turbulent velocity u from the scales relevant to both the medium stratification and the oscillating mode (respectively, the pressure scale height and the mode wavelength). The scale separation approximation is not realistic in the subsurface layers (in particular, the mode wavelength is comparable to the typical correlation length associated with turbulence); however, for want of a better alternative, we are led to use this approximation in the following.

Since the integral defining $\widehat{{\mathit{v}}_{\mathrm{osc}}}(\omega )$ in Eq. (12) contains the turbulent velocity fluctuations squared, expanding the square of its modulus will raise these fluctuations to the fourth. The contribution of turbulence to the expression of v_osc thus takes the form of two-points, fourth-order correlation moments of the turbulent velocity. We use the closure relation presented and detailed in Sect. 2.3.1 to express them as a function of second-order moments.

We then use analytical expressions for the second-order moments of the turbulent velocity. We describe the second-order moment of the ith and jth component of the turbulent velocity in terms of its spatial and temporal Fourier transform ϕ_ij(k, ω). For isotropic turbulence, it reads (Batchelor 1953):

$$\begin{array}{c}\hfill {\varphi }_{\mathit{ij}}(\mathit{k},\omega )={\displaystyle \frac{E(k,\omega )}{4\pi k^2}}({\delta }_{\mathit{ij}}-{\displaystyle \frac{k_i{k}_j}{k^2}}),\end{array}$$(17)

where E(k, ω) is the specific turbulent kinetic energy spectrum, k is the norm of the wavevector k, k_i and k_j are its ith and jth component, and δ_ij is the Kronecker symbol. The integration over the solid angle of wave vectors k is straightforward, and only an integral over the norm of k remains. However, solar turbulence close to the photosphere is known to be highly anisotropic. To take this anisotropy into account, we follow the formalism developed by Gough (1977). In this formalism, the integral over the solid angle of k is simply readjusted by adding an anisotropy factor given by Eq. (B.10) (see Appendix B in Samadi & Goupil 2001).

Following Stein (1967), we then decompose E(k, ω) into a spatial part E(k), which describes how the turbulent kinetic energy is distributed among modes of different wave numbers, and a temporal part χ_k(ω), which describes the statistical distribution of the life-time of eddies of wavenumber k

$$\begin{array}{c}\hfill E(k,\omega )=E(k){\chi }_k(\omega ).\end{array}$$(18)

In order to model the spatial and temporal part of the spectrum of turbulent kinetic energy, we followed two different approaches, described in the following.

The “theoretical spectrum” model. We use theoretical prescriptions to model both the spatial spectrum E(k) and the temporal spectrum χ_k(ω) of turbulent velocity. Based on the assumption that turbulent flows are self-similar, Kolmogorov’s theory of turbulence leads to a spatial spectrum E(k)∝k^−5/3 in the inertial range, between k = k₀ (where k₀ is the scale at which the kinetic energy is injected in the turbulent cascade, and is henceforth referred to as the injection scale) and the dissipation scale (at which the turbulent kinetic energy is converted into heat). Given the very high Reynolds number characterising solar turbulence (Re ∼ 10¹⁴), we cast the dissipation scale to infinity. Then, following Musielak et al. (1994), we extend the turbulent spectrum below the injection scale by considering that E(k) takes a constant value for k <  k₀. This extended spectrum, referred to as the broadened Kolmogorov spectrum (BKS) was introduced to account for the broadness of the maximum of E(k). The BKS can be written as

$$\begin{array}{c}\hfill E(k)=\{\begin{array}{cc}0.652{\displaystyle \frac{u_0^2}{k_0}}\hfill & \mathrm{if}0.2k_0<k<k_0\hfill \\ 0.652{\displaystyle \frac{u_0^2}{k_0}}{\left({\displaystyle \frac{k}{k_0}}\right)}^{-5/3}\hfill & \mathrm{if}{k}_0<k,\hfill \end{array}\end{array}$$(19)

where $u_0^2\equiv \langle {\mathit{u}}^2(\mathit{r})\rangle /3$ and the 0.652 factor is introduced so that the total specific kinetic energy of the turbulent spectrum matches $u_0^2/2$. Therefore, the spatial spectrum is parametrised solely by the injection scale k₀. However, the injection scale varies significantly between the sub-surface layers and the atmosphere (Samadi et al. 2003), so that we keep it free in our model and allow for it to depend on the radial coordinate r.

Following Samadi et al. (2003), we consider a Lorentzian shape for the temporal spectrum χ_k(ω), which is supported both by numerical simulations (Samadi et al. 2003) and by theoretical arguments. Indeed, a noise described by a stationary, Gaussian Markov process in time is expected to relax exponentially, meaning that the resulting eddy-time correlation is expected to be a decreasing exponential, and its Fourier transform a Lorentzian function (Belkacem et al. 2011). The width ω_k associated to eddy-time correlation is linked to the life-time of the eddies of wavenumber k. Dimensional arguments would suggest that ω_k ∝ ku_k, where u_k is the typical velocity associated to the eddies of wavenumber k. However, there remains a substantial indetermination on the actual value of ω_k, so that, following Balmforth (1992), we consider:

$$\begin{array}{c}\hfill {\omega }_k=2k{u}_k/\lambda ,\end{array}$$(20)

where λ is a dimensionless, constant parameter. Overall, the only input parameters of this model are k₀(r) and λ.

The “numerical spectrum” model. In the second model, we extract the spatial spectrum E(k) from a 3D hydrodynamic simulation of the solar atmosphere, provided by the CO⁵BOLD code (see Sect. 3 for details). This simulation gives us access to the velocity field as a function of all three spatial coordinates and time. In order to extract the turbulent spectrum E(k), we average the velocity field temporally, then isolate each horizontal slice in the simulated cube and perform a 2D Fourier transform of each slice separately, thereof which we only retain the radial part. This gives us a spectrum E(k) for each vertical point in the simulation. Finally, we renormalise each spectrum so that

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{+\infty }\mathrm{d}kE(k)={\displaystyle \frac{u_0^2}{2}},}\end{array}$$(21)

where u₀ is also extracted from the 3D atmospheric simulation, by averaging the fluid velocity squared temporally and horizontally, and using the definition $u_0^2=\langle {\mathit{u}}^2(\mathit{r})\rangle /3$.

The temporal spectrum χ_k(ω) is also treated in a slightly different manner than in the “theoretical model” above. Indeed, the arguments invoked above to justify the Lorentzian shape of the spectrum, while valid for most of the relevant time scales associated to turbulent eddies, are no longer valid for shorter time scales, that is, for higher angular frequencies. Belkacem et al. (2010) argued that if the time correlation associated to small eddies indeed originates from their advection by larger, energy-bearing eddies – a hypothesis referred to as the sweeping assumption – one recovers a Gaussian spectrum instead of a Lorentzian one. The transition between a Lorentzian spectrum, valid for low angular frequencies, and a Gaussian spectrum, valid for high angular frequencies, occurs at the cut-off angular frequency ω_E, which is given by the curvature of the eddy-time correlation function at τ = 0 (Belkacem et al. 2010):

$$\begin{array}{c}\hfill {\omega }_E=k{u}_0.\end{array}$$(22)

Since a Gaussian spectrum would fall off much more rapidly than a Lorentzian spectrum, we simply consider that χ_k vanishes entirely for ω >  ω_E,

$$\begin{array}{c}\hfill {\chi }_k(\omega )=\{\begin{array}{cc}{\displaystyle \frac{1}{2arctan({\omega }_E/{\omega }_k){\omega }_k}}{\displaystyle \frac{1}{1+{(\omega /{\omega }_k)}^2}}\hfill & \mathrm{if}\omega <{\omega }_E\hfill \\ 0\hfill & \mathrm{if}{\omega }_E<\omega .\hfill \end{array}\end{array}$$(23)

We modified the prefactor so that χ_k meets the normalisation condition. The typical life time of eddies of wavenumber k, parametrised by ω_k, is still given by Eq. (20). We note that the convolution of the function χ_k(ω) with itself must be computed to evaluate the leading term (see Eq. (B.13)). While the convolution of a Lorentzian function with itself straightforwardly yields a Lorentzian function with a width twice as large, the convolution of the modified spectrum above with itself is slightly different, but can be obtained analytically as

$$\begin{array}{cc}\hfill ({\chi }_k\ast {\chi }_k)(\omega )=& {\displaystyle \frac{1}{2\pi {\omega }_k}}{\displaystyle \frac{1}{1+{(\omega /2{\omega }_k)}^2}}\hfill \\ \hfill & \times {\displaystyle \frac{\pi (arctan\left({\displaystyle \frac{{\omega }_E}{{\omega }_k}}\right)-arctan\left({\displaystyle \frac{\omega -{\omega }_E}{{\omega }_k}}\right))}{4{arctan}^2\left({\displaystyle \frac{{\omega }_E}{{\omega }_k}}\right)}}·\hfill \end{array}$$(24)

Physically, taking the cut-off frequency into account significantly decreases the predicted amplitudes for high frequency modes. As far as mode asymmetry is concerned, we found that it did not have a significant impact in the “theoretical spectrum” model. In contrast, it substantially changes mode asymmetry at high frequency in the “numerical spectrum” model, which is why we only introduce it in the latter. The reason is the following: the spatial spectrum of turbulent energy falls off much more rapidly with k in the “theoretical spectrum” than in the “numerical spectrum”. Therefore, small spatial scales play a more important role in the latter. Since the cut-off frequency only impacts these small scales, it is natural that taking it into account should only have a significant impact on the “numerical spectrum” model.

2.3.3. The cross term

The derivation of the cross term follows essentially the same steps as for the leading term. The main difference is that the turbulent velocity correlation moments that appear are now two-point, third-order moments. We use the closure relation presented in Sect. 2.3.1 to express them, as we did with the fourth-order moments, as a function of two-point, second-order moments of the turbulent velocity. We then use the same analytical description of second-order moments as the one we used for the leading term. The rest of the calculations is very similar to those described in Sect. 2.3.2 and leads to the second term in Eq. (13).

These two models – the “theoretical spectrum” and “numerical spectrum” models – are complementary in the sense that the first one allows us to explore the impact of the properties of turbulence on mode asymmetries and gives physical insight into this problem, whereas the second one relies on fewer input parameters and, therefore, has more predictive capability (we recall here that the former requires the parameters λ and k₀(r) to be set, while the latter only requires λ). Consequently, in the following, we present and develop the results yielded by both.

3. Numerical implementation

In this section, we detail how we numerically implemented the model presented in Sect. 2. We describe how we obtained the solar equilibrium state in which the acoustic modes develop and how we integrated the inhomogeneous wave equation given by Eq. (8). Having obtained the solar radial p-mode line profiles, we then detail how we extracted their asymmetries and perform several tests to validate our model and its numerical implementation.

3.1. The solar equilibrium state

The acoustic potential given by Eq. (7) depends only on the equilibrium structure of the Sun. We extracted the potential from a 1D solar model provided by the evolutionary code CESTAM (Morel 1997; Marques et al. 2013). The 1D model includes treatment of the convective flux (using standard mixing-length theory with no overshoot) and of the radiative flux (using the Eddington approximation). On the other hand, turbulent pressure, rotation, and diffusion processes are ignored.

However, 1D stellar models do not fully take into account the more complex physical phenomena taking place in the uppermost layers of a star; in particular, the rapid transition between the convective zone and the superficial radiative region (Kupka & Muthsam 2017). This leads to significant biases in the equilibrium structure. Since the excitation of solar oscillations precisely takes place in these layers, it is essential that we model them more accurately. To do so, we use a 3D hydrodynamic simulation of the solar atmosphere computed using the CO⁵BOLD code (Freytag et al. 2012). The modelled region includes the super-adiabatic peak just below the photosphere and goes up to the lower atmosphere of the star.

It is now possible to construct a “patched” model of the solar interior. We use the solar patched model computed by Manchon et al. (2018). The process of constructing patched models has been extensively discussed (e.g. Trampedach 1997; Samadi et al. 2008) and the particular case of the patched model used in this paper is described in much detail in Manchon et al. (2018). The basic idea is to transform the 3D atmosphere into a 1D atmosphere through temporal and horizontal averaging and then to replace the surface layers of the 1D stellar interior with this 1D atmosphere. We note that the input parameters of the CESTAM model used to describe the solar interior (age, total stellar mass, mixing-length parameter α_MLT, and helium abundance) are chosen so that the top layers match the bottom layers of the CO⁵BOLD atmosphere. Here the model was computed with the mixing-length parameter α_MLT = 1.65, an initial helium abundance of Y_init = 0.249, and an initial metallicity of Z_init = 0.0135. Figure 1 shows the acoustic potential profile V(r) given by this solar patched model, computed using Eq. (7).

Finally, we use the same simulation of the solar atmosphere to extract the various parameters appearing in the analytical description of the source term; in particular, the standard deviation u₀ associated to turbulent velocities, the anisotropy factor G given by Eq. (B.10), as well as the parameters used in the Plume closure model (see Eq. (15)). Specifically, the ensemble average appearing in the definition of u₀ was computed by performing a temporal and horizontal average of the norm of the velocity squared in the 3D simulation of the solar atmosphere. We also used this same simulation to extract the spatial spectrum of turbulent kinetic energy in the “numerical spectrum” model (see Sect. 2).

3.2. Integration of the inhomogeneous wave equation

To compute one value of P(ω) for one value of the angular frequency ω (i.e. one point in the velocity power spectrum), we convolve the Green’s function G_ω(r_s) associated to Eq. (8) with the stochastic source term S(r_s) (see Sect. 2). It is then possible to reconstruct the velocity power spectral density, and in particular the line profile of the resonant modes, point by point (typically, we only need 10 points regularly spaced between ω₀ − Γ_ω0 and ω₀ + Γ_ω0, where ω₀ is the angular eigenfrequency of the mode and Γ_ω0 its linewidth). In the following, we describe how the wave equation was integrated and how we extracted its Green’s function.

For a given angular frequency ω, we carry out the integration using a 4th-order Runge-Kutta scheme (Press et al. 1986) with the acoustic potential V(r) given by the solar equilibrium state (see Sect. 3.1). Given that radial modes develop in the entire solar volume, we perform this integration over the entire solar radius, between r = 0 and r = r_max. We note that r_max refers not to the photospheric radius, but to the maximum radial extent of the solar model described in Sect. 3.1, so that r_max >  r_photosphere.

We imposed Dirichlet boundary conditions on the wave equation. The condition at the centre is straightforward: by definition, Ψ_ω(r = 0) = 0. At r_max, we impose a vanishing Lagrangian pressure perturbation (which physically means that the atmosphere of the Sun is force-free). The continuity equation and the equation of state allow us to rewrite this latter condition in terms of Ψ_ω:

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_{\omega }}{\mathrm{d}r}}+{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}[ln\left({\displaystyle \frac{r}{c\sqrt{{\rho }_0}}}\right)]{\mathrm{\Psi }}_{\omega }=0.\end{array}$$(25)

The use of Dirichlet boundary conditions leads us to implement a shooting method: we perform the integration with Ψ_ω(r = 0) = 0, and tune the initial slope (i.e. the value of dΨ_ω/dr at r = 0) until the other boundary condition is met. Note that this method is not similar to the shooting method usually implemented to solve the eigenvalue problem associated to the determination of acoustic mode eigenfrequencies: here, the pulsation ω is fixed, and it is the initial slope that is tuned to meet the boundary condition at the surface. The difference between these two methods is that in the inhomogeneous problem, the initial slope (or, alternatively, the amplitude of the mode) is fixed by the amplitude of the source of excitation. The shape of the eigenfunction, however, remains the same as in the homogeneous problem.

This method enables us to extract the Green’s function associated to the wave equation (Eq. (8)). To obtain one value of the Green’s function G_ω(r_s), for one value of the pulsation ω and one value of the source position r_s, we carry out the integration of the inhomogeneous equation as described above, adding a point-like source term to the numerical scheme. The source is normalised in such a way that the right-hand side equals 1/h when the source falls within the integration radial step, and 0 otherwise (h is the radial step of the integration).

This integration gives us the radial oscillation profile Ψ_ω(r), and we simply extract its value at a fixed radius r_o, which corresponds to the height in the atmosphere at which the spectrum is measured. We note that the presence of damping in the wave equation implies that it is complex-valued. As such, the Green’s function is complex-valued as well.

Finally, to calculate the integrals over source positions which appear in Eqs. (B.19) and (B.28), we compute the Green’s function using the above method for a grid of source positions r_s, while ω is kept constant. This grid corresponds to the radial grid provided by the 3D atmospheric model described in Sect. 3.1.

We also use the aforementioned 3D model to extract the physical quantities appearing in both the leading term of Eq. (13) (the anisotropy factor G, the turbulent velocity fluctuations u₀), and the cross term (the parameters a and δw in the PCM; see Sect. 2 for a definition of these parameters).

Using Eq. (13) provided by the model developed in Sect. 2 and the radial grid of Green’s functions computed using the above method finally allows us to extract the value of P(ω) for each value of ω.

3.3. Fitting of the mode asymmetries

We fit the line profile of the modes following Nigam & Kosovichev (1998) with the formula

$$\begin{array}{c}\hfill P(\omega )=H_0{\displaystyle \frac{{(1+Bx)}^2+B^2}{1+x^2}},\end{array}$$(26)

where x = 2(ω − ω₀)/Γ_ω0 is the reduced pulsation frequency. The fit contains three free parameters (H₀, ν₀ and B), the last of which is defined as the asymmetry parameter. We illustrate the dependence of the line profile on B in Fig. 2. In particular, B <  0 means that the peak contains more power in the low-frequency side (that corresponds to negative asymmetry), B >  0 means that the high-frequency side contains more power (that corresponds to positive asymmetry), while with B = 0 we recover a classic, Lorentzian profile. Figure 2 also shows that the mode does not peak exactly at the eigenfrequency, but rather at a slightly higher (for B >  0) or lower value (for B <  0). This can have important repercussions for the determination of the mode eigenfrequencies from observations, as we discuss in Sect. 6.3.

Note that several other definitions of the asymmetry parameters can be found in the literature. Korzennik (2005) prefers to adjust the mode line profiles with

$$\begin{array}{c}\hfill P(\omega )\propto {\displaystyle \frac{1+\alpha (x-\alpha /2)}{x^2+1}},\end{array}$$(27)

and defines the asymmetry parameter as α_n, l. Meanwhile, Vorontsov & Jefferies (2013) use the following formula:

$$\begin{array}{c}\hfill P=H[{\left({\displaystyle \frac{Acos(\varphi -S)}{1-R^2}}\right)}^2+B^2],\end{array}$$(28)

where the frequency variable is ϕ, and the asymmetry parameter is defined as S. While these three formulas have been derived in different ways, they are perfectly equivalent close to the eigenfrequency (for x ≪ 1, or equivalently for ϕ ≡ 0 (mod π)), with S ∼ B ∼ α/2.

Finally, Gizon (2006) provides the asymmetry parameter defined as (see also Benomar et al. 2018):

$$\begin{array}{c}\hfill \chi =2B{\omega }_0/{\mathrm{\Gamma }}_{\omega }.\end{array}$$(29)

The author quantified the mode asymmetry by means of the relative positions of the local maxima (peaks) and minima (troughs) in the power spectral density: minima located half-way between the neighbouring maxima lead to symmetric line profiles, while a deviation from this behaviour leads to asymmetric line profiles. The parameter χ derived from these considerations is independent from both the amplitude and the line-widths of the modes.

The formulas presented above only lead to different line profiles far from the central frequency, whereas they are equivalent in our range of interest. We opted for the definition given by Nigam & Kosovichev (1998) (Eq. (26)) because it is the most commonly used.

To ensure the significance of fitting an asymmetric profile to the mode obtained through our model, we compared the results produced by the fitting formula Eq. (26) and by a symmetric, Lorentzian profile (that is, imposing B = 0 in Eq. (26)). The asymmetric fits led to excellent agreement with the modelled line profiles; however, the symmetric fits led to substantial discrepancies, with one wing consistently falling off more rapidly than the numerical line profile and the other too slowly. Finally, it should be noted that the excellent fit given by Eq. (26) to the numerical line profile is independent from the number of points used for the adjustment; we have indeed performed a similar fit with thrice the number of points, without any loss of accuracy and the resulting asymmetry parameter B was the same to an excellent approximation.

3.4. Validation of the method

Using the method presented above, we extracted solar radial modes of radial orders n = 6 to n = 30. Indeed, the formula used for the fit and given by Eq. (26) does not converge properly for higher-order modes (because the increasing linewidths lead to mode overlapping), while we did not have access to observed linewidths for lower-order modes. In addition to their line profile asymmetries, we also extract other fundamental properties, namely their eigenfrequencies, amplitudes, and eigenfunctions. In the following, to support the validity of our model, we compare the mode properties we obtained with similar properties obtained through other methods.

First, we compare the eigenfrequencies obtained through our model to the eigenfrequencies of the 1D adiabatic oscillations calculated using the ADIPLS code (Christensen-Dalsgaard 2011). For this validation, we did not make use of the patched model described in Sect. 3.1 but, rather, the corresponding unpatched model. The reason is that the patching procedure produces a small discontinuity of the physical quantities at the patching point, which can affect the eigenfrequencies calculated by ADIPLS. We recover the correct eigenfrequencies, with errors not exceeding ∼0.1%. Since mode asymmetry is only expected to vary on the scale of ∼mHz, modelled asymmetries will not be significantly affected by such small discrepancies of the eigenfrequencies.

Our numerical method also allows us to extract the radial profile Ψ(r) of the eigenmodes. We compare them in Fig. 3 to the eigenfunctions calculated using the same 1D adiabatic oscillations obtained through the ADIPLS code and presented above. The figure shows that the modes obtained through our model have eigenfunctions that are very similar to those obtained through this dedicated code, which further supports the validity of the model we have used.

Finally, we compare the mode amplitudes obtained through our model to the observed ones. To that end, we estimated the velocity power spectrum at an observation height of 340 km, which corresponds to the observation height of the GOLF instrument as estimated by Baudin et al. (2005), following Bruls et al. (1992). By definition, the velocity amplitude squared is the total area under the mode peak, so that it depends both on its maximum H and on its width Γ,

$$\begin{array}{c}\hfill v_{\mathrm{osc}}=\sqrt{\pi H\mathrm{\Gamma }}.\end{array}$$(30)

We note that when it is used to treat observational data, this formula also contains a geometric factor to account for instrumental effects, including mode visibility. This factor is, however, irrelevant in our case.

We show in Fig. 4 the comparison between the mode amplitudes v_osc obtained through our “numerical spectrum” model and the mode amplitudes inferred from observations performed by the GOLF instrument (Baudin et al. 2005). The free parameter λ of our model (cf. Sect. 2) has been adjusted so as to obtain the best possible agreement. As a consequence, our model does not hold any predictive power when it comes to mode amplitudes. However, the fact that we manage to retrieve a very good agreement with observational data by using reasonable values of the input parameters is still a solid sign that our model is valid. In particular, we correctly recover the frequency at maximum amplitude ν_max, as well as the slopes on both the low-frequency and the high frequency limit. To conclude on the matter, we emphasise that the asymmetry parameter B is independent from the mode amplitude, so that potential discrepancies concerning the latter should not affect the former.

4. Results for the asymmetry profile B(ν)

Using the model presented in Sect. 2, numerically implemented using the method presented in Sect. 3, we extract the solar p-modes line profile asymmetries B(ν) throughout a large part of the spectrum, between n  =  6 (ν  ∼  1 mHz) and n  =  30 (ν  ∼  4.3 mHz). In this section, we present the results yielded by our model, focusing on the dependence of the asymmetry parameter B on frequency (which we hereafter shorten to the expression “asymmetry profile”) and on the impact of our different input parameters on the asymmetry profile.

As we detailed in Sect. 2, we followed two different approaches to model the turbulent kinetic energy spectrum. The first one, which we refer to as the “theoretical spectrum” model, uses the prescription given by Kolmogorov’s theory of turbulence and which we have described in detail in Sect. 2. The second approach, which we refer to as the “numerical spectrum” model, uses the turbulent spectrum extracted from the 3D hydrodynamic simulation of the solar atmosphere described in Sect. 3. In this section, we present separately the asymmetry profile B(ν) yielded by both models.

4.1. The “theoretical spectrum” model

This model relies on a prescription for the properties of turbulence. It contains the following input parameters: the temporal spectrum of turbulent kinetic energy, parametrised by the dimensionless quantity λ, which is defined by Eq. (20), and its spatial spectrum, parametrised by k₀(r), which is defined as the (radius-dependent) injection wavenumber of turbulent kinetic energy. We let the latter depend on r in order to account for the fact that the typical size of turbulent eddies drastically depends on where they are located with respect to the photosphere. It is known that the size of the energy-bearing eddies increases with height, so that the injection scale k₀ decreases with r (Samadi et al. 2003). We simplify the situation by considering that the injection rate only takes two values: k₀(r) = k_0, int below the photosphere, and k₀(r) = k_0, atm above the photosphere. This picture crudely corresponds to what is observed in 3D atmospheric simulations (Samadi et al. 2003). In the following, we denote the ratio between the two as R_k ≡ k_0, int/k_0, atm. This leaves only three input parameters in our model: λ, k_0, int, and k_0, atm; or equivalently λ, k_0, int, and R_k.

In Fig. 5, we keep λ and R_k constant, and we show the asymmetry profile B(ν) for several values of k_0, int. Despite the fact that we vary k_0, int across almost one order of magnitude, the asymmetry profile B(ν) does not depend significantly on the absolute value of k₀, except close to ν ∼ 1.7 mHz. By comparison, its dependence on both R_k and λ is more substantial, especially at high frequencies (cf. Figs. 6 and 7). Since k_0, int does not seem to play an important role, we keep it fixed in the following, and focus on the impact of the other two input parameters, λ and R_k.

Figure 5 illustrates the main qualitative features of the asymmetry profile B(ν). In fact, together with Figs. 6 and 7, it shows that the qualitative behaviour of the asymmetry profile is largely model-independent. Thus the asymmetries of the solar radial p-mode line profiles are negative across a large part of the spectrum, in agreement with solar observations (see for instance Duvall et al. 1993). Furthermore, the asymmetry profile B(ν) exhibits three distinct local extrema: the absolute value of B increases below ∼1.7 mHz, decreases between ∼1.7 and ∼3 mHz, increases again between ∼3 mHz and ∼4 mHz, and finally decreases again above ∼4 mHz. Note, however, that this last extremum is, unlike the other ones, somewhat impacted by the values given to the different input parameters of the model.

The first two local extrema (∼1.7 and ∼3 mHz) correspond respectively to the beginning and end of the damping rate plateau. Indeed, the asymmetry parameter B depends on the linewidth of the modes, so it is natural that a sudden change in the behaviour of the latter should reflect on the former. The third extrema is not so easily explained and it will be discussed in Sect. 5.

Figure 6 shows how the asymmetry profile B(ν) depends on R_k. An increase of this parameter attenuates low-frequency mode asymmetries (below ν_max ∼ 3 mHz), while on the contrary, the high-frequency modes (above ν_max) become more asymmetric. The effect is significantly more substantial for the latter than for the former. Asymmetries close to ν_max, however, are not affected by the parameter R_k whatsoever.

Figure 7 shows that the impact of λ on the asymmetry profile B(ν) is quite similar, albeit inverted, in the sense that |B| increases with λ for low-frequency modes and decreases for high-frequency modes. Similarly, B is barely impacted by a change of λ close to ν_max. Furthermore, the asymmetry profile B(ν) undergoes saturation, in the sense that it ceases to depend on λ when it is increased above a certain value. In the following, we denote this threshold as λ_sat. Figure 7 shows that λ_sat ∼ 1. This dichotomy between λ ≲ 1 and λ ≳ 1 originates in the Lorentzian nature of the temporal turbulent spectrum: depending on the value of λ, the angular frequencies relevant to solar p-modes are either in the low frequency part or in the high frequency part of the spectrum. We do not go into too much detail here as we discuss this matter further in Sect. 5.

4.2. The “numerical spectrum” model

In the “numerical spectrum” model, which describes the properties of turbulence more realistically, there is only one input parameter left, λ. In this sense, it has a greater predictive power than the previous model. The qualitative behaviour of the asymmetry profile B(ν) and, in particular, the positions of the different local extrema featured by B(ν), are, in this model, rather independent from λ and in agreement with what we observed in the scope of the previous model.

However, the input parameter λ does have an impact on the quantitative behaviour of the asymmetry profile B(ν). We show in Fig. 8 the asymmetry profile B(ν) obtained with the “numerical spectrum” model (see Sect. 2) for several values of λ. As for the “theoretical spectrum” model, B is always negative and features several local extrema at ν ∼ 1.7, 3 and 4 mHz.

As for the dependence of B(ν) on λ, two distinct regimes can be separated. Below ν_max ∼ 3 mHz, we recover the same dependence of the asymmetry parameter B with λ as we obtained in the scope of the “theoretical spectrum” model, with absolute values of B increasing with λ. The picture at frequencies higher than ν_max is, however, somewhat different. The asymmetry profile B(ν) features a local minimum at ν ∼ 4 mHz; the curve inflexion grows sharper as λ increases up to λ ∼ 1, after which this part of the asymmetry profile does not significantly depend on λ. In this sense, the asymmetry profile B(ν) seems to undergo the same saturation behaviour as in the “theoretical spectrum” model (see Sect. 4.1), for the same value λ_sat ∼ 1. The fact that we recover approximately the same threshold gives us confidence that this particular feature of the asymmetry profile B(ν) is not a mere artefact of one model or the other but, rather, it is a genuine effect based on a physical origin. Again, we postpone the discussion of the physical origin of this behaviour to Sect. 5.

5. Impact of the properties of turbulence on mode asymmetry

Line profile asymmetry of solar-like oscillations have two main causes: localisation of the source of excitation (see for instance Duvall et al. 1993) and correlation with the turbulent perturbations (see for instance Nigam et al. 1998). In the following, we investigate both contributions in light of the results yielded by the “theoretical spectrum” model and presented in Sect. 4.1. With its various input parameters, the “theoretical spectrum” model allows us to understand the physical origin of mode asymmetry. In this section, then we only consider this model, although the conclusions are valid for the “numerical spectrum” model as well. We first discuss how source localisation and correlated turbulent perturbations can skew the mode line profiles. In particular, we support the discussion concerning source localisation with a simplified toy-model of mode excitation, which we describe in Appendix C. We then use this discussion to interpret the results yielded by our model. Additionally, we show that the contribution of the correlated turbulent perturbations to the mode asymmetries is negligible in the velocity power spectrum.

5.1. Origin of mode asymmetry

5.1.1. Effect of source localisation on mode asymmetry: generic arguments

The fact that the source of excitation of a mode is spatially localised can affect the skewness of the mode line profile in Fourier space. There are several ways of describing the impact of source localisation on mode asymmetry.

One way is to make use of the analogy between the development of acoustic modes in the stellar cavity and the phenomenon of optical interference in a Fabry-Pérot cavity. This analogy was used to account for the acoustic mode asymmetry in the Sun by Gabriel (1992), Duvall et al. (1993), among others. The idea is that acoustic, stationary modes in the Sun can be described by means of two progressing waves, propagating in opposite directions. Each of these waves follows the same cycle: they propagate one way, get refracted on the lower turning point of the acoustic cavity, then propagate backwards, get reflected on the upper turning point, and so forth. As a result of these multiple reflections and refractions on both turning points, the acoustic waves pass multiple times through the same regions and, therefore, interfere with each other (and with themselves). This interference pattern leads to the development of resonant modes in the cavity. What we observe then is the evanescent tail of these modes in the atmosphere, which lies outside the resonant cavity.

Let us now consider that the source of the waves is located at a certain point within the cavity. The waves propagating outwards and inwards will have travelled over different distances before interfering with one another and this difference of travel times will depend on the location of the source. The shape of the mode line profile is directly related to the dependence of the phase difference between the outwards and inwards interfering waves on frequency. Since this phase difference is not exactly symmetric about the mode eigenfrequency, neither is its line profile; and given that it depends on the source location, mode asymmetry is indeed a marker of source localisation.

Another physical interpretation of how source localisation can bring about mode asymmetry has been proposed by Rast & Bogdan (1998), and later refined by Rosenthal (1998). They remarked that mode asymmetry could be mathematically described by the relative position of local maxima (or peaks) and local minima (or troughs) in the power spectrum. Peaks located exactly halfway between their neighbouring troughs feature symmetric, Lorentzian line profiles. However, if one of the neighbouring trough is closer than the other, the peak in question appears skewed and, depending on which trough is closest, its asymmetry parameter is either positive or negative.

The position of the peaks are simply related to the eigenmodes of the solar acoustic cavity. As for the position of the troughs, in the special case of a point-like source of excitation, with a given multipolar decomposition, the authors showed that it is related to the eigenmodes of the atmosphere truncated at the source position, with a vanishing external boundary condition depending on the multipolar nature of the source. In that interpretation, the position of the troughs thus depends on both the position and the multipolar decomposition of the source.

Yet another way to describe the impact of source localisation on mode asymmetry is to consider the eigenfunction of the mode. In order to illustrate this, we present in Appendix C a very simplified toy-model of mode excitation, where the source is considered point-like and the acoustic cavity is simplified to a square well potential. From this toy-model we draw the following conclusion: for a given frequency, the amplitude of the wave is proportional to the eigenfunction associated with the wave at the source of excitation. In particular, excitation at a mode’s antinode is much more efficient than at a mode’s node.

With this conclusion in mind, let us consider the situation illustrated by Fig. 9. The blue and red curves represent the radial profile of the acoustic wave for two different angular frequencies. It can be seen that an increase of ω causes the radial profile of the oscillation to “shrink” radially. Therefore, the amplitude of the oscillation as seen by the source will either increase or decrease with ω, depending on its position. More specifically, a source at r = r₁ (see illustration in Fig. 9) will see the amplitude of the oscillation increase with ω, and a source at r = r₂ will see it decrease. In light of the conclusion presented in the previous paragraph, it can be deduced that if the source is located at r = r₁, the right wing of the mode line profile will be slightly elevated compared to the left wing, thus leading to positive asymmetry. Likewise, the asymmetry generated by a source at r = r₂ will be negative.

From the illustration in Fig. 9, it is straightforward to see that the dichotomy between the r = r₁ case and the r = r₂ case is based on the relative position of the source and the nodes and antinodes of the mode, or, in other words, on the sign of the derivative of the absolute value of the eigenfunction. To be more specific, one has to separate the case of a source inside and outside the acoustic cavity. If the source is inside the cavity, the r = r₁ case (i.e. case where source localisation entails positive asymmetry) corresponds to any source position located above a node and below an antinode of the oscillation profile, whereas the r = r₂ case (i.e. the case where source localisation entails negative asymmetry) corresponds to any source position located above an antinode and below a node. Here we recall that a node is a point at which the wave amplitude is zero and an antinode is a point at which it is maximal. If the source is outside the cavity, however, it is always as in the r = r₂ case and, thus, it always generates negative asymmetry: indeed, the outside of the cavity corresponds to an evanescent zone for the acoustic waves so that the absolute value of the eigenfunction always decreases in this region.

It should be noted that we only consider this toy-model in the present subsection. In the following sections, we return to the discussion of our model, simply using the conclusions drawn above to interpret the results which it yields.

5.1.2. Correlated turbulent fluctuations

Acoustic modes in the Sun are excited by fluctuations of turbulent nature – more specifically by turbulent fluctuations of the Reynolds stress or non-adiabatic pressure perturbations. It is therefore natural that a part of the turbulent fluctuations should be not only coherent, but statistically correlated with the oscillating mode.

The resulting interference between the mode and the turbulent fluctuations leads, in turn, to mode asymmetry. In order to illustrate this, let us consider a mode whose line profile is intrinsically Lorentzian and turbulent fluctuations whose power spectral density is constant over the width of the mode under consideration. We then have

$$\begin{array}{c}\hfill P(x)={|{\displaystyle \frac{A_m}{x+j}}+A_n{e}^{j{\varphi }_n}|}^2,\end{array}$$(31)

where P is the total power spectral density, x = 2(ω − ω₀)/Γ_ω is the reduced frequency (ω₀ is the angular eigenfrequency of the mode, and Γ_ω its linewidth), A_m and A_n are the (real) amplitudes associated to the mode and the noise respectively, ϕ_n is the phase difference between the mode and the noise, and j is the imaginary unit. Expanding the module squared, we obtain

$$\begin{array}{c}\hfill P(x)={\displaystyle \frac{A_m^2}{1+x^2}}+{\displaystyle \frac{2A_m{A}_n}{\sqrt{1+x^2}}}sin(arctanx+{\varphi }_n)+A_n^2.\end{array}$$(32)

The first term of the right-hand side of Eq. (32) corresponds to a Lorentzian profile and is symmetric about x = 0. The third term simply acts as an offset and does not introduce any mode asymmetry. The second term, however, is clearly not symmetric at x = 0, unless ϕ_n = ±π/2. For instance, if ϕ_n = 0, this term is even antisymmetric. In other words, the interference between the mode and the noise is destructive in the left wing of the mode and constructive in its right wing. As such, the power spectral density P(x) is higher than the Lorentzian profile in the right wing and lower in the left wing, thus entailing positive mode asymmetry. The sign and magnitude of the mode asymmetry depends on the amplitude A_n, that is, on the degree of correlation between the mode and the turbulent fluctuations, as well as on the phase difference ϕ_n, both of which are included in the model we developed in Sect. 2.

5.2. Contribution of source localisation to B(ν)

In the previous subsection, we summarised the impact of source localisation on mode asymmetry by stating the following: a source within the resonant cavity of a mode entails negative asymmetry if it is located above an antinode and below a node of the associated eigenfunction and positive asymmetry otherwise; a source outside the resonant cavity always entails negative mode asymmetry. With this in mind, we set out to interpret the results obtained in Sect. 4 in the scope of the “theoretical spectrum” model.

Once applied to the case of solar p-mode excitation, this rule can be rephrased in the following way. There is a dichotomy between the effect of the turbulent eddies located below the upper turning point of the mode and those located above. The former skew the mode line profile one way or the other depending on their height relative to the nodes and antinodes of the mode eigenfunction. The latter always skew the mode line profile so that it feature negative asymmetry. Until now, we have only discussed the case of a point-like source of excitation. However, the driving region of the solar p-modes, while localised around the super-adiabatic peak just below the photosphere, has a certain spatial extent. As such, driving turbulent eddies can be found both below and above the upper turning point of the modes and the observed mode asymmetry is due to the combination of both.

As we mention above, the sense of asymmetry created by turbulent eddies below the upper turning point depends on their position with respect to the nodes and antinodes of the modes. Since the wavelength of the modes is much larger than the spatial extent of the driving region, it is sufficient to study the position of the super-adiabatic peak with respect to the nodes and antinodes of the mode eigenfunction. We illustrate this in Fig. 10, which shows the position of the node and antinode of the radial modes obtained through our model, with respect to the super-adiabatic peak, where the stochastic excitation is mainly located. Below ν ∼ 3.5 mHz, the super-adiabatic peak lies above the closest antinode (dashed blue line in Fig. 10), thus generating negative asymmetries. Above the aforementioned threshold, the closest antinode is above the super-adiabatic peak and it generates positive asymmetry. Close to 3.5 mHz, the super-adiabatic peak coincides with an antinode, so that the asymmetry is very low.

The dichotomy between turbulent eddies below and above the upper turning point stands thus: excitation localised below the upper turning point makes modes with ν  ≲  ν_max negatively asymmetric and modes with ν  ≳  ν_max positively asymmetric; excitation localised above the upper turning point makes all modes negatively asymmetric; since the source of excitation has a certain spatial extent, the total asymmetry is a combination of both cases. This is in perfect accordance with the dashed blue curve of Fig. 6 which shows the asymmetry profile B(ν) with R_k  =  λ  =  1 (i.e. imposing the same turbulent spectrum everywhere). Indeed, this curve shows that B is negative for low frequencies, and positive for high frequencies. If the source of excitation was only located below the upper turning point of the modes, the threshold between the two regimes would be ∼3.5 mHz; however, turbulent eddies above the upper turning point generate additional negative asymmetry, thus shifting the curve downwards and increasing the threshold between the B <  0 and the B >  0 regimes.

The results presented in Figs. 5 and 6 are also easily interpreted. Indeed, Eqs. (B.19) and (B.28) show that the efficiency of stochastic excitation scales as $k_0^{-4}$, where k₀ is the injection wavenumber of turbulent kinetic energy. Therefore, keeping R_k constant does not change the contribution of atmospheric turbulence relatively to the contribution of turbulence below the upper turning point, and consequently only impacts the mode amplitude, not its asymmetry.

However, decreasing k_0, atm with respect to k_0, int increases the contribution of atmospheric eddies relatively to eddies below the upper turning point. Therefore, increasing R_k makes asymmetries at high frequencies decrease, and the frequency above which B >  0 increases. The ratio R_k needed for the asymmetry profile to be negative throughout the entire spectrum is only R_k ∼ 2, which is explained by the high sensitivity of the excitation efficiency on k₀ (since it scales to $k_0^{-4}$).

The physical interpretation of the influence of λ on the asymmetry profile B(ν) is not as straightforward. It cannot be interpreted in the same way as the influence of k₀(r), since we consider λ to be uniform throughout the region of excitation. Furthermore, it cannot be interpreted in terms of the relative contribution of the leading and cross term (see Eq. (13)), because both terms scale to λ. The existence of the threshold λ_sat can be explained as follows: depending on λ, the typical period of solar p-modes compares differently to the typical turbulent eddy-time correlation, that is, ω compares differently to ω_k, which indeed depends on λ. Depending on whether ω <  ω_k or ω >  ω_k, the temporal turbulent spectrum χ_k(ω) vary differently with ω: indeed, χ_k is almost flat for very low frequencies, whereas it decreases as ω⁻² for high frequencies. As such, as λ is increased, it is expected that the qualitative behaviour of the mode properties – including mode asymmetry – changes when ω ∼ ω_k. This explanation is supported by the value found for λ_sat. Indeed, if we take k ∼ 10⁻⁶ m⁻¹, u_k ∼ 10³ m s⁻¹ and ω ∼ 10⁻³ rad s⁻¹, we obtain ω/ω_k ∼ λ. Consequently, λ ∼ 1 does correspond to the threshold at which ω and ω_k have the same order of magnitude.

5.3. Contribution of correlated turbulent perturbations to B(ν)

Earlier in this paper, we illustrate how a certain degree of correlation between the oscillating modes and the turbulent fluctuations can create mode asymmetry. Furthermore, it has been claimed (Nigam et al. 1998) that correlation between pulsational velocity and acoustic turbulent perturbations is at the root of the inversion of the sign of mode asymmetry between spectrometric and photometric measurements. This suggests that this correlation plays a crucial role when it comes to interpreting photometric data, although other explanations exist (Duvall et al. 1993; Georgobiani et al. 2003). However, determining whether or not this role can be disregarded in the velocity spectrum or if it must be taken into account at all (even if it is not so significant as to change the sign of the mode asymmetries) remains an open question.

Our model allows us to shed some light upon this issue. In the following, we compute the asymmetry profile B(ν) alternatively with and without the cross term C(ω) in Eq. (13). In terms of amplitude, the leading term will, unsurprisingly, dominate over the correlation term C(ω); however, one must keep in mind that the asymmetry of the mode line profile is a subtle effect, and it is possible that, albeit negligible in amplitude, C(ω) impacts the asymmetry as much as the leading term.

Figure 11 shows the relative difference between the asymmetry profile B(ν) when the cross term is respectively taken into account and discarded. This relative difference is highest at low frequency, and almost vanishes when ν >  ν_max. This can be easily explained by the fact that low-frequency modes have the smallest power spectral density, while power spectral density associated to the noise is higher: in contrast, the effect of correlated turbulent perturbations is at its most substantial at the lowest frequencies. Nevertheless, even in the lowest part of the spectrum, the relative difference in asymmetry between a model with and without the cross term does not exceed 3%, which is much smaller than the dispersion characterising observed asymmetries. We therefore conclude that the dominant source of asymmetry in velocity data is the source localisation and that the effect of correlated turbulent perturbations can be disregarded.

It is possible to support this conclusion with a simple order of magnitude estimation. Indeed, B is of order 10⁻², which means it is necessary for the cross term to represent at least 1% of the leading term to have a significant impact on the asymmetry profile. The Cauchy-Schwarz inequality provides an upper bound to the cross term: $C(\omega )⩽{\widehat{\mathit{v}}}_{\mathrm{rms}}(\omega ){\widehat{u}}_{\mathrm{rms}}(\omega )$, where the case of equality happens when the correlation between the mode v and the turbulence u is optimal. Therefore

$$\begin{array}{c}\hfill {\displaystyle \frac{\langle {\left|\widehat{v}(\omega )\right|}^2\rangle }{C(\omega )}}\gtrsim {\displaystyle \frac{{\widehat{v}}_{\mathrm{rms}}(\omega )}{{\widehat{u}}_{\mathrm{rms}}(\omega )}},\end{array}$$(33)

where C(ω) corresponds to the second term in Eq. (13), and is defined by Eq. (B.28).

If we estimate the power spectral density of the mode as ${\widehat{\mathit{v}}}_{\mathrm{rms}}^2\sim {10}^5$ m² s⁻² Hz⁻¹ and that of the turbulent fluctuations as ${\widehat{u}}_{\mathrm{rms}}^2\sim 10$ m² s⁻² Hz⁻¹ (Fig. 2 Turck-Chièze et al. 2004), the above ratio becomes $\sqrt{{10}^5/10}\sim {10}^2$. We note that at this stage, the correlation term is about 1% of the leading term, which is barely enough to impact the asymmetries significantly. However, by considering the limiting case in the Cauchy-Schwarz inequality, we assumed that the mode velocity and the turbulent velocity were both optimally coherent and completely independent. This is not, however, the case; consequently we probably overestimated the importance of C(ω) by at least an order of magnitude. Therefore, this crude order of magnitude estimation indeed tends to support the conclusion that correlated turbulent perturbations can be disregarded when interpreting mode asymmetries in the velocity spectrum.

6. Comparison with observations

Observed properties of solar-like oscillations depend not only on the observable (velocity or intensity), but also on the specifics of each instruments. As far as velocity measurements are concerned, all instruments do not perform their Doppler observations on the same spectral line. Since different spectral lines form at different altitudes in the atmosphere, and since the properties of turbulence change throughout the atmosphere, mode properties, and especially mode asymmetries, may depend on the instrument.

In the following, we compare the results of our model to observations performed with the GONG network. We then focus on the dependence of the asymmetry profile B(ν) on the height at which the velocity spectrum is observed. We use the “numerical spectrum” model (see Sect. 2) to compare the asymmetry profiles B(ν) as observed by several instruments performing solar velocity spectrum measurements. Finally, we focus on the bias in the determination of mode eigenfrequencies entailed by mode asymmetry, whose understanding is of primary importance for accurate inference of mode properties.

6.1. Comparison with GONG observations

First – and in order to support the validity of our model in terms of mode asymmetry – we compare the asymmetries yielded by our model to those inferred from observations. We use the data points extracted from the spectrum analysis of Barban et al. (2004). The authors chose an equivalent, albeit different, set of parameters to fit the acoustic modes observed by the GONG network. Therefore, we reconstructed the shape of the modes point by point using the parameters extracted from their fit and fitted them again with the formula given by Eq. (26). We note that the modes analysed in their study are not radial, but have angular degrees ranging from l = 15 to l = 50. It is known, however, that the dependence of mode asymmetry on the angular degree is very weak (see for instance Duvall et al. 1993), who studied modes of angular degree up to l = 170) because the eigenfunctions associated to the acoustic modes are very weakly dependent on l close to the photosphere, as long as l is not too high. Thus, comparing our radial study to their non-radial observations remains relevant.

We showcase this comparison in Fig. 12. The parameter λ was adjusted for a better agreement with the observations. We obtained λ = 0.5 (see Eq. (20) for a definition of this parameter), which is approximately the value obtained by constraining mode amplitudes (Samadi & Goupil 2001). It is rather clear that we reproduce the main features of the asymmetry profile B(ν), especially its sense of variation on the different intervals (when the frequency range of the observations overlaps ours), as well as the positions and values of its local extrema. Alternatively, we show the same comparison in terms of the asymmetry parameter χ in the bottom panel of Fig. 12. This parameter is defined by Eq. (29), and is more robust when it comes to comparing theory and observations because it does not depend on the determination of the mode linewidths, which may introduce extra uncertainty in the determination of the observed asymmetries. In conclusion, the good agreement obtained between our model and observations show that the model developed in this paper is relevant to account for acoustic mode asymmetry quantitatively.

6.2. Dependence of asymmetry profile B(ν) on observation height

We consider three different observation heights in this paper associated, respectively, with the Michelson Doppler Imager (MDI), the GOLF instrument (both onboard the SOHO spacecraft), the Global Oscillation Network Group (ground-based), and the Helioseismic and Magnetic Imager (HMI; onboard the SDO spacecraft). We provide further information on the observation height of each of these instruments in Table 2. However, determining the formation height of a given absorption line is extremely difficult in that it does not depend only on the nature of the line (Fleck et al. 2011). Therefore, the values given in Table 2 are not to be considered as precise estimates but, rather, as approximate figures.

In Fig. 13, we compare the asymmetry profile B(ν) as it would be observed by the various instruments listed in Table 2. We control the height by tuning the radial coordinate r_o at which the Green’s function of the homogeneous wave equation is calculated (see Sect. 2 for more details). Mode asymmetry is only slightly dependent on the observation height. This is not so surprising since mode asymmetry is a global property of the modes. It is noticeable, however, that the difference between these asymmetry profiles is most prominent at high frequency (≳3.5 mHz). This is because by observing the modes at a lower height in the atmosphere, we effectively increase the contribution of atmospheric turbulence compared to the contribution of sub-photospheric turbulence on mode asymmetry. As we discussed in Sect. 5, it is at high frequency that the atmospheric turbulence has a more significant impact on mode asymmetry. Therefore, changing the observation height has more impact at high frequency than at low frequency.

6.3. Bias in eigenfrequency determination

In order to infer the stellar internal structure from asteroseismic measurements, it is important to determine the eigenfrequencies of the observed modes not only precisely, but also accurately. As was noted early on (Duvall et al. 1993; Abrams & Kumar 1996; Chaplin et al. 1999; Thiery et al. 2000; Toutain et al. 1998), using a symmetric, Lorentzian model to fit asymmetric line profiles results in an appreciable bias in the best-fit parameters, thus rendering inaccurate the frequency determination.

We illustrate this bias in the following. We fit the line profiles we obtained numerically with a formula given either by Eq. (26) or by the symmetric version (i.e. fixing B = 0). Following Abrams & Kumar (1996), we denote the reduced frequency bias as

$$\begin{array}{c}\hfill \delta \nu ={\displaystyle \frac{{\nu }_0^{B=0}-{\nu }_0^{B\ne 0}}{{\mathrm{\Gamma }}_{{\nu }_0}}}·\end{array}$$(34)

Figure 14 shows the results obtained for all radial modes described in Sect. 4. Given that the asymmetry of solar p-modes in the velocity spectrum are all negative, using a symmetric fit introduces an underestimation of ν₀. It is very clear that this underestimation grows linearly with the asymmetry B of the modes. We superimpose on Fig. 14 the best linear fit to the numerical data, given by: δν = 0.463 B + 0.000296. Given the typical values of δν, it is safe to say that the intercept of this linear regression is negligible.

This result can be easily interpreted by taking a closer look to the asymmetric fitting formula given by Eq. (26). Indeed, after applying elementary algebra, it becomes clear that for leading order in B, this expression reaches its maximum at x ∼ B (or ν_max = ν₀ + BΓ_ν0/2). Therefore, while the asymmetric fit accurately finds the eigenfrequency at ${\nu }_0^{B\ne 0}={\nu }_0$, the symmetric fit finds it at ${\nu }_0^{B=0}={\nu }_{\text{max}}$. One can therefore derive the simple expression δν = B/2, which is in accordance with the linear fit shown in Fig. 14.

The frequency bias introduced when not taking the line profile asymmetry into account can easily reach several percents of the mode linewidth. This bias largely exceeds the frequency resolution achievable in current solar measurements, especially for high-frequency modes, which are the widest.

7. Summary and conclusion

In this paper, we detail the development of a realistic and predictive model for the asymmetry displayed by solar radial p-modes in the velocity spectrum. The basic idea behind this model is to compute the Green’s function associated to the radial acoustic wave equation, as well as its inhomogeneous part (which corresponds to the source of excitation) and to convolve the two to reconstruct the velocity power spectral density point by point. Once the power spectral density is reconstructed, we extract its resonant modes and study their asymmetry. In particular, and unlike previous attempts to such modelling, we included in our model the correlation of the oscillating modes with the fluctuations associated to turbulent velocity.

First, the Green’s function associated with the wave equation was computed numerically. We put the wave equation in the form of a 1D stationary Schrödinger equation, whose potential only depends on the equilibrium structure of the Sun. We extracted the acoustic potential from a solar patched model: the solar interior is calculated using the 1D evolutionary code CESTAM and the solar atmosphere is calculated using the 3D hydrodynamic code CO⁵BOLD and horizontally averaged. We integrated the wave equation along the solar radius, and added a point-like, normalised source to the integration scheme in order to compute the Green’s function.

Secondly, the source term of the wave equation being of stochastic nature, we modelled the statistical properties of the source by means of theoretical developments. We made use of the adequate closure relation to express the third and fourth-order correlation products of the turbulent velocity as functions of second-order products; more specifically on their spatial and temporal Fourier transform. We developed two distinct models: one is based on theoretical prescriptions for the spatial and temporal spectrum of turbulent kinetic energy; the other is based on theoretical modelling of the temporal spectrum only, whereas the spatial spectrum is extracted from a 3D hydrodynamic simulation of the solar atmosphere. We refer to the former as the “theoretical spectrum” model, and to the latter as the “numerical spectrum” model.

The asymmetry B displayed by the modes in our model drastically depends on their frequency ν. This is because the shape of the eigenfunctions close to the photosphere is very dependent on ν. We find that B is negative throughout the p-mode spectrum, and that its behaviour weakly depends on the input parameters of our model. It drops from −0.01 to −0.05 between 1 mHz and 1.7 mHz, then rises to 0.015 at 3 mHz, and decreases again from 3 mHz to 4 mHz. Above 4 mHz, the behaviour of B(ν) is much more dependent on the value given to our input parameters and, in particular, on the injection scale associated to the turbulent cascade above the photosphere, compared to below the photosphere. This is related to the fact that the contribution of atmospheric turbulence to mode excitation only becomes significant at high frequency, so that only in this part of the spectrum it may have an impact on mode asymmetry.

The asymmetry of the modes can have two different origins: localisation of their source of excitation within a region of lesser spatial extent than the mode wavelength and correlation between the oscillating modes and the fluctuations associated to turbulent velocity. Formally, these two phenomena have the same impact on mode asymmetry, so that they cannot be separated using observational data only. Our model allows us to make this distinction and to study their relative weight in the total mode asymmetry. We find that the correlation with turbulent fluctuations is negligible in the velocity spectrum, and that the observed asymmetries are exclusively due to source localisation. More precisely, we interpret the results of our model in terms of the source position with respect to the various nodes and antinodes featured by the eigenfunctions of the modes. In the case of a point-like source of excitation, mode asymmetry drastically depends on whether it is located within or outside the mode acoustic cavity. In our model, however, the source of excitation has a certain spatial extent, so that the total asymmetry is a combination of the contributions from the source outside and inside the mode acoustic cavity.

We find that it is impossible to interpret even the qualitative behaviour of the asymmetry profile B(ν) by considering that the source of excitation is point-like (either outside or inside the modes cavity). On the contrary, taking into account the spatial extent of the source allows us to reproduce the observed asymmetries, not only qualitatively, but also quantitatively. This positive result shows that our model is indeed relevant to describe – and, more importantly, to predict – acoustic mode asymmetry in solar-like oscillators. It also shows that any model that assumes a point-like source of excitation cannot give reliable results as far as mode asymmetry is concerned. In particular, such a model would predict positive asymmetries for high-frequency modes, whereas observations show that all asymmetries are negative when measured in terms of velocity power spectral density.

Finally, we study the eigenfrequency bias entailed by neglecting to fit observations with an asymmetric profile. We find that for the most asymmetric modes, this bias can reach several percent of the mode linewidth. Therefore, this bias is higher for high frequency modes, which are the widest. In particular, for ν ∼ 4 mHz, the asymmetry parameter is of order B ∼ −0.04, and the linewidth is of order Γ ∼ 10 μHz, so that the eigenfrequency bias is of order δν ∼ 0.2 μHz. This is in perfect accordance with actual biases obtained from observation fit of the solar spectrum (see Benomar et al. 2018, Fig. 6, topmost panel). Since the eigenfrequency bias is most pronounced for higher frequency (because it is proportional to the mode line-width, which is widest at high frequency), it is likely to have a non-negligible impact on inversion methods, especially those based on asymptotic formulae. One must keep in mind, however, that the deviation of the modelled eigenfrequencies from the observed ones, induced by surface effects, largely dominates the eigenfrequency bias entailed by symmetric fits.

In this paper, we have restricted ourselves to the study of solar radial p-modes. Our formalism can be easily adapted to the study of non-radial modes simply by using a non-radial wave equation instead of the radial one. However, since the eigenfunction associated to p-modes in solar-like oscillators are very weakly independent on angular degree l close to the photosphere, which is precisely where the excitation takes place, the mode asymmetry is not expected to vary significantly with l, at least as long as l remains reasonably small. Observational data tend to confirm this (see e.g. Vorontsov & Jefferies 2013, who report that the spectral parameters of individual modes collapse to slowly varying functions of frequency only for modes with l ≲ 100).

We only considered one type of acoustic source in this study, that is, the turbulent fluctuations of the Reynolds stress. Indeed, it has been shown by Stein & Nordlund (2001) that this is the dominant source of excitation of solar acoustic modes (see also Chaplin et al. 2005; Samadi et al. 2007; Nordlund et al. 2009). Therefore, our objective was to start by considering only this source. However, further refinements of the model will have to include other sources of excitation, in the form of non-adiabatic, turbulent pressure fluctuations.

Our formalism can also be easily applied to other solar-like oscillators. Comparing the asymmetries featured by the velocity spectra of several solar-like oscillators as modelled by the method presented in this paper and, in particular, the trend followed by mode asymmetry with stellar parameters such as effective temperature or surface gravity, undoubtedly constitutes the next step of this study. In the long run, mode asymmetry may serve as a useful tool for seismic diagnoses of solar-like oscillators. However, the one major difference that remains between the solar case and other stars is that the Sun is the only solar-like oscillator for which spectra obtained by spectrometric measurements are sufficiently resolved to allow for a determination of their mode asymmetry. The asymmetry of acoustic modes of all other stars can only be observed in intensity spectra. As has been reported numerous times (see e.g. Duvall et al. 1993), asymmetry in intensity and in velocity spectra are drastically different. It is, therefore, necessary to adapt our formalism to the intensity spectrum, which is another key element of any further considerations on the matter treated here.

Acknowledgments

J.P. wishes to warmly thank Louis Manchon for having provided us with the solar patched model we used in this paper. J.P. and K.B. would also like to thank Marie-Jo Goupil for her thorough reading of the manuscript, as well as John Leibacher for useful discussions. J.P. and K.B. acknowledge financial support from the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU. H.G.L. acknowledges financial support by the Sonderforschungsbereich SFB 881 “The Milky Way System” (subprojects A4) of the German Research Foundation (DFG).

References

Appendix A: The inhomogeneous wave equation (Eq. (8))

A.1. Hydrodynamic equations and their linearisation

We linearise the governing, hydrodynamic equations in order to derive the wave equation with its source term. We consider that the mode velocity and the turbulent velocity obey separately their own continuity equation. Furthermore, we only consider radial modes, such that the mode velocity may we written in terms of the radial fluid displacement as v_osc = dξ_r/dte_r. In this context, the governing equations are as follows:

• the continuity equation associated to the mode velocity can be written as

  $$\begin{array}{c}\hfill {\displaystyle \frac{\partial \rho }{\partial t}}+\mathbf{\nabla }(\rho {\mathit{v}}_{\mathrm{osc}})=0.\end{array}$$(A.1)

  Writing ρ = ρ₀ + ρ′ (where ρ′ is the Eulerian density perturbation corresponding to the mode), linearising this equation around a motionless state, and integrating with respect to time yields

  $$\begin{array}{c}\hfill {\rho }^{\prime }+\mathbf{\nabla }({\rho }_0\mathit{\xi })=0.\end{array}$$(A.2)

  We then introduce the Lagrangian density perturbation δρ = ρ′+(ξ.∇)ρ₀, which allows us to write the linearised continuity equation in its final form:

  $$\begin{array}{c}\hfill {\displaystyle \frac{\delta \rho }{{\rho }_0}}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\mathrm{d}(r^2{\xi }_r)}{\mathrm{d}r}}=0;\end{array}$$(A.3)

• the Euler equation:

  $$\begin{array}{c}\hfill {\displaystyle \frac{\partial \rho \mathit{v}}{\partial t}}+\mathbf{\nabla }:(\rho \mathit{v}\mathit{v})=-\mathbf{\nabla }P+\rho \mathit{g}.\end{array}$$(A.4)

  Unlike what we did for the continuity equation, the velocity v now includes the mode velocity v_osc as well as the turbulent velocity u_turb. We further decompose the latter into a mean value U ≡ ⟨u_turb⟩ (where the notation ⟨.⟩ refers to an ensemble average) and fluctuations around this mean value u ≡ u_turb − U. As such, we have

  $$\begin{array}{c}\hfill \mathit{v}=\mathit{U}+{\mathit{v}}_{\mathrm{osc}}+\mathit{u}=\mathit{U}+\mathrm{d}{\xi }_r/\mathrm{d}t{\mathit{e}}_{\mathit{r}}+\mathit{u}.\end{array}$$(A.5)

  The last two terms are treated as small perturbations compared to the first one. In the term ∂ρv/∂t, the contribution of U vanishes because we consider that U is independent of time (in other words, we consider a stationary turbulence), and the contribution of u vanishes after ensemble averaging. Concerning the advection term, among the 9 terms of its development, only 2 survive after the linearisation and ensemble averaging, namely ∇  :  (ρUU) and ∇  :  (ρuu). The first one can be rewritten as ∇p_t, where p_t is the turbulent pressure, and is of order zero, so that it will only impact the equilibrium structure. The second one can be equivalently rewritten as $\mathit{\nabla }{p}_t^{\prime }$, where $p_t^{\prime }$ refers to the perturbation of the turbulent pressure. Finally, performing a Fourier transform with respect to time, the radial component of the Euler equation reads:

  $$\begin{array}{c}\hfill -{\omega }^2{\xi }_r+{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\mathrm{d}{p}^{\prime }}{\mathrm{d}r}}+{\displaystyle \frac{{\rho }^{\prime }}{{\rho }_0}}{g}_0-g^{\prime }=-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\mathrm{d}{p}_t^{\prime }}{\mathrm{d}r}},\end{array}$$(A.6)

  where p′ is the Eulerian pressure perturbation, g₀ is the mean gravitational acceleration, g′ is its Eulerian perturbation and $\text{d}{p}_t^{\prime }/\text{d}r$ refers to the turbulent fluctuations of the Reynolds stress around its mean value. Since we only aim at modelling radial modes, using the Cowling approximation to eliminate g′ would not reduce the order of the final wave equation, and is therefore of no particular use. Instead, we follow Unno et al. (1989) and express g′ as a function of the radial fluid displacement (see their Eq. (14.36)):

  $$\begin{array}{c}\hfill g^{\prime }=-{\displaystyle \frac{\mathrm{d}{\varphi }^{\prime }}{\mathrm{d}r}}=4\pi G{\rho }_0{\xi }_r·\end{array}$$(A.7)

  One can note that this is equivalent to saying that the Lagrangian perturbation of the gravitational potential is zero.

• the equation of state we will use to close the system: after some algebra, a linearised version of the equation of state in terms of the Lagrangian perturbations can be derived:

  $$\begin{array}{c}\hfill {\displaystyle \frac{\delta \rho }{{\rho }_0}}={\displaystyle \frac{1}{{\mathrm{\Gamma }}_1}}{\displaystyle \frac{\delta p}{p_0}}-{\displaystyle \frac{{\rho }_0{T}_0}{p_0}}{\mathrm{\nabla }}_{\mathrm{ad}}\delta s,\end{array}$$(A.8)

  where δs corresponds to the turbulent fluctuation of the specific entropy of the fluid and we define the various thermodynamic coefficients as

  $$\begin{array}{cc}& {\mathrm{\Gamma }}_1\equiv {\left({\displaystyle \frac{\partial lnp}{\partial ln\rho }}\right)}_s\hfill \\ \hfill & {\mathrm{\nabla }}_{\mathrm{ad}}\equiv {\left({\displaystyle \frac{\partial lnT}{\partial lnp}}\right)}_s·\hfill \end{array}$$(A.9)

  In order to facilitate the following calculations, we replace the Lagrangian pressure perturbation δp with the Eulerian one p′, and we derive two versions of the linearised equation of state, one with the Lagrangian density perturbation, one with the Eulerian one. Noting that the hydrostatic equilibrium gives us

  $$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}{p}_0}{\mathrm{d}r}}=-{\rho }_0{g}_0,\end{array}$$(A.10)

  and that by definition of the Brunt-Väisälä frequency, we have

  $$\begin{array}{c}\hfill {\displaystyle \frac{N^2}{g_0}}={\displaystyle \frac{1}{{\mathrm{\Gamma }}_1}}{\displaystyle \frac{\mathrm{d}ln{p}_0}{\mathrm{d}r}}-{\displaystyle \frac{\mathrm{d}ln{\rho }_0}{\mathrm{d}r}},\end{array}$$(A.11)

  and we finally obtain

  $$\begin{array}{l}\frac{\delta \rho }{{\rho }_0}=\frac{1}{{\Gamma }_1}\frac{p^{\prime }}{p_0}-\frac{g_0{\rho }_0}{{\Gamma }_1{p}_0}{\xi }_r-\frac{{\rho }_0{T}_0}{p_0}{\nabla }_{\text{ad}}\delta s\\ \frac{{\rho }^{\prime }}{{\rho }_0}=\frac{1}{{\Gamma }_1}\frac{p^{\prime }}{p_0}+\frac{N^2}{g_0}{\xi }_r-\frac{{\rho }_0{T}_0}{p_0}{\nabla }_{\text{ad}}\delta s.\end{array}$$(A.12)

A.2. Changing variables

The two variables that we wish to keep in these equations are ξ_r and p′. We first make use of Eq. (A.12) to eliminate the density fluctuations. Noting that c² = Γ₁p₀/ρ₀ (where c is the sound speed), the continuity and Euler equations then yield:

$$\begin{array}{cc}& {\displaystyle \frac{\mathrm{d}(r^2{\xi }_r)}{\mathrm{d}r}}-g_0{\displaystyle \frac{r^2}{c^2}}{\xi }_r+{\displaystyle \frac{r^2}{{\rho }_0{c}^2}}{p}^{\prime }=r^2{\mathrm{\nabla }}_{\mathrm{ad}}{\displaystyle \frac{{\rho }_0{T}_0}{p_0}}\delta s\hfill \\ \hfill & {\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\mathrm{d}{p}^{\prime }}{\mathrm{d}r}}+{\displaystyle \frac{g_0}{c^2}}{\displaystyle \frac{p^{\prime }}{{\rho }_0}}+(N^2-{\omega }^2-4\pi G{\rho }_0){\xi }_r={\displaystyle \frac{{\rho }_0{g}_0{T}_0}{p_0}}{\mathrm{\nabla }}_{\mathrm{ad}}\delta s-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\mathrm{d}{p}_t^{\prime }}{\mathrm{d}r}}·\hfill \end{array}$$(A.13)

In order to remove ξ_r from the 0th order term in the ξ_r equation, and same for p′, the required variable change is then (Unno et al. 1989):

$$\begin{array}{c}\hfill r^2{\xi }_r(r)=\stackrel{\sim }{\xi }(r)exp\left({\displaystyle {\int }_0^r{\displaystyle \frac{g_0}{c^2}}\mathrm{d}{r}^{\prime }}\right)\\ \hfill p^{\prime }={\rho }_0\stackrel{\sim }{\eta }(r)exp\left({\displaystyle {\int }_0^r{\displaystyle \frac{N^2}{g_0}}\mathrm{d}{r}^{\prime }}\right)·\end{array}$$(A.14)

Plugging this into Eq. ((A.13)), we obtain

$$\begin{array}{cc}& {\displaystyle \frac{\mathrm{d}\stackrel{\sim }{\xi }}{\mathrm{d}r}}+{\displaystyle \frac{r^2}{c^2}}exp\left({\displaystyle {\int }_0^r{\displaystyle \frac{N^2}{g_0}}-{\displaystyle \frac{g_0}{c^2}}\mathrm{d}{r}^{\prime }}\right)\stackrel{\sim }{\eta }\hfill \\ \hfill & =r^{2}exp\left({\displaystyle -{\int }_0^r{\displaystyle \frac{g_0}{c^2}}\mathrm{d}{r}^{\prime }}\right){\mathrm{\nabla }}_{\mathrm{ad}}{\displaystyle \frac{{\rho }_0{T}_0}{p_0}}\delta s,\hfill \end{array}$$(A.15)

and

$$\begin{array}{cc}& {\displaystyle \frac{\mathrm{d}\stackrel{\sim }{\eta }}{\mathrm{d}r}}+{\displaystyle \frac{1}{r^2}}exp\left({\displaystyle {\int }_0^r{\displaystyle \frac{g_0}{c^2}}-{\displaystyle \frac{N^2}{g_0}}\mathrm{d}{r}^{\prime }}\right)(N^2-{\omega }^2-4\pi G{\rho }_0)\stackrel{\sim }{\xi }\hfill \\ \hfill & =exp\left({\displaystyle -{\int }_0^r{\displaystyle \frac{N^2}{g_0}}\mathrm{d}{r}^{\prime }}\right)[{\mathrm{\nabla }}_{\mathrm{ad}}{\displaystyle \frac{{\rho }_0{g}_0{T}_0}{p_0}}\delta s-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\mathrm{d}{p}_t^{\prime }}{\mathrm{d}r}}],\hfill \end{array}$$(A.16)

where we denote the right-hand side terms of Eqs. (A.15) and (A.16) as S₀ and S₁ respectively in the following. We also define

$$\begin{array}{cc}& I(r)\equiv exp\left({\displaystyle {\int }_0^r{\displaystyle \frac{N^2}{g_0}}-{\displaystyle \frac{g_0}{c^2}}\mathrm{d}{r}^{\prime }}\right)\hfill \\ \hfill & x(r)\equiv {\displaystyle \frac{r\sqrt{I}}{c}}\hfill \\ \hfill & k^2\equiv {\displaystyle \frac{{\omega }^2-N^2+4\pi G{\rho }_0}{c^2}}·\hfill \end{array}$$(A.17)

The above set of equations can be rewritten as

$$\begin{array}{l}\frac{\text{d}\stackrel{\sim }{\xi }}{\text{d}r}+x^2\stackrel{\sim }{\eta }=S_0\\ \frac{\text{d}\stackrel{\sim }{\eta }}{\text{d}r}-\frac{k^2}{x^2}\stackrel{\sim }{\xi }=S_1.\end{array}$$(A.18)

We can now eliminate $\stackrel{\sim }{\eta }$ to get a single second-order wave equation. Using the first of Eq. (A.18) to express $\stackrel{\sim }{\eta }$ as a function of $\stackrel{\sim }{\xi }$, and plugging it in the second equation, we get the following equation:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^2\stackrel{\sim }{\xi }}{\mathrm{d}{r}^2}}-{\displaystyle \frac{2}{x}}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}r}}{\displaystyle \frac{\mathrm{d}\stackrel{\sim }{\xi }}{\mathrm{d}r}}+k^2\stackrel{\sim }{\xi }={\displaystyle \frac{\mathrm{d}{S}_0}{\mathrm{d}r}}-{\displaystyle \frac{2}{x}}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}r}}{S}_0-x^2{S}_1·\end{array}$$(A.19)

Similarly to what has been done for the first change of variables, we wish for the left-hand side to contain no first-order term, but only second-order and 0th-order ones. Thus we introduce yet another variable: $\mathrm{\Psi }(r)\equiv \stackrel{\sim }{\xi }/x$. Plugging this new variable into Eq. (A.19), we easily obtain a wave equation that assumes the form of a 1D stationary Schrödinger equation

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^2\mathrm{\Psi }}{\mathrm{d}{r}^2}}+({\displaystyle \frac{{\omega }^2}{c^2}}-V(r))\mathrm{\Psi }={\displaystyle \frac{1}{x}}({\displaystyle \frac{\mathrm{d}{S}_0}{\mathrm{d}r}}-{\displaystyle \frac{2}{x}}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}r}}{S}_0-x^2{S}_1),\end{array}$$(A.20)

with an acoustic potential V(r) that only depends on the star’s equilibrium state:

$$\begin{array}{c}\hfill V(r)={\displaystyle \frac{N^2-4\pi G{\rho }_0}{c^2}}+{\displaystyle \frac{2}{x^2}}{\left({\displaystyle \frac{\mathrm{d}x}{\mathrm{d}r}}\right)}^2-{\displaystyle \frac{1}{x}}{\displaystyle \frac{{\mathrm{d}}^{2}x}{\mathrm{d}{r}^2}}·\end{array}$$(A.21)

A.3. The source term

With the above notations, the parameters intervening in the source term of Eq. (A.20) have the following expressions:

$$\begin{array}{cc}& S_0(r)=r^{2}exp\left({\displaystyle -{\int }_0^r{\displaystyle \frac{g_0}{c^2}}\mathrm{d}{r}^{\prime }}\right){\mathrm{\nabla }}_{\mathrm{ad}}{\displaystyle \frac{{\rho }_0{T}_0}{p_0}}\delta s\hfill \\ \hfill & S_1(r)=exp\left({\displaystyle -{\int }_0^r{\displaystyle \frac{N^2}{g_0}}\mathrm{d}{r}^{\prime }}\right)[{\mathrm{\nabla }}_{\mathrm{ad}}{\displaystyle \frac{{\rho }_0{g}_0{T}_0}{p_0}}\delta s-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\mathrm{d}{p}_t^{\prime }}{\mathrm{d}r}}]\hfill \\ \hfill & x(r)={\displaystyle \frac{r}{c}}exp\left({\displaystyle {\displaystyle \frac{1}{2}}{\int }_0^r{\displaystyle \frac{N^2}{g_0}}-{\displaystyle \frac{g_0}{c^2}}\mathrm{d}{r}^{\prime }}\right)·\hfill \end{array}$$(A.22)

Furthermore, one can easily derive the following relationship between ∇_ad and α_s ≡ (∂P/∂s)_ρ by means of the adequate Schwarz relation:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mathrm{ad}}{\displaystyle \frac{{\rho }_0{T}_0}{p_0}}={\displaystyle \frac{{\alpha }_s}{{\rho }_0{c}^2}}·\end{array}$$(A.23)

After some manipulations, one finally obtain the source term in the form:

$$\begin{array}{cc}\hfill S(r)=& {\displaystyle \frac{r}{c{\rho }_0}}exp\left({\displaystyle -{\displaystyle \frac{1}{2}}{\int }_0^r{\displaystyle \frac{N^2}{g_0}}+{\displaystyle \frac{g_0}{c^2}}\mathrm{d}{r}^{\prime }}\right)\hfill \\ \hfill & \times [{\alpha }_s\delta s{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}ln\left({\displaystyle \frac{{\alpha }_s\delta s}{{\rho }_0}}\right)-{\alpha }_s\delta s({\displaystyle \frac{N^2}{g_0}}+{\displaystyle \frac{g_0}{c^2}})+{\displaystyle \frac{\mathrm{d}{p}_t^{\prime }}{\mathrm{d}r}}]·\hfill \end{array}$$(A.24)

Finally, since

$$\begin{array}{c}\hfill {\displaystyle \frac{N^2}{g_0}}+{\displaystyle \frac{g_0}{c^2}}=-{\displaystyle \frac{\mathrm{d}ln{\rho }_0}{\mathrm{d}r}},\end{array}$$(A.25)

this expression can be drastically simplified to

$$\begin{array}{c}\hfill S(r)={\displaystyle \frac{r}{c\sqrt{{\rho }_0(r=0){\rho }_0(r)}}}({\displaystyle \frac{\mathrm{d}({\alpha }_s\delta s)}{\mathrm{d}r}}+{\displaystyle \frac{\mathrm{d}{p}_t^{\prime }}{\mathrm{d}r}})·\end{array}$$(A.26)

This form clearly shows that the source term can be split three ways: a monopolar source term (proportional to δsdα_s/dr) due to non-adiabatic pressure fluctuations in a stratified environment, a dipolar term (proportional to α_sdδs/dr) due to a stratification in the non-adiabatic pressure fluctuations themselves, and a quadripolar term (proportional to $\text{d}{p}_t^{\prime }/\text{d}r$) due to Reynolds stress fluctuations. In the following, we only consider this last term but we also show here how the effect of non-adiabaticity can be introduced as well.

To conclude, note that the value of the fluid density at the centre of the star ρ₀(r = 0) appears both in the definition of the variable Ψ and in the source term S(r). This is due to the particular change of variable we have performed, and it can be factored out of the wave equation. Finally, we can put the wave equation in the following form:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^2\mathrm{\Psi }}{\mathrm{d}{r}^2}}+({\displaystyle \frac{{\omega }^2}{c^2}}-V(r))\mathrm{\Psi }=S(r),\end{array}$$(A.27)

with V(r) given by Eq. (A.21), and the source term and wave variable are given by

$$\begin{array}{cc}& S(r)={\displaystyle \frac{r}{c\sqrt{{\rho }_0(r)}}}({\displaystyle \frac{\mathrm{d}({\alpha }_s\delta s)}{\mathrm{d}r}}+{\displaystyle \frac{\mathrm{d}{p}_t^{\prime }}{\mathrm{d}r}})\hfill \\ \hfill & \mathrm{\Psi }(r)=rc(r)\sqrt{{\rho }_0(r)}{\xi }_r(r).\hfill \end{array}$$(A.28)

Appendix B: From the Green’s function to the power spectral density

Here we detail the calculations carried out to obtain the expression of the velocity power spectral density (Eq. (13)) as a function of the Green’s function associated with the homogeneous wave Eq. (8). Note that these calculations correspond to the “theoretical spectrum” model described in Sect. 2.3.2. The calculations in the “numerical spectrum” model being fairly similar, we do not detail it. We start with the development given by Eq. (4), with the expression of $\widehat{{\mathit{v}}_{\mathrm{osc}}}$ given by Eq. (12). We detail the treatment of both terms in the development (4) (leading term and cross term) separately.

B.1. The leading term

For more clarity, in the following, we introduce

$$\begin{array}{c}\hfill X_{\omega }(\mathit{r})\equiv G_{\omega }(\mathit{r}){\displaystyle \frac{||\mathit{r}||}{c(\mathit{r})\sqrt{{\rho }_0(\mathit{r})}}}·\end{array}$$(B.1)

We then have

$$\begin{array}{cc}& \langle {\left|\widehat{v_{\mathrm{osc}}}(\omega )\right|}^2\rangle \hfill \\ \hfill & =\langle {\displaystyle {\displaystyle \frac{{\omega }^2}{r_{\mathrm{o}}^2{c}_{\mathrm{o}}^2{\rho }_0(r_{\mathrm{o}})}}\times \int \int {\mathrm{d}}^3{\mathit{r}}_{\mathit{s}\mathbf{1}}{\mathrm{d}}^3{\mathit{r}}_{\mathit{s}\mathbf{2}}(\mathbf{\nabla }{X}_{\omega }.{\rho }_0\widehat{u_r\mathit{u}})({\mathit{r}}_{\mathit{s}\mathbf{1}})}\hfill \\ \hfill & \times (\mathbf{\nabla }{X}_{\omega }^{\star }.{\rho }_0{\widehat{u_r\mathit{u}}}^{\star })({\mathit{r}}_{\mathit{s}\mathbf{2}})\rangle ,\hfill \end{array}$$(B.2)

where r_o is the radius at which the spectrum is observed, c_o is the speed of sound at that radius, and the notation ⋆ refers to the complex conjugate.

We then perform the following change of variable: R = (r_s1 + r_s2)/2 and r = (r_s1 − r_s2)/2, the former being a “slow” variable, and the latter a “fast” variable. This allows us to separate the scales relevant to the turbulent velocity u from the scales relevant to the medium stratification and the mode wavelength, with turbulent quantities only relevant in the r scale and the stratification and Green function only relevant in the R scale. The scale separation approximation is not realistic in the subsurface layers (in particular, the mode wavelength is comparable to the typical correlation length associated with turbulence); however, for want of a better alternative, we are led to use this approximation in the following.

Therefore, we make the assumption that the second-order correlation product of the turbulent velocity vanishes for lengths much shorter than the scale associated to the variations of the equilibrium structure. Being able to separate the two scales, as well as the fact that, for radial modes, X_ω only depends on the radial coordinate, allows us to rewrite the leading term as

$$\begin{array}{cc}& \langle {\left|\widehat{v_{\mathrm{osc}}}(\omega )\right|}^2\rangle ={\displaystyle \frac{{\omega }^2}{r_{\mathrm{o}}^2{c}_{\mathrm{o}}^2{\rho }_0(r_{\mathrm{o}})}}\hfill \\ \hfill & {\displaystyle \times \int \mathrm{d}m{\left|{\displaystyle \frac{\mathrm{d}{X}_{\omega }}{\mathrm{d}R}}\right|}^2{\rho }_0(R)\int {\mathrm{d}}^3\mathit{r}\langle \widehat{u_r^2}(\mathbf{0},\omega ){\widehat{u_r^2}}^{\star }(\mathit{r},\omega )\rangle ,}\hfill \end{array}$$(B.3)

where we have dropped the variable R in favor of the more practical mass variable m. We note that we can only perform this change of variable because the wave equation is radial so that the function X_ω(r) only depends on the radial coordinate.

In the following, we focus on establishing the expression of the integral over the fast variable r. By definition of the temporal Fourier transform appearing in said integral, we have

$$\begin{array}{cc}& {\displaystyle \int {\mathrm{d}}^3\mathit{r}\langle \widehat{u_r^2}(\mathbf{0},\omega ){\widehat{u_r^2}}^{\star }(\mathit{r},\omega )\rangle }\hfill \\ \hfill & {\displaystyle ={\displaystyle \frac{1}{{(2\pi )}^2}}\int \int {\mathrm{d}}^3\mathit{r}\mathrm{d}\tau e^{-j\omega \tau }\langle u_r^2(\mathbf{0},0)u_r^2(\mathit{r},\tau )\rangle .}\hfill \end{array}$$(B.4)

We then use the Quasi-Normal Approximation (hereby abbreviated QNA), under which any fourth-order correlation product can be decomposed into a sum of three second-order correlation products, so that (Lesieur 2008)

$$\begin{array}{c}\hfill \langle u_r^2(\mathbf{0},0)u_r^2(\mathit{r},\tau )\rangle =2{\langle u_r(\mathbf{0},0)u_r(\mathit{r},\tau )\rangle }^2+{\langle u_r(\mathbf{0},0)\rangle }^2{\langle u_r(\mathit{r},\tau )\rangle }^2.\end{array}$$(B.5)

The last term does not depend on τ or r if the turbulence is homogeneous and uniform, and thus yields zero when the Fourier transform is performed. We can then write

$$\begin{array}{cc}& {\displaystyle \int {\mathrm{d}}^3\mathit{r}\langle \widehat{u_r^2}(\mathbf{0},\omega ){\widehat{u_r^2}}^{\star }(\mathit{r},\omega )\rangle }\hfill \\ \hfill & {\displaystyle ={\displaystyle \frac{2}{{(2\pi )}^2}}\int \int {\mathrm{d}}^3\mathit{r}\mathrm{d}\tau {\mathrm{e}}^{-j\omega \tau }{\langle u_r(\mathbf{0},0)u_r(\mathit{r},\tau )\rangle }^2.}\hfill \end{array}$$(B.6)

Using the Parseval identity, we can express this as an integral over wave vectors k and angular frequencies ω

$$\begin{array}{cc}& {\displaystyle \int {\mathrm{d}}^3\mathit{r}\langle \widehat{u_r^2}(\mathbf{0},\omega ){\widehat{u_r^2}}^{\star }(\mathit{r},\omega )\rangle }\hfill \\ \hfill & =2\times {(2\pi )}^2\hfill \\ \hfill & {\displaystyle \times \int \int {\mathrm{d}}^3\mathit{k}\mathrm{d}{\omega }^{\prime }\mathrm{TF}\left[{\mathrm{e}}^{-j\omega \tau }\langle u_r(\mathbf{0},0)u_r(\mathit{r},\tau )\rangle \right]}\hfill \\ \hfill & \times \mathrm{TF}\left[\langle u_r(\mathbf{0},0)u_r(\mathit{r},\tau )\rangle \right],\hfill \end{array}$$(B.7)

where the notation TF[.] refers to temporal and spatial Fourier transform.

We then proceed to describe the second-order correlation product not in terms of time and space increments, but in terms of angular frequencies ω and spatial modes k. We denote the temporal and spatial Fourier transform of the second-order correlation product of the ith and jth component of the turbulent velocity as ϕ_ij(k, ω), so that

$$\begin{array}{cc}& {\displaystyle \int {\mathrm{d}}^3\mathit{r}\langle \widehat{u_r^2}(\mathbf{0},\omega ){\widehat{u_r^2}}^{\star }(\mathit{r},\omega )\rangle }\hfill \\ \hfill & {\displaystyle =8{\pi }^2\int \int {\mathrm{d}}^3\mathit{k}\mathrm{d}{\omega }^{\prime }{\varphi }_{\mathit{rr}}(\mathit{k},{\omega }^{\prime }-\omega ){\varphi }_{\mathit{rr}}(\mathit{k},{\omega }^{\prime }).}\hfill \end{array}$$(B.8)

For isotropic turbulence, ϕ_ij can be expressed analytically (Batchelor 1953) as

$$\begin{array}{c}\hfill {\varphi }_{\mathit{ij}}={\displaystyle \frac{E(k,\omega )}{4\pi k^2}}({\delta }_{\mathit{ij}}-{\displaystyle \frac{k_i{k}_j}{k^2}}),\end{array}$$(B.9)

where E(k, ω) is the specific turbulent kinetic energy spectrum, k, k_i and k_j are the norm, ith component and jth component of the wave vector k, and δ_ij is the Kronecker symbol. The integration over the solid angle of k is straightforward and only an integral over its norm remains. However, solar turbulence close to the photosphere is known to be highly anisotropic. To take this anisotropy into account, we follow the formalism developed by Gough (1977). In this formalism, the integral over the solid angle of k is simply readjusted by adding an anisotropy factor G, given by (see Appendix B in Samadi & Goupil 2001)

$$\begin{array}{c}\hfill {\displaystyle G={\int }_{-1}^1\mathrm{d}\mu {(1-{\displaystyle \frac{Q^2{\mu }^2}{(Q^2-1){\mu }^2+1}})}^2,}\end{array}$$(B.10)

where

$$\begin{array}{c}\hfill Q^2={\displaystyle \frac{\langle u_x^2\rangle }{\langle u_r^2\rangle }}={\displaystyle \frac{\langle u_y^2\rangle }{\langle u_r^2\rangle }},\end{array}$$(B.11)

u_x and u_y referring to the two horizontal components of the turbulent velocity.

This anisotropy factor depends on the ratio between horizontal and vertical turbulent velocities, and therefore depends on the slow R variable – or equivalently, on the mass variable m. Under this formalism, the integral over k and ω remains the same as in the isotropic case.

Following Stein (1967), we decompose E(k, ω) into a spatial part E(k), which describes how the turbulent kinetic energy is distributed among modes of different wave numbers, and a temporal part χ_k(ω), which describes the statistical life-time distribution of eddies of wavenumber k

$$\begin{array}{c}\hfill E(k,\omega )=E(k){\chi }_k(\omega ).\end{array}$$(B.12)

Finally, the integral given by Eq. (B.8) can be rewritten as

$$\begin{array}{cc}& {\displaystyle \int {\mathrm{d}}^3\mathit{r}\langle \widehat{u_r^2}(\mathbf{0},\omega ){\widehat{u_r^2}}^{\star }(\mathit{r},\omega )\rangle }\hfill \\ \hfill & {\displaystyle =2\pi G\int \mathrm{d}k{\displaystyle \frac{E{(k)}^2}{k^2}}\int \mathrm{d}{\omega }^{\prime }{\chi }_k({\omega }^{\prime }-\omega ){\chi }_k({\omega }^{\prime }).}\hfill \end{array}$$(B.13)

As mentioned in the main body of the paper, we have followed two different leads to model the functions E(k) and χ_k(ω) in this study. In the following, we only detail what we refer to as the “theoretical spectrum” model, which is based on theoretical prescriptions.

Based on the assumption that turbulent flows are self-similar, Kolmogorov’s theory of turbulence leads to a spatial spectrum E(k)∝k^−5/3 in the inertial range, between k = k₀ (where k₀ is the scale at which the kinetic energy is injected in the turbulent cascade, and is henceforth referred to as the injection scale) and the dissipation scale (at which the turbulent kinetic energy is converted into heat). Given the very high Reynolds number characterising solar turbulence (Re ∼ 10¹⁴), we cast the dissipation scale to infinity. Then, following Musielak et al. (1994), we extend the turbulent spectrum below the injection scale by considering that E(k) takes a constant value for k <  k₀. This extended spectrum, referred to as the BKS was introduced to account for the broadness of the maximum of E(k). Finally, the BKS can be written thus:

$$\begin{array}{c}\hfill E(k)=\{\begin{array}{cc}0.652{\displaystyle \frac{u_0^2}{k_0}}\hfill & \mathrm{if}0.2k_0<k<k_0\hfill \\ 0.652{\displaystyle \frac{u_0^2}{k_0}}{\left({\displaystyle \frac{k}{k_0}}\right)}^{-5/3}\hfill & \mathrm{if}{k}_0<k,\hfill \end{array}\end{array}$$(B.14)

where $u_0^2\equiv \langle {\mathit{u}}^2(\mathit{r})\rangle /3$ and the 0.652 factor is introduced so that the total specific kinetic energy of the turbulent spectrum matches $u_0^2/2$.

Following Samadi et al. (2003), we consider a Lorentzian shape for the temporal spectrum χ_k(ω), which is supported both by numerical simulations (Samadi et al. 2003) and by theoretical arguments (if the noise is characterised by a time-correlation function which decays exponentially, its spectrum is Lorentzian). Thus:

$$\begin{array}{c}\hfill {\chi }_k(\omega )={\displaystyle \frac{1}{\pi {\omega }_k}}{\displaystyle \frac{1}{1+{(\omega /{\omega }_k)}^2}}·\end{array}$$(B.15)

The width of the Lorentzian is the inverse of the typical correlation time-scale and by dimensional arguments, it is proportional to ku_k, where u_k is the typical velocity of eddies of wavenumber k. However, there remains a substantial indetermination on the actual value of ω_k. To account for this indetermination, we follow Balmforth (1992) and introduce the dimensionless parameter λ, so that

$$\begin{array}{c}\hfill {\omega }_k=2k{u}_k/\lambda .\end{array}$$(B.16)

For a Kolmogorov spectrum, u_k scales as k^−1/3, which means that we have

$$\begin{array}{c}\hfill {\omega }_k={\omega }_{k_0}{\left({\displaystyle \frac{k}{k_0}}\right)}^{2/3}\equiv {\displaystyle \frac{2k_0{u}_{k_0}}{\lambda }}{\left({\displaystyle \frac{k}{k_0}}\right)}^{2/3},\end{array}$$(B.17)

There is a temptation to approximate the typical velocities of eddies of wavenumber k₀ with u₀. This assumption requires some discussion, however. Indeed, Stein (1967) has pointed out that eddies of all sizes have the same Eulerian velocity fluctuations u₀. As far as Lagrangian fluctuations go, the fluctuations u_k can be expressed as (Stein 1967)

$$\begin{array}{c}\hfill {\displaystyle u_k^2={\int }_k^{2k}\mathrm{d}kE(k).}\end{array}$$(B.18)

Using the expression of E(k) given in Eq. (B.14) and applying it to k = k₀, we finally find u_k0 = 0.602u₀. Under all these assumptions, all further calculations being carried out, we ultimately obtain the leading term,

$$\begin{array}{cc}& {\displaystyle \langle {\left|\widehat{v_{\mathrm{osc}}}(\omega )\right|}^2\rangle ={\displaystyle \frac{0.353\lambda {\omega }^2}{r_{\mathrm{o}}^2{c}_{\mathrm{o}}^2{\rho }_0(r_{\mathrm{o}})}}\int \mathrm{d}m[{\rho }_{0}G{\left|{\displaystyle \frac{\mathrm{d}{X}_{\omega }}{\mathrm{d}r}}\right|}^2{\displaystyle \frac{u_0^3}{k_0^4}}}\hfill \\ \hfill & \left({\displaystyle {\int }_{0.2}^1{f}_1(K)\mathrm{d}K+{\int }_1^{\infty }{f}_2(K)\mathrm{d}K}\right)],\hfill \end{array}$$(B.19)

with

$$\begin{array}{c}\hfill \begin{array}{c}{f}_1(K)={\displaystyle \frac{K^{-8/3}}{1+{\left({\displaystyle \frac{\lambda \omega }{2.408u_0{k}_0}}\right)}^2{K}^{-4/3}}}\hfill \\ f_2(K)={\displaystyle \frac{K^{-6}}{1+{\left({\displaystyle \frac{\lambda \omega }{2.408u_0{k}_0}}\right)}^2{K}^{-4/3}}},\hfill \end{array}\end{array}$$(B.20)

and where K is the reduced inverse eddy scale (K ≡ k/k₀). We note that all the terms appearing in the integrand depend on the mass variable m, particularly u₀ and k₀, even when the dependence does not appear explicitly. Free parameters are left in this description of the leading term in the form of λ and k₀(m); solar turbulence close to the photosphere being as poorly constrained as it is today, we cannot hope to achieve a non-parametrised description of stochastically excited modes of oscillation.

B.2. The cross term

Similarly to the leading term, the cross term in the development of P(ω) shown in Eq. (4) can be written as

$$\begin{array}{cc}& \mathrm{Re}\left({\displaystyle \int \mathrm{d}\mathrm{\Omega }\stackrel{\sim }{h}(\mu )\langle \widehat{v_{\mathrm{osc}}}(\omega ){\widehat{u_n}}^{\star }(\omega )\rangle }\right)\hfill \\ \hfill & =-{\displaystyle \frac{\omega }{r_{\mathrm{o}}{c}_{\mathrm{o}}\sqrt{{\rho }_0(r_{\mathrm{o}})}}}\hfill \\ \hfill & {\displaystyle \times \int \mathrm{d}\mathrm{\Omega }\stackrel{\sim }{h}(\mu )\mathrm{Re}\left({\displaystyle j\int {\mathrm{d}}^3{\mathit{r}}_{\mathrm{s}}{\rho }_0{\displaystyle \frac{\mathrm{d}{X}_{\omega }}{\mathrm{d}r}}\langle \widehat{u_r^2({\mathit{r}}_{\mathrm{s}})}{\widehat{u_n({\mathit{r}}_{\mathrm{o}})}}^{\star }\rangle }\right),}\hfill \end{array}$$(B.21)

where μ is the cosine of the angle between the local vertical direction and the direction of the line of sight.

We note that this time, one of the velocities in the correlation product is estimated at a fixed location corresponding to the observation height, so that only one variable is left. Expressing the line-of-sight component of u as u_n = u_r cos θ − u_θ sin θ, this transforms into

$$\begin{array}{cc}& \mathrm{Re}\left({\displaystyle \int \mathrm{d}\mathrm{\Omega }\stackrel{\sim }{h}(\mu )\langle \widehat{v_{\mathrm{osc}}}(\omega ){\widehat{u_n}}^{\star }(\omega )\rangle }\right)\hfill \\ \hfill & =-{\displaystyle \frac{\omega }{r_{\mathrm{o}}{c}_{\mathrm{o}}\sqrt{{\rho }_0(r_{\mathrm{o}})}}}\times [{\displaystyle \mathrm{Re}\left({\displaystyle j\int \mathrm{d}m{\displaystyle \frac{\mathrm{d}{X}_{\omega }}{\mathrm{d}r}}\langle \widehat{u_r^2({\mathit{r}}_{\mathrm{s}})}{\widehat{u_r({\mathit{r}}_{\mathrm{o}})}}^{\star }\rangle }\right)\int \mathrm{d}\mathrm{\Omega }\stackrel{\sim }{h}(\mu )\mu }\hfill \\ \hfill & {\displaystyle +\mathrm{Re}\left({\displaystyle j\int \mathrm{d}m{\displaystyle \frac{\mathrm{d}{X}_{\omega }}{\mathrm{d}r}}\langle \widehat{u_r^2({\mathit{r}}_{\mathrm{s}})}{\widehat{u_{\theta }({\mathit{r}}_{\mathrm{o}})}}^{\star }\rangle }\right)\int \mathrm{d}\mathrm{\Omega }\stackrel{\sim }{h}(\mu )\sqrt{1-{\mu }^2}}].\hfill \end{array}$$(B.22)

Since u_θ is a horizontal component of the turbulent velocity, if we consider there is no preferential horizontal direction as far as turbulence goes, the third-order correlation product appearing in the second integral cancels out, so that we are left with the first integral only. The latter can be rewritten thus:

$$\begin{array}{cc}& {\displaystyle \int \mathrm{d}m{\displaystyle \frac{\mathrm{d}X}{\mathrm{d}r}}\langle \widehat{u_r^2({\mathit{r}}_{\mathrm{s}})}{\widehat{u_r({\mathit{r}}_{\mathrm{o}})}}^{\star }\rangle }\hfill \\ \hfill & {\displaystyle =\int {\mathrm{d}}^3{\mathit{r}}_{\mathrm{s}\mathbf{1}}{\mathrm{d}}^3{\mathit{r}}_{\mathrm{s}\mathbf{2}}{\rho }_0({\mathit{r}}_{\mathrm{s}\mathbf{1}}){\displaystyle \frac{\mathrm{d}X}{\mathrm{d}r}}({\mathit{r}}_{\mathrm{s}\mathbf{1}})\langle \widehat{u_r^2({\mathit{r}}_{\mathrm{s}\mathbf{1}})}{\widehat{u_r({\mathit{r}}_{\mathrm{s}\mathbf{2}})}}^{\star }\rangle \delta ({\mathit{r}}_{\mathrm{s}\mathbf{2}}-{\mathit{r}}_{\mathrm{o}}),}\hfill \end{array}$$(B.23)

where we have artificially introduced a second spatial variable r_s2, so as to get an expression formally similar to that of the leading term above. Performing the same change of variables, and plugging the definition of the temporal Fourier transform, we have

$$\begin{array}{cc}& {\displaystyle \int \mathrm{d}m{\displaystyle \frac{\mathrm{d}X}{\mathrm{d}r}}\langle \widehat{u_r^2({\mathit{r}}_{\mathrm{s}})}{\widehat{u_r({\mathit{r}}_{\mathrm{o}})}}^{\star }\rangle }\hfill \\ \hfill & {\displaystyle ={\displaystyle \frac{1}{{(2\pi )}^2}}{\rho }_0({\mathit{r}}_{\mathrm{o}}){\displaystyle \frac{\mathrm{d}X}{\mathrm{d}r}}({\mathit{r}}_{\mathrm{o}})\int \mathrm{d}\tau {\mathrm{d}}^3\mathit{r}{e}^{-j\omega \tau }\langle u_r^2(0,\mathbf{0})u_r(\tau ,\mathit{r})\rangle .}\hfill \end{array}$$(B.24)

The challenge in estimating the contribution of correlated turbulent perturbations fundamentally lies in the correct determination of non-local, third-order correlation products. While the QNA provides an adequate closure relationship for fourth-order moments, it yields vanishing third-order moments, which is known to lead to serious violations of the energy conservation principle, as well as an impossibility for the turbulent cascade to develop (see Kraichnan 1957).

In order to estimate the third-order moments, we follow Belkacem et al. (2006b) and use the Plume closure model, which consists of separating the flow into upward flows and downward plumes, each normally distributed, with different mean values and standard deviations. In addition, we consider that the downflows are much more turbulent than the upflows (which is supported by Goode et al. 1998, according to whom the intergranular lanes harbour stronger turbulence than the granules themselves at the Sun’s surface), and that the two types of flows considered separately have zero third-order correlation products. In Belkacem et al. (2006b), the authors use the same approximations but focused on one-point correlation products; however, the calculations can be easily extended to the two-point correlation products that we need and the model yields

$$\begin{array}{cc}\hfill ⟨u_r{(\mathit{R},t)}^2{u}_r(\mathit{R}+\mathit{r},t+\tau )⟩=& [a{(1-a)}^3-a^3(1-a)]\delta u^3\hfill \\ \hfill & -a(1-a)[2⟨\stackrel{\sim }{u_d}(\mathit{R},t)\stackrel{\sim }{u_d}(\mathit{R}+\mathit{r},t+\tau )⟩\hfill \\ \hfill & +⟨\stackrel{\sim }{u_d}{(\mathit{R},t)}^2⟩]\delta u,\hfill \end{array}$$(B.25)

where a is the relative horizontal section of the upflows, δu is the difference between the mean velocity of the two types of flows (considering their respective signs, it actually is the sum of their absolute values), and $\stackrel{\sim }{u_d}$ is the fluctuation of the downflow velocity around its mean value. The only additional approximation we make to adapt these calculations to two-point correlation products is that the parameters of the model a and δu vary on scales much greater than the typical correlation length, which is another illustration of the scale separation approximation, which we have already used for the leading term (see above).

Note that strictly speaking, the third-order moment appearing in Eq. (B.25) and yielded by the PCM are centred. However, we consider that the mean value of the overall vertical velocity of the flow is sufficiently low (compared to its standard deviation for instance) to be neglected. Therefore, the moment described by Eq. (B.25) may interchangeably refer either to a centred or non-centred moment.

Also note that this closure relation is written here in terms of $\stackrel{\sim }{u_d}$ (i.e. the turbulent fluctuations in the downflows only). It would be more practical to rewrite it in terms of u_r (i.e. the total turbulent fluctuations). The two are related through

$$\begin{array}{c}\hfill ⟨\stackrel{\sim }{u_d}(\mathit{R},t)\stackrel{\sim }{u_d}(\mathit{R}+\mathit{r},t+\tau )⟩={\displaystyle \frac{1}{1-a}}⟨u_r(\mathit{R},t)u_r(\mathit{R}+\mathit{r},t+\tau )⟩-a\delta u^2.\end{array}$$(B.26)

Plugging Eq. (B.26) into Eq. (B.25), we obtain a closure model that allows us to write the third-order moments as a function of second-order moments only, after which we can use the same prescriptions for turbulence as we did for the leading term. We note that while this closure relation contains many terms, not many survive the Fourier transform in Eq. (B.24) as all terms not depending on r or τ will not contribute to the Fourier transform. The integral over m appearing in Eq. (B.21) becomes

$$\begin{array}{c}\hfill {\displaystyle \int \mathrm{d}m{\displaystyle \frac{\mathrm{d}X}{\mathrm{d}r}}\langle \widehat{u_r^2({\mathit{r}}_{\mathit{s}})}{\widehat{u_r({\mathit{r}}_{\mathit{o}})}}^{\star }\rangle =-8{\pi }^{2}a\delta u{\rho }_0(r_o){\displaystyle \frac{\mathrm{d}X}{\mathrm{d}r}}(r_o){\varphi }_{\mathit{rr}}(\mathbf{0},\omega ).}\end{array}$$(B.27)

Formally, the integration over r makes the value of ϕ_rr at k = 0 appear. Physically, that means only the largest eddies have a significant impact on the correlation between the mode and the turbulent perturbations. Since considering eddies characterised by k = 0 does not actually make physical sense, we considered a comprise by assuming that the largest contributing eddies are those at the injection scale k₀, so that the correlated turbulent perturbations term finally becomes

$$\begin{array}{cc}& {\displaystyle \mathrm{Re}\left({\displaystyle \int \mathrm{d}\mathrm{\Omega }\stackrel{\sim }{h}(\mu )\langle \widehat{v_{\mathrm{osc}}}(\omega ){\widehat{u_n}}^{\star }(\omega )\rangle }\right)=C(\omega )\int \mathrm{d}\mathrm{\Omega }\stackrel{\sim }{h}(\mu )\mu }\hfill \\ \hfill & C(\omega )=-{\displaystyle \frac{1.083\lambda a\omega u_0\delta u}{r_{\mathrm{o}}{c}_{\mathrm{o}}{k}_0^4}}\hfill \\ \hfill & \times \sqrt{{\rho }_0(r_{\mathrm{o}})}{{\displaystyle \frac{\mathrm{d}\mathrm{Im}{X}_{\omega }}{\mathrm{d}r}}|}_{r_o}\times {(1+{\left({\displaystyle \frac{\lambda \omega }{1.204k_0{u}_0}}\right)}^2)}^{-1},\hfill \end{array}$$(B.28)

where every parameter is estimated at the observation height r_o, even when not specified explicitly. Note that we introduce no new free parameter compared to those used for the leading term, as the parameters a and δu appearing in the PCM are extracted from numerical hydrodynamic simulations of the stellar atmosphere. Plugging Eqs. (B.19) and (B.28) in the expression of the velocity power spectral density given by Eq. (4), we obtain Eq. (13).

Appendix C: A simplified toy model

Here we consider a simplified model of solar acoustic mode excitation, which has already been developed, used, and analysed in previous works (see e.g. Abrams & Kumar 1996; Chaplin & Appourchaux 1999). In the scope of this toy-model, the acoustic potential appearing in Eq. (8) takes the form of a square well, and the sound speed c is considered constant throughout the entire stellar radius. The latter approximation allows us to substitute the radial coordinate r in the wave equation for the acoustic depth τ, defined such that dτ = dr/c. In this approximation, the wave equation simply yields

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^2\mathrm{\Psi }}{\mathrm{d}{\tau }^2}}+({\omega }^2-V(\tau )+j\omega \gamma )\mathrm{\Psi }=\delta (\tau -{\tau }_{\mathrm{s}}),\end{array}$$(C.1)

where we have introduced a point-like source at acoustic depth τ_s, and the acoustic potential is

$$\begin{array}{c}\hfill V(\tau )=\{\begin{array}{cc}+\infty \hfill & \mathrm{if}\tau <0\hfill \\ 0\hfill & \mathrm{if}0\le \tau \le a\hfill \\ {\alpha }^2\hfill & \mathrm{if}a<\tau .\hfill \end{array}\end{array}$$(C.2)

In this model, a represents the acoustic length of the cavity (for radial modes, it corresponds to the time it takes for sound waves to propagate throughout the entire stellar radius) and α is the acoustic cut-off angular frequency above which waves are no longer confined. We added an infinite step at τ = 0 to force the wave variable Ψ to vanish at the centre.

C.1. Analytic solution of Eq. (C.1)

In order to solve the wave equation, it should be rewritten in terms of matrices:

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}X}{\mathrm{d}\tau }}=AX+B,\end{array}$$(C.3)

where

$$\begin{array}{c}\hfill \begin{array}{c}X=\left[\begin{array}{c}\mathrm{\Psi }\\ \mathrm{d}\mathrm{\Psi }/\mathrm{d}\tau \end{array}\right]\hfill \\ A=\left[\begin{array}{cc}0& 1\\ V(\tau )-{\omega }^2-j\omega \gamma & 0\end{array}\right]\hfill \\ B=\left[\begin{array}{c}0\\ \delta (\tau -{\tau }_{\mathrm{s}})\end{array}\right].\hfill \end{array}\end{array}$$(C.4)

The general solution to the homogeneous equation is X_h(τ) = exp(Aτ)C, where C is a constant vector. A particular solution to the inhomogeneous equation can then be sought in the form X_p(τ) = exp(Aτ)C(τ). For each domain in which the matrix A is constant, injecting this form in Eq. (C.3) yields

$$\begin{array}{c}\hfill C(\tau )=\{\begin{array}{cc}\left[\begin{array}{c}0\\ 0\end{array}\right]\hfill & \mathrm{if}\tau <{\tau }_{\mathrm{s}}\hfill \\ \hfill \\ \left[\begin{array}{c}0\\ exp(-A{\tau }_{\mathrm{s}})\end{array}\right]\hfill & \mathrm{if}{\tau }_{\mathrm{s}}⩽\tau .\hfill \end{array}\end{array}$$(C.5)

A being piecewise constant, we can thus write the general solution to Eq. (C.3) as

$$\begin{array}{c}\hfill X(\tau )=\{\begin{array}{cc}exp(A\tau )C\hfill & \mathrm{if}\tau <{\tau }_{\mathrm{s}}\hfill \\ exp(A\tau )(C+\left[\begin{array}{c}0\\ exp(-A{\tau }_{\mathrm{s}})\end{array}\right])\hfill & \mathrm{if}{\tau }_{\mathrm{s}}⩽\tau ,\hfill \end{array}\end{array}$$(C.6)

where the integration constant C is constant on whichever domain A is; in other words, C is piecewise constant on the same domains as A, that is, between 0 and a, and above a separately. In the following, we denote the column vector C as [A_i B_i] in the former domain, and [A_o B_o] in the latter, the indices i and o referring to the inner and outer regions, respectively.

The simple form of A allows for a straightforward computation of its exponential. We obtain

$$\begin{array}{c}\hfill exp(A\tau )=\{\begin{array}{cc}\left[\begin{array}{cc}cos{\omega }_{\mathrm{i}}\tau & {\displaystyle \frac{sin{\omega }_{\mathrm{i}}\tau }{{\omega }_{\mathrm{i}}}}\\ -{\omega }_{\mathrm{i}}sin{\omega }_{\mathrm{i}}\tau & cos{\omega }_{\mathrm{i}}\tau \end{array}\right]\hfill & \mathrm{if}0<\tau <a\hfill \\ \left[\begin{array}{cc}cosh{\omega }_{\mathrm{o}}\tau & {\displaystyle \frac{sinh{\omega }_{\mathrm{o}}\tau }{{\omega }_{\mathrm{o}}}}\\ {\omega }_{\mathrm{o}}sinh{\omega }_{\mathrm{o}}\tau & cosh{\omega }_{\mathrm{o}}\tau \end{array}\right]\hfill & \mathrm{if}a<\tau ,\hfill \end{array}\end{array}$$(C.7)

where ${\omega }_{\text{i}}^2={\omega }^2+j\omega \gamma$ and ${\omega }_{\text{o}}^2={\alpha }^2-{\omega }^2-j\omega \gamma$ (with the understanding that 0 <  ω <  α).

Finally, injecting this into the general solution (C.6) and only keeping the first line (at this point, the second one only gives the derivative of the solution and is redundant with the solution itself) yields the following expression, depending on whether the source is inside the cavity or not:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathrm{i}}(\tau )=\{\begin{array}{cc}{A}_{\mathrm{i}}cos{\omega }_{\mathrm{i}}\tau +{\displaystyle \frac{B_{\mathrm{i}}}{{\omega }_{\mathrm{i}}}}sin{\omega }_{\mathrm{i}}\tau \hfill & \mathrm{if}\tau <{\tau }_{\mathrm{s}}\hfill \\ A_{\mathrm{i}}cos{\omega }_{\mathrm{i}}\tau +{\displaystyle \frac{B_{\mathrm{i}}}{{\omega }_{\mathrm{i}}}}sin{\omega }_{\mathrm{i}}\tau \hfill & \\ +{\displaystyle \frac{1}{{\omega }_{\mathrm{i}}}}sin{\omega }_{\mathrm{i}}(\tau -{\tau }_{\mathrm{s}})\hfill & \mathrm{if}{\tau }_{\mathrm{s}}<\tau <a\hfill \\ A_{\mathrm{o}}cosh{\omega }_{\mathrm{o}}\tau +{\displaystyle \frac{B_{\mathrm{o}}}{{\omega }_{\mathrm{o}}}}sinh{\omega }_{\mathrm{o}}\tau \hfill & \mathrm{if}a<{\tau }_{\mathrm{s}},\hfill \end{array}\end{array}$$(C.8)

and

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathrm{o}}(\tau )=\{\begin{array}{cc}{A}_{\mathrm{i}}cos{\omega }_{\mathrm{i}}\tau +{\displaystyle \frac{B_{\mathrm{i}}}{{\omega }_{\mathrm{i}}}}sin{\omega }_{\mathrm{i}}\tau \hfill & \mathrm{if}\tau <a\hfill \\ A_{\mathrm{o}}cosh{\omega }_{\mathrm{o}}\tau +{\displaystyle \frac{B_{\mathrm{o}}}{{\omega }_{\mathrm{o}}}}sinh{\omega }_{\mathrm{o}}\tau \hfill & \mathrm{if}a<\tau <{\tau }_{\mathrm{s}}\hfill \\ A_{\mathrm{o}}cosh{\omega }_{\mathrm{o}}\tau +{\displaystyle \frac{B_{\mathrm{o}}}{{\omega }_{\mathrm{o}}}}sinh{\omega }_{\mathrm{o}}\tau \hfill & \\ +{\displaystyle \frac{1}{{\omega }_{\mathrm{o}}}}sinh{\omega }_{\mathrm{o}}(\tau -{\tau }_{\mathrm{s}})\hfill & \mathrm{if}{\tau }_{\mathrm{s}}<\tau ,\hfill \end{array}\end{array}$$(C.9)

where Ψ_i is the solution if the source is inside the cavity (τ_s <  a) and Ψ_o is the solution if the source is outside (a <  τ_s). Predictably, this general solution contains 4 constants of integration (in the form of A_i, B_i, A_o and B_o), and we therefore need 4 boundary conditions to find Ψ_i and Ψ_o. We impose that Ψ(τ = 0) = 0 and that the solution do not diverge when τ → +∞; furthermore, we impose that both Ψ and its derivative be continuous at τ = a. With all calculations having been carried out, this set of boundary conditions finally gives us:

$$\begin{array}{cc}& {\mathrm{\Psi }}_{\mathrm{i},\mathrm{o}}(\tau )=-{\displaystyle \frac{f_{\mathrm{i},\mathrm{o}}({\tau }_{\mathrm{s}})}{{\omega }_{\mathrm{i}}cos{\omega }_{\mathrm{i}}a+{\omega }_{\mathrm{o}}sin{\omega }_{\mathrm{i}}a}}{exp}^{-{\omega }_{\mathrm{o}}(\tau -a)}\hfill \\ \hfill & f_{\mathrm{i}}({\tau }_{\mathrm{s}})=sin{\omega }_{\mathrm{i}}{\tau }_{\mathrm{s}}\hfill \\ \hfill & f_{\mathrm{o}}({\tau }_{\mathrm{s}})={\omega }_{\mathrm{i}}cos{\omega }_{\mathrm{i}}asinh{\omega }_{\mathrm{o}}({\tau }_{\mathrm{s}}-a)\hfill \\ \hfill & +{\omega }_{\mathrm{o}}sin{\omega }_{\mathrm{i}}acosh{\omega }_{\mathrm{o}}({\tau }_{\mathrm{s}}-a).\hfill \end{array}$$(C.10)

We note that the above expression is only valid if τ >  a and τ >  τ_s. Since the first condition means that we observe the mode of oscillation above the upper turning point (which is always the case in practice) and since the second condition means that the excitation of the mode occurs in layers located deeper than the observation height in the atmosphere (which is also always the case in practice, at least for spectrometric measurements), these are not restrictive conditions.

C.2. Discussion

In each case (source inside or outside the cavity), the solution Ψ_i, o can be decomposed into three parts:

• the denominator corresponds to the Wronskian W(ω) of the wave equation. |1/W(ω)|² peaks at the eigenfrequencies associated to the acoustic cavity and is responsible for the presence of resonant modes in the spectrum. The line profile it generates have a Lorentzian profile, so long as the damping rate of the modes are much smaller than their frequency;

• the exponential part accounts for the fact that the modes are evanescent outside the acoustic cavity, so that the higher in the atmosphere the mode is observed, the lower its amplitude as we observe it. Therefore, it only affects the observed amplitude of the mode, not its line profile;

• the numerator f_o, i(τ_s) is responsible for the mode line profile asymmetry. Because it is the only term that depends on the source position τ_s, it is commonly said that mode asymmetry is caused by source localisation.

Regardless of whether the source is located inside or outside the cavity, it can be seen from Eqs. (C.8) and (C.9) that the third term f_o, i(τ_s) actually corresponds to the amplitude of the solution Ψ_i, o at τ = τ_s. This leads us to the following conclusion, which we reproduce in the main body of the paper: for a given frequency, the amplitude of the excited wave is proportional to the radial profile associated with the wave at the source of excitation.

All Tables

All Figures