Modelling the asymmetries of the Sun’s radial p-mode line profiles

Context. The advent of space-borne missions has substantially increased the number and quality of the measured power spectrum of solar-like oscillators. It now allows for the p-mode 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, point-like 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 p-mode 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 p-modes 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 p-mode 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.


Introduction
Solar-like 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 Lorentzian-shaped peak centred around their eigenfrequencies. This idealised line profile has been extensively used to fit observations (see e.g. Jefferies et al. 1991). However, as the resolution reached in helioseismic measurements (both ground-based and space-borne) has increased, it has become apparent that the observed line profiles feature a certain degree of asymmetry (see e.g. Duvall et al. 1993 for observations made at the geographic South Pole; Toutain et al. 1998 for data from the MDI and SPM instruments aboard the SOHO spacecraft).
Since the discovery of this skew symmetry in solar p-mode line profiles, several studies have been devoted to explaining this feature. In particular, it had been recognised early on that a source of excitation that is highly localised compared to the mode wavelength (which we refer in the rest of the paper as 'source localisation') could lead to a certain degree of mode asymmetry, depending on the position of the source (Gabriel 1992;Duvall et al. 1993;Gabriel 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;. Furthermore, Duvall et al. (1993) noticed an inversion of the sense of asymmetry between spectrometric and photometric measurements, with line profiles in the velocity spectrum featuring more power in their low-frequency wing than in their high-frequency wing and vice-versa for line profiles in the intensity spectrum. Since intensity perturbations were expected to be proportional to velocity perturbations, one would have expected the asymmetries to be the same. Many hypotheses were posited to explain this puzzling result. Duvall et al. (1993) suggested that it was due to non-adiabatic effects lifting the proportionality relationship between the two kinds of perturbations (fluid displacement and temperature) but this hypothesis was later contradicted by Rast & Bogdan (1998). Non-adiabaticity was brought up again later on by Georgobiani et al. (2003) who suggested that the explanation resided in radiative transfer between the mode and the medium. Indeed, the observed radiation temper-Article number, page 1 of 26 arXiv:2001.10271v1 [astro-ph.SR] 28 Jan 2020 A&A proofs: manuscript no. 36847corr ature corresponds to the gas temperature at local optical depth τ = 1. But optical depth depends on opacity, which non-linearly depends on temperature. Therefore, the temperature fluctuations due to the oscillating mode entails opacity fluctuations, which in turn impacts the 'observed' radiation temperature. Given the non-linear nature of the κ − T relation, this modulation decreases the observed temperature fluctuations more significantly in the low-frequency wing of the mode than in its high-frequency wing. Since this radiative transfer does not impact the velocity measurements, this could explain the asymmetry reversal between velocity and intensity spectra. Using 3D simulations, Georgobiani et al. (2003) computed mode line profiles in both the velocity and the intensity power spectrum alternatively at mean unity optical depth and instantaneous unity optical depth. Their results indeed show that the modulation of the 'observed' intensity fluctuations due to radiative transfer could be significant enough to reverse the sense of mode asymmetry. One of the hypothesis enjoying the most support for asymmetry reversal, however, is based on the effect of turbulent perturbations partially correlated with the mode, which thus impact its line profile Roxburgh & Vorontsov 1997;Rast & Bogdan 1998;Kumar & Basu 1999). Indeed, a part of these perturbations is coherent with the mode and, thus, leads to interference. This interference term may be constructive or destructive, depending on the phase difference between the mode and the coherent turbulent perturbations. For frequencies at which the interference is constructive, the power spectral density is slightly elevated, whereas it drops slightly for frequencies at which it is destructive. Typically, in the vicinity of a resonant mode, the dependence of the phase difference between mode and turbulent perturbation is such that the interference term is constructive for frequencies located in one wing of the mode and destructive in the other. Therefore, as a result of this interference behaviour, one of the wings falls off more slowly and the other more rapidly, leading to mode asymmetry. It has been suggested that the degree of correlation between the turbulent perturbations and the oscillation it excites is higher in intensity than in velocity, so that it changes the sign of mode asymmetry only in the intensity spectrum. While it is widely accepted that correlated turbulent perturbations must be taken into account to explain asymmetries in the intensity spectrum, the question of whether it has a significant impact on the velocity spectrum remains an open issue (see e.g. Jefferies et al. 2003).
The possibility that correlated turbulent fluctuations have an affect on mode asymmetry has led many authors to include them in their models for the power spectrum. Even though correlated noise was introduced to explain the particular puzzle of asymmetry reversal between velocity and intensity measurements, several models include correlated noise in the velocity spectrum as well as in the intensity spectrum. This is the case for the model developed by Severino et al. (2001) and later used, for instance, by Barban et al. (2004), which includes three types of noise (coherent-correlated, coherent-uncorrelated and incoherent) in both the velocity spectrum, the intensity spectrum, and the velocity-intensity cross-spectrum. They considered, however, that the 'pure oscillation' (without the noise) has a Lorentzian shape, thus discarding the contribution of source localisation. This model was later refined by Wachter & Kosovichev (2005) to take this contribution into account.
These prior studies have one thing in common, however, and that is that they all treat the various sources of asymmetry (mainly source localisation and correlated noise) in a simplified, parametrised fashion. Indeed, the excitation is consistently modelled as a point-like source, with radial position and multipolar development left as free parameters. This prescription remains somewhat unsatisfactory in the sense that it does not take into account the finer properties of the source of excitation, such as its spatial extent or its dependence on frequency, for instance. As such, these prior models lack a realistic description of the source of excitation. Likewise, for models including the effect of noise on the power spectrum, the various relative amplitudes and phase differences between modal oscillation and correlated noise in both spectra are also left as free parameters. The approach followed by these studies is to find best-fit values for all their free parameters by fitting their model to observations in order to localise the source.
In contrast, in the present paper, we follow a different approach: we model both the source of excitation and the correlated background by constraining their properties using an analytical model of stochastic excitation, coupled with a 3D simulation of the solar atmosphere. The novelty of our approach lies in the fact that we do not fit a parametrised model to the observations but, instead, we predict the dependence of mode asymmetry on frequency, which we then compare to observations in order to validate our model. Our model of mode asymmetry is, therefore, both more realistic (in its description of the source of excitation) and more complete (in its lack of freely adjustable parameters). It can then be used to deepen our understanding of the underlying physical mechanisms behind mode asymmetry by exploring how varying a given physical constraint impacts the results yielded by our model. Finally, our model allows for a much higher predictability of mode asymmetry, which is essential when it comes to applying these results to other solar-like oscillators. We note that this paper is devoted to the modelling of the velocity power spectrum only and, as a result, we do not address the problem of asymmetry reversal, which is a separate challenge altogether. These efforts to model the line profiles of solar-like oscillations are also necessary in order to correctly infer mode properties from observations. Indeed, it was discovered early on that using a Lorentzian shape to fit skew symmetric line profiles led to a significant bias in the eigenfrequency determination, which may be higher than the frequency resolution achieved by helioseismic measurements (Duvall et al. 1993;Abrams & Kumar 1996;Thiery et al. 2000;Toutain et al. 1998). Such eigenfrequency determination bias has also been revealed for solar-like oscillations in stars other than the Sun by Benomar et al. (2018). Inversion methods used to infer the internal structure of solar-like oscillators, whether they be analytical or numerical, require a very accurate determination of the mode eigenfrequencies. For spectra extracted from very long time series, the resolution is high enough that this bias in eigenfrequencies impacts the results obtained by inversion methods (see e.g. Toutain et al. 1998, who show that the difference between the sound speed squared inferred from symmetric and asymmetric fits can reach 0.3% in the core). When fitting these observations, mode asymmetry must, therefore, be taken into account. Since it has proven very difficult to determine accurate mode eigenfrequency without prior knowledge on its line shape, obtaining an a priori model of p-mode line profiles is of primary importance.
In this paper, we present a predictive model of solar radial 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 non-parametrised way, unlike what was done in previous works (see e.g. Severino et al. 2001). The paper is structured as follows: we present the analytical model of the Sun's velocity power spectral density in Section 2 and its numerical implemen-tation in Section 3. We then present the results yielded by our model concerning mode asymmetry in Section 4. In Section 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 Section 6 and discuss the issue of eigenfrequency determination bias entailed by the skewness of the mode line profiles.

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

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 Sun-as-a-star measurements. We note, however, that other instrumental effects -in particular mode leakage -are not accounted for.
To derive an expression for the observed power spectral density, we separate the total surface velocity into an oscillatory part v osc and a turbulent part u, where it is understood that the modes are described by the oscillatory part. The observations made for the Sun-as-a-star are obtained by integrating the velocities over the entire solar disk. Neither the mode velocity (for radial modes), nor the turbulent perturbations depend on the point of the disk at which it is estimated; however, the projection on the line of sight n does. This integration over the solar disk is affected by limb-darkening h(µ) (where µ refers to the cosine of the angle between the local radial direction and the line of sight). Furthermore, since it is the turbulent perturbations that excite the mode, a certain fraction of the former must be correlated with the latter, so that the contribution of turbulent perturbations to the velocity spectrum must be considered.
With these considerations, the observed velocity power spectral density can be expressed as 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 limb-darkening h(µ) 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 | u n | 2 is negligible compared to the terms that contain | v osc | 2 and Re u n v osc , 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. Turck-Chiè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. Roxburgh & Vorontsov 1997;Kumar & Basu 1999), although its importance in velocity measurements is not as clear.

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 Article number, page 3 of 26 A&A proofs: manuscript no. 36847corr and the acoustic potential and source term are given by S (r) = r c ρ 0 (r) dp t dr , where r is the radial coordinate, ρ 0 is the density, N is the Brunt-Väisälä frequency, g 0 is the gravitational acceleration, G is the gravitational constant, and p t refers to the fluctuations of the Reynolds stress around its mean value. Indeed, only the fluctuating part of the Reynolds stress contributes to the source term S (r) and its mean value only modifies the equilibrium structure. The subscript 0 refers to the equilibrium structure and all the above quantities are dependent on the radius at which they are estimated, even when not explicitly specified. We note that we only model radial modes in this paper, so that the wave equation (Eq. 5) is of the second order despite the fact that we did not use the Cowling approximation. Mode damping is not included in Eq. (5). Indeed, we did not take into account the feedback of modal oscillations on the equilibrium state through modulations in the fluid density, pressure, opacity, etc. Such feedback allows mechanical work and thermal transfer to occur from the mode to the medium in which it develops; depending on the phase-lag between these different modulations energy can be exchanged with the surrounding medium. However, the modelling of damping rates of solar-like oscillations is extremely difficult (Samadi et al. 2015). Thus, we directly introduce damping in the wave equation in the form of a mode lifetime, or, equivalently, by a (frequency-dependent) linewidth Γ ω , so that the wave equation takes the following form where j denotes the imaginary unit and the linewidths Γ ω are inferred from observations. We used the line-widths presented in Baudin et al. (2005) (see their Table 2), which were inferred from GOLF data. Note, however, that our definition of Γ ω corresponds to their Γ multiplied by 2π, or equivalently to twice their damping rate η. We completed these data with low-frequency 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.

Expression of the velocity power spectral density
By definition, the Green's function G ω (r o , r s ) is the value taken by the function Ψ ω at the radius r = r o (the variable r o refers to the height in the atmosphere at which the spectrum is observed and the variable r s refers to the position of the point- 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 v osc (r o ) = jω Using the source term given by Eq. (7) in Eq. (10), and Eq. (11) and after finally performing an integration by part, we write the velocity Fourier transform at angular frequency ω as v osc (ω, 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 v osc (ω) 2 and C(ω) are given, respectively, by Eqs.
(B.19), and (B.28). We note that the effects of limb-darkening and disk integration are now contained in a single factor and, thus, these will only have an effect on mode amplitude. Since the asymmetry of a mode does not depend on its amplitude, it is not impacted by such a factor. The calculations leading from Eq. (4) to Eq. (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.

Closure models
The calculations leading from Eq. (4) to Eq. (13) involve the evaluation of fourth-order and third-order two-point correlation moments of the turbulent velocity. Therefore, an appropriate closure model is needed to express these high-order moments as a function of second-order moments. We devote the following subsection to presenting and developing these closure models.
Fourth-order moments To describe the fourth-order correlation moments of the turbulent velocity, we make use of the Quasi-Normal Approximation (QNA hereafter). This closure model consists in considering that all turbulent quantities are normally distributed, in which case fourth-order moments can be analytically expressed as a combination of second-order moments (Lesieur 2008): where a, b, c , and d refer to any turbulent scalar quantity. Applying the QNA to isotropic, homogeneous turbulence inhibits energy transfers among modes of different wave numbers, thus leading to violations of the energy conservation principle (Kraichnan 1957). This is due to the fact the QNA entails vanishing third-order correlation moments. When it comes to estimating the fourth-order moments, however, the picture is different. Belkacem et al. (2006b) have studied the validity of the QNA for two-points, fourth-order correlation moments of the vertical turbulent velocity, in the form of u 2 r,1 u 2 r,2 (where the indices 1 and 2 refer to two different points in space), using 3D simulations of the solar atmosphere. They found that the dependence of this correlation moment on the distance ∆X between the two points is correctly estimated by the QNA but that its absolute value (which can be taken as the corresponding one-point moment) is not. Consequently, the amplitude of the modes are largely underestimated when the QNA is used. However, the asymmetry of the modes does not depend on their amplitude, so that mode asymmetry will be unaffected by a discrepancy in the absolute value of the two-points, fourth-order moments. As such, the decomposition given by Eq. (14) can be considered valid when it comes to studying mode asymmetry.
Third-order moments While the QNA provides an adequate closure relation for fourth-order moments, as mentioned above, it assumes vanishing third-order moments. Therefore, in order to estimate these third-order moments, we make use of another closure model, the Plume closure model (PCM hereafter), which was developed by Belkacem et al. (2006a). The idea behind this closure model is to separate the flows directed upwards from those directed downwards (the latter being referred to as plumes) and to apply the QNA to both separately. The anisotropy between the two types of flow -in particular, turbulence is more prominent in the downwards plumes (e.g. Goode et al. 1998) -yields non-vanishing third-order correlation moments: 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 u d is the fluctuation of the vertical velocity around its mean value in the downflows.
We note that, strictly speaking, the third-order moment given by Eq. (15) and yielded by the PCM are centred. However, we consider that the mean value of the overall vertical velocity of the flow is sufficiently low (compared to its standard deviation for instance) to be neglected. Therefore, these moments may interchangeably refer here either to centred or non-centred moments.
We also note that this closure relation is written here in terms of u d (i.e. the turbulent fluctuations in the downflows only). It would be more practical to rewrite it in terms of u r (i.e. the total turbulent fluctuations). The two are related through u d (R, t) u d (R+r, t+τ) = 1 1 − a u r (R, t)u r (R+r, t+τ) −aδu 2 . (16)

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 v osc (ω) in Eq. (12) contains the turbulent velocity fluctuations squared, expanding the square of its modulus will raise these fluctuations to the fourth. The contribution of turbulence to the expression of v osc thus takes the form of two-points, fourth-order correlation moments of the turbulent velocity. We use the closure relation presented and detailed in Subsection 2.3.1 to express them as a function of second-order moments.
We then use analytical expressions for the second-order moments of the turbulent velocity. We describe the second-order moment of the i-th and j-th component of the turbulent velocity in terms of its spatial and temporal Fourier transform φ i j (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 i-th and j-th component, and δ i j is the Kronecker symbol. The integration over the solid angle of wave vectors k is straightforward, and only an integral over the norm of k remains. However, solar turbulence close to the photosphere is known to be highly anisotropic. To take this anisotropy into account, we follow the formalism developed by Gough (1977). In this formalism, the integral over the solid angle of k is simply readjusted by adding an anisotropy factor given by Eq. (B.10) (see Appendix B in Samadi & Goupil 2001). Following Stein (1967), we then decompose E(k, ω) into a spatial part E(k), which describes how the turbulent kinetic energy is distributed among modes of different wave numbers, and a temporal part χ k (ω), which describes the statistical distribution of the life-time of eddies of wavenumber k In order to model the spatial and temporal part of the spectrum of turbulent kinetic energy, we followed two different approaches, described in the following.
The 'theoretical spectrum' model We use theoretical prescriptions to model both the spatial spectrum E(k) and the temporal spectrum χ k (ω) of turbulent velocity. Based on the assumption that turbulent flows are self-similar, Kolmogorov's theory of turbulence leads to a spatial spectrum E(k) ∝ k −5/3 in the inertial range, between k = k 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 hereafter) was introduced to account for the broadness of the maximum of E(k). The BKS can be written as where u 2 0 ≡ u 2 (r) /3 and the 0.652 factor is introduced so that the total specific kinetic energy of the turbulent spectrum matches u 2 0 /2. Therefore, the spatial spectrum is parametrised solely by the injection scale k 0 . However, the injection scale varies significantly between the sub-surface layers and the atmosphere (Samadi et al. 2003), so that we keep it free in our model and allow for it to depend on the radial coordinate r.
Following Samadi et al. (2003), we consider a Lorentzian shape for the temporal spectrum χ k (ω), which is supported both by numerical simulations (Samadi et al. 2003) and by theoretical arguments. Indeed, a noise described by a stationary, Gaussian Markov process in time is expected to relax exponentially, meaning that the resulting eddy-time correlation is expected to be a decreasing exponential, and its Fourier transform a Lorentzian function (Belkacem et al. 2011). The width ω k associated to eddy-time correlation is linked to the life-time of the eddies of wavenumber k. Dimensional arguments would suggest that ω k ∝ ku k , where u k is the typical velocity associated to the eddies of wavenumber k. However, there remains a substantial indetermination on the actual value of ω k , so that, following Balmforth (1992), we consider: 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 u 2 0 = u 2 (r) /3. The temporal spectrum χ k (ω) is also treated in a slightly different manner than in the 'theoretical model' above. Indeed, the arguments invoked above to justify the Lorentzian shape of the spectrum, while valid for most of the relevant time scales associated to turbulent eddies, are no longer valid for shorter time scales, that is, for higher angular frequencies. Belkacem et al. (2010) argued that if the time correlation associated to small eddies indeed originates from their advection by larger, energy-bearing eddies -a hypothesis referred to as the sweeping assumption -one recovers a Gaussian spectrum instead of a Lorentzian one. The transition between a Lorentzian spectrum, valid for low angular frequencies, and a Gaussian spectrum, valid for high angular frequencies, occurs at the cut-off angular frequency ω E , which is given by the curvature of the eddy-time correlation function at τ = 0 (Belkacem et al. 2010): 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 cut-off frequency into account significantly decreases the predicted amplitudes for high frequency modes. As far as mode asymmetry is concerned, we found that it did not have a significant impact in the 'theoretical spectrum' model. In contrast, it substantially changes mode asymmetry at high frequency in the 'numerical spectrum' model, which is why we only introduce it in the latter. The reason is the following: the spatial spectrum of turbulent energy falls off much more rapidly with k in the 'theoretical spectrum' than in the 'numerical spectrum'. Therefore, small spatial scales play a more important role in the latter. Since the cut-off frequency only impacts these small scales, it is natural that taking it into account should only have a significant impact on the 'numerical spectrum' model.

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, third-order 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 two-point, second-order moments of the turbulent velocity. We then use the same analytical description of second-order moments as the one we used for the leading term. The rest of the calculations is very similar to those described in Sect. 2.3.2 and leads to the second term in Eq. (13).
These two models -the 'theoretical spectrum' and 'numerical spectrum' models -are complementary in the sense that the first one allows us to explore the impact of the properties of turbulence on mode asymmetries and gives physical insight into this problem, whereas the second one relies on fewer input parameters and, therefore, has more predictive capability (we recall here that the former requires the parameters λ and k 0 (r) to be set, while the latter only requires λ). Consequently, in the following, we present and develop the results yielded by both.

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

The solar equilibrium state
The acoustic potential given by Eq. (7) depends only on the equilibrium structure of the Sun. We extracted the potential from a 1D solar model provided by the evolutionary code CESTAM (Morel 1997;Marques et al. 2013). The 1D model includes treatment of the convective flux (using standard mixing-length theory with no overshoot) and of the radiative flux (using the Eddington approximation). On the other hand, turbulent pressure, rotation, and diffusion processes are ignored.
However, 1D stellar models do not fully take into account the more complex physical phenomena taking place in the uppermost layers of a star; in particular, the rapid transition between the convective zone and the superficial radiative region (Kupka & Muthsam 2017). This leads to significant biases in the equilibrium structure. Since the excitation of solar oscillations precisely takes place in these layers, it is essential that we model them more accurately. To do so, we use a 3D hydrodynamic simulation of the solar atmosphere computed using the CO 5 BOLD code (Freytag et al. 2012). The modelled region includes the super-adiabatic peak just below the photosphere and goes up to the lower atmosphere of the star.
It is now possible to construct a 'patched' model of the solar interior. We use the solar patched model computed by Manchon et al. (2018). The process of constructing patched models has been extensively discussed (e.g. Trampedach 1997;Samadi et al. 2008) and the particular case of the patched model used in this paper is described in much detail in Manchon et al. (2018). The basic idea is to transform the 3D atmosphere into a 1D atmosphere through temporal and horizontal averaging and then to replace the surface layers of the 1D stellar interior with this 1D atmosphere. We note that the input parameters of the CES-TAM model used to describe the solar interior (age, total stellar mass, mixing-length parameter α MLT , and helium abundance) are chosen so that the top layers match the bottom layers of the CO 5 BOLD atmosphere. Here the model was computed with the mixing-length parameter α MLT = 1.65, an initial helium abundance of Y init = 0.249, and an initial metallicity of Z init = 0.0135. Figure 1 shows the acoustic potential profile V(r) given by this solar patched model, computed using Eq. (7).
Finally, we use the same simulation of the solar atmosphere to extract the various parameters appearing in the analytical description of the source term; in particular, the standard deviation u 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).

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  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 R −2 here (where R is the radius of the solar photosphere) so that it is dimensionless.
wave equation was integrated and how we extracted its Green's function.
For a given angular frequency ω, we carry out the integration using a 4th-order Runge-Kutta scheme (Press et al. 1986) with the acoustic potential V(r) given by the solar equilibrium state (see Sect. 3.1). Given that radial modes develop in the entire solar volume, we perform this integration over the entire solar radius, between r = 0 and r = r max . We note that r max refers not to the photospheric radius, but to the maximum radial extent of the solar model described in Sect. 3.1, so that r max > r photosphere .
We imposed Dirichlet boundary conditions on the wave equation. The condition at the centre is straightforward: by definition, Ψ ω (r = 0) = 0. At r max , we impose a vanishing Lagrangian pressure perturbation (which physically means that the atmosphere of the Sun is force-free). The continuity equation and the equation of state allow us to rewrite this latter condition in terms of Ψ ω : The use of Dirichlet boundary conditions leads us to implement a shooting method: we perform the integration with Ψ ω (r = 0) = 0, and tune the initial slope (i.e. the value of dΨ ω / dr at r = 0) until the other boundary condition is met. Note that this method is not similar to the shooting method usually implemented to solve the eigenvalue problem associated to the determination of acoustic mode eigenfrequencies: here, the pulsation ω is fixed, and it is the initial slope that is tuned to meet the boundary condition at the surface. The difference between these two methods is that in the inhomogeneous problem, the initial slope (or, alternatively, the amplitude of the mode) is fixed by the amplitude of the source of excitation. The shape of the eigenfunction, however, remains the same as in the homogeneous problem.
This method enables us to extract the Green's function associated to the wave equation (Eq. 8). To obtain one value of the Green's function G ω (r s ), for one value of the pulsation ω and one value of the source position r s , we carry out the integration of the inhomogeneous equation as described above, adding a point-like source term to the numerical scheme. The source is normalised in such a way that the right-hand side equals 1/h when the source falls within the integration radial step, and 0 otherwise (h is the radial step of the integration).
This integration gives us the radial oscillation profile Ψ ω (r), and we simply extract its value at a fixed radius r o , which corresponds to the height in the atmosphere at which the spectrum is measured. We note that the presence of damping in the wave equation implies that it is complex-valued. As such, the Green's function is complex-valued as well.
Finally, to calculate the integrals over source positions which appear in Eqs. (B.19) and (B.28), we compute the Green's function using the above method for a grid of source positions r s , while ω is kept constant. This grid corresponds to the radial grid provided by the 3D atmospheric model described in Sect. 3.1.
We also use the aforementioned 3D model to extract the physical quantities appearing in both the leading term of Eq. (13) (the anisotropy factor G, the turbulent velocity fluctuations u 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 ω.

Fitting of the mode asymmetries
We fit the line profile of the modes following  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 low-frequency side (that corresponds to negative asymmetry), B > 0 means that the high-frequency side contains more power (that corresponds to positive asymmetry), while with B = 0 we recover a classic, Lorentzian profile. Figure 2 also shows that the mode does not peak exactly at the eigenfrequency, but rather at a slightly higher (for B > 0) or lower value (for B < 0). This can have important repercussions for the determination of the mode eigenfrequencies from observations, as we discuss in Section 6.3. Note that several other definitions of the asymmetry parameters can be found in the literature. Korzennik (2005) prefers to adjust the mode line profiles with 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 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 half-way between the neighbouring maxima lead to symmetric line profiles, while a deviation from this behaviour leads to asymmetric line profiles. The parameter χ derived from these considerations is independent from both the amplitude and the line-widths of the modes.
The formulas presented above only lead to different line profiles far from the central frequency, whereas they are equivalent in our range of interest. We opted for the definition given by 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. 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.

Validation of the method
Using the method presented above, we extracted solar radial modes of radial orders n = 6 to n = 30. Indeed, the formula used for the fit and given by Eq. (26) does not converge properly for higher-order modes (because the increasing linewidths lead to mode overlapping), while we did not have access to observed linewidths for lower-order modes. In addition to their line profile asymmetries, we also extract other fundamental properties, namely their eigenfrequencies, amplitudes, and eigenfunctions. In the following, to support the validity of our model, we compare the mode properties we obtained with similar properties obtained through other methods.
First, we compare the eigenfrequencies obtained through our model to the eigenfrequencies of the 1D adiabatic oscillations calculated using the ADIPLS code (Christensen-Dalsgaard 2011). For this validation, we did not make use of the patched model described in Sec. (3.1) but, rather, the corresponding unpatched model. The reason is that the patching procedure produces a small discontinuity of the physical quantities at the patching point, which can affect the eigenfrequencies calculated by ADIPLS. We recover the correct eigenfrequencies, with errors not exceeding ∼ 0.1%. Since mode asymmetry is only expected to vary on the scale of ∼ mHz, modelled asymmetries will not be significantly affected by such small discrepancies of the eigenfrequencies.
Our numerical method also allows us to extract the radial profile Ψ(r) of the eigenmodes. We compare them in Fig. 3 to the eigenfunctions calculated using the same 1D adiabatic oscillations obtained through the ADIPLS code and presented above. The figure shows that the modes obtained through our model have eigenfunctions that are very similar to those obtained through this dedicated code, which further supports the validity of the model we have used.
Finally, we compare the mode amplitudes obtained through our model to the observed ones. To that end, we estimated the ve- locity 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 low-frequency and the high frequency limit. To conclude on the matter, we emphasise that the asymmetry parameter B is independent from the mode amplitude, so that potential discrepancies concerning the latter should not affect the former.

Results for the asymmetry profile B(ν)
Using the model presented in Sect. 2, numerically implemented using the method presented in Sect. 3, we extract the solar p-modes line profile asymmetries B(ν) throughout a large part of the spectrum, between n = 6 (ν ∼ 1 mHz) and n = 30 (ν ∼ 4.3 mHz). In this section, we present the results yielded by our model, focusing on the dependence of the asymmetry parameter B on frequency (which we hereafter shorten to the expression 'asymmetry profile') and on the impact of our different input parameters on the asymmetry profile.
As we detailed in Sect. 2, we followed two different approaches to model the turbulent kinetic energy spectrum. The first one, which we refer to as the 'theoretical spectrum' model, uses the prescription given by Kolmogorov's theory of turbulence and which we have described in detail in Sect. 2. The second approach, which we refer to as the 'numerical spectrum' model, uses the turbulent spectrum extracted from the 3D hydrodynamic simulation of the solar atmosphere described in Sect. 3. In this section, we present separately the asymmetry profile B(ν) yielded by both models.

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 energy-bearing 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 illustrates the main qualitative features of the asymmetry profile B(ν). In fact, together with Figs. 6 and 7, it shows that the qualitative behaviour of the asymmetry profile is largely model-independent. Thus the asymmetries of the solar radial 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. Fig. 6 shows how the asymmetry profile B(ν) depends on R k . An increase of this parameter attenuates low-frequency mode asymmetries (below ν max ∼ 3 mHz), while on the contrary, the high-frequency modes (above ν max ) become more asymmetric. The effect is significantly more substantial for the latter than for the former. Asymmetries close to ν max , however, are not affected by the parameter R k whatsoever. Figure 7 shows that the impact of λ on the asymmetry profile B(ν) is quite similar, albeit inverted, in the sense that |B| increases with λ for low-frequency modes and decreases for 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 p-modes are either in the low frequency part or in the high frequency part of the spectrum. We do not go into too much detail here as we discuss this matter further in Sect. 5.

The 'numerical spectrum' model
In the 'numerical spectrum' model, which describes the properties of turbulence more realistically, there is only one input parameter left, λ. In this sense, it has a greater predictive power than the previous model. The qualitative behaviour of the asymmetry profile B(ν) and, in particular, the positions of the different local extrema featured by B(ν), are, in this model, rather independent from λ and in agreement with what we observed in the scope of the previous model. However, the input parameter λ does have an impact on the quantitative behaviour of the asymmetry profile B(ν). We show in Fig.8 the asymmetry profile B(ν) obtained with the 'numerical spectrum' model (see Sect. 2) for several values of λ. As for the 'theoretical spectrum' model, B is always negative and features several local extrema at ν ∼ 1.7, 3 and 4 mHz.
As for the dependence of B(ν) on λ, two distinct regimes can be separated. Below ν max ∼ 3 mHz, we recover the same dependence of the asymmetry parameter B with λ as we obtained in the scope of the 'theoretical spectrum' model, with absolute values of B increasing with λ. The picture at frequencies higher than ν max is, however, somewhat different. The asymmetry profile B(ν) features a local minimum at ν ∼ 4 mHz; the curve inflexion grows sharper as λ increases up to λ ∼ 1, after which this part of the asymmetry profile does not significantly depend on λ. In this sense, the asymmetry profile B(ν) seems to undergo the same saturation behaviour as in the 'theoretical spectrum' model (see Sect. 4.1), for the same value λ sat ∼ 1. The fact that we recover approximately the same threshold gives us confidence that this particular feature of the asymmetry profile B(ν) is not a mere artefact of one model or the other but, rather, it is a genuine effect based on a physical origin. Again, we postpone the discussion of the physical origin of this behaviour to Sect. 5. Article number, page 11 of 26 A&A proofs: manuscript no. 36847corr

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

Effect of source localisation on mode asymmetry: generic arguments
The fact that the source of excitation of a mode is spatially localised can affect the skewness of the mode line profile in Fourier space. There are several ways of describing the impact of source localisation on mode asymmetry.
One way is to make use of the analogy between the development of acoustic modes in the stellar cavity and the phenomenon of optical interference in a Fabry-Pérot cavity. This analogy was used to account for the acoustic mode asymmetry in the Sun by Gabriel (1992); Duvall et al. (1993), among others. The idea is that acoustic, stationary modes in the Sun can be described by means of two progressing waves, propagating in opposite directions. Each of these waves follows the same cycle: they propagate one way, get refracted on the lower turning point of the acoustic cavity, then propagate backwards, get reflected on the upper turning point, and so forth. As a result of these multiple reflections and refractions on both turning points, the acoustic waves pass multiple times through the same regions and, therefore, interfere with each other (and with themselves). This interference pattern leads to the development of resonant modes in the cavity. What we observe then is the evanescent tail of these modes in the atmosphere, which lies outside the resonant cavity.
Let us now consider that the source of the waves is located at a certain point within the cavity. The waves propagating outwards and inwards will have travelled over different distances before interfering with one another and this difference of travel times will depend on the location of the source. The shape of the mode line profile is directly related to the dependence of the phase difference between the outwards and inwards interfering waves on frequency. Since this phase difference is not exactly symmetric about the mode eigenfrequency, neither is its line profile; and given that it depends on the source location, mode asymmetry is indeed a marker of source localisation.
Another physical interpretation of how source localisation can bring about mode asymmetry has been proposed by Rast & Bogdan (1998), and later refined by Rosenthal (1998). They remarked that mode asymmetry could be mathematically described by the relative position of local maxima (or peaks) and local minima (or troughs) in the power spectrum. Peaks located exactly halfway between their neighbouring troughs feature symmetric, Lorentzian line profiles. However, if one of the neighbouring trough is closer than the other, the peak in question appears skewed and, depending on which trough is closest, its asymmetry parameter is either positive or negative.
The position of the peaks are simply related to the eigenmodes of the solar acoustic cavity. As for the position of the troughs, in the special case of a point-like source of excitation, with a given multipolar decomposition, the authors showed that it is related to the eigenmodes of the atmosphere truncated at the source position, with a vanishing external boundary condition depending on the multipolar nature of the source. In that interpretation, the position of the troughs thus depends on both the position and the multipolar decomposition of the source.
Yet another way to describe the impact of source localisation on mode asymmetry is to consider the eigenfunction of the mode. In order to illustrate this, we present in Appendix C a very simplified toy-model of mode excitation, where the source is considered point-like and the acoustic cavity is simplified to a square well potential. From this toy-model we draw the following conclusion: for a given frequency, the amplitude of the wave is proportional to the eigenfunction associated with the wave at the source of excitation. In particular, excitation at a mode's antinode is much more efficient than at a mode's node.
With this conclusion in mind, let us consider the situation illustrated by Fig. 9. The blue and red curves represent the radial profile of the acoustic wave for two different angular frequencies. It can be seen that an increase of ω causes the radial profile of the oscillation to 'shrink' radially. Therefore, the amplitude of the oscillation as seen by the source will either increase or decrease with ω, depending on its position. More specifically, a source at r = r 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.
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 toy-model in the present subsection. In the following sections, we return to the discussion of our model, simply using the conclusions drawn above to interpret the results which it yields.

Correlated turbulent fluctuations
Acoustic modes in the Sun are excited by fluctuations of turbulent nature -more specifically by turbulent fluctuations of the Reynolds stress or non-adiabatic pressure perturbations. It is therefore natural that a part of the turbulent fluctuations should be not only coherent, but statistically correlated with the oscillating mode.
The resulting interference between the mode and the turbulent fluctuations leads, in turn, to mode asymmetry. In order to illustrate this, let us consider a mode whose line profile is intrinsically Lorentzian and turbulent fluctuations whose power spec-tral 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 right-hand side of Eq. (32) corresponds to a Lorentzian profile and is symmetric about x = 0. The third term simply acts as an offset and does not introduce any mode asymmetry. The second term, however, is clearly not symmetric at x = 0, unless φ n = ±π/2. For instance, if φ n = 0, this term is even antisymmetric. In other words, the interference between the mode and the noise is destructive in the left wing of the mode and constructive in its right wing. As such, the power spectral density P(x) is higher than the Lorentzian profile in the right wing and lower in the left wing, thus entailing positive mode asymmetry. The sign and magnitude of the mode asymmetry depends on the amplitude A n , that is, on the degree of correlation between the mode and the turbulent fluctuations, as well as on the phase difference φ n , both of which are included in the model we developed in Sect. 2.

Contribution of source localisation to B(ν)
In the previous subsection, we summarised the impact of source localisation on mode asymmetry by stating the following: a source within the resonant cavity of a mode entails negative asymmetry if it is located above an antinode and below a node of the associated eigenfunction and positive asymmetry otherwise; a source outside the resonant cavity always entails negative mode asymmetry. With this in mind, we set out to interpret the results obtained in Sect. 4 in the scope of the 'theoretical spectrum' model.
Once applied to the case of solar p-mode excitation, this rule can be rephrased in the following way. There is a dichotomy between the effect of the turbulent eddies located below the upper turning point of the mode and those located above. The former skew the mode line profile one way or the other depending on their height relative to the nodes and antinodes of the mode eigenfunction. The latter always skew the mode line profile so that it feature negative asymmetry. Until now, we have only discussed the case of a point-like source of excitation. However, the driving region of the solar p-modes, while localised around the super-adiabatic peak just below the photosphere, has a certain spatial extent. As such, driving turbulent eddies can be found both below and above the upper turning point of the modes and the observed mode asymmetry is due to the combination of both.
As we mention above, the sense of asymmetry created by turbulent eddies below the upper turning point depends on their position with respect to the nodes and antinodes of the modes. Since the wavelength of the modes is much larger than the spatial extent of the driving region, it is sufficient to study the position of the super-adiabatic peak with respect to the nodes and antinodes of the mode eigenfunction. We illustrate this in Fig.  10, which shows the position of the node and antinode of the radial modes obtained through our model, with respect to the super-adiabatic peak, where the stochastic excitation is mainly located. Below ν ∼ 3.5 mHz, the super-adiabatic peak lies above the closest antinode (dashed blue line in Fig. 10), thus generating negative asymmetries. Above the aforementioned threshold, the closest antinode is above the super-adiabatic peak and it generates positive asymmetry. Close to 3.5 mHz, the super-adiabatic peak coincides with an antinode, so that the asymmetry is very low.
The dichotomy between turbulent eddies below and above the upper turning point stands thus: excitation localised below the upper turning point makes modes with ν ν max negatively asymmetric and modes with ν ν max positively asymmetric; excitation localised above the upper turning point makes all modes negatively asymmetric; since the source of excitation has a certain spatial extent, the total asymmetry is a combination of both cases. This is in perfect accordance with the dashed blue curve of Fig. 6 which shows the asymmetry profile B(ν) with R k = λ = 1 (i.e. imposing the same turbulent spectrum everywhere). Indeed, this curve shows that B is negative for low frequencies, and positive for high frequencies. If the source of excitation was only located below the upper turning point of the modes, the threshold between the two regimes would be ∼ 3.5 mHz; however, turbulent eddies above the upper turning point generate additional negative asymmetry, thus shifting the curve downwards and increasing the threshold between the B < 0 and the B > 0 regimes.
The results presented in Figs. 5 and 6 are also easily interpreted. Indeed,Eqs. (B.19) and (B.28) show that the efficiency of stochastic excitation scales as k −4 0 , 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 k −4 0 ). 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 p-modes compares differently to the typical turbulent eddy-time correlation, that is, ω compares differently to ω k , which indeed depends on λ. Depending on whether ω < ω k or ω > ω k , the temporal turbulent spectrum χ k (ω) vary differently with ω: indeed, χ k is almost flat for very low frequencies, whereas it decreases as ω −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.

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 ) 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. Fig. 11 shows the relative difference between the asymmetry profile B(ν) when the cross term is respectively taken into account and discarded. This relative difference is highest at low frequency, and almost vanishes when ν > ν max . This can be easily explained by the fact that low-frequency modes have the smallest power spectral density, while power spectral density associated to the noise is higher: in contrast, the effect of correlated turbulent perturbations is at its most substantial at the lowest frequencies. Nevertheless, even in the lowest part of the spectrum, the relative difference in asymmetry between a model with and without the cross term does not exceed 3%, which is much smaller than the dispersion characterising observed asymmetries. We therefore conclude that the dominant source of asymmetry in ve- 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. locity data is the source localisation and that the effect of correlated turbulent perturbations can be disregarded.
It is possible to support this conclusion with a simple order of magnitude estimation. Indeed, B is of order 10 −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 Cauchy-Schwarz inequality provides an upper bound to the cross term: C(ω) v rms (ω) u rms (ω), where the case of equality happens when the correlation between the mode v and the turbulence u is optimal. Therefore v(ω) 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 v 2 rms ∼ 10 5 m 2 .s −2 .Hz −1 and that of the turbulent fluctuations as u 2 rms ∼ 10 m 2 .s −2 .Hz −1 (Turck-Chièze et al. 2004, fig. 2), the above ratio becomes 10 5 /10 ∼ 10 2 . We note that at this stage, the correlation term is about 1% of the leading term, which is barely enough to impact the asymmetries significantly. However, by considering the limiting case in the Cauchy-Schwarz inequality, we assumed that the mode velocity and the turbulent velocity were both optimally coherent and completely independent. This is not, however, the case; consequently we probably overestimated the importance of C(ω) by at least an order of magnitude. Therefore, this crude order of magnitude estimation indeed tends to support the conclusion that correlated turbulent perturbations can be disregarded when interpreting mode asymmetries in the velocity spectrum.

Comparison with observations
Observed properties of solar-like oscillations depend not only on the observable (velocity or intensity), but also on the specifics of each instruments. As far as velocity measurements are concerned, all instruments do not perform their Doppler observations on the same spectral line. Since different spectral lines form at different altitudes in the atmosphere, and since the properties of turbulence change throughout the atmosphere, mode properties, and especially mode asymmetries, may depend on the instrument.
In the following, we compare the results of our model to observations performed with the GONG network. We then focus on the dependence of the asymmetry profile B(ν) on the height at which the velocity spectrum is observed. We use the 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.

Comparison with GONG observations
First -and in order to support the validity of our model in terms of mode asymmetry -we compare the asymmetries yielded by our model to those inferred from observations. We use the data points extracted from the spectrum analysis of Barban et al. (2004). The authors chose an equivalent, albeit different, set of parameters to fit the acoustic modes observed by the GONG network. Therefore, we reconstructed the shape of the modes point by point using the parameters extracted from their fit and fitted them again with the formula given by Eq. (26). We note that the modes analysed in their study are not radial, but have angular degrees ranging from l = 15 to l = 50. It is known, however, that the dependence of mode asymmetry on the angular degree is very weak (see for instance Duvall et al. 1993, who studied modes of angular degree up to l = 170) because the eigenfunctions associated to the acoustic modes are very weakly dependent on l close to the photosphere, as long as l is not too high. Thus, comparing our radial study to their non-radial observations remains relevant.
We showcase this comparison in Fig. 12. The parameter λ was adjusted for a better agreement with the observations. We obtained λ = 0.5 (see Eq. 20 for a definition of this parameter), which is approximately the value obtained by constraining mode amplitudes (Samadi & Goupil 2001). It is rather clear that we reproduce the main features of the asymmetry profile B(ν), especially its sense of variation on the different intervals (when the frequency range of the observations overlaps ours), as well as the positions and values of its local extrema. Alternatively, we show the same comparison in terms of the asymmetry parameter χ in the bottom panel of Fig. 12. This parameter is defined by Eq. (29), and is more robust when it comes to comparing theory and observations because it does not depend on the determination of the mode linewidths, which may introduce extra uncertainty in the determination of the observed asymmetries. In conclusion, the good agreement obtained between our model and observations show that the model developed in this paper is relevant to account for acoustic mode asymmetry quantitatively.

Dependence of asymmetry profile B(ν) on observation height
We consider three different observation heights in this paper associated, respectively, with the Michelson Doppler Imager (MDI), the GOLF instrument (both onboard the SOHO spacecraft), the Global Oscillation Network Group (ground-based), 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. and the Helioseismic and Magnetic Imager (HMI; onboard the SDO spacecraft). We provide further information on the observation height of each of these instruments in Table 2. However, determining the formation height of a given absorption line is extremely difficult in that it does not depend only on the nature of the line (Fleck et al. 2011). Therefore, the values given in Table 2 are not to be considered as precise estimates but, rather, as approximate figures.
In Fig. 13, we compare the asymmetry profile B(ν) as it would be observed by the various instruments listed in Table  2. We control the height by tuning the radial coordinate r o at which the Green's function of the homogeneous wave equation is calculated (see Sect. 2 for more details). Mode asymmetry is only slightly dependent on the observation height. This is not so surprising since mode asymmetry is a global property of the modes. It is noticeable, however, that the difference between these asymmetry profiles is most prominent at high frequency ( 3.5 mHz). This is because by observing the modes at a lower   height in the atmosphere, we effectively increase the contribution of atmospheric turbulence compared to the contribution of 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.

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;Thiery et al. 2000;Toutain et al. 1998), using a symmetric, Lorentzian model to fit asymmetric line profiles results in an appreciable bias in the best-fit parameters, thus rendering inaccurate the frequency determination.
We illustrate this bias in the following. We fit the line profiles we obtained numerically with a formula given either by Eq. (26) or by the symmetric version (i.e. fixing B = 0). Following Abrams & Kumar (1996), we denote the reduced frequency bias as Figure 14 shows the results obtained for all radial modes described in Sect. 4. Given that the asymmetry of solar p-modes in the velocity spectrum are all negative, using a symmetric fit introduces an underestimation of ν 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.
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 ν B 0 0 = ν 0 , the symmetric fit finds it at ν B=0 0 = ν max . One can therefore derive the simple expression δν = B/2, which is in accordance with the linear fit shown in Fig. 14.
The frequency bias introduced when not taking the line profile asymmetry into account can easily reach several percents of the mode linewidth. This bias largely exceeds the frequency resolution achievable in current solar measurements, especially for high-frequency modes, which are the widest.

Summary and conclusion
In this paper, we detail the development of a realistic and predictive model for the asymmetry displayed by solar radial p-modes in the velocity spectrum. The basic idea behind this model is to compute the Green's function associated to the radial acoustic wave equation, as well as its inhomogeneous part (which corresponds to the source of excitation) and to convolve the two to reconstruct the velocity power spectral density point by point. Once the power spectral density is reconstructed, we extract its resonant modes and study their asymmetry. In particular, and unlike previous attempts to such modelling, we included in our model the correlation of the oscillating modes with the fluctuations associated to turbulent velocity.
First, the Green's function associated with the wave equation was computed numerically. We put the wave equation in the form of a 1D stationary Schrödinger equation, whose potential only depends on the equilibrium structure of the Sun. We extracted the acoustic potential from a solar patched model: the solar interior is calculated using the 1D evolutionary code CES-TAM 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 point-like, normalised source to the integration scheme in order to compute the Green's function.
Secondly, the source term of the wave equation being of stochastic nature, we modelled the statistical properties of the source by means of theoretical developments. We made use of the adequate closure relation to express the third and fourthorder correlation products of the turbulent velocity as functions of second-order products; more specifically on their spatial and temporal Fourier transform. We developed two distinct models: one is based on theoretical prescriptions for the spatial and temporal spectrum of turbulent kinetic energy; the other is based on theoretical modelling of the temporal spectrum only, whereas the spatial spectrum is extracted from a 3D hydrodynamic simulation of the solar atmosphere. We refer to the former as the 'theoretical spectrum' model, and to the latter as the 'numerical spectrum' model.
The asymmetry B displayed by the modes in our model drastically depends on their frequency ν. This is because the shape of the eigenfunctions close to the photosphere is very dependent on ν. We find that B is negative throughout the p-mode spectrum, and that its behaviour weakly depends on the input parameters of our model. It drops from −0.01 to −0.05 between 1 mHz and 1.7 mHz, then rises to 0.015 at 3 mHz, and decreases again from 3 mHz to 4 mHz. Above 4 mHz, the behaviour of B(ν) is much more dependent on the value given to our input parameters and, in particular, on the injection scale associated to the turbulent cascade above the photosphere, compared to below the photosphere. This is related to the fact that the contribution of atmospheric turbulence to mode excitation only becomes significant at high frequency, so that only in this part of the spectrum it may have an impact on mode asymmetry.
The asymmetry of the modes can have two different origins: localisation of their source of excitation within a region of lesser spatial extent than the mode wavelength and correlation between the oscillating modes and the fluctuations associated to turbulent velocity. Formally, these two phenomena have the same impact on mode asymmetry, so that they cannot be separated using observational data only. Our model allows us to make this distinction and to study their relative weight in the total mode asymmetry. We find that the correlation with turbulent fluctuations is negligible in the velocity spectrum, and that the observed asymmetries are exclusively due to source localisation. More precisely, we interpret the results of our model in terms of the source position with respect to the various nodes and antinodes featured by the eigenfunctions of the modes. In the case of a point-like source of excitation, mode asymmetry drastically depends on whether it is located within or outside the mode acoustic cavity. In our model, however, the source of excitation has a certain spatial extent, so that the total asymmetry is a combination of the contributions from the source outside and inside the mode acoustic cavity.
We find that it is impossible to interpret even the qualitative behaviour of the asymmetry profile B(ν) by considering that the source of excitation is point-like (either outside or inside the modes cavity). On the contrary, taking into account the spatial extent of the source allows us to reproduce the observed asymmetries, not only qualitatively, but also quantitatively. This positive result shows that our model is indeed relevant to describe -and, more importantly, to predict -acoustic mode asymmetry in solar-like oscillators. It also shows that any model that assumes a point-like source of excitation cannot give reliable results as far as mode asymmetry is concerned. In particular, such a model would predict positive asymmetries for high-frequency modes, whereas observations show that all asymmetries are negative when measured in terms of velocity power spectral density.
Finally, we study the eigenfrequency bias entailed by neglecting to fit observations with an asymmetric profile. We find that for the most asymmetric modes, this bias can reach several percent of the mode linewidth. Therefore, this bias is higher for high frequency modes, which are the widest. In particular, for ν ∼ 4 mHz, the asymmetry parameter is of order B ∼ −0.04, and the linewidth is of order Γ ∼ 10 µHz, so that the eigenfrequency bias is of order δν ∼ 0.2 µHz. This is in perfect accordance with actual biases obtained from observation fit of the solar spectrum (see Benomar et al. 2018, Fig. 6, topmost panel). Since the eigenfrequency bias is most pronounced for higher frequency (because it is proportional to the mode line-width, which is widest at high frequency), it is likely to have a non-negligible impact on inversion methods, especially those based on asymptotic formulae. One must keep in mind, however, that the deviation of the modelled eigenfrequencies from the observed ones, induced by surface effects, largely dominates the eigenfrequency bias entailed by symmetric fits.
In this paper, we have restricted ourselves to the study of solar radial p-modes. Our formalism can be easily adapted to the study of non-radial modes simply by using a non-radial wave equation instead of the radial one. However, since the eigenfunction associated to p-modes in solar-like oscillators are very weakly independent on angular degree l close to the photosphere, which is precisely where the excitation takes place, the mode asymmetry is not expected to vary significantly with l, at least as long as l remains reasonably small. Observational data tend to confirm this (see e.g. Vorontsov & Jefferies 2013, who report that the spectral parameters of individual modes collapse to slowly varying functions of frequency only for modes with l 100).
We only considered one type of acoustic source in this study, that is, the turbulent fluctuations of the Reynolds stress. Indeed, it has been shown by Stein & Nordlund (2001) that this is the dominant source of excitation of solar acoustic modes (see also Chaplin et al. 2005;Samadi et al. 2007;Nordlund et al. 2009). Therefore, our objective was to start by considering only this source. However, further refinements of the model will have to include other sources of excitation, in the form of non-adiabatic, turbulent pressure fluctuations.
Our formalism can also be easily applied to other solar-like oscillators. Comparing the asymmetries featured by the velocity spectra of several solar-like oscillators as modelled by the method presented in this paper and, in particular, the trend followed by mode asymmetry with stellar parameters such as effective temperature or surface gravity, undoubtedly constitutes the next step of this study. In the long run, mode asymmetry may serve as a useful tool for seismic diagnoses of solar-like oscil-lators. However, the one major difference that remains between the solar case and other stars is that the Sun is the only solar-like oscillator for which spectra obtained by spectrometric measurements are sufficiently resolved to allow for a determination of their mode asymmetry. The asymmetry of acoustic modes of all other stars can only be observed in intensity spectra. As has been reported numerous times (see e.g. Duvall et al. 1993), asymmetry in intensity and in velocity spectra are drastically different. It is, therefore, necessary to adapt our formalism to the intensity spectrum, which is another key element of any further considerations on the matter treated here.
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 ) u x and u y referring to the two horizontal components of the turbulent velocity. This anisotropy factor depends on the ratio between horizontal and vertical turbulent velocities, and therefore depends on the slow R variable -or equivalently, on the mass variable m. Under this formalism, the integral over k and ω remains the same as in the isotropic case.
Following Stein (1967), we decompose E(k, ω) into a spatial part E(k), which describes how the turbulent kinetic energy is distributed among modes of different wave numbers, and a temporal part χ k (ω), which describes the statistical life-time distribution of eddies of wavenumber k E(k, ω) = E(k)χ k (ω) . (B.12) 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 self-similar, 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 thereafter) was introduced to account for the broadness of the maximum of E(k). Finally, the BKS can be written thus: where u 2 0 ≡ u 2 (r) /3 and the 0.652 factor is introduced so that the total specific kinetic energy of the turbulent spectrum matches u 2 0 /2.
Following Samadi et al. (2003), we consider a Lorentzian shape for the temporal spectrum χ k (ω), which is supported both by numerical simulations (Samadi et al. 2003) and by theoretical arguments (if the noise is characterised by a timecorrelation function which decays exponentially, its spectrum is Lorentzian). Thus: The width of the Lorentzian is the inverse of the typical correlation time-scale and by dimensional arguments, it is proportional to ku k , where u k is the typical velocity of eddies of wavenumber k. However, there remains a substantial indetermination on the actual value of ω k . To account for this indetermination, we follow Balmforth (1992) and introduce the dimensionless parameter λ, so that 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 k 0 = 0.602u 0 . Under all these assumptions, all further calculations being carried out, we ultimately obtain the leading term, v osc (ω) 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 non-parametrised description of stochastically excited modes of oscillation.
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 non-local, third-order correlation products. While the QNA provides an adequate closure relationship for fourthorder moments, it yields vanishing third-order moments, which is known to lead to serious violations of the energy conservation principle, as well as an impossibility for the turbulent cascade to develop (see Kraichnan 1957).
In order to estimate the third-order moments, we follow Belkacem et al. (2006a) and use the Plume closure model, which consists of separating the flow into upward flows and downward plumes, each normally distributed, with different mean values and standard deviations. In addition, we consider that the downflows are much more turbulent than the upflows (which is supported by Goode et al. 1998, according to whom the intergranular lanes harbour stronger turbulence than the granules themselves at the Sun's surface), and that the two types of flows considered separately have zero third-order correlation products. In Belkacem et al. (2006a), the authors use the same approximations but focused on one-point correlation products; however, the calculations can be easily extended to the two-point correlation products that we need and the model yields 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 u d is the fluctuation of the downflow velocity around its mean value. The only additional approximation we make to adapt these calculations to two-point correlation products is that the parameters of the model a and δu vary on scales much greater than the typical correlation length, which is another illustration of the scale separation approximation, which we have already used for the leading term (see above). Note that strictly speaking, the third-order moment appearing in Eq. (B.25) and yielded by the PCM are centred. However, we consider that the mean value of the overall vertical velocity of the flow is sufficiently low (compared to its standard deviation for instance) to be neglected. Therefore, the moment described by Eq. (B.25) may interchangeably refer either to a centred or non-centred moment.
Also note that this closure relation is written here in terms of u d (i.e. the turbulent fluctuations in the downflows only). It would be more practical to rewrite it in terms of u r (i.e. the total turbulent fluctuations). The two are related through u d (R, t) u d (R+r, t+τ) = 1 1 − a u r (R, t)u r (R+r, t+τ) −aδu 2 . (B.26) Plugging Eq. (B.26) into Eq. (B.25), we obtain a closure model that allows us to write the third-order moments as a function of second-order moments only, after which we can use the same prescriptions for turbulence as we did for the leading term. We note that while this closure relation contains many terms, not many survive the Fourier transform in Eq. (B.24) as all terms not depending on r or τ will not contribute to the Fourier transform. The integral over m appearing in (B.21) becomes dm dX dr u 2 r (r s ) u r (r o ) = −8π 2 aδuρ 0 (r o ) dX dr (r o )φ rr (0, ω) .

(B.27)
Formally, the integration over r makes the value of φ rr at k = 0 appear. Physically, that means only the largest eddies have a significant impact on the correlation between the mode and the turbulent perturbations. Since considering eddies characterised by k = 0 does not actually make physical sense, we considered a comprise by assuming that the largest contributing eddies are those at the injection scale k 0 , so that the correlated turbulent perturbations term finally becomes Re dΩ h(µ) v osc (ω) u n (ω) = C(ω) dΩ h(µ)µ where every parameter is estimated at the observation height r o , even when not specified explicitly. Note that we introduce no new free parameter compared to those used for the leading term, as the parameters a and δu appearing in the PCM are extracted from numerical hydrodynamic simulations of the stellar atmosphere. Plugging Eqs. (B.19) and (B.28) in the expression of the velocity power spectral density given by Eq. (4), we obtain Eq. (13).

Appendix C: A simplified toy model
Here we consider a simplified model of solar acoustic mode excitation, which has already been developed, used, and analysed in previous works (see e.g. Abrams & Kumar 1996;Chaplin & Appourchaux 1999). In the scope of this toy-model, the acoustic potential appearing in Eq. (8) takes the form of a square well, and the sound speed c is considered constant throughout the entire stellar radius. The latter approximation allows us to substitute the radial coordinate r in the wave equation for the acoustic depth τ, defined such that dτ = dr/c. In this approximation, the wave equation simply yields d 2 Ψ dτ 2 + (ω 2 − V(τ) + jωγ)Ψ = δ(τ − τ s ) , (C.1) where we have introduced a point-like source at acoustic depth τ s , and the acoustic potential is if a < τ .

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

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

(C.5)
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 ω 2 i = ω 2 + jωγ and ω 2 o = α 2 − ω 2 − jωγ (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: 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: (C.10) We note that the above expression is only valid if τ > a and τ > τ s . Since the first condition means that we observe the mode of oscillation above the upper turning point (which is always the case in practice) and since the second condition means that the excitation of the mode occurs in layers located deeper than the observation height in the atmosphere (which is also always the case in practice, at least for spectrometric measurements), these are not restrictive conditions.

Appendix 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.