Signature of solar g modes in first-order p-mode frequency shifts

Context. Solar gravity modes (g modes) are buoyancy waves trapped in the solar radiative zone that have been very difficult to detect at the surface. Solar g modes would complement solar pressure modes (p modes) in probing the central regions of the Sun, for example the core rotation rate. Aims. A detection of g modes using changes in the large frequency separation of p modes has recently been reported. However, it is unclear how p and g modes interact. The aim of this study is to evaluate to what extent g modes can perturb the frequencies of p modes. Methods. We computed the first-order perturbation to global p-mode frequencies due to a flow field and perturbations to solar structure caused by a g mode. We focused on long-period g modes and assumed that the g-mode perturbations are constant in time. The surface amplitude of g modes is assumed to be $1$ mm s$^{-1}$, which is close to the observational limit. Results. Gravity modes do perturb p-mode frequencies to first order if the harmonic degree of the g mode is even and if its azimuthal order is zero. The effect is extremely small. For dipole and quadrupole p modes, all frequency shifts are smaller than $0.1$ nHz, or $2\times10^{-8}$ in relative numbers. This is because the relative perturbation to solar structure quantities caused by a g mode of realistic amplitude is of the order of $10^{-6}$ to $10^{-5}$. We find that structural changes dominate over advection. Surprisingly, the interaction of g and p modes takes place to a large part near the surface, where p modes spend most of their propagation times and g modes generate the largest relative changes to solar structure. Conclusions. It appears to be impossible to detect g modes solely through their signature in p-mode frequency shifts. Whether g modes leave a detectable signature in p-mode travel times under a given observational setup remains an open question.


Introduction
Using observations of pressure modes (p modes), helioseismology has been very successful in determining solar interior structure and differential rotation in the convection zone and in most of the radiative zone (see, e.g. Howe 2009;Aerts et al. 2010;Basu 2016). However, a number of outstanding questions about the deep solar interior remain, such as the rotation of the solar core. Solar internal gravity modes (g modes), if detected, would provide complementary information to the p modes. This is because g modes have most of their kinetic energy in the central regions of the Sun. In addition, they have a smaller radial wavelength in this region (e.g. Berthomieu & Provost 1990, which would give a higher spatial resolution than the p modes. The quest for detecting solar g modes has, therefore, attracted considerable attention in the past with several reported detections (see Appourchaux et al. 2010;Appourchaux & Pallé 2013, for reviews). Unfortunately, none of the reported detections have been confirmed so far.
In line with Kennedy et al. (1993) and Duvall (2004), Fossat et al. (2017) proposed a new method for detecting g modes in the first-arrival travel time of low-degree p modes using diskintegrated observations. Under the simplest interpretation, this p-mode travel time is inversely proportional to the large frequency separation of p modes, which is the frequency difference between p modes with identical harmonic degrees whose radial orders differ by one. The p-mode travel time may be affected by a g-mode oscillation through changes in sound speed, for example. Fossat et al. (2017) measured the p-mode travel time using eighthour intervals of Global Oscillations at Low Frequency (GOLF, Gabriel et al. 1995) data and a sampling rate of four hours. The idea is that such a time series can sample changes in solar structure and dynamics due to high-order g modes that have periods longer than eight hours. Using this method, Fossat et al. (2017) and Fossat & Schmider (2018) reported periodic signals in the power spectrum of the data that were interpreted in terms of the rotational splittings of the g modes, implying a solar core rotating nearly four times faster than the radiative zone. Subsequently, Schunker et al. (2018) reproduced the signal of Fossat et al. (2017), but warned that it is very sensitive to a number of choices and parameters in the analysis procedure. Very recently, Scherrer & Gough (2019) also showed that the statistical significance of the signal depends on the helioseismic instrument used (from SOHO, SDO, or from the ground).
The motivation for the present study is to understand to what extent g modes can affect p-mode frequencies and whether and Article number, page 1 of 16 arXiv:1907.02379v2 [astro-ph.SR] 1 Aug 2019 A&A proofs: manuscript no. article how such an effect would be measurable. We do not address the question of the signature of g modes in the observed p-mode travel times, which will require a separate study in order to include the details of the measurement procedure.
One may first ask the general question whether p and g modes should interact at all. In a linear wave equation, the superposition principle holds. It ensures that any two wave fields can be added and still fulfil the wave equation. As a consequence, the presence of one wave does not change the properties of a second wave at all. The situation is different, however, when nonlinear effects are taken into account. Non-linearities in the wave equation include non-linearities in the Navier-Stokes equations. For example, this is the case for advection, temporal changes in the background (sound speed, density, pressure), dissipation, and non-adiabatic effects (e.g. Chapters 2 and 3 in Hamilton & Blackstock 1998). In a non-linear wave equation, two solutions to the wave equation do not, in general, add up to a third solution because the waves interact. This effect may in principle be exploited. A successful example of detecting the presence of waves in the Sun using other waves is given by the recent detection of solar equatorial Rossby waves, which were detected because they advect p modes (Löptien et al. 2018;Liang et al. 2019;Hanasoge & Mandal 2019).
In classical acoustics, the non-linear scattering of sound by sound has been observed experimentally (e.g. Thuras et al. 1935) and modelled theoretically (e.g. Westervelt 1957Westervelt , 1963Hamilton & Blackstock 1998, and references therein). For example, Westervelt (1957) obtained an approximative non-linear acoustic wave equation by inserting a solution to the linear wave equation into the full non-linear Navier-Stokes equations. In addition, the non-linearity introduced by a second order term in the equation of state was taken into account. In the presence of sound waves of two different frequencies, the generation of waves with frequencies equal to the sum and difference of the original frequencies is observed (e.g. Thuras et al. 1935).
The non-linear coupling of stellar oscillation eigenmodes has been studied by several authors. Perdang & Blacher (1982) used a Hamiltonian formalism and showed that a non-linear coupling can lead to chaotic motions with an approximate periodicity. Dziembowski (1982) studied non-linear interactions of two and three modes and found that two-mode interaction has an amplitude-limiting effect on stellar pulsations (see also Kumar & Goldreich 1989). Based on Dziembowski (1982), Wentzel (1987) showed that two p modes can generate a g mode by a non-linear interaction. Wentzel (1987) argued that the evanescent behaviour of the g modes in the convection zone is overcompensated by the increase of p-mode sensitivity towards the surface. As a consequence, the coupling would take place predominantly in the convection zone. The non-linear coupling of p and g modes was also studied for tidally interacting neutron star binaries by Weinberg et al. (2013), where the coupling is found to alter the phase of the emitted gravitational wave and thereby its detection signal. Dziembowski (1982) considered modes that are in a resonance condition, that is one frequency is the sum or difference of the other two frequencies. In Weinberg et al. (2013), the mode frequencies are non-resonant, but the radial wave numbers of the p and g modes are similar in a region in the stellar interior, which means that the radial wave numbers fulfil a resonance condition.
A possible signature of g modes in the p-mode frequency signal has been studied previously by Kennedy et al. (1993), Lou (2001), and Scherrer & Gough (2019). In these studies, the g-mode period was assumed to be much longer than the pmode period, and hence the perturbation to solar structure caused by a g mode was assumed to be time-independent on shorter temporal scales. Kennedy et al. (1993) used perturbation theory on the asymptotic p-mode dispersion relation and found that g modes can produce a frequency modulation of the p modes. This would create sidelobes near each p-mode frequency, which was also found by Lou (2001). Using degenerate perturbation theory, Scherrer & Gough (2019) gave an expression for the first-order perturbation of p modes by g modes and estimated the order of magnitude of the signal. They found that the signal reported by Fossat & Schmider (2018) would be due to relative frequency shifts of the order of 10 −5 , which would correspond to a g-mode amplitude of 50 cm s −1 and which is much higher than the observational upper limit on g-mode amplitudes (for a review, see Appourchaux et al. 2010). Both Kennedy et al. (1993) and Scherrer & Gough (2019) found that only g modes with even harmonic degrees can couple to p modes to leading order.
As the non-linearity in the wave equation is certainly a small effect, the interaction of p modes with g modes can be studied using perturbation theory. The simplest approach is to model the effect of a time-independent perturbation on mode frequencies (Kennedy et al. 1993;Scherrer & Gough 2019). In addition, it is possible to consider the change in eigenfunctions (e.g. Lavely & Ritzwoller 1992), changes in mode amplitude ratios (e.g. Schad et al. 2011;, or the coupling of oscillations at different frequencies (e.g. Woodard 2007Woodard , 2016Hanasoge et al. 2017).
In this work, we expand on the approach taken by Kennedy et al. (1993) and Scherrer & Gough (2019), and study the firstorder perturbation of p-mode frequencies by a frozen g mode, that is we assume that the perturbation caused by the g mode is constant in time. Our aim is to obtain quantitative estimates of the induced frequency shifts, and to understand in a more detailed manner where the interaction is taking place and which of the perturbative effects of the g mode is causing the interaction.
The effect of aspherical perturbations on p-mode frequencies has been studied by Gough (1993), on which Kennedy et al. (1993) and Scherrer & Gough (2019) are based. The aforementioned work was done under the assumption that the perturbation to solar structure preserves hydrostatic equilibrium. In this paper, we include the effect of a departure from hydrostatic equilibrium as it is the case for a g mode. In addition, we provide detailed quantitative estimates of the induced first-order frequency shift and we provide a detailed analysis of the resulting selection rules. To do that, we will include a number of effects such as changes to solar structure through changes in sound speed, pressure, density, gravitational potential, and advection due to the flow field from a g mode.
In addition, we need numerical solar p and g mode eigenfunctions which can be computed using stellar oscillation codes like ADIPLS (Christensen-Dalsgaard 2008) or GYRE (Townsend & Teitler 2013). Such a computation is based on a stellar model. In this work, we use model S (Christensen-Dalsgaard et al. 1996) and a result from the MESA stellar evolution code (Paxton et al. 2011). This paper is structured as follows. We first introduce the oscillation equations (Section 2) and the changes in structure and dynamics due to a g mode (Section 3). In Section 4 we setup the perturbation theory and we derive the equations for the firstorder frequency shift of a p mode in the presence of a g mode, including selection rules. We then present numerical results in Section 5 and conclusions in Section 6.

Linearised equations of stellar oscillations
To model the effect of solar g on p modes, we briefly review the equations of stellar oscillations following Aerts et al. (2010, Chapter 3) and Basu (2016). Solar g and p modes are characterised by an eigenfunction ξ(r) which describes the displacement of a parcel of gas during an oscillation, depending on the location r in the solar interior. The displacement eigenfunction is a solution to the linear wave equation where ω is the cyclic eigenfrequency of the mode and the wave operator L can be written Here, ρ, p, and g denote density, pressure, and gravitational acceleration, and u 0 is a background flow, which we assume to be equal to differential rotation (e.g. Schou et al. 1998), Quantities with a zero subscript are values in hydrostatic equilibrium of the Sun. Quantities with a one as subscript denote Eulerian first-order perturbations to hydrostatic equilibrium and their explicit dependence on the displacement is given in Appendix A. It can be shown (Lynden-Bell & Ostriker 1967) that the wave operator L is self-adjoint under the scalar product for a suitable boundary condition, for example a vacuum boundary with ∇ · ξ = 0. As a consequence, two eigenfunctions ξ and ξ with eigenvalues ω 2 ω 2 are orthogonal with regard to the above scalar product. For comparison purposes, we use four different boundary conditions in this work, see Section 3 and Appendix G.
In hydrostatic equilibrium and in absence of a background flow (u 0 = 0), the Sun is assumed to be a spherically symmetric star. In this case, the eigenfunctions can be written in terms of spherical harmonics, where Ψ lm is a vector spherical harmonic (Barrera et al. 1985), the quantitiesr,θ,φ are unit vectors in the directions of r, θ, φ, the functions ξ r , ξ h are radial and horizontal displacement eigenfunctions, and the quantum numbers l, m, n denote harmonic degree, azimuthal order, and radial order, respectively. We adopt the usual convention that the radial order n for a p mode is positive and for a g mode is negative (e.g. Aerts et al. 2010, Chapter 3) and normalise the eigenfunctions so that ξ lmn , ξ l m n = δ ll δ mm δ nn .
The linearity of the wave equation implies that any mode, for example a p mode, is not affected to first order by the restoring force of a different mode, for example a g mode. By first order, we mean here the linearization of the Navier-Stokes equations that yields the linear wave equation. However, as g modes do perturb the hydrostatic equilibrium of the Sun, they may have some affect on p modes at a higher order, if non-linear effects are taken into account.

Perturbation to solar structure and flows due to g modes
In order to determine the effect of a g mode on p-mode frequencies, we need to determine the perturbation of a g mode to solar structure. To do so, we computed eigenfunctions for g and p modes using the codes GYRE (Townsend & Teitler Unno et al. 1989), as well as the atmospheric boundary condition from Dziembowski (1971). The Eulerian perturbations to pressure, density, and gravitational potential for a particular eigenmode are a direct output from the GYRE code while they also can be obtained from the displacement eigenfunctions directly (GYRE and ADIPLS).
As an example, we show eigenfunctions for l g = 2 g modes with n g = −10 and n g = −36 in Figures 1 and 2. The eigenfunctions were computed using the slightly smoothed version of model S, the GYRE code, and the boundary condition of an isothermal atmosphere from Christensen-Dalsgaard (2008). The panels in the left columns of Figures 1 and 2 show that the kinetic energy density of the g modes is concentrated in the radiative zone. The g modes are trapped in this region, the eigenfunctions are oscillatory and the absolute changes to solar structure variables caused by the g mode, i.e. Eulerian perturbations to density, pressure, and squared sound speed c 2 , reach their maximum amplitude. In the convection zone, the amplitude of the kinetic energy density is lower and decays. The mode is evanescent in this region. This is well known (e.g. Christensen-Dalsgaard 2003).
On the other hand, the panels in the right columns of Figures 1 and 2 show that relative changes to density, pressure, and sound speed are very small, namely at the order of 10 −6 − 10 −5 for a g mode normalised to a surface amplitude of 1 mm s −1 . This is close to the upper limit for the g-mode amplitude from observations (e.g. Appourchaux et al. 2010, and references therein). Even more interestingly, the relative changes to solar structure caused by the g mode reach local maxima close to the solar surface and their amplitude is larger than in the core, at least for g modes of lower radial order. For modes with high radial order, the relative changes to solar structure are larger but of similar order of magnitude in the core compared to the surface. This concentration of relative changes to solar structure near the surface is due to the sharp decrease of density and pressure in this region and takes place despite the concentration of kinetic energy near the core and despite the evanescent character of the g modes in the convection zone. To guarantee the correctness of the eigenfunction computation, we have performed a number of tests, see Appendix B. These tests include a comparison of eigenfunctions to asymptotic formulae, see Figures B.4 and B.5.
The question is now, how p modes are sensitive to the perturbation to solar structure from a g mode. It is well known that p modes are trapped between their lower turning point, which mostly depends on their horizontal surface phase speed ωR / √ l(l + 1), and their upper turning point close to the photosphere (e.g. Christensen-Dalsgaard 2003). The kinetic energy of p modes is concentrated in this region and it is increasing to- c 2 1 /c 2 0 l g = 2, n g = 10, g = 103.2 Hz . The mode was normalised to a surface amplitude of 1 mm s −1 , which is at the order of observational upper limits (e.g. Appourchaux et al. 2010). Displayed are (from top to bottom) radial and horizontal eigenfunctions, and the Eulerian perturbations to density, pressure, and squared sound speed. The left column shows that the kinetic energy density and absolute changes to solar structure caused by the mode are centred near the solar core. The right column shows that relative changes to solar structure are very small, and in the case of density, pressure, and sound speed, reach the largest values near the surface. The exact behaviour of the relative Eulerian changes to solar structure at the surface depends on the choice of boundary condition used in the eigenfunction computation.
wards the surface. In addition, p modes spend most of their lifetime close to the surface due to the strong increase of the sound speed in the solar interior. As a consequence, the sensitivity of p modes to relative changes in solar structure increases towards the surface (e.g. Figures 10 and 11 in Basu 2016). At the same time, relative changes to solar structure due to g modes reach considerable amplitudes near the surface as we have illustrated in Figures 1 and 2.
We thus conclude that the coupling of a g mode with a p mode takes place to a large extent near the solar surface, if there is any such coupling present in the Sun. The hypothesis that the increase of p-mode sensitivity in the convection zone would overcompensate the decrease in amplitude of the g mode was already put forward by Wentzel (1987).

P-mode frequency perturbation due to a frozen g mode
In the following, we derive equations for the first-order perturbation to a p-mode eigenfrequency in the presence of a g mode.
To do that, we assume that the period of the g mode is much longer than that of the p mode. Furthermore, we assume that the g-mode period is long compared to the observational time scale.
To simplify the analysis, we therefore assume that the perturbation to solar structure caused by a g mode is constant in time (see also Kennedy et al. 1993;Scherrer & Gough 2019). In this picture, the time-dependence of the g mode may be added after the frequency shift has been calculated. The perturbed wave operator for changes to solar structure has been derived before (see, e.g. Basu 2016, and references therein). However, these derivations were made for a perturbation to solar structure that preserves hydrostatic equilibrium. As the perturbation in our case is due to a g-mode oscillation, this assumption is no longer valid. We therefore rederive the perturbed wave operator for a perturbation to solar structure that does not necessarily conserve hydrostatic equilibrium, see Appendix C.
The perturbed wave operator L + δL acts on a p-mode with eigenfunction ξ p and frequency ω p , where the index p stands for a p mode with p = (l p , m p , n p ) = (l, m, n), n p > 0.
On the other hand, the perturbation is caused by a g mode with eigenfunction ξ g and frequency ω g with index g = (l g , m g , n g ) = (l , m , n ), n g < 0.
Consequently, we can express changes to all structure quantities, namely density, pressure, gravitational acceleration, and squared sound-speed, as well as the flow field associated to the g-mode oscillation in terms of the g-mode eigenfunction or Eulerian per- c 2 1 /c 2 0 l g = 2, n g = 36, g = 31.2 Hz where Γ 1,0 (r) is the first adiabatic exponent of the equilibrium solar model (Christensen-Dalsgaard 2008; Aerts et al. 2010). In Equation (14), we have assumed that the Lagrangian perturbation to the first adiabatic exponent during the g-mode oscillation is zero. The adiabatic exponent is a material property and reflects the chemical composition and the ionization state of the plasma (see, e.g. Kippenhahn & Weigert 1990, chapter 14). During the oscillation, the chemical composition does not change when moving with the plasma. The degree of ionization, however, is a function of the gas pressure and temperature of the plasma. As temperature and pressure do change during the oscillation, the degree of ionization can change as well. This may have a small effect on the adiabatic exponent in the ionization zones of the Sun, which we neglect in this analysis.

Frequency shifts in the rotating Sun
The p-mode eigenfunctionsξ p of the wave equation (1) for the rotating Sun can be expressed in terms of the solutions of the eigenfunctions in the non-rotating case, ξ p , to first order in quasi-degenerate perturbation theory (Lavely & Ritzwoller 1992), according tô where we note that in Equation (17) only terms of the same m as in ξ p enter. The perturbation to the eigenfunction,ξ (1) p , has a much smaller amplitude than the unperturbed eigenfunction if the harmonic degree l p of the unperturbed mode is small (compare to Woodard 1989). Similarly, the perturbed eigenfrequencies can be written, to first order, as (e.g. Ritzwoller & Lavely 1991) We assume that rotation does not couple p modes and g modes since their frequencies differ largely. For g modes, the same expansion as in Equations (16)-(19) is possible, and we note that the perturbation to the eigenfunction,ξ (1) g , has a much smaller amplitude than the eigenfunction of the non-rotating model also for g modes in the case of small harmonic degrees as considered in this paper.
Given the eigenfunctions and mode frequencies for the rotating solar model, we can now evaluate the frequency perturbation of the p mode due to a g mode using non-degenerate perturbation theory. The use of non-degenerate perturbation theory will be justified later because the frequency shift from the Article number, page 5 of 16 A&A proofs: manuscript no. article g modes is much smaller than any frequency difference between the eigenfrequencies of the rotating eigenmodes. In this case, the frequency shift from a g mode is simply given by where δL = δL + δL (1) indicates that the perturbation is caused by a g mode in a rotating Sun,ξ g = ξ g +ξ (1) g . Using this and ξ p ,ξ p = 1 + ξ (1) p ,ξ (1) p (see Equations [7] and [16]), Equation (20) can be written to leading order as Here, all terms involvingξ (1) p or δL (1) are of higher than first order and have therefore been omitted. Hence, we evaluate δ(ω 2 p ) and the frequency shift as 2π δν p = δω p = δ(ω 2 p )/(2ω p ) in the following.
The scalar product in Equation (21) leads to a number of integrals over one g-mode and two p-mode eigenfunctions, which are given in detail in Appendix D for the terms involved in δL. Similar integrals have been found by other authors for the coupling of g and p modes (e.g. Dziembowski 1982;Kennedy et al. 1993;Scherrer & Gough 2019), as well as for flows and aspherical perturbations to solar structure preserving hydrostatic equilibrium that are decomposed in spherical harmonics (e.g. Lavely & Ritzwoller 1992;Gough 1993). We note here that we omitted one term (δL 4 ) for simplicity. This term models the change of the Eulerian perturbation to the gravitational potential during the p mode oscillation, or in other words changes to Φ 1,p , in the presence of the g mode. This term is expected to have a small contribution to the wave equation.

Selection rules
Our derivation of the first-order frequency shift gives rise to selection rules for mode coupling, which are derived in Appendix E. The horizontal integrals over three spherical harmonics appearing in Appendix D can be evaluated using Gaunt's formula (e.g, Dahlen & Tromp 1998, see also Hanasoge et al. 2017;Kiefer et al. 2017;Fournier et al. 2018), see Appendices E and F.
We find that a first-order frequency shift of a p mode by a g mode is possible only if all of the following conditions are satisfied.
I.) The harmonic degree l g of the g mode is even. II.) The azimuthal order m g of the g mode is zero. III.) The harmonic degree l g of the g mode is not larger than twice the harmonic degree of the p mode (l g ≤ 2l p ). (2019) have found previously that only g modes of even harmonic degree produce a p-mode frequency shift at first order, see also Gough (1993) for perturbations to solar structure that preserve hydrostatic equilibrium. Although Scherrer & Gough (2019) do not note a selection rule regarding m g , it can be derived using their framework of degenerate perturbation theory, where both the g mode and rotation where treated as perturbations to the solar model. This can be seen by considering the coupling matrix that leads to the perturbed p-mode eigenfunctions and perturbed pmode frequencies at first order. If m g 0, the diagonal of that matrix is equal to the frequency perturbation from rotation. The coupling due to g modes appears at the off-diagonal matrix elements where m p −m p = m g . In the notation of Scherrer & Gough (2019), this corresponds to µ − µ = m, see Equations (25) -(27) in their work. All other matrix elements are equal to zero. As the off-diagonal matrix elements due to the g mode are much smaller than the diagonal matrix elements due to rotation, the off-diagonal matrix elements only lead to a second order perturbation in the eigenfrequency and the first-order frequency shift from the g mode is zero.

Kennedy et al. (1993) and Scherrer & Gough
Somewhat similarly, the selection rules I-II can be obtained using the framework of quasi-degenerate perturbation theory outlined in Lavely & Ritzwoller (1992, see also Schad 2013. Here, we have modelled the first-order frequency shift in the p modes but we did not model the change in p-mode travel time as measured by Fossat et al. (2017). Given the above selection rules, global p-mode frequency shifts can be ruled out as a cause for the observations of Fossat et al. (2017). However, in our eyes, it is not straightforward to make a direct connection between frequency shifts and the measurement of the round-trip travel time by Fossat et al. (2017). For example, it is in principle possible that perturbations to the p-mode eigenfunctions cause a signature of g modes in this measurement.

Numerical results
In the following, we numerically evaluate the analytical expression for the first-order perturbation to the p mode eigenfrequency that we derived in the preceding section. Here, we present results for g modes with harmonic degree l g = 2 and p modes with l p = 2. Results for l p = 1 are very similar. The g mode eigenfunctions are normalised to a surface amplitude of 1 mm s −1 , which is at the order of observational upper limits for the amplitude of solar g modes (e.g. Appourchaux et al. 2010). The frequency shifts computed in this section can therefore either be seen as approximative upper limits, or they can be rescaled to different g-mode surface amplitudes.
Inspired by the analysis of Fossat et al. (2017), we assume that the observational time scale is eight hours and therefore apply our results to g modes with periods longer than that, which means frequencies ν g < 34.7 µHz. Moreover, we obtain results also for low-order g modes with frequencies above 34.7 µHz for information purposes, even though our assumptions of a clear separation of observational time scale and g-mode period may break down in this range. We do not consider results obtained for n g < −200 (ν g < 5.67 µHz) in the following. We checked however, that this limitation does not effect the major conclusions presented here. In addition to the p-mode range considered by Fossat et al. (2017), we obtain results for p modes with frequencies ν p < 5 mHz. The results presented in this section were obtained from eigenmodes computed with GYRE and solar model S (Christensen-Dalsgaard et al. 1996), and we assumed an isothermal atmosphere as surface boundary condition (Christensen-Dalsgaard 2008), see also Appendix G.

The effect is tiny in magnitude
As a first result, we find that the magnitude of the frequency shift of a p mode caused by a frozen g mode is extremely small, see Figure 3. All frequency shifts δν p computed are found to be below 0.1 nHz. The relative frequency shifts δν p /ν p are all below an order of magnitude of 2 × 10 −8 when converted to the g-mode amplitude assumed in this work.
Generally, we find the frequency shifts to be greatest for large p-mode and large g-mode frequencies. The study by Fossat et al. (2017) considered a range of g mode frequencies, which does not cover the g modes with the largest frequencies, see Figure 3 (right panel and dashed region in left panel). For these modes, the p mode frequency shifts are smaller by about one order of magnitude and do not exceed 0.01 nHz, or relative shifts of 10 −9 .
The observational upper limit on the g-mode amplitude depends on mode frequency, see for instance Figure 15 in Appourchaux et al. (2010). While it is around 1 mm s −1 for the g modes with the lowest radial orders, it is at the order of 1 cm s −1 for modes with moderately higher radial orders, for example ν g ≈ 35 µHz, n g ≈ −33. This is where the analysis of Fossat et al. (2017) is sensitive to. Even if we rescale the frequency shifts for all modes with ν g < 35 µHz to an amplitude of 1 cm s −1 , the maximum relative frequency shift among all modes does not change. This can also be seen by inspecting the relevant region in Figure 3.

Contributions from density, pressure, sound speed, and flow
The question naturally arises, which of the effects of the g mode, the induced changes to density, pressure, sound speed, gravitional acceleration, or its flow field, is the dominant contributor to the frequency shift.
In order to answer this question, we sum the contributions from the individual terms derived in Appendix D for each of these quantities. The results are shown in Figure 4. In addition, we also show the contribution of an example term that usually is of large amplitude (δL 3−2 , see Eq. [D.33]), and which is due to changes in density. It can be seen that this term is usually about one order of magnitude larger than the total density contribution. Similarly, the total contribution from relevant structural changes, that is from density, pressure, and sound-speed, is again about an order of magnitude smaller than their individual contribution. This is because these terms tend to cancel each other to a large extent. The contribution from changes to gravity is negligible.
Even though the contribution of the flow is about two orders of magnitude smaller than from individual density terms, it still has a relevant contribution to the total frequency shift of order 10 %. In some cases, this number is larger, if the total frequency shift is close to zero.
In general, however, the effect of structural changes dominates over the effect of a flow caused by the g mode.

Regions of interaction
We further investigate in which region in the Sun the interaction predominantly occurs. To do that, we rewrite the total integral leading to the frequency shift from Equation (21) as one radial integral. In addition, we decompose this integrand in the contribution from changes to solar structure and from the flow, see Figure 5. We find that the total integrand takes the largest values very close to the surface and that it is nearly identical to the integrand from the structure terms. This is in agreement with our finding from Section 3 that the relative changes to solar structure from a g mode are concentrated near the surface. However, the integrands change sign near the surface, and the total con-tribution from the near-surface layers is thus smaller than the magnitude of the integrands suggest.
To investigate this further, we computed the contributions to the total integral from the convection zone and from the radiative zone, which are shown in Figure 6. It can be seen that the contribution from the radiative zone is about an order of magnitude smaller than from the convection zone (n g = −36 in Fig. 6). A similar conclusion was reached by Wentzel (1987) by the simple argument that the evanescent behaviour of the g modes in the convection zone is overcompensated by the increase in sensitivity of the p modes. However, the contribution from the radiative zone can reach a comparable magnitude as from the convection zone for higher-order modes, for example n g ≤ −100, see Figure 6. Following Christensen-Dalsgaard & Thompson (1997), one may ask whether the above conclusions stay the same if we use Lagrangian instead of Eulerian perturbations to solar structure. These authors found that contributions to the frequency shift from, e.g, Eulerian changes to sound speed and density nearly cancel as it is the case in our study. If one uses Lagrangian perturbations, such a near cancellation does not occur as strongly in their computation.
We therefore recomputed our main results using Lagrangian perturbations to solar structure for a selected number of mode pairs as test cases. As expected, the resulting total frequency shifts agree with the ones from the Eulerian perturbations (successful for low radial order g-modes |n g | 40).
We further considered the example of quadrupolar modes with n g = −36 and n p = 20 in more detail. For this pair, the radial integrand from the Lagrangian variables has the same order of magnitude as the Eulerian one and a qualitatively similar behaviour. Both assume maximal values in the same region close to the surface. Their actual functional dependence on depth is different, which is consistent with Christensen- Dalsgaard & Thompson (1997). The contributions from the radiative zone and the convection zone are of comparable amplitude in the Lagrangian case, but of opposite sign. These contributions have a larger amplitude as in the Eulerian case, which leads to additional near-cancellation. This is found for all pairs of modes we considered. Thus, Lagrangian perturbations do not seem to us as a better choice for studying this problem.
For both Eulerian and Lagrangian variables, we find in general that the contribution from near-surface layers is of similar magnitude as the contribution from the convection zone. We thus conclude that the interaction of g and p modes is taking place to a large part near the surface, whether solar structure perturbations are computed in Eulerian or Lagrangian variables.

Dependence on solar model and surface boundary conditions
Due to the large magnitude of the integrands near the surface, one may of course ask whether and how the numerical results depend on the surface boundary condition applied in the eigenfunction computation. We therefore computed eigenfunctions using a number of different boundary conditions, which are implemented in the GYRE code and which are detailed in Appendix G.
In addition, we checked whether our results depend on the solar model used to compute the eigenfunctions. The second derivative of the equilibrium density enters in our results, although it is not quite smooth in the original version of model S. This leads to spikes in some of the integrands near the core and wiggles near the surface, see for example the right panel in Figure 5. We therefore also computed eigenfunctions using a version of model S, where the density was slightly smoothed, see details in Appendix G. In order to get a rough idea about the magnitude of possible systematic errors, we also obtained results for a 1 M stellar model obtained using MESA (Paxton et al. 2011). We find that the numerical results depend to some extent on the boundary condition and solar model used. Qualitative results, however, are unchanged, like the overall magnitude of the effect, selection rules, the interaction taking place to a large part near the surface, that the effects of density, pressure and sound speed nearly cancel, and that the changes to solar structure have a larger effect than the flow.
For all boundary conditions and models, the general pattern and the order of magnitude is the same as in Figure 3. For some boundary conditions, models and mode combinations, however, parts of the pattern visible in the right panel in Figure 3 are shifted by about 0.5 − 1.0 × 10 −12 Hz.
There are a number of good reasons why the observed sensitivity to solar model and surface boundary conditions is not surprising. Firstly, we found that the interaction of p and g modes takes place to a large part near the surface. This layer has thus to be accurately modelled to guarantee an accurate answer. Secondly, p modes are sensitive to relative perturbations to solar structure. Close to the surface, the relative Eulerian perturbations from a g mode depend on the surface boundary conditions used in the eigenfunction computation. Thirdly, we found that the effects of pressure, density, and sound speed almost cancel, with the final frequency shift being nearly two orders of magnitude smaller than the largest term involved in the sum. Substracting terms of similar magnitude may result in an amplification of systematic errors.
As a consequence, all effects have to be modelled in a very accurate way if quantitative predictions from this method are to be used.
Finally, it is unknown whether convection has effects on g modes that are not included in the eigenfunction computation such as by turbulence. For example, if the amplitude of the g-mode eigenfunction was significantly reduced at the surface compared to the core, the contribution of the radiative zone to the frequency shifts would be enhanced. Attenuation of p and g modes has not been included in our model.

Conclusion
In this paper, we have developed a model for the first-order frequency shift of solar p modes caused by solar g modes in the framework of non-degenerate time-independent perturbation theory. To do so, the time-dependent perturbation to solar structure caused by a g mode is assumed to be constant at the observational time scale and the degeneracy in mode frequency is assumed to be lifted by rotation.
Our study was motivated by a recent study by Fossat et al. (2017), who proposed to detect g modes by measuring their impact on the wave travel time across the Sun, which is inversely proportional to the large frequency separation of p-mode frequencies.
Firstly, we find that indeed g modes can perturb global pmode eigenfrequencies. However, the effect is extremely small. For a g-mode surface amplitude of A g = 1 mm s −1 , all frequency shifts computed with our model are smaller than 0.  5. Radial integrands for computing the frequency shift of a p mode with l p = 2, m p = 0, n p = 20 caused by a g mode with l g = 2, m g = 0, n g = −36. When integrating over radius, the total integrand yields the total frequency shift δν p (both panels). The contributions to the total integrand resulting from the structural changes and the flow are also shown (right panel). The spike at the surface boundary is due to numerical difficulties in computing second order derivatives at the boundary. Therefore, the three outermost radial grid points were discarded. relative numbers where results for p modes with l p = 1, 2 and g modes with l g = 2 are taken into account. Measuring such a small effect would require an extremely long observation time. Our results are in broad agreement with the estimates obtained by Scherrer & Gough (2019).
Secondly, we find that the interaction of g and p modes takes place to a large part in near-surface layers. This is due to two reasons. On the one hand, p modes are most sensitive to relative changes to solar structure close to the surface. On the other hand, relative changes caused by g modes are relatively large in this region, too. This is the case despite the large amplitude of g modes in the solar core, and caused by the sharp decrease of pressure and density near the surface. Furthermore, the effect from structural changes dominates over the effect from the flow caused by the g mode.
Thirdly, we find that a first-order frequency shift of a p mode by a g mode is only possible if all of the following selection rules are satisfied. The harmonic degree l g of the g mode has to be even, the azimuthal order m g of the g mode has to be zero, and the harmonic degree l g of the g mode has to be smaller than or equal to twice the harmonic degree of the p mode (l g ≤ 2l p ). The condition that the harmonic degree of the g mode has to be even has been previously obtained by Kennedy et al. (1993) and Scherrer & Gough (2019).
We caution that we only modelled the first-order frequency shift in the p modes but not the change in p-mode travel time as measured by Fossat et al. (2017). For example, we did not model the effect of g-modes on the p-mode eigenfunctions, which may result in different selection rules.
While challenging, the exploration of various methods for detecting g modes certainly remains a worthwhile endeavour.
As a consequence, the wave operator L in Eq.
(2) can be written as

Appendix B: Checking the eigenfunction computation
To check the correctness of the eigenfunction computations, we have performed a number of tests. Firstly, we reproduced plots of eigenfunctions in the literature, namely Berthomieu & Provost (1990, Figs. 6, 7a, 7b), Provost et al. (2000, Fig . 1 Secondly, we compared our results for g-mode eigenfunctions to asymptotic relations, both for the trapping region (Eqs. 7.129 and 7.131 in Christensen-Dalsgaard 2003) and the evanescent region, which mostly coincides with the convection zone (Eqs. 7.137 and 7.142 in Christensen-Dalsgaard 2003). We find that the asymptotic relations qualitatively well reproduce the displacement eigenfunctions, and in most cases, a reasonable approximation to the displacement eigenfunction can be obtained by glueing together eigenfunctions obtained from the asymptotic relations for the trapping and evanescent regions near the bottom of the convection zone, see Figure B.4. From the asymptotic displacement eigenfunctions, asymptotic expressions for the structure eigenfunctions can be obtained using Eqs. (A.1), (A.2), and (D.7), which agree well with the full structure eigenfunctions near the surface, see Figure B.5.
Thirdly, we checked for modes of low order that eigenfunctions obtained using GYRE and ADIPLS agree closely provided the same solar model and the same boundary conditions are used.

Appendix C: The perturbed wave operator without hydrostatic equilibrium
In the following, we compute the perturbation to the wave operator δL in the presence of perturbations to solar structure, δρ 0 , δp 0 , δ(c 2 0 ), δg 0 , which potentially do not fulfil the hydrostatic equilibrium condition δ∇p 0 = δ(ρ 0 g 0 ). Omitting the zero Kinetic energy density (ρr 2 (ξ r (r) 2 + l(l + 1)ξ h (r) 2 ) [g cm]) as a function of fractional solar radius r/R for modes with degree l p = l g = 1. The radial orders and frequencies of the modes are given in each panel. Compare to Figure 1 from Provost et al. (2000). We note that the mode frequencies are slightly different compared to Provost et al. (2000) due to slight differences in the solar model.

Appendix D: Perturbation to p-mode frequencies
We will now assume that the perturbation to solar structure is given by a g mode as outlined in Eqs. (10)-(15). Given these perturbations, we will determine the perturbation to p-mode frequencies using Eq. (21) and Eqs. (C.1)-(C.5).

Appendix D.1: Some useful relations
To ease the necessary derivations, we assume that in addition to the eigenfunctions ξ r,g , ξ h,g , the Eulerian perturbations p 1,g , ρ 1,g , Φ 1,g are given and we separate their radial and horizontal dependence (e.g. Aerts et al. 2010), In addition, we need and an expression for the divergence of the displacement eigenfunction, ∇ · ξ p = ∇ · ξ r,p Y lmr + ∇ · ξ h,p Ψ lm (D.10) = dξ r,p dr + 2 r ξ r,p − l(l + 1) r ξ h,p Y lm (D.11) = Div ξ p (r)Y lm , (D.12) where the function Div ξ p (r) is defined as the term in the large brackets in Eq. (D.11).

Appendix D.2: Solar structure terms
We first separately look at the perturbations to solar structure, which means contributions to δ(ω 2 p ) arising from δL 1 till δL 3 . We omit one term (δL 4 ) for simplicity. This term is expected to have a small contribution to the wave equation. It models the change in Eulerian perturbation to the gravitational potential Φ 1,p due to the presence of the g mode. To simplify our derivations, we write δL i− j for the jth term in δL i in Eqs. (C.1)-(C.3) and we use the relations given in Appendix D.1.