© V. G. A. Böning et al. 2019

1. 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, 1991), 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 et al. (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 disk-integrated 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 eight-hour 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 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 non-linear 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. Chaps. 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 1957, 1963; Hamilton & 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 p-mode 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⁻⁵, which would correspond to a g-mode amplitude of 50 cm s⁻¹ 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, 2013), or the coupling of oscillations at different frequencies (e.g. Woodard 2007, 2016; Hanasoge et al. 2017).

In this work, we expand on the approach taken by Kennedy et al. (1993) and Scherrer & Gough (2019), and study the first-order 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 (Sect. 2) and the changes in structure and dynamics due to a g mode (Sect. 3). In Sect. 4 we setup the perturbation theory and we derive the equations for the first-order frequency shift of a p mode in the presence of a g mode, including selection rules. We then present numerical results in Sect. 5 and conclusions in Sect. 6.

2. 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, Chap. 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

$$\begin{array}{c}\hfill {\omega }^2\mathit{\xi }=\mathcal{L}[\mathit{\xi }],\end{array}$$(1)

where ω is the cyclic eigenfrequency of the mode and the wave operator ℒ can be written

$$\begin{array}{c}\hfill \mathcal{L}[\mathit{\xi }]=\frac{1}{{\rho }_0}(\mathrm{\nabla }{p}_1-{\rho }_1{\mathit{g}}_0-{\rho }_0{\mathit{g}}_1)-2i\omega ({\mathit{u}}_0·\mathrm{\nabla })\mathit{\xi }.\end{array}$$(2)

Here, ρ, p, and g denote density, pressure, and gravitational acceleration, and u₀ is a background flow, which we assume to be equal to differential rotation (e.g. Schou et al. 1998),

$$\begin{array}{c}\hfill {\mathit{u}}_0=rsin\theta \mathrm{\Omega }(r,\theta )\widehat{\mathit{\varphi }}.\end{array}$$(3)

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 ℒ is self-adjoint under the scalar product

$$\begin{array}{c}\hfill ⟨\mathit{\xi },{\mathit{\xi }}^{\prime }⟩={{\displaystyle \int \odot }}_{}{\mathit{\xi }}^{\ast }·{\mathit{\xi }}^{\prime }\rho {\mathrm{d}}^3\mathit{r}\end{array}$$(4)

for a suitable boundary condition, for example a vacuum boundary with ∇ ⋅ ξ = 0. As a consequence, two eigenfunctions ξ and ξ′ with eigenvalues ω² ≠ ω′² are orthogonal with regard to the above scalar product. For comparison purposes, we use four different boundary conditions in this work, see Sect. 3 and Appendix G.

In hydrostatic equilibrium and in absence of a background flow (u₀ = 0), the Sun is assumed to be a spherically symmetric star. In this case, the eigenfunctions can be written in terms of spherical harmonics,

$$\begin{array}{cc}\hfill {\mathit{\xi }}_{\mathit{lmn}}(r,\theta ,\varphi )& ={\xi }_{r,ln}(r)Y_{\mathit{lm}}(\theta ,\varphi )\widehat{\mathit{r}}+{\xi }_{h,ln}(r)(\frac{\partial Y_{\mathit{lm}}}{\partial \theta }\widehat{\mathit{\theta }}+\frac{1}{sin\theta }\frac{\partial Y_{\mathit{lm}}}{\partial \varphi }\widehat{\mathit{\varphi }})\hfill \\ \hfill & ={\xi }_{r,ln}(r)Y_{\mathit{lm}}(\theta ,\varphi )\widehat{\mathit{r}}+{\xi }_{h,ln}{\mathbf{\Psi }}_{\mathit{lm}},\hfill \end{array}$$(5)

where Ψ_lm is a vector spherical harmonic (Barrera et al. 1985), the quantities $\widehat{\mathit{r}},\widehat{\mathit{\theta }},\widehat{\mathit{\varphi }}$ 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, Chap. 3) and normalise the eigenfunctions so that

$$\begin{array}{c}\hfill ⟨{\mathit{\xi }}_{\mathit{lmn}},{\mathit{\xi }}_{l^{\prime }{m}^{\prime }{n}^{\prime }}⟩={\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}{\delta }_{n{n}^{\prime }}.\end{array}$$(6)

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.

3. 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 2013; Townsend et al. 2018) and ADIPLS (Christensen-Dalsgaard 2008). The computations are based on three non-rotating solar equilibrium models, namely solar model S (Christensen-Dalsgaard et al. 1996), a version of model S where the density was smoothed slightly, and a solar model obtained from MESA (Paxton et al. 2011), see Appendix G. Furthermore, we tested the effect of several surface boundary conditions like zero Lagrangian pressure perturbation (vacuum boundary) at the surface and two conditions for an isothermal atmosphere (Christensen-Dalsgaard 2008; 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 Figs. 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 Figs. 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², 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).

Fig. 1. Eigenfunctions and Eulerian perturbations to solar structure for a typical g mode with lg = 2, ng = −10, νg = 103.2 μHz (compare to Fig. 5.10 in Christensen-Dalsgaard 2003). 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. Left column: kinetic energy density and absolute changes to solar structure caused by the mode are centred near the solar core. Right column: 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.

Fig. 2. Same as Fig. 1, but for a g mode with ng = −36, νg = 31.2 μHz. See also Figs. B.4 and B.5.

On the other hand, the panels in the right columns of Figs. 1 and 2 show that relative changes to density, pressure, and sound speed are very small, namely at the order of 10⁻⁶ − 10⁻⁵ for a g mode normalised to a surface amplitude of 1 mm s⁻¹. 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 Figs. 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 $\omega R_{\odot }/\sqrt{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 towards the surface. In addition, p modes spend most of their life-time 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. Figs. 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 Figs. 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).

4. 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 ℒ + δℒ acts on a p-mode with eigenfunction ξ_p and frequency ω_p, where the index p stands for a p mode with

$$\begin{array}{c}\hfill p=(l_p,m_p,n_p)=(l,m,n),n_p>0.\end{array}$$(7)

On the other hand, the perturbation is caused by a g mode with eigenfunction ξ_g and frequency ω_g with index

$$\begin{array}{c}\hfill g=(l_g,m_g,n_g)=(l^{\prime },m^{\prime },n^{\prime }),n_g<0.\end{array}$$(8)

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 perturbations (see also Eqs. (A.1)–(A.3)),

$$\begin{array}{cc}& \delta p_0=p_{1,g},\hfill \end{array}$$(9)

$$\begin{array}{cc}& \delta {\rho }_0={\rho }_{1,g},\hfill \end{array}$$(10)

$$\begin{array}{cc}& \delta {\mathit{g}}_0={\mathit{g}}_{1,g},\hfill \end{array}$$(11)

$$\begin{array}{cc}& \delta (c_0^2)=-\frac{{\mathrm{\Gamma }}_{1,0}{p}_0}{{\rho }_0^2}\delta {\rho }_0+\frac{{\mathrm{\Gamma }}_{1,0}}{{\rho }_0}\delta p_0+\frac{p_0}{{\rho }_0}\delta {\mathrm{\Gamma }}_{1,0},\hfill \end{array}$$(12)

$$\begin{array}{cc}& \delta {\mathrm{\Gamma }}_{1,0}=-\mathrm{\nabla }{\mathrm{\Gamma }}_{1,0}·{\mathit{\xi }}_g,\hfill \end{array}$$(13)

$$\begin{array}{cc}& \delta {\mathit{u}}_0=-i\omega {\mathit{\xi }}_g,\hfill \end{array}$$(14)

where Γ_1, 0(r) is the first adiabatic exponent of the equilibrium solar model (Christensen-Dalsgaard 2008; Aerts et al. 2010). In Eq. (13), 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, Chap. 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.

4.1. Frequency shifts in the rotating Sun

The p-mode eigenfunctions ${\widehat{\mathit{\xi }}}_p$ of the wave equation (Eq. (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 to

$$\begin{array}{cc}& {\widehat{\mathit{\xi }}}_p={\mathit{\xi }}_p+{\widehat{\mathit{\xi }}}_p^{(1)},\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\widehat{\mathit{\xi }}}_p^{(1)}={\displaystyle \sum _{\begin{array}{c}{p}^{\prime }\ne p\\ m_{p^{{}^{\prime }}}=m_p\end{array}}}{c}_{p,p^{\prime }}{\mathit{\xi }}_{p^{\prime }},\hfill \end{array}$$(16)

$$\begin{array}{cc}& c_{p,p^{\prime }}=\frac{-2i{\omega }_p⟨{\mathit{\xi }}_{p^{{}^{\prime }}},({\mathit{u}}_0·\mathrm{\nabla }){\mathit{\xi }}_p⟩}{{\omega }_p^2-{\omega }_{p^{\prime }}^2},\hfill \end{array}$$(17)

where we note that in Eq. (16) only terms of the same m as in ξ_p enter. The perturbation to the eigenfunction, ${\widehat{\mathit{\xi }}}_p^{(1)}$, 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)

$$\begin{array}{c}\hfill {\widehat{\omega }}_p={\omega }_p+{\widehat{\omega }}_p^{(1)}.\end{array}$$(18)

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 Eqs. (15)–(18) is possible, and we note that the perturbation to the eigenfunction, ${\widehat{\mathit{\xi }}}_g^{(1)}$, 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 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

$$\begin{array}{c}\hfill \delta ({\widehat{\omega }}_p^2)=\frac{⟨{\widehat{\mathit{\xi }}}_p,\delta \widehat{\mathcal{L}}[{\widehat{\mathit{\xi }}}_p]⟩}{⟨{\widehat{\mathit{\xi }}}_p,{\widehat{\mathit{\xi }}}_p⟩},\end{array}$$(19)

where $\delta \widehat{\mathcal{L}}=\delta \mathcal{L}+\delta {\widehat{\mathcal{L}}}^{(1)}$ indicates that the perturbation is caused by a g mode in a rotating Sun, ${\widehat{\mathit{\xi }}}_g={\mathit{\xi }}_g+{\widehat{\mathit{\xi }}}_g^{(1)}$. Using this and $⟨{\widehat{\mathit{\xi }}}_p,{\widehat{\mathit{\xi }}}_p⟩=1+⟨{\widehat{\mathit{\xi }}}_p^{(1)},{\widehat{\mathit{\xi }}}_p^{(1)}⟩$ (see Eqs. (6) and (15)), Eq. (19) can be written to leading order as

$$\begin{array}{c}\hfill \delta ({\widehat{\omega }}_p^2)=⟨{\mathit{\xi }}_p,\delta \mathcal{L}[{\mathit{\xi }}_p]⟩=\delta ({\omega }_p^2).\end{array}$$(20)

Here, all terms involving ${\widehat{\mathit{\xi }}}_p^{(1)}$ or $\delta {\widehat{\mathcal{L}}}^{(1)}$ are of higher than first order and have therefore been omitted. Hence, we evaluate δ(ω_p²) and the frequency shift as 2π δν_p = δω_p = δ(ω_p²)/(2ω_p) in the following.

The scalar product in Eq. (20) 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 δℒ. 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 (δℒ₄) 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.

4.2. 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).

Kennedy et al. (1993) and Scherrer & Gough (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 p-mode 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 Eqs. (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.

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.

5. 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⁻¹, 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.

5.1. 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 Fig. 3. All frequency shifts δν_p computed are found to be below 0.1 nHz. The relative frequency shifts δν_p/ν_p are all below 2 × 10⁻⁸. These numbers are in very broad agreement with the estimates presented in Scherrer & Gough (2019), who reported an order of magnitude of 2 × 10⁻⁸ when converted to the g-mode amplitude assumed in this work.

Fig. 3. First-order perturbations to the p-mode frequency due to a g-mode oscillation assumed to be frozen over the observation period. Results are shown for g modes and p modes with harmonic degrees lg = lp = 2 and azimuthal orders mg = mp = 0. Left: frequency shift as a function of radial order for mode frequencies νp <  5 mHz and νg >  5.67 μHz. The range of modes studied by Fossat et al. (2017), Fossat & Schmider (2018), and Schunker et al. (2018) is inside the dashed rectangle, extending outwards of the plotted region at |ng| > 200. Right: cuts through the results shown in the dashed region in the left panel for a selection of g-mode radial orders.

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 Fig. 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⁻⁹.

The observational upper limit on the g-mode amplitude depends on mode frequency, see for instance Fig. 15 in Appourchaux et al. (2010). While it is around 1 mm s⁻¹ for the g modes with the lowest radial orders, it is at the order of 1 cm s⁻¹ 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⁻¹, the maximum relative frequency shift among all modes does not change. This can also be seen by inspecting the relevant region in Fig. 3.

5.2. 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 Fig. 4. In addition, we also show the contribution of an example term that usually is of large amplitude (δℒ₃₋₂, see Eq. (D.21)), 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.

Fig. 4. Contribution to the frequency shift of p modes by changes to solar structure variables ρ, p, c2, g, and the flow field u caused by a g mode in absolute values. The total frequency shift |δνp| and the contribution from an example term (δℒ3−2, see Eq. (D.21)) are also given for comparison. The azimuthal order of the modes is mp = mg = 0.

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.

5.3. 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 Eq. (20) as one radial integral. In addition, we decompose this integrand in the contribution from changes to solar structure and from the flow, see Fig. 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 Sect. 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 contribution from the near-surface layers is thus smaller than the magnitude of the integrands suggest.

Fig. 5. Radial integrands for computing the frequency shift of a p mode with lp = 2, mp = 0, np = 20 caused by a g mode with lg = 2, mg = 0, ng = −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.

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 Fig. 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 Fig. 6.

Fig. 6. Contribution to the frequency shift from the convection zone and from the radiative zone as a function of p-mode frequency for several typical g modes considered in this study. Results are shown for lp = lg = 2, mp = mg = 0.

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.

5.4. 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 Fig. 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 Fig. 3. For some boundary conditions, models and mode combinations, however, parts of the pattern visible in the right panel in Fig. 3 are shifted by about 0.5−1.0 × 10⁻¹² 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.

6. 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 p-mode eigenfrequencies. However, the effect is extremely small. For a g-mode surface amplitude of A_g = 1 mm s⁻¹, all frequency shifts computed with our model are smaller than 0.1 nHz, or in relative numbers

$$\begin{array}{c}\hfill \frac{|\delta {\nu }_p|}{{\nu }_p}\le 2\times {10}^{-8}\frac{A_g}{1\mathrm{mm}{\mathrm{s}}^{-1}},\end{array}$$(21)

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.

Acknowledgments

We thank Jesper Schou, Aaron Birch, and Damien Fournier for scientific discussions, Earl Bellinger and Warrick Ball for help with installing and running GYRE, and Richard Townsend for answering questions on eigenfunction computations via the GYRE forum. We thank the anonymous referee for a number of useful suggestions that improved the manuscript. We used the tomso python package written by Warrick Ball for reading ADIPLS eigenfunctions and solar models. This work was supported in part by the Max Planck Society through a grant on PLATO Science. The computational resources were provided by the German Data Center for SDO through a grant from the German Aerospace Center (DLR). We acknowledge partial support from the European Research Council Synergy Grant WHOLE SUN #810218. HH acknowledges the support of the Erasmus Mundus Joint Masters Programme, AstroMundus. Parts of this work were submitted as a MSc. thesis by HH.

References

Appendix A: Wave operator

The Eulerian perturbations to the solar structure variables are (e.g. Christensen-Dalsgaard 2003; Basu 2016):

$$\begin{array}{cc}\hfill p_1& =-c_0^2{\rho }_0\mathrm{\nabla }·\mathit{\xi }-\mathrm{\nabla }{p}_0·\mathit{\xi },\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill {\rho }_1& =-\mathrm{\nabla }·({\rho }_0\mathit{\xi })\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill {\mathit{g}}_1& =-\mathrm{\nabla }{\mathrm{\Phi }}_1=-G{\mathrm{\nabla }}_{\mathit{r}}{\displaystyle {\int }_V}\frac{{\mathrm{\nabla }}_{{\mathit{r}}^{\prime }}·({\rho }_0({\mathit{r}}^{\prime })\mathit{\xi }({\mathit{r}}^{\prime }))}{|\mathit{r}-{\mathit{r}}^{\prime }|}{\mathrm{d}}^3{\mathit{r}}^{\prime }.\hfill \end{array}$$(A.3)

As a consequence, the wave operator ℒ in Eq. (2) can be written as

$$\begin{array}{cc}& \mathcal{L}[\mathit{\xi }]={\mathcal{L}}_1[\mathit{\xi }]+{\mathcal{L}}_2[\mathit{\xi }]+{\mathcal{L}}_3[\mathit{\xi }]+{\mathcal{L}}_4[\mathit{\xi }]\hfill \end{array}$$(A.4)

$$\begin{array}{cc}& {\mathcal{L}}_1[\mathit{\xi }]=-\frac{1}{{\rho }_0}\mathrm{\nabla }(c_0^2{\rho }_0\mathrm{\nabla }·\mathit{\xi })\hfill \\ \hfill & =-\mathrm{\nabla }(c_0^2\mathrm{\nabla }·\mathit{\xi })-\frac{\mathrm{\nabla }{\rho }_0}{{\rho }_0}{c}_0^2\mathrm{\nabla }·\mathit{\xi }\hfill \end{array}$$(A.5)

$$\begin{array}{cc}& {\mathcal{L}}_2[\mathit{\xi }]=-\frac{1}{{\rho }_0}\mathrm{\nabla }(\mathrm{\nabla }{p}_0·\mathit{\xi })\hfill \end{array}$$(A.6)

$$\begin{array}{cc}& {\mathcal{L}}_3[\mathit{\xi }]=\frac{1}{{\rho }_0}\mathrm{\nabla }·({\rho }_0\mathit{\xi }){\mathit{g}}_0=(\mathrm{\nabla }·\mathit{\xi }){\mathit{g}}_0+\frac{1}{{\rho }_0}(\mathrm{\nabla }{\rho }_0·\mathit{\xi }){\mathit{g}}_0\hfill \end{array}$$(A.7)

$$\begin{array}{cc}& {\mathcal{L}}_4[\mathit{\xi }]=G{\mathrm{\nabla }}_{\mathit{r}}{\displaystyle {\int }_V}\frac{{\mathrm{\nabla }}_{{\mathit{r}}^{\prime }}·({\rho }_0({\mathit{r}}^{\prime })\mathit{\xi }({\mathit{r}}^{\prime }))}{|\mathit{r}-{\mathit{r}}^{\prime }|}{\mathrm{d}}^3{\mathit{r}}^{\prime }.\hfill \end{array}$$(A.8)

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), and Christensen-Dalsgaard (2003, Figs. 5.10 and 7.4) with the results from our eigenfunction computation, see Figs. B.1–B.3. The results look the same as in the literature.

Fig. B.1. Kinetic energy density (ρr2(ξr(r)2 + l(l + 1)ξh(r)2) [g cm]) as a function of fractional solar radius r/R⊙ for modes with degree lp = lg = 1. The radial orders and frequencies of the modes are given in each panel. Compare to Fig. 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.

Fig. B.2. Scaled radial g-mode eigenfunction. Compare to Fig. 7.4 in Christensen-Dalsgaard (2003).

Fig. B.3. Scaled radial g-mode eigenfunction. Compare to panel a in Fig. 5.10 from Christensen-Dalsgaard (2003).

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 Fig. 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 Fig. B.5.

Fig. B.4. Comparison of radial (top) and horizontal (bottom) g-mode eigenfunctions (computed with GYRE) with asymptotic formulae. For the convection zone (CZ) and radiative zone (RZ), different approximations were used, see Appendix B. The same mode as in Fig. 2 is shown.

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.

Fig. B.5. Comparison of g-mode solar structure eigenfunctions (computed with GYRE) with asymptotic formulae for the convection zone (CZ). The same mode as in Fig. 2 is shown. See Appendix B for details.

Appendix C: The perturbed wave operator without hydrostatic equilibrium

In the following, we compute the perturbation to the wave operator δℒ in the presence of perturbations to solar structure, δρ₀, δp₀, δ(c₀²),δg₀, which potentially do not fulfil the hydrostatic equilibrium condition δ∇p₀ = δ(ρ₀g₀). Omitting the zero subscript in perturbed quantities, for example δρ = δρ₀, and perturbing Eqs. (A.4)–(A.8), we obtain

$$\begin{array}{cc}\hfill \delta {\mathcal{L}}_1[\mathit{\xi }]& =-\mathrm{\nabla }(\delta (c^2)\mathrm{\nabla }·\mathit{\xi })-\frac{\mathrm{\nabla }\delta \rho }{{\rho }_0}{c}_0^2\mathrm{\nabla }·\mathit{\xi }\hfill \\ \hfill & +\frac{\mathrm{\nabla }{\rho }_0}{{\rho }_0^2}\delta \rho c_0^2\mathrm{\nabla }·\mathit{\xi }-\frac{\mathrm{\nabla }{\rho }_0}{{\rho }_0}\delta (c^2)\mathrm{\nabla }·\mathit{\xi }\hfill \end{array}$$(C.1)

$$\begin{array}{c}\hfill \delta {\mathcal{L}}_2[\mathit{\xi }]=\frac{\delta \rho }{{\rho }_0^2}\mathrm{\nabla }(\mathrm{\nabla }{p}_0·\mathit{\xi })-\frac{1}{{\rho }_0}\mathrm{\nabla }(\mathrm{\nabla }\delta p·\mathit{\xi })\end{array}$$(C.2)

$$\begin{array}{cc}\hfill \delta {\mathcal{L}}_3[\mathit{\xi }]& =(\mathrm{\nabla }·\mathit{\xi })\delta \mathit{g}-\frac{\delta \rho }{{\rho }_0^2}(\mathrm{\nabla }{\rho }_0·\mathit{\xi }){\mathit{g}}_0\hfill \\ \hfill & +\frac{1}{{\rho }_0}(\mathrm{\nabla }\delta \rho ·\mathit{\xi }){\mathit{g}}_0+\frac{1}{{\rho }_0}(\mathrm{\nabla }{\rho }_0·\mathit{\xi })\delta \mathit{g}\hfill \end{array}$$(C.3)

$$\begin{array}{c}\hfill \delta {\mathcal{L}}_4[\mathit{\xi }]=G{\mathrm{\nabla }}_{\mathit{r}}{\displaystyle {\int }_V}\frac{{\mathrm{\nabla }}_{{\mathit{r}}^{\prime }}·(\delta \rho ({\mathit{r}}^{\prime })\mathit{\xi }({\mathit{r}}^{\prime }))}{|\mathit{r}-{\mathit{r}}^{\prime }|}{\mathrm{d}}^3{\mathit{r}}^{\prime }.\end{array}$$(C.4)

In addition, in the presence of a perturbation to the flow field δu(r), this flow gives rise to a fifth perturbation term (e.g. Christensen-Dalsgaard 2003)

$$\begin{array}{c}\hfill \delta {\mathcal{L}}_5[\mathit{\xi }]=-2i\omega (\delta \mathit{u}·\mathrm{\nabla })\mathit{\xi },\end{array}$$(C.5)

such that $\delta \mathcal{L}={\displaystyle \sum _{i=1}^5}\delta {\mathcal{L}}_i$.

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. (9)–(14). Given these perturbations, we will determine the perturbation to p-mode frequencies using Eq. (20) and Eqs. (C.1)–(C.5).

For brevity in notation, we use

$$\begin{array}{cc}& (l,m,n)=(l_p,m_p,n_p),\hfill \end{array}$$(D.1)

$$\begin{array}{cc}& (l^{\prime },m^{\prime },n^{\prime })=(l_g,m_g,n_g)\hfill \end{array}$$(D.2)

in Appendices D and E.

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),

$$\begin{array}{cc}\hfill p_{1,g}& ={\stackrel{\sim }{p}}_{1,g}(r)Y_{l^{\prime }{m}^{\prime }},\hfill \end{array}$$(D.3)

$$\begin{array}{cc}\hfill {\rho }_{1,g}& ={\stackrel{\sim }{\rho }}_{1,g}(r)Y_{l^{\prime }{m}^{\prime }},\hfill \end{array}$$(D.4)

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{1,g}& ={\stackrel{\sim }{\mathrm{\Phi }}}_{1,g}(r)Y_{l^{\prime }{m}^{\prime }}.\hfill \end{array}$$(D.5)

Similarly, we obtain using Eqs. (9)–(13) and Eqs. (D.3)–(D.4),

$$\begin{array}{cc}& \delta (c^2)={\stackrel{\sim }{c}}_{1,g}^2(r)Y_{l^{\prime }{m}^{\prime }},\hfill \end{array}$$(D.6)

$$\begin{array}{cc}& {\stackrel{\sim }{c}}_{1,g}^2=-\frac{{\mathrm{\Gamma }}_{1,0}{p}_0}{{\rho }_0^2}{\stackrel{\sim }{\rho }}_{1,g}+\frac{{\mathrm{\Gamma }}_{1,0}}{{\rho }_0}{\stackrel{\sim }{p}}_{1,g}-\frac{p_0}{{\rho }_0}\frac{\mathrm{d}{\mathrm{\Gamma }}_{1,0}}{\mathrm{d}r}{\xi }_{r,g}.\hfill \end{array}$$(D.7)

In addition, we need

$$\begin{array}{cc}\hfill {\mathit{g}}_{1,g}& =-\mathrm{\nabla }{\mathrm{\Phi }}_{1,g}=-\mathrm{\nabla }({\stackrel{\sim }{\mathrm{\Phi }}}_{1,g}{Y}_{l^{\prime }{m}^{\prime }})\hfill \\ \hfill & =-\frac{\mathrm{d}{\stackrel{\sim }{\mathrm{\Phi }}}_{1,g}}{\mathrm{d}r}{Y}_{l^{\prime }{m}^{\prime }}\widehat{\mathit{r}}-\frac{{\stackrel{\sim }{\mathrm{\Phi }}}_{1,g}}{r}{\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}\hfill \end{array}$$(D.8)

and an expression for the divergence of the displacement eigenfunction,

$$\begin{array}{cc}\hfill \mathrm{\nabla }·{\mathit{\xi }}_p& =\mathrm{\nabla }·\left({\xi }_{r,p}{Y}_{\mathit{lm}}\widehat{\mathit{r}}\right)+\mathrm{\nabla }·\left({\xi }_{h,p}{\mathbf{\Psi }}_{\mathit{lm}}\right)\hfill \\ \hfill & =(\frac{\mathrm{d}{\xi }_{r,p}}{\mathrm{d}r}+\frac{2}{r}{\xi }_{r,p}-\frac{l(l+1)}{r}{\xi }_{h,p})Y_{\mathit{lm}}\hfill \\ \hfill & ={\mathrm{Div}}_{{\mathit{\xi }}_p}(r)Y_{\mathit{lm}},\hfill \end{array}$$(D.9)

where the function Div_ξp(r) is defined as the term in the large brackets in Eq. (D.9).

D.2. Solar structure terms

We first separately look at the perturbations to solar structure, which means contributions to δ(ω_p²) arising from δℒ₁ till δℒ₃. We omit one term (δℒ₄) 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 δℒ_i − j for the jth term in δℒ_i in Eqs. (C.1)–(C.3) and we use the relations given in Appendix D.1.

We start by the first term in δℒ₁ (Eq. (C.1)),

$$\begin{array}{cc}\hfill \delta {\mathcal{L}}_{1-1}[{\mathit{\xi }}_p]& =-\mathrm{\nabla }(\delta (c^2)\mathrm{\nabla }·{\mathit{\xi }}_p)=-\mathrm{\nabla }\left({\stackrel{\sim }{c}}_{1,g}^2{Y}_{l^{\prime }{m}^{\prime }}{\mathrm{Div}}_{{\mathit{\xi }}_p}{Y}_{\mathit{lm}}\right)\hfill \\ \hfill & =-(\frac{\mathrm{d}{\stackrel{\sim }{c}}_{1,g}^2}{\mathrm{d}r}{Y}_{l^{\prime }{m}^{\prime }}\widehat{\mathit{r}}+\frac{{\stackrel{\sim }{c}}_{1,g}^2}{r}{\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}){\mathrm{Div}}_{{\mathit{\xi }}_p}{Y}_{\mathit{lm}}\hfill \\ \hfill & -{\stackrel{\sim }{c}}_{1,g}^2{Y}_{l^{\prime }{m}^{\prime }}(\frac{\mathrm{d}{\mathrm{Div}}_{{\mathit{\xi }}_p}}{\mathrm{d}r}{Y}_{\mathit{lm}}\widehat{\mathit{r}}+\frac{{\mathrm{Div}}_{{\mathit{\xi }}_p}}{r}{\mathbf{\Psi }}_{\mathit{lm}}).\hfill \end{array}$$(D.10)

The frequency shift arising from this term can then be evaluated using Eq. (20) and Eq. (5) applied to the p-mode eigenfunction,

$$\begin{array}{cc}\hfill ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{1-1}[{\mathit{\xi }}_p]⟩& =-{\displaystyle \int {\xi }_{r,p}}(\frac{\mathrm{d}{\stackrel{\sim }{c}}_{1,g}^2}{\mathrm{d}r}{\mathrm{Div}}_{{\mathit{\xi }}_p}+{\stackrel{\sim }{c}}_{1,g}^2\frac{\mathrm{d}{\mathrm{Div}}_{{\mathit{\xi }}_p}}{\mathrm{d}r}){\rho }_0{r}^2\mathrm{d}r\hfill \\ \hfill & \times {\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & -{\displaystyle \int {\xi }_{h,p}}{\stackrel{\sim }{c}}_{1,g}^2{\mathrm{Div}}_{{\mathit{\xi }}_p}{\rho }_{0}r\mathrm{d}r({\displaystyle \int {\mathbf{\Psi }}_{\mathit{lm}}^{\ast }}·{\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & +{\displaystyle \int {\mathbf{\Psi }}_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }).\hfill \end{array}$$(D.11)

Here, dΩ = sin θ d θ dϕ is the surface element on the unit sphere.

For the following terms, we proceed in a similar way and obtain

$$\begin{array}{cc}& \delta {\mathcal{L}}_{1-2}[{\mathit{\xi }}_p]=-\frac{\mathrm{\nabla }\delta \rho }{{\rho }_0}{c}_0^2\mathrm{\nabla }·{\mathit{\xi }}_p=-\frac{c_0^2}{{\rho }_0}\mathrm{\nabla }({\stackrel{\sim }{\rho }}_{1,g}{Y}_{l^{\prime }{m}^{\prime }})\mathrm{\nabla }·{\mathit{\xi }}_p\hfill \\ \hfill & =-\frac{c_0^2}{{\rho }_0}(\frac{\mathrm{d}{\stackrel{\sim }{\rho }}_{1,g}}{\mathrm{d}r}{Y}_{l^{\prime }{m}^{\prime }}\widehat{\mathit{r}}+{\stackrel{\sim }{\rho }}_{1,g}\frac{1}{r}{\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}){\mathrm{Div}}_{{\mathit{\xi }}_p}{Y}_{\mathit{lm}},\hfill \end{array}$$(D.12)

$$\begin{array}{cc}& ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{1-2}[{\mathit{\xi }}_p]⟩=-{\displaystyle \int {\xi }_{r,p}}{c}_0^2\frac{\mathrm{d}{\stackrel{\sim }{\rho }}_{1,g}}{\mathrm{d}r}{\mathrm{Div}}_{{\mathit{\xi }}_p}{r}^2\mathrm{d}r{\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & -{\displaystyle \int {\xi }_{h,p}}{c}_0^2{\stackrel{\sim }{\rho }}_{1,g}{\mathrm{Div}}_{{\mathit{\xi }}_p}r\mathrm{d}r{\displaystyle \int {\mathbf{\Psi }}_{\mathit{lm}}^{\ast }}·{\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }.\hfill \end{array}$$(D.13)

Using similar intermediate steps, we obtain the results for the following terms,

$$\begin{array}{cc}\hfill ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{1-3}[{\mathit{\xi }}_p]⟩& ={\displaystyle \int {\xi }_{r,p}}\frac{{\stackrel{\sim }{\rho }}_{1,g}}{{\rho }_0}\frac{\mathrm{d}{\rho }_0}{\mathrm{d}r}{c}_0^2{\mathrm{Div}}_{{\mathit{\xi }}_p}{r}^2\mathrm{d}r{\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega },\hfill \end{array}$$(D.14)

$$\begin{array}{cc}\hfill ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{1-4}[{\mathit{\xi }}_p]⟩& =-{\displaystyle \int {\xi }_{r,p}}\frac{\mathrm{d}{\rho }_0}{\mathrm{d}r}{\stackrel{\sim }{c}}_{1,g}^2{\mathrm{Div}}_{{\mathit{\xi }}_p}{r}^2\mathrm{d}r{\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega },\hfill \end{array}$$(D.15)

$$\begin{array}{cc}\hfill ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{2-1}[{\mathit{\xi }}_p]⟩& ={\displaystyle \int {\xi }_{r,p}}\frac{1}{{\rho }_0}{\stackrel{\sim }{\rho }}_{1,g}(\frac{{\mathrm{d}}^2{p}_0}{\mathrm{d}{r}^2}{\xi }_{r,p}+\frac{\mathrm{d}{p}_0}{\mathrm{d}r}\frac{\mathrm{d}{\xi }_{r,p}}{\mathrm{d}r})r^2\mathrm{d}r\hfill \\ \hfill & \times {\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & +{\displaystyle \int {\xi }_{h,p}}\frac{1}{{\rho }_0}{\stackrel{\sim }{\rho }}_{1,g}\frac{\mathrm{d}{p}_0}{\mathrm{d}r}{\xi }_{r,p}r\mathrm{d}r\hfill \\ \hfill & \times {\displaystyle \int {\mathbf{\Psi }}_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }.\hfill \end{array}$$(D.16)

The next term requires a bit more work. We start with

$$\begin{array}{cc}\hfill \delta {\mathcal{L}}_{2-2}[{\mathit{\xi }}_p]& =-\frac{1}{{\rho }_0}\mathrm{\nabla }(\frac{\mathrm{d}{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}r}{Y}_{l^{\prime }{m}^{\prime }}{\xi }_{r,p}{Y}_{\mathit{lm}}+{\stackrel{\sim }{p}}_{1,g}\frac{1}{r}{\xi }_{h,p}{\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}}),\hfill \\ \hfill & =-\frac{1}{{\rho }_0}(\mathrm{\nabla }(\frac{\mathrm{d}{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}r}{\xi }_{r,p})Y_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}+\frac{\mathrm{d}{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}r}{\xi }_{r,p}\mathrm{\nabla }(Y_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}})\hfill \\ \hfill & +\mathrm{\nabla }({\stackrel{\sim }{p}}_{1,g}\frac{1}{r}{\xi }_{h,p}){\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}}\hfill \\ \hfill & +{\stackrel{\sim }{p}}_{1,g}\frac{1}{r}{\xi }_{h,p}\mathrm{\nabla }({\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}})),\hfill \\ \hfill & =-\frac{1}{{\rho }_0}((\frac{{\mathrm{d}}^2{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}{r}^2}{\xi }_{r,p}+\frac{\mathrm{d}{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}r}\frac{\mathrm{d}{\xi }_{r,p}}{\mathrm{d}r})Y_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\widehat{\mathit{r}}\hfill \\ \hfill & +\frac{\mathrm{d}{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}r}{\xi }_{r,p}\frac{1}{r}({\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}+Y_{l^{\prime }{m}^{\prime }}{\mathbf{\Psi }}_{\mathit{lm}})\hfill \\ \hfill & +(\frac{\mathrm{d}{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}r}\frac{1}{r}{\xi }_{h,p}-{\stackrel{\sim }{p}}_{1,g}\frac{1}{r^2}{\xi }_{h,p}+{\stackrel{\sim }{p}}_{1,g}\frac{1}{r}\frac{\mathrm{d}{\xi }_{h,p}}{\mathrm{d}r})\hfill \\ \hfill & \times ({\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}})\widehat{\mathit{r}}\hfill \\ \hfill & +{\stackrel{\sim }{p}}_{1,g}\frac{1}{r}{\xi }_{h,p}\mathrm{\nabla }({\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}})),\hfill \end{array}$$(D.17)

To evaluate the last term in ⟨ξ_p, δℒ₂₋₂[ξ_p]⟩ arising from Eq. (D.17), we first note that ∇(Ψ_l′m′ ⋅ Ψ_lm) is purely horizontal. The horizontal gradient ${\mathrm{\nabla }}_h=\widehat{\mathit{\theta }}\frac{\partial }{\partial \theta }+\widehat{\mathit{\varphi }}\frac{1}{sin\theta }\frac{\partial }{\partial \varphi }$ applied to a function f(θ, ϕ) that is independent of radius, satisfies r∇f = ∇_hf. It can then be shown using Gauss’s divergence theorem and the property $\mathrm{\nabla }·(g(r){\mathbf{\Psi }}_{\mathit{lm}})=-\frac{l(l+1)}{r}g(r)Y_{\mathit{lm}}$ that

$$\begin{array}{c}\hfill {\displaystyle \int {\mathbf{\Psi }}_{\mathit{lm}}^{\ast }}·{\mathrm{\nabla }}_h({\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}})\mathrm{d}\mathrm{\Omega }=l(l+1){\displaystyle \int Y_{\mathit{lm}}^{\ast }}({\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}})\mathrm{d}\mathrm{\Omega }.\end{array}$$(D.18)

Using this, we can now evaluate

$$\begin{array}{cc}\hfill ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{2-2}[{\mathit{\xi }}_p]⟩& =-{\displaystyle \int r^2}{\xi }_{r,p}(\frac{{\mathrm{d}}^2{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}{r}^2}{\xi }_{r,p}+\frac{\mathrm{d}{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}r}\frac{\mathrm{d}{\xi }_{r,p}}{\mathrm{d}r})\mathrm{d}r\hfill \\ \hfill & \times {\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & -{\displaystyle \int r^2}{\xi }_{h,p}\frac{\mathrm{d}{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}r}{\xi }_{r,p}\frac{1}{r}\mathrm{d}r\hfill \\ \hfill & \times {\displaystyle \int {\mathbf{\Psi }}_{\mathit{lm}}^{\ast }}·({\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}+Y_{l^{\prime }{m}^{\prime }}{\mathbf{\Psi }}_{\mathit{lm}})\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & -{\displaystyle \int (\frac{\mathrm{d}{\stackrel{\sim }{p}}_{1,g}}{\mathrm{d}r}\frac{1}{r}{\xi }_{h,p}-{\stackrel{\sim }{p}}_{1,g}\frac{1}{r^2}{\xi }_{h,p}+{\stackrel{\sim }{p}}_{1,g}\frac{1}{r}\frac{\mathrm{d}{\xi }_{h,p}}{\mathrm{d}r})}\hfill \\ \hfill & \times r^2{\xi }_{r,p}\mathrm{d}r{\displaystyle \int Y_{\mathit{lm}}^{\ast }}({\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}})\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & -{\displaystyle \int r^2}{\xi }_{h,p}{\stackrel{\sim }{p}}_{1,g}\frac{1}{r^2}{\xi }_{h,p}\mathrm{d}r\hfill \\ \hfill & \times l(l+1){\displaystyle \int Y_{\mathit{lm}}^{\ast }}({\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}})\mathrm{d}\mathrm{\Omega }.\hfill \end{array}$$(D.19)

The remaining terms can again be evaluated in a similar way as in Eqs. (D.10)–(D.11),

$$\begin{array}{cc}\hfill ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{3-1}[{\mathit{\xi }}_p]⟩& =-{\displaystyle \int {\xi }_{r,p}}\frac{\mathrm{d}{\stackrel{\sim }{\mathrm{\Phi }}}_{1,g}}{\mathrm{d}r}{\mathrm{Div}}_{{\mathit{\xi }}_p}{\rho }_0{r}^2\mathrm{d}r{\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & -{\displaystyle \int {\xi }_{h,p}}{\stackrel{\sim }{\mathrm{\Phi }}}_{1,g}{\mathrm{Div}}_{{\mathit{\xi }}_p}{\rho }_{0}r\mathrm{d}r\hfill \\ \hfill & \times {\displaystyle \int {\mathbf{\Psi }}_{\mathit{lm}}^{\ast }}·{\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega },\hfill \end{array}$$(D.20)

$$\begin{array}{cc}\hfill ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{3-2}[{\mathit{\xi }}_p]⟩& =-{\displaystyle \int {({\xi }_{r,p})}^2}\frac{1}{{\rho }_0^2}{\stackrel{\sim }{\rho }}_{1,g}\frac{\mathrm{d}{\rho }_0}{\mathrm{d}r}\frac{\mathrm{d}{p}_0}{\mathrm{d}r}{r}^2\mathrm{d}r\hfill \\ \hfill & \times {\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega },\hfill \end{array}$$(D.21)

$$\begin{array}{cc}\hfill ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{3-3}[{\mathit{\xi }}_p]⟩& ={\displaystyle \int {({\xi }_{r,p})}^2}\frac{\mathrm{d}{\stackrel{\sim }{\rho }}_{1,g}}{\mathrm{d}r}\frac{1}{{\rho }_0}\frac{\mathrm{d}{p}_0}{\mathrm{d}r}{r}^2\mathrm{d}r{\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & +{\displaystyle \int {\xi }_{r,p}}{\stackrel{\sim }{\rho }}_{1,g}{\xi }_{h,p}\frac{1}{{\rho }_0}\frac{\mathrm{d}{p}_0}{\mathrm{d}r}r\mathrm{d}r\hfill \\ \hfill & \times {\displaystyle \int Y_{\mathit{lm}}^{\ast }}{\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}·{\mathbf{\Psi }}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega },\hfill \end{array}$$(D.22)

$$\begin{array}{cc}\hfill ⟨{\mathit{\xi }}_p,\delta {\mathcal{L}}_{3-4}[{\mathit{\xi }}_p]⟩& =-{\displaystyle \int {({\xi }_{r,p})}^2}\frac{\mathrm{d}{\rho }_0}{\mathrm{d}r}\frac{\mathrm{d}{\stackrel{\sim }{\mathrm{\Phi }}}_{1,g}}{\mathrm{d}r}{r}^2\mathrm{d}r{\displaystyle \int Y_{\mathit{lm}}^{\ast }}{Y}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill & -{\displaystyle \int {\xi }_{h,p}}\frac{\mathrm{d}{\rho }_0}{\mathrm{d}r}{\xi }_{r,p}{\stackrel{\sim }{\mathrm{\Phi }}}_{1,g}r\mathrm{d}r{\displaystyle \int {\mathbf{\Psi }}_{\mathit{lm}}^{\ast }}·{\mathbf{\Psi }}_{l^{\prime }{m}^{\prime }}{Y}_{\mathit{lm}}\mathrm{d}\mathrm{\Omega }.\hfill \end{array}$$(D.23)

D.3. Flow terms

We now compute the terms due to the flow field caused by the g-mode that are involved in δℒ₅. The perturbation to the solar model due to flows is given by Eq. (14), so that, from Eq. (C.5),

$$\begin{array}{c}\hfill \delta {\mathcal{L}}_5[{\mathit{\xi }}_p]=-2{\omega }_p{\omega }_g({\mathit{\xi }}_g·\mathrm{\nabla }){\mathit{\xi }}_p.\end{array}$$(D.24)

The computation of ⟨ξ_p, δℒ₅[ξ_p]⟩ has essentially been done by Lavely & Ritzwoller (1992). This is because a g mode produces a poloidal flow field which can be expanded in vector spherical harmonics as (Lavely & Ritzwoller 1992, Eqs. (1)–(3), (C14))

$$\begin{array}{cc}& \mathit{u}(\mathit{r})={\displaystyle \sum _{s,t}}{u}_{\mathit{st}}(r)Y_{\mathit{st}}\widehat{\mathit{r}}+v_{\mathit{st}}(r){\mathbf{\Psi }}_{\mathit{st}},\hfill \end{array}$$(D.25)

$$\begin{array}{cc}& u_{\mathit{st}}(r)=-{\delta }_{s,l^{\prime }}{\delta }_{t,m^{\prime }}i{\omega }_g{\xi }_{r,g}(r),\hfill \end{array}$$(D.26)