Issue 
A&A
Volume 635, March 2020



Article Number  A81  
Number of page(s)  25  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201936847  
Published online  13 March 2020 
Modelling the asymmetries of the Sun’s radial pmode line profiles
^{1}
LESIA, Observatoire de Paris, PSL Research University, CNRS, Université Pierre et Marie Curie, Université Paris Diderot, 92195 Meudon, France
email: jordan.philidet@obspm.fr
^{2}
Zentrum für Astronomie der Universität Heidelberg, Landessternwarte, Königstuhl 12, 69117 Heidelberg, Germany
^{3}
GEPI, Observatoire de Paris, PSL Research University, CNRS, Université Pierre et Marie Curie, Université Paris Diderot, 92195 Meudon, France
Received:
5
October
2019
Accepted:
16
January
2020
Context. The advent of spaceborne missions has substantially increased the number and quality of the measured power spectrum of solarlike oscillators. It now allows for the pmode line profiles to be resolved and facilitates an estimation of their asymmetry. The fact that this asymmetry can be measured for a variety of stars other than the Sun calls for a revisiting of acoustic mode asymmetry modelling. This asymmetry has been shown to be related to a highly localised source of stochastic driving in layers just beneath the surface. However, existing models assume a very simplified, pointlike source of excitation. Furthermore, mode asymmetry could also be impacted by a correlation between the acoustic noise and the oscillating mode. Prior studies have modelled this impact, but only in a parametrised fashion, which deprives them of their predictive power.
Aims. In this paper, we aim to develop a predictive model for solar radial pmode line profiles in the velocity spectrum. Unlike the approach favoured by prior studies, this model is not described by free parameters and we do not use fitting procedures to match the observations. Instead, we use an analytical turbulence model coupled with constraints extracted from a 3D hydrodynamic simulation of the solar atmosphere. We then compare the resulting asymmetries with their observationally derived counterpart.
Methods. We model the velocity power spectral density by convolving a realistic stochastic source term with the Green’s function associated with the radial homogeneous wave equation. We compute the Green’s function by numerically integrating the wave equation and we use theoretical considerations to model the source term. We reconstruct the velocity power spectral density and extract the line profile of radial pmodes as well as their asymmetry.
Results. We find that stochastic excitation localised beneath the mode upper turning point generates negative asymmetry for ν < ν_{max} and positive asymmetry for ν > ν_{max}. On the other hand, stochastic excitation localised above this limit generates negative asymmetry throughout the pmode spectrum. As a result of the spatial extent of the source of excitation, both cases play a role in the total observed asymmetries. By taking this spatial extent into account and using a realistic description of the spectrum of turbulent kinetic energy, both a qualitative and quantitative agreement can be found with solar observations performed by the GONG network. We also find that the impact of the correlation between acoustic noise and oscillation is negligible for mode asymmetry in the velocity spectrum.
Key words: methods: numerical / turbulence / Sun: helioseismology / Sun: oscillations / line: profiles
© J. Philidet et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Solarlike oscillations are known to be stochastically excited and damped by turbulence occurring close to the surface of lowmass 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 Lorentzianshaped 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 groundbased and spaceborne) 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 pmode 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 lowfrequency wing than in their highfrequency wing and viceversa 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 nonadiabatic 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). Nonadiabaticity 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 nonlinearly 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 nonlinear nature of the κ − T relation, this modulation decreases the observed temperature fluctuations more significantly in the lowfrequency wing of the mode than in its highfrequency 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 (coherentcorrelated, coherentuncorrelated and incoherent) in both the velocity spectrum, the intensity spectrum, and the velocityintensity crossspectrum. 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 pointlike 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 bestfit 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 solarlike 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 solarlike 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 solarlike oscillations in stars other than the Sun by Benomar et al. (2018). Inversion methods used to infer the internal structure of solarlike 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 pmode line profiles is of primary importance.
In this paper, we present a predictive model of solar radial pmode 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 nonparametrised 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 pmode line profiles
To extract the asymmetries of solar radial pmodes, 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 discintegrated 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 limbdarkening and of disk integration that affect the Sunasastar 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 Sunasastar 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 limbdarkening 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
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 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
where u_{n} is the component of the turbulent velocity along the line of sight. We introduced the reduced limbdarkening 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 is negligible compared to the terms that contain and , 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^{5} m^{2} s^{−2} Hz^{−1}, while the latter is of order 10 m^{2} s^{−2} Hz^{−1}, e.g. TurckChièze et al. 2004, Fig. 2), so that
where the notation Re refers to the real part of a complex quantity, and ⋆ refers to its complex conjugate. Finally, we obtain
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 nonadiabatic 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:
where c is the sound speed, the wave variable Ψ_{ω}(r) is related to the radial fluid displacement ξ_{r}(r) through
and the acoustic potential and source term are given by
where r is the radial coordinate, ρ_{0} is the density, N is the BruntVäisälä frequency, g_{0} is the gravitational acceleration, G is the gravitational constant, and 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 phaselag between these different modulations energy can be exchanged with the surrounding medium. However, the modelling of damping rates of solarlike 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 (frequencydependent) linewidth Γ_{ω}, so that the wave equation takes the following form
where j denotes the imaginary unit and the linewidths Γ_{ω} are inferred from observations. We used the linewidths 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 lowfrequency linewidths 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 linewidths 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 pointlike source term), where Ψ_{ω} is the solution to the inhomogeneous wave equation,
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
where the source term is given by Eq. (7). The pulsational velocity v_{osc} is related to the variable Ψ_{ω} through
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
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}):
where and C(ω) are given, respectively, by Eqs. (B.19) and (B.28). We note that the effects of limbdarkening 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 fourthorder and thirdorder twopoint correlation moments of the turbulent velocity. Therefore, an appropriate closure model is needed to express these highorder moments as a function of secondorder moments. We devote the following subsection to presenting and developing these closure models.
Fourthorder moments. To describe the fourthorder correlation moments of the turbulent velocity, we make use of the quasinormal approximation (QNA). This closure model consists in considering that all turbulent quantities are normally distributed, in which case fourthorder moments can be analytically expressed as a combination of secondorder moments (Lesieur 2008):
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 thirdorder correlation moments. When it comes to estimating the fourthorder moments, however, the picture is different. Belkacem et al. (2006a) have studied the validity of the QNA for twopoints, fourthorder correlation moments of the vertical turbulent velocity, in the form of (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 onepoint 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 twopoints, fourthorder moments. As such, the decomposition given by Eq. (14) can be considered valid when it comes to studying mode asymmetry.
Thirdorder moments. While the QNA provides an adequate closure relation for fourthorder moments, as mentioned above, it assumes vanishing thirdorder moments. Therefore, in order to estimate these thirdorder 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 nonvanishing thirdorder correlation moments:
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 is the fluctuation of the vertical velocity around its mean value in the downflows.
We note that, strictly speaking, the thirdorder 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 noncentred moments.
We also note that this closure relation is written here in terms of (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
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 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 twopoints, fourthorder 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 secondorder moments.
We then use analytical expressions for the secondorder moments of the turbulent velocity. We describe the secondorder 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):
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 lifetime of eddies of wavenumber k
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 selfsimilar, Kolmogorov’s theory of turbulence leads to a spatial spectrum E(k)∝k^{−5/3} in the inertial range, between k = k_{0} (where k_{0} 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^{14}), 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_{0}. 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
where and the 0.652 factor is introduced so that the total specific kinetic energy of the turbulent spectrum matches . Therefore, the spatial spectrum is parametrised solely by the injection scale k_{0}. However, the injection scale varies significantly between the subsurface 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 eddytime correlation is expected to be a decreasing exponential, and its Fourier transform a Lorentzian function (Belkacem et al. 2011). The width ω_{k} associated to eddytime correlation is linked to the lifetime 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:
where λ is a dimensionless, constant parameter. Overall, the only input parameters of this model are k_{0}(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^{5}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
where u_{0} is also extracted from the 3D atmospheric simulation, by averaging the fluid velocity squared temporally and horizontally, and using the definition .
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, energybearing 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 cutoff angular frequency ω_{E}, which is given by the curvature of the eddytime correlation function at τ = 0 (Belkacem et al. 2010):
Since a Gaussian spectrum would fall off much more rapidly than a Lorentzian spectrum, we simply consider that χ_{k} vanishes entirely for ω > ω_{E},
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
Physically, taking the cutoff 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 cutoff 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 twopoint, thirdorder moments. We use the closure relation presented in Sect. 2.3.1 to express them, as we did with the fourthorder moments, as a function of twopoint, secondorder moments of the turbulent velocity. We then use the same analytical description of secondorder 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_{0}(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 pmode 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 mixinglength 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^{5}BOLD code (Freytag et al. 2012). The modelled region includes the superadiabatic 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, mixinglength parameter α_{MLT}, and helium abundance) are chosen so that the top layers match the bottom layers of the CO^{5}BOLD atmosphere. Here the model was computed with the mixinglength 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).
Fig. 1. Acoustic potential V(r) used in Eq. (8), calculated using Eq. (7) and the equilibrium state of the Sun given by the solar patched model described in the text. The radius is normalised by the photospheric radius R_{⊙}, and only the outermost region is shown. We note that the acoustic potential is normalised by here (where R_{⊙} is the radius of the solar photosphere) so that it is dimensionless. 

Open with DEXTER 
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_{0} 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_{0} 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} − Γ_{ω0} and ω_{0} + Γ_{ω0}, where ω_{0} 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 4thorder RungeKutta 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 forcefree). The continuity equation and the equation of state allow us to rewrite this latter condition in terms of Ψ_{ω}:
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 pointlike source term to the numerical scheme. The source is normalised in such a way that the righthand 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 complexvalued. As such, the Green’s function is complexvalued 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_{0}), 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
where x = 2(ω − ω_{0})/Γ_{ω0} is the reduced pulsation frequency. The fit contains three free parameters (H_{0}, ν_{0} 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 lowfrequency side (that corresponds to negative asymmetry), B > 0 means that the highfrequency 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.
Fig. 2. Dependence of the line profile given by Eq. (26) on the asymmetry parameter B. The line profiles are normalised with H_{0} = 1. 

Open with DEXTER 
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
and defines the asymmetry parameter as α_{n, l}. Meanwhile, Vorontsov & Jefferies (2013) use the following formula:
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):
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 halfway 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 linewidths 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 higherorder modes (because the increasing linewidths lead to mode overlapping), while we did not have access to observed linewidths for lowerorder 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 (ChristensenDalsgaard 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.
Fig. 3. Eigenfunction Ψ(r) of several radial acoustic modes (black: n = 10 (ν_{0} = 1.580 mHz); red: n = 20 (ν_{0} = 2.912 mHz); blue: n = 30 (ν_{0} = 4.267 mHz)) as computed by our model (solid lines), and calculated using the ADIPLS code (dashed lines). The radial axis is zoomed to show only the outermost region. The eigenfunctions have been normalised so that their maximum over the entire solar interior equals 1. 

Open with DEXTER 
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 Γ,
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 lowfrequency 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.
Fig. 4. Velocity amplitudes of radial acoustic modes as computed by our model, using Eq. (30) (black dashed line), and as observed by GOLF (black points). The data points are taken from Baudin et al. (2005). The free parameters in the model have been tuned to obtain the best possible agreement with observational data. 

Open with DEXTER 
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 pmodes 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_{0}(r), which is defined as the (radiusdependent) 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 energybearing eddies increases with height, so that the injection scale k_{0} decreases with r (Samadi et al. 2003). We simplify the situation by considering that the injection rate only takes two values: k_{0}(r) = k_{0, int} below the photosphere, and k_{0}(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_{0}, 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}.
Fig. 5. Asymmetry parameter B as a function of frequency obtained by the “theoretical spectrum” model, for different values of k_{0, int}, with λ = 1 and R_{k} = 2 fixed. The subphotospheric injection scale k_{0, int} is expressed in Mm^{−1}. 

Open with DEXTER 
Fig. 6. Same as Fig. 5, but only R_{k} varies, λ = 1 and k_{0, int} = 2 Mm^{−1}. 

Open with DEXTER 
Fig. 7. Same as Fig. 5, but only λ varies, R_{k} = 2 and k_{0, int} = 2 Mm^{−1}. 

Open with DEXTER 
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 modelindependent. Thus the asymmetries of the solar radial pmode 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 lowfrequency mode asymmetries (below ν_{max} ∼ 3 mHz), while on the contrary, the highfrequency 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 lowfrequency modes and decreases for highfrequency 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 pmodes 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.
Fig. 8. Asymmetry parameter B as a function of frequency obtained by the “numerical spectrum” model, for several values of λ. 

Open with DEXTER 
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 solarlike 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 toymodel 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 FabryPé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 pointlike 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 toymodel of mode excitation, where the source is considered pointlike and the acoustic cavity is simplified to a square well potential. From this toymodel 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_{1} (see illustration in Fig. 9) will see the amplitude of the oscillation increase with ω, and a source at r = r_{2} 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_{1}, 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_{2} will be negative.
Fig. 9. Illustration of the importance of source position with respect to nodes and antinodes of the eigenfunction associated to a mode to explain its asymmetry. The blue and red lines show the radial profile of the oscillation for two angular frequencies very close to one another (ω_{red} < ω_{blue}). The bold vertical dashed lines show two source positions generating opposite mode asymmetries: positive for r_{1}, negative for r_{2}. The third vertical dashed line marks the edge of the acoustic cavity r = a. 

Open with DEXTER 
From the illustration in Fig. 9, it is straightforward to see that the dichotomy between the r = r_{1} case and the r = r_{2} 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_{1} 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_{2} 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_{2} 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 toymodel 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 nonadiabatic 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
where P is the total power spectral density, x = 2(ω − ω_{0})/Γ_{ω} is the reduced frequency (ω_{0} 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
The first term of the righthand 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 pmode 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 pointlike source of excitation. However, the driving region of the solar pmodes, while localised around the superadiabatic 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 superadiabatic 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 superadiabatic peak, where the stochastic excitation is mainly located. Below ν ∼ 3.5 mHz, the superadiabatic 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 superadiabatic peak and it generates positive asymmetry. Close to 3.5 mHz, the superadiabatic peak coincides with an antinode, so that the asymmetry is very low.
Fig. 10. Radial location of the nodes (red symbols connected by a dashed line) and antinodes (blue symbols connected by a dashed line) that are closest to the superadiabatic peak, for each radial mode between n = 16 and n = 30. The vertical black line represents the maximum of the superadiabatic peak, where the excitation is most efficient. Horizontal black dotted lines are added for readability only. 

Open with DEXTER 
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 , where k_{0} 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_{0} (since it scales to ).
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_{0}(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 pmodes compares differently to the typical turbulent eddytime 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 ω^{−2} 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^{−6} m^{−1}, u_{k} ∼ 10^{3} m s^{−1} and ω ∼ 10^{−3} rad s^{−1}, 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 lowfrequency 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.
Fig. 11. Relative difference (in percentages) between the asymmetry parameter B obtained by respectively taking into account and discarding the cross term C(ω) (k_{0, int} = 2 Mm^{−1}, R_{k} = 2 and λ = 1), as a function of frequency. 

Open with DEXTER 
It is possible to support this conclusion with a simple order of magnitude estimation. Indeed, B is of order 10^{−2}, 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 CauchySchwarz inequality provides an upper bound to the cross term: , where the case of equality happens when the correlation between the mode v and the turbulence u is optimal. Therefore
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 m^{2} s^{−2} Hz^{−1} and that of the turbulent fluctuations as m^{2} s^{−2} Hz^{−1} (Fig. 2 TurckChièze et al. 2004), the above ratio becomes . 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 CauchySchwarz 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 solarlike 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 nonradial 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.
Fig. 12. Asymmetry profile obtained by the “numerical spectrum” model (solid red line), compared to the observed asymmetry profile (blue crosses). The data points are taken from Barban et al. (2004). For more readability, only data points corresponding to modes with angular degrees 15 ≤ l ≤ 20 have been retained. The asymmetry profile is given in terms of the parameter B (top panel), and alternatively in terms of the parameter χ (bottom panel), which is defined by Eq. (29). The error bars in the top panel correspond to the uncertainty on the observed values of the mode linewidths, which propagates to the asymmetry parameter B. 

Open with DEXTER 
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 (groundbased), 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.
Summary of the nature, wavelength, and formation height of the absorption line used by the instruments considered in this paper.
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 subphotospheric 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.
Fig. 13. Asymmetry profile B(ν) obtained by the “numerical spectrum” model (with λ = 0.5), as would be observed at three different heights in the atmosphere, corresponding respectively to the observation heights of GOLF (dashed red line), MDI/GONG (dashed blue line) and HMI (dashed black line). 

Open with DEXTER 
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 bestfit 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
Figure 14 shows the results obtained for all radial modes described in Sect. 4. Given that the asymmetry of solar pmodes in the velocity spectrum are all negative, using a symmetric fit introduces an underestimation of ν_{0}. 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.
Fig. 14. Reduced frequency bias δν, defined by Eq. (34), as a function of the asymmetry parameter B for each radial pmode between n = 6 and n = 30. The blue dashed line shows the best linear fit. 

Open with DEXTER 
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} = ν_{0} + BΓ_{ν0}/2). Therefore, while the asymmetric fit accurately finds the eigenfrequency at , the symmetric fit finds it at . 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 highfrequency 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 pmodes 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^{5}BOLD and horizontally averaged. We integrated the wave equation along the solar radius, and added a pointlike, 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 fourthorder correlation products of the turbulent velocity as functions of secondorder 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 pmode 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 pointlike 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 pointlike (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 solarlike oscillators. It also shows that any model that assumes a pointlike source of excitation cannot give reliable results as far as mode asymmetry is concerned. In particular, such a model would predict positive asymmetries for highfrequency 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 linewidth, which is widest at high frequency), it is likely to have a nonnegligible 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 pmodes. Our formalism can be easily adapted to the study of nonradial modes simply by using a nonradial wave equation instead of the radial one. However, since the eigenfunction associated to pmodes in solarlike 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 nonadiabatic, turbulent pressure fluctuations.
Our formalism can also be easily applied to other solarlike oscillators. Comparing the asymmetries featured by the velocity spectra of several solarlike 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 solarlike oscillators. However, the one major difference that remains between the solar case and other stars is that the Sun is the only solarlike 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 MarieJo 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
 Abrams, D., & Kumar, P. 1996, ApJ, 472, 882 [NASA ADS] [CrossRef] [Google Scholar]
 Balmforth, N. J. 1992, MNRAS, 255, 639 [NASA ADS] [Google Scholar]
 Barban, C., Hill, F., & Kras, S. 2004, ApJ, 602, 516 [NASA ADS] [CrossRef] [Google Scholar]
 Batchelor, G. K. 1953, The Theory of Homogeneous Turbulence (Cambridge: Cambridge University Press) [Google Scholar]
 Baudin, F., Samadi, R., Goupil, M.J., et al. 2005, A&A, 433, 349 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., Kupka, F., & Baudin, F. 2006a, A&A, 460, 183 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., & Kupka, F. 2006b, A&A, 460, 173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M.J., & Dupret, M.A. 2008, A&A, 478, 163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M. J., et al. 2010, A&A, 522, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., & Goupil, M. J. 2011, GONGSoHO 24: A New Era of Seismology of the Sun and SolarLike Stars, 271, 012047 [Google Scholar]
 Benomar, O., Goupil, M., Belkacem, K., et al. 2018, ApJ, 857, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Bruls, J. H. M. J., Rutten, R. J., & Shchukina, N. G. 1992, A&A, 265, 237 [NASA ADS] [Google Scholar]
 Chaplin, W. J., & Appourchaux, T. 1999, MNRAS, 309, 761 [NASA ADS] [CrossRef] [Google Scholar]
 Chaplin, W. J., Elsworth, Y., Isaak, G. R., Miller, B. A., & New, R. 1999, MNRAS, 308, 424 [NASA ADS] [CrossRef] [Google Scholar]
 Chaplin, W. J., Houdek, G., Elsworth, Y., et al. 2005, MNRAS, 360, 859 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J. 2011, Astrophysics Source Code Library [record ascl:1109.002] [Google Scholar]
 Davies, G. R., Broomhall, A. M., Chaplin, W. J., Elsworth, Y., & Hale, S. J. 2014, MNRAS, 439, 2025 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L., Jefferies, S. M., Harvey, J. W., Osaki, Y., & Pomerantz, M. A. 1993, ApJ, 410, 829 [NASA ADS] [CrossRef] [Google Scholar]
 Fleck, B., Couvidat, S., & Straus, T. 2011, Sol. Phys., 271, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H.G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
 Gabriel, M. 1992, A&A, 265, 771 [NASA ADS] [Google Scholar]
 Gabriel, M. 1993, A&A, 274, 935 [NASA ADS] [Google Scholar]
 Georgobiani, D., Stein, R. F., & Nordlund, Å. 2003, ApJ, 596, 698 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L. 2006, Cent. Eur. Astrophys. Bull., 30, 1 [NASA ADS] [Google Scholar]
 Goldreich, P., & Keeley, D. A. 1977a, ApJ, 211, 934 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Keeley, D. A. 1977b, ApJ, 212, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Goode, P. R., Strous, L. H., Rimmele, T. R., & Stebbins, R. T. 1998, ApJ, 495, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Gough, D. O. 1977, ApJ, 214, 196 [NASA ADS] [CrossRef] [Google Scholar]
 Jefferies, S. M., Duvall, Jr., T. L., Harvey, J. W., Osaki, Y., & Pomerantz, M. A. 1991, ApJ, 377, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Jefferies, S. M., Severino, G., Moretti, P.F., Oliviero, M., & Giebink, C. 2003, ApJ, 596, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Korzennik, S. G. 2005, ApJ, 626, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Kraichnan, R. H. 1957, Phys. Rev., 107, 1485 [NASA ADS] [CrossRef] [Google Scholar]
 Kumar, P., & Basu, S. 1999, ApJ, 519, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2017, Liv. Rev. Comput. Astrophys., 3, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Lesieur, M. 2008, Turbulence in Fluids (Berlin: Springer) [CrossRef] [Google Scholar]
 Manchon, L., Belkacem, K., Samadi, R., et al. 2018, A&A, 620, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morel, P. 1997, A&AS, 124, 597 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Musielak, Z. E., Rosner, R., Stein, R. F., & Ulmschneider, P. 1994, ApJ, 423, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Nigam, R., & Kosovichev, A. G. 1998, ApJ, 505, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Nigam, R., Kosovichev, A. G., Scherrer, P. H., & Schou, J. 1998, ApJ, 495, L115 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 2 [Google Scholar]
 Press, W. H., Flannery, B. P., & Teukolsky, S. A. 1986, Numerical Recipes. The Art of Scientific Computing (Cambridge: University Press) [Google Scholar]
 Rast, M. P., & Bogdan, T. J. 1998, ApJ, 496, 527 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenthal, C. S. 1998, ApJ, 508, 864 [NASA ADS] [CrossRef] [Google Scholar]
 Roxburgh, I. W., & Vorontsov, S. V. 1995, MNRAS, 272, 850 [NASA ADS] [Google Scholar]
 Roxburgh, I. W., & Vorontsov, S. V. 1997, MNRAS, 292, L33 [NASA ADS] [Google Scholar]
 Samadi, R., & Goupil, M.J. 2001, A&A, 370, 136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Nordlund, Å., Stein, R. F., Goupil, M. J., & Roxburgh, I. 2003, A&A, 403, 303 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Georgobiani, D., Trampedach, R., et al. 2007, A&A, 463, 297 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Belkacem, K., Goupil, M. J., Dupret, M.A., & Kupka, F. 2008, A&A, 489, 291 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Belkacem, K., & Sonoi, T. 2015, EAS Publ. Ser., 73, 111 [CrossRef] [Google Scholar]
 Severino, G., Magrì, M., Oliviero, M., Straus, T., & Jefferies, S. M. 2001, ApJ, 561, 444 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, R. F. 1967, Sol. Phys., 2, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, R. F., & Nordlund, Å. 2001, ApJ, 546, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Thiery, S., Boumier, P., Gabriel, A. H., et al. 2000, A&A, 355, 743 [NASA ADS] [Google Scholar]
 Toutain, T., Appourchaux, T., Fröhlich, C., et al. 1998, ApJ, 506, L147 [NASA ADS] [CrossRef] [Google Scholar]
 Trampedach, R. 1997, Master’s Thesis, Aarhus University, Denmark [Google Scholar]
 TurckChièze, S., García, R. A., Couvidat, S., et al. 2004, ApJ, 604, 455 [NASA ADS] [CrossRef] [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (Tokyo: University of Tokyo Press) [Google Scholar]
 Vorontsov, S. V., & Jefferies, S. M. 2013, ApJ, 778, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Wachter, R., & Kosovichev, A. G. 2005, ApJ, 627, 550 [NASA ADS] [CrossRef] [Google Scholar]
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
Writing ρ = ρ_{0} + ρ′ (where ρ′ is the Eulerian density perturbation corresponding to the mode), linearising this equation around a motionless state, and integrating with respect to time yields
We then introduce the Lagrangian density perturbation δρ = ρ′+(ξ.∇)ρ_{0}, which allows us to write the linearised continuity equation in its final form:

the Euler equation:
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
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 , where 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:
where p′ is the Eulerian pressure perturbation, g_{0} is the mean gravitational acceleration, g′ is its Eulerian perturbation and 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)):
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:
where δs corresponds to the turbulent fluctuation of the specific entropy of the fluid and we define the various thermodynamic coefficients as
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
and that by definition of the BruntVäisälä frequency, we have
and we finally obtain
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^{2} = Γ_{1}p_{0}/ρ_{0} (where c is the sound speed), the continuity and Euler equations then yield:
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):
Plugging this into Eq. ((A.13)), we obtain
and
where we denote the righthand side terms of Eqs. (A.15) and (A.16) as S_{0} and S_{1} respectively in the following. We also define
The above set of equations can be rewritten as
We can now eliminate to get a single secondorder wave equation. Using the first of Eq. (A.18) to express as a function of , and plugging it in the second equation, we get the following equation:
Similarly to what has been done for the first change of variables, we wish for the lefthand side to contain no firstorder term, but only secondorder and 0thorder ones. Thus we introduce yet another variable: . 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
with an acoustic potential V(r) that only depends on the star’s equilibrium state:
A.3. The source term
With the above notations, the parameters intervening in the source term of Eq. (A.20) have the following expressions:
Furthermore, one can easily derive the following relationship between ∇_{ad} and α_{s} ≡ (∂P/∂s)_{ρ} by means of the adequate Schwarz relation:
After some manipulations, one finally obtain the source term in the form:
Finally, since
this expression can be drastically simplified to
This form clearly shows that the source term can be split three ways: a monopolar source term (proportional to δsdα_{s}/dr) due to nonadiabatic pressure fluctuations in a stratified environment, a dipolar term (proportional to α_{s}dδs/dr) due to a stratification in the nonadiabatic pressure fluctuations themselves, and a quadripolar term (proportional to ) due to Reynolds stress fluctuations. In the following, we only consider this last term but we also show here how the effect of nonadiabaticity can be introduced as well.
To conclude, note that the value of the fluid density at the centre of the star ρ_{0}(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:
with V(r) given by Eq. (A.21), and the source term and wave variable are given by
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 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
We then have
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 secondorder 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
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
We then use the QuasiNormal Approximation (hereby abbreviated QNA), under which any fourthorder correlation product can be decomposed into a sum of three secondorder correlation products, so that (Lesieur 2008)
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
Using the Parseval identity, we can express this as an integral over wave vectors k and angular frequencies ω
where the notation TF[.] refers to temporal and spatial Fourier transform.
We then proceed to describe the secondorder 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 secondorder correlation product of the ith and jth component of the turbulent velocity as ϕ_{ij}(k, ω), so that
For isotropic turbulence, ϕ_{ij} can be expressed analytically (Batchelor 1953) as
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)
where
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 lifetime distribution of eddies of wavenumber k
Finally, the integral given by Eq. (B.8) can be rewritten as
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 selfsimilar, Kolmogorov’s theory of turbulence leads to a spatial spectrum E(k)∝k^{−5/3} in the inertial range, between k = k_{0} (where k_{0} 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^{14}), 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_{0}. 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:
where and the 0.652 factor is introduced so that the total specific kinetic energy of the turbulent spectrum matches .
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 timecorrelation function which decays exponentially, its spectrum is Lorentzian). Thus:
The width of the Lorentzian is the inverse of the typical correlation timescale 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
For a Kolmogorov spectrum, u_{k} scales as k^{−1/3}, which means that we have
There is a temptation to approximate the typical velocities of eddies of wavenumber k_{0} with u_{0}. This assumption requires some discussion, however. Indeed, Stein (1967) has pointed out that eddies of all sizes have the same Eulerian velocity fluctuations u_{0}. As far as Lagrangian fluctuations go, the fluctuations u_{k} can be expressed as (Stein 1967)
Using the expression of E(k) given in Eq. (B.14) and applying it to k = k_{0}, we finally find u_{k0} = 0.602u_{0}. Under all these assumptions, all further calculations being carried out, we ultimately obtain the leading term,
with
and where K is the reduced inverse eddy scale (K ≡ k/k_{0}). We note that all the terms appearing in the integrand depend on the mass variable m, particularly u_{0} and k_{0}, 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_{0}(m); solar turbulence close to the photosphere being as poorly constrained as it is today, we cannot hope to achieve a nonparametrised 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
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 lineofsight component of u as u_{n} = u_{r} cos θ − u_{θ} sin θ, this transforms into
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 thirdorder 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:
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
The challenge in estimating the contribution of correlated turbulent perturbations fundamentally lies in the correct determination of nonlocal, thirdorder correlation products. While the QNA provides an adequate closure relationship for fourthorder moments, it yields vanishing thirdorder 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 thirdorder 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 thirdorder correlation products. In Belkacem et al. (2006b), the authors use the same approximations but focused on onepoint correlation products; however, the calculations can be easily extended to the twopoint correlation products that we need and the model yields
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 is the fluctuation of the downflow velocity around its mean value. The only additional approximation we make to adapt these calculations to twopoint 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 thirdorder 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 noncentred moment.
Also note that this closure relation is written here in terms of (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
Plugging Eq. (B.26) into Eq. (B.25), we obtain a closure model that allows us to write the thirdorder moments as a function of secondorder 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
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_{0}, so that the correlated turbulent perturbations term finally becomes
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 toymodel, 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
where we have introduced a pointlike source at acoustic depth τ_{s}, and the acoustic potential is
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 cutoff 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:
where
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
A being piecewise constant, we can thus write the general solution to Eq. (C.3) as
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
where and (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:
and
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:
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(ω)^{2} 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
Summary of the nature, wavelength, and formation height of the absorption line used by the instruments considered in this paper.
All Figures
Fig. 1. Acoustic potential V(r) used in Eq. (8), calculated using Eq. (7) and the equilibrium state of the Sun given by the solar patched model described in the text. The radius is normalised by the photospheric radius R_{⊙}, and only the outermost region is shown. We note that the acoustic potential is normalised by here (where R_{⊙} is the radius of the solar photosphere) so that it is dimensionless. 

Open with DEXTER  
In the text 
Fig. 2. Dependence of the line profile given by Eq. (26) on the asymmetry parameter B. The line profiles are normalised with H_{0} = 1. 

Open with DEXTER  
In the text 
Fig. 3. Eigenfunction Ψ(r) of several radial acoustic modes (black: n = 10 (ν_{0} = 1.580 mHz); red: n = 20 (ν_{0} = 2.912 mHz); blue: n = 30 (ν_{0} = 4.267 mHz)) as computed by our model (solid lines), and calculated using the ADIPLS code (dashed lines). The radial axis is zoomed to show only the outermost region. The eigenfunctions have been normalised so that their maximum over the entire solar interior equals 1. 

Open with DEXTER  
In the text 
Fig. 4. Velocity amplitudes of radial acoustic modes as computed by our model, using Eq. (30) (black dashed line), and as observed by GOLF (black points). The data points are taken from Baudin et al. (2005). The free parameters in the model have been tuned to obtain the best possible agreement with observational data. 

Open with DEXTER  
In the text 
Fig. 5. Asymmetry parameter B as a function of frequency obtained by the “theoretical spectrum” model, for different values of k_{0, int}, with λ = 1 and R_{k} = 2 fixed. The subphotospheric injection scale k_{0, int} is expressed in Mm^{−1}. 

Open with DEXTER  
In the text 
Fig. 6. Same as Fig. 5, but only R_{k} varies, λ = 1 and k_{0, int} = 2 Mm^{−1}. 

Open with DEXTER  
In the text 
Fig. 7. Same as Fig. 5, but only λ varies, R_{k} = 2 and k_{0, int} = 2 Mm^{−1}. 

Open with DEXTER  
In the text 
Fig. 8. Asymmetry parameter B as a function of frequency obtained by the “numerical spectrum” model, for several values of λ. 

Open with DEXTER  
In the text 
Fig. 9. Illustration of the importance of source position with respect to nodes and antinodes of the eigenfunction associated to a mode to explain its asymmetry. The blue and red lines show the radial profile of the oscillation for two angular frequencies very close to one another (ω_{red} < ω_{blue}). The bold vertical dashed lines show two source positions generating opposite mode asymmetries: positive for r_{1}, negative for r_{2}. The third vertical dashed line marks the edge of the acoustic cavity r = a. 

Open with DEXTER  
In the text 
Fig. 10. Radial location of the nodes (red symbols connected by a dashed line) and antinodes (blue symbols connected by a dashed line) that are closest to the superadiabatic peak, for each radial mode between n = 16 and n = 30. The vertical black line represents the maximum of the superadiabatic peak, where the excitation is most efficient. Horizontal black dotted lines are added for readability only. 

Open with DEXTER  
In the text 
Fig. 11. Relative difference (in percentages) between the asymmetry parameter B obtained by respectively taking into account and discarding the cross term C(ω) (k_{0, int} = 2 Mm^{−1}, R_{k} = 2 and λ = 1), as a function of frequency. 

Open with DEXTER  
In the text 
Fig. 12. Asymmetry profile obtained by the “numerical spectrum” model (solid red line), compared to the observed asymmetry profile (blue crosses). The data points are taken from Barban et al. (2004). For more readability, only data points corresponding to modes with angular degrees 15 ≤ l ≤ 20 have been retained. The asymmetry profile is given in terms of the parameter B (top panel), and alternatively in terms of the parameter χ (bottom panel), which is defined by Eq. (29). The error bars in the top panel correspond to the uncertainty on the observed values of the mode linewidths, which propagates to the asymmetry parameter B. 

Open with DEXTER  
In the text 
Fig. 13. Asymmetry profile B(ν) obtained by the “numerical spectrum” model (with λ = 0.5), as would be observed at three different heights in the atmosphere, corresponding respectively to the observation heights of GOLF (dashed red line), MDI/GONG (dashed blue line) and HMI (dashed black line). 

Open with DEXTER  
In the text 
Fig. 14. Reduced frequency bias δν, defined by Eq. (34), as a function of the asymmetry parameter B for each radial pmode between n = 6 and n = 30. The blue dashed line shows the best linear fit. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.