GREEN'S FUNCTIONS FOR FAR-SIDE SEISMIC IMAGES: A POLAR-EXPANSION APPROACH

The idea that the observed wave field at the solar surface can be used to estimate the wave field at any other location in the Sun has resulted in the development of several local helioseismology methods such as the seismic-holography technique (Lindsey & Braun 1997). Although the technique has been used to characterize flows in the solar subphotosphere (Braun et al. 2004) and to investigate seismic sources from flares (Donea et al. 2006), its application to map areas of strong magnetic field (Lindsey & Braun 1998) has remained a major focus. In particular, the mapping of active regions in the non-visible (or far side) of the Sun (Lindsey & Braun 1990, 2000a; González Hernández et al. 2009) has proven to be a powerful tool for space weather forecasting. A general review on local helioseismology can be found in Gizon & Birch (2005).

The use of phase-sensitive seismic holography to map areas of strong magnetic field is based on the fact that there is a travel delay or phase shift between waves entering and exiting an active region (Braun et al. 1992; Duvall et al. 1996). The technique compares theoretical waves going in and out of a particular point in the Sun (the focus) with the observed wave field, $\tilde{\Psi }_{\mathrm{obs}} (\Omega, t)$, in some portion of the solar surface (the pupil $\mathcal {P}_\pm$), by computing the ingression and egression functions given by

$$\begin{equation} \tilde{H}_{\pm }({\bf r}_0, t_0) = \int {d}t \int _{\mathcal {P}_\pm } {d}\Omega \, \tilde{\mathcal {G}}_{\pm }({\bf r}_0, t_0 ; R,\Omega, t) \, \tilde{\Psi }_{\mathrm{obs}} (\Omega, t)\,, \end{equation} \tag{ 1 }$$

where r₀ is the location of the focus, R is the solar radius, Ω = (θ, φ), $\tilde{\mathcal {G}}_+ ({\bf r}_0, t_0 ; R,\Omega, t)$ is a Green's function that represents the theoretical disturbance at (R, Ω, t) resulting from a unit impulse originating at (r₀, t₀), and $\tilde{\mathcal {G}}_-$ is the time-reverse counterpart (Lindsey & Braun 2000b, 2004).

In the quiet Sun, the wave field can be well represented by a superposition of gravito-acoustic Green's functions of a standard reference model and the ingression and egression functions are very similar. However, when the focus is located in an area of high magnetic activity, the measured signal deviates from that, effectively introducing a phase shift between the waves going in and out. This can be detected by computing the correlation that in the Fourier domain becomes the product of the Fourier transforms, H₊(r₀, ω) and H*₋(r₀, ω),

$$\begin{equation} C({\bf r}_0) \equiv \int _{\omega _1}^{\omega _2} H_+({\bf r}_0, \omega) H_-^*({\bf r}_0,\omega) {d}\omega \,. \end{equation} \tag{ 2 }$$

The phase of the correlation, $\phi ({\bf r}_0) \equiv \arg C({\bf r}_0)$, is related to the perturbed travel time by ϕ ≃ ωΔτ. The seismic-holography maps are representations of ϕ(r₀) as a measurement of the disturbance of the wave field due to the presence of magnetic regions.

Up to the present, far-side seismic holography has used Green's functions calculated using the (acoustic) ray theory, e.g., Lindsey & Braun (2000b). Formally, this approximation is only valid for high frequencies and high angular degrees. In this paper, we develop the formalism to determine Green's functions through their polar expansion. Although the method allows a full numerical computation both in the radial and angular variables, in the present work we have used the asymptotic approximation for the radial dependences of Green's function. On the other hand, the angular factors are computed without any asymptotic approximation and hence Green's functions are expected to be more accurate for low and intermediate degrees.

The outline of the paper is as follows. In Section 2, we describe the basic wave equations. Section 3 presents the calculation of Green's functions using the gravito-acoustic ray-path approximation, and Section 4 presents Green's functions using the spherical polar expansion. The phase shift introduced by the non-adiabatic sub-surface layers of the Sun is reviewed in Section 5. In Section 6, we present an example of the improvement achieved by using the new Green's functions to calculate far-side maps of solar activity.

We start the analysis with the adiabatic wave equation in the Cowling and local plane–parallel approximations. In terms of the scalar field $\tilde{\Psi }(t,{\bf r})= - \rho ^{-1/2} \delta p$, where ρ is the density of the background state and δp(t, r) is the Lagrangian pressure fluctuation, this equation is given by (see Gough 1993)

$$\begin{equation} c^{-2} \left(\frac{\partial ^2}{\partial t^2} + \omega _c^2 \right) \frac{\partial ^2 \tilde{\Psi }}{\partial t^2} - \frac{\partial ^2}{\partial t^2} \nabla ^2 \tilde{\Psi }- N^2 \nabla _{\mathrm{h}}^2 \tilde{\Psi }= 0, \end{equation} \tag{ 3 }$$

where c and N are the adiabatic sound speed and buoyancy frequency, defined as usual, ∇²_h represent the terms involving horizontal derivatives in the ∇² operator and the critical acoustic frequency is given by

$$\begin{equation} \omega _c = \frac{c}{2H} (1 - 2 \hat{\bf r}\cdot \nabla H)^{1/2} \,. \end{equation} \tag{ 4 }$$

Here H is the density-scale height, and $\hat{\bf r}$ is an upward directed unit vector.

If the background state does not depend on time, as we will assume, the scalar field can be Fourier decomposed such that the wave field is a superposition of monochromatic waves given by $\tilde{\Psi }({\bf r}, t) = \Psi ({\bf r}, \omega) \exp {(-i\omega t),}$ and Equation (3) becomes

$$\begin{equation} \nabla ^2 \Psi - \frac{N^2}{\omega ^2} \nabla _{\rm h}^2 \Psi + \frac{1}{c^2} \left(\omega ^2 - \omega _c^2 \right) \Psi = 0 \,. \end{equation} \tag{ 5 }$$

We define the operator

$$\begin{equation} \mathcal {L}_0 = \nabla ^2 - \frac{N_0^2}{\omega ^2} \nabla _{\rm h}^2 + \frac{1}{c_0^2} \left(\omega ^2 - \omega _{c0}^2 \right) \,, \end{equation} \tag{ 6 }$$

where N₀, c₀, and ω_c0 are the buoyancy frequency, sound speed, and cut-off frequency of a spherically symmetric reference model. In this work, we will use Model S from Christensen-Dalsgaard et al. (1996). We assume that the actual solar oscillations satisfy a wave equation of the form

$$\begin{equation} \mathcal {L} \Psi = \mathcal {L}_0 \Psi + \mathcal {L}^{\prime } \Psi = 0 \,, \end{equation} \tag{ 7 }$$

with $\mathcal {L}^{\prime }\Psi$ accounting for all departures from Equation (5) (e.g., non-adiabatic effects) and from the reference model (e.g., non-spherical perturbations to the sound speed).

2.1. Integral Representation

Let us take a target point r₀, the focus, anywhere in the Sun. We introduce Green's functions, G(r|r₀; ω) that satisfy the source-point wave equation for the standard solar model,

$$\begin{equation} \mathcal {L}_0 G({\bf r}| {\bf r}_0; \omega)= - \delta ({\bf r}- {\bf r}_0) \,. \end{equation} \tag{ 8 }$$

Multiplying Equation (7) by G, Equation (8) by Ψ, subtracting both, integrating over a volume $\mathcal {V}$ that containing the focus r₀ and, finally, using Green's theorem, the above equation can be written as

$$\begin{equation} \int _{S} \left[ G {\bf V} \Psi - \Psi {\bf V} G \right] \cdot {d}{\bf S} = \Psi ({\bf r}_0,\omega) + \int _{\mathcal {V}} G \mathcal {L}^{\prime } \Psi {d}\mathcal {V} \,, \end{equation} \tag{ 9 }$$

where S is the surface enclosing $\mathcal {V}$, dS is directed outward and V = ∇ − (N²₀/ω²)∇_h . If S is the solar surface with radius R, the expression is simplified to

$$\begin{equation} R^2 \oint _{S} \left(G \frac{\partial \Psi }{\partial r} - \Psi \frac{\partial G}{\partial r} \right) {d}\Omega = \Psi ({\bf r}_0, \omega) + \int _{\mathcal {V}} G \mathcal {L}^{\prime } \Psi {d}\mathcal {V} \,. \end{equation} \tag{ 10 }$$

A filtered version of the left-hand side of Equation (10) can be obtained from the observed wave field and a reference model. Moreover, if we compare the left-hand side of Equation (10) obtained from two sets of Green's functions, the differences will be due to the term with $\mathcal {L}^{\prime }$. This includes what can be called the "anomalies" (e.g., magnetic activity), but also systematic errors in the wave operator $\mathcal {L}_0$.

As indicated in the introduction we will consider Green's functions G_± corresponding to outgoing and ingoing waves from the focus. In this case, Equation (10) is a Kirchhoff representation for the wave field Ψ in the frequency domain. We compute the egression and ingression functions as

$$\begin{equation} H_{\pm }({\bf r}_0, \omega) = \int _{\mathcal {P}_{\pm }} \Psi _{\mathrm{obs}} (\Omega, \omega)\frac{\partial G_{\pm }}{\partial r} (\Omega,\omega ; {\bf r}_0) {d}\Omega \;, \end{equation} \tag{ 11 }$$

where Ψ_obs(Ω, ω) is the observed signal (Doppler shift velocity in the frequency domain), and the integral is over a limited region of the solar surface, the pupil ($\mathcal {P}_{\pm }$). Correlating the ingression and egression functions one can expect to obtain a filtered version of the second term on the right-hand side of Equation (10).

We will compare our results with those obtained from the ray theory. Rather than the acoustic approximation used by Lindsey & Braun (2000b), we will consider a more general approximation.

For high frequencies, the general solution of $\mathcal {L}_0 \Psi =0$ can be expressed as a superposition of "rays" of the form

$$\begin{equation} \Psi ({\bf r}) = A({\bf r}) \exp {\left(\pm i\int _l k \, {d}l \right),} \end{equation} \tag{ 12 }$$

where we have introduced the local wavenumber, k. The integration is along a ray path and we have explicitly indicated the two possible directions for each path. In this equation and hereafter we will remove the sub-index 0 for the reference model.

In the ray theory, the ray path and the local wavevector, ${\bf k}({\bf r})=k_r \hat{\bf r}+ {\bf k}_{\rm h}$, are determined from a first-order three-dimensional Liouville–Green expansion (a "WKB asymptotic approximation"). Specifically the dispersion relation is found to be given by

$$\begin{equation} \frac{1}{c^2} \left(\omega ^2 -\omega _c^2 \right) - k^2 + \frac{N^2}{\omega ^2} k_{\rm h}^2 = 0 \,. \end{equation} \tag{ 13 }$$

If the background state is spherically symmetric, k²_h = L²/r², where L is a constant.

The ray path is determined from the group velocity. For the spherically symmetric case, the rays are contained in a plane with a path given by

$$\begin{equation} \frac{{d}\theta }{{d}r} = \frac{v_{\theta }}{r v_r} \;. \end{equation} \tag{ 14 }$$

Here r and θ are the usual polar coordinates, and the group velocity ${\bf v} = v_r \hat{\bf r}+ {\bf v}_{\theta } \hat{\theta }$ has components

$$\begin{equation} v_r = \frac{\partial \omega }{\partial k_r} = \frac{k_r \omega ^3 c^2}{\omega ^4 - k_h^2 c^2 N^2}, \end{equation} \tag{ 15 }$$

$$\begin{equation} v_{\theta } = \frac{\partial \omega }{\partial k_{\rm h}} = k_{\rm h} \, \omega \, c^2 \left(\frac{\omega ^2 - N^2}{\omega ^4 -k_{\rm h}^2 c^2 N^2} \right) \,. \end{equation} \tag{ 16 }$$

The complex amplitude A(r) is obtained from a second-order Liouville–Green expansion called the transport equation. In our case, this equation reduces to

$$\begin{equation} \nabla \cdot \left(A^2 \left\lbrace {\bf k}- \frac{N^2}{\omega ^2} {\bf k}_{\rm h} \right\rbrace \right) =0 \,. \end{equation} \tag{ 17 }$$

Integrating this equation over the volume of a ray tube, applying the Gauss theorem, and taking ΔS → 0 for the tube section, one gets

$$\begin{equation} A^2 \left({\bf k}- \frac{N^2}{\omega ^2} {\bf k}_{\rm h} \right) \Delta S \cdot {\bf n}\simeq {\rm constant} \,, \end{equation} \tag{ 18 }$$

where n is a unit vector normal to the section. A representation of a ray tube similar to that considered here can be found in Figure 7 of Lindsey & Braun (2000b). If we take the section such that n is in the radial direction, the amplitude can be determined by

$$\begin{equation} A \, r \left(k_r \sin \theta \, \Delta \theta \right)^{1/2}= {\rm constant} \,, \end{equation} \tag{ 19 }$$

where Δθ is the angular size of the tube of section ΔS → 0.

3.1. Green's Functions

Green's functions defined by Equation (8) have also a solution of the form (12) but with prescribed initial conditions, corresponding to a Dirac-delta impulse at the focus. Thus, the (complex) constant in Equation (19) is fixed. Specifically, the amplitude for a Green's function at any point r far from the focus is given by

$$\begin{equation} A_{\pm }({\bf r}|{\bf r}_0) = \frac{\exp {(\pm i k_0 |{\bf r}- {\bf r}_{\epsilon }|)}}{4\pi |{\bf r}_{\epsilon } - {\bf r}_0|} \frac{r_{\epsilon }}{r} \left| \frac{(k_r \sin \theta)_{{\bf r}_{\epsilon }}}{k_r \sin \theta } \right|^{1/2} \left| \frac{\Delta \theta _{\epsilon }}{\Delta \theta } \right|^{1/2} \,. \end{equation} \tag{ 20 }$$

Here r is a point close to r₀ where the integration along the ray path in Equation (12) starts up. For Green's functions, here and in Equation (12), the + sign corresponds to waves moving out of the focus r₀ (amplitude A₊) and the − sign to waves moving into the focus (amplitude A₋).

We are interested in computing outward and inward Green's functions but, to be more general, after a given number of bounces at the surface and the interior. At the turning points, k_r = 0 and the amplitude can no longer be computed by means of Equation (19). In fact, we know that at the turning points A(r) is just retarded by π/2. Since these Green's functions are defined only within the resonance cavity, both r and r₀ must be inside the cavity. In particular, at least formally, Green's functions cannot be computed in the photosphere or above. Thus, for a ray that initially goes inward and then it is observed near the surface at its first arrival, the phase shift is −π/2. For any additional pair of inner and outer reflections, the phase shift is −π. So after s inner and s − 1 surface bounces (hereafter s skips), the phase shift is −(2s − 1)π/2 and the inward and outward Green's functions at a point r are given by

$$\begin{eqnarray} && G_{\pm }({\bf r}|{\bf r}_0) = A_{\pm }({\bf r}|{\bf r}_0) \exp { \left\lbrace \mp i \frac{(2s-1)\pi }{2} \right\rbrace } \nonumber \\ &&\;\;\; \times \exp { \left(\! \pm i \left\lbrace \!\int _{{\bf r}_0\rightarrow {\bf r}_1}\! k {d}l + 2 (s-1) \!\int _{{\bf r}_1 \rightarrow {\bf r}_2}\! k {d}l + \!\int _{{\bf r}_1 \rightarrow {\bf r}}\! k {d}l \right\rbrace \right) }, \nonumber\\ \end{eqnarray} \tag{ 21 }$$

where r₁ and r₂ are the inner and upper turning points, and A_±(r|r₀) is given by (20). The limits on the integrals indicate the interval over the ray path to be considered.

For a spherically symmetric reference model, Equation (8) can be solved by expanding the solution in terms of spherical harmonics in the form

$$\begin{eqnarray} G ({\bf r}| {\bf r}_0) &=& \sum _{\ell =0}^{\infty } g^{\ell } (r | r_0) \sum _{m=-\ell }^{\ell } Y_{\ell m}^* (\Omega _0) Y_{\ell m} (\Omega) \nonumber \\ &=&\frac{1}{4\pi } \sum _{\ell =0}^{\infty } g^{\ell } (r | r_0) (2\ell +1) P_{\ell }(\mu), \end{eqnarray} \tag{ 22 }$$

where P_ℓ(μ) is a Legendre polynomial and μ = cos θcos θ₀ + sin θsin θ₀cos(φ − φ₀). Here r, θ, and φ are standard spherical coordinates and Ω = (θ, φ).

Since we are considering a spherically symmetric reference model, we can choose, without losing generality, the polar axis in a such a way that θ₀ = 0 and φ₀ = 0 for the focus. In this way μ = cos θ and Green's functions have not dependence on φ.

Inserting the decomposition given by Equation (22) into Equation (8), it is found that the radial factors g_ell satisfy the equation

$$\begin{equation} \frac{{d}^2 \tilde{g}_{\ell }}{{d}r} + k_r^2 \tilde{g}_{\ell } = 0 \qquad {\rm at} \quad {\bf r}\ne {\bf r}_0, \end{equation} \tag{ 23 }$$

where $\tilde{g}_{\ell }= r g_{\ell }$ and k_r is given by Equation (13) with L² = l(l + 1). The Dirac-delta term is transformed into the following inhomogeneous initial conditions:

$$\begin{equation} \tilde{g}_{\ell } (r_0^- |r_0)= \tilde{g}_{\ell } (r_0^+| r_0) \end{equation} \tag{ 24 }$$

and

$$\begin{equation} \frac{{d}\tilde{g}_{\ell }}{{d}r} (r_0^+| r_0)- \frac{{d}\tilde{g}_{\ell }}{{d}r} (r_0^- |r_0)= -\frac{1}{r_0} \;. \end{equation} \tag{ 25 }$$

We also recall that the angular factors, given by the Legendre polynomials, satisfy the equation

$$\begin{equation} \frac{{d}^2 Q_{\ell }}{{d}\theta ^2} + \kappa _{\ell }^2 Q_{\ell } = 0, \end{equation} \tag{ 26 }$$

where $P_{\ell } (\cos \theta)= (\csc ^{1/2} \theta) Q_{\ell } (\theta)$ and

$$\begin{equation} \kappa _{\ell }^2 = \left(\ell + \frac{1}{2} \right)^2 + \frac{1}{4} \frac{1}{\sin ^2\theta } \,. \end{equation} \tag{ 27 }$$

4.1. Asymptotic Inward and Outward Green's Functions

We are interested in keeping the same kind of ingoing and outgoing Green's functions than in the ray theory, since this will allow us to do a straightforward comparison. Since we have done a separation of radial and angular variables, we need to search for inward and outward radial solutions and prograde and retrograde angular solutions and combine both in a suitable way. These kinds of "travel-wave" solutions are well defined in the asymptotic limit and we start with that.

In a first-order Liouville–Green expansion the general solution of the angular factors can be written as Q_ℓ(θ) = β_ℓ+B_ℓ+(θ) + β_ℓ−B_ℓ−(θ), where

$$\begin{equation} B_{\ell \pm } (\theta) = \kappa _{\ell }^{-1/2} \exp {\left(\pm i \int _{\theta _0}^{\theta } \kappa _{\ell } \, {d}\theta \right)} \end{equation} \tag{ 28 }$$

and β_ℓ+ and β_ℓ− are constants. The value θ₀ will be taken at the focus. In a similar way, the solution of the radial factors can be written as $\tilde{g}_{\ell } (r) = \alpha _{\ell +} C_{\ell +}(r) + \alpha _{\ell -} C_{\ell -}(r),$ where α_ℓ+ and α_ℓ− are constants and

$$\begin{equation} C_{\ell \pm } (r)= k_r^{-1/2} \exp {\left(\pm i \int _{r_1}^r k_r {d}r \right)} \,. \end{equation} \tag{ 29 }$$

Here r₁ is the inner turning point. The constants α_ℓ+ and α_ℓ− can be obtained by introducing the equation $\tilde{g}_{\ell } (r) = \alpha _{\ell +} C_{\ell +}(r) + \alpha _{\ell -} C_{\ell -}(r)$ into the initial conditions, Equations (24) and (25). One gets

$$\begin{equation} \alpha _{\ell \pm } = \frac{i}{2r_0} \frac{1}{\sqrt{k_r(r_0)}} \exp {\left(\mp i \int _{r_1}^{r_0} k_r {d}r \right)} \;. \end{equation} \tag{ 30 }$$

In principle, we can take the outward prograde part of the solution, proportional to α_ℓ+C_ℓ+(r)β_ℓ+B_ℓ+(θ), as Green's function going out of the focus. However, the solution α_ℓ+C_ℓ+ corresponds to an outgoing wave emanating from the focus r₀ and reaching the surface without any reflection. Since we are interested in the outgoing Green's function after s skips, we need to define new inward and outward asymptotic solutions. Let us call C^s_ℓ+ to the wave emanating inward from the focus and now going outward after s skips. This function and his inward counterparts are

$$\begin{equation} C_{\ell \pm }^s = C_{\ell \pm } \exp { \left[ \pm i \left\lbrace 2 (s-1) \left(\int _{r_1}^{r_2} k_r {d}r - \frac{\pi }{2} \right) -\frac{\pi }{2} \right\rbrace \right] } \;, \end{equation} \tag{ 31 }$$

where we have taken into account a phase jump of −π/2 every time the wave reach a turning point.

The angular factors have no turning points, so we can use the asymptotic functions defined in Equation (28). Then, Green's function for a given degree l and frequency ω going out of the focus, proportional to α_ℓ+C^s_ℓ+(r)β_ℓ+B_ℓ+(θ), is

$$\begin{eqnarray} &&G_ {\ell +}^s (\vec{r} | \vec{r}_0, \omega) = \frac{i(2\ell +1)}{8 \pi r_0 \sqrt{k_r (r_0})} r^{-1} k_r^{-1/2} \kappa _{\ell }^{-1/2} \csc ^{1/2}\theta \nonumber \\ &&\quad\times \exp { \left\lbrace i \left(\int ^{\theta }_{\theta _0} \kappa _l {d}\theta + \int _{r_0}^{r_1} k_r {d}r \right) \right\rbrace } \nonumber \\ &&\quad\times \exp { \left\lbrace i \left(2 (s-1) \left(\int _{r_1}^{r_2} k_r {d}r - \frac{\pi }{2} \right) - \frac{\pi }{2} + \int _{r_1}^r k_r {d}r \right) \right\rbrace }.\nonumber\\ \end{eqnarray} \tag{ 32 }$$

The ingoing Green's function is the complex conjugate, G^s_ℓ− = G^s*_ℓ+ .

4.2. Non-asymptotic Inward and Outward Green's Functions

Let us start with the angular factors. This corresponds to Equation (26) that of course has the solution $Q_{\ell }=\sqrt{\sin \theta } P_{\ell }(\cos \theta)$. To split it in prograde and retrograde parts, we written the solution as

$$\begin{equation} Q_{\ell } (\theta) = \beta _{\ell +}(\theta) B_{\ell +}(\theta) + \beta _{\ell -}(\theta) B_{\ell -}(\theta) \;, \end{equation} \tag{ 33 }$$

where B_ℓ± are still given by Equation (28). In the asymptotic approximation, β_ℓ± were taken to be constants but now they are functions that change smoothly, at least where the asymptotic approximation is expected to work. Thus, we replace the constants by osculating parameters that satisfies the additional condition

$$\begin{equation} Q_{\ell }^{\prime } (\theta)= \beta _{\ell +} (\theta) B_{\ell +}^{\prime }(\theta) + \beta _{\ell -}(\theta) B_{\ell -}^{\prime } (\theta), \end{equation} \tag{ 34 }$$

where the primes mean derivatives respect to θ. With this definition, it follows that

$$\begin{equation} \beta _{\ell +} = \frac{i}{2} \big(Q_{\ell } B_{\ell -}^{\prime } - Q_{\ell }^{\prime } B_{\ell -}\big), \qquad \beta _{\ell -} = -\frac{i}{2} \big(Q_{\ell } B_{\ell +}^{\prime } - Q_{\ell }^{\prime } B_{\ell +}\big). \end{equation} \tag{ 35 }$$

Since Q_ℓ is known, β_ℓ±, and hence the prograde β_ℓ+B_ℓ+ and retrograde parts β_ℓ−B_ℓ−, are obtained.

For reasons explained later, for the applications shown in this work we will use the asymptotic approximation for the radial factors. However, for completeness, we explain here how to obtain the non-asymptotic radial Green's functions. The general solution of Equation (23) can be written as

$$\begin{equation} \tilde{g}_{\ell } (r|r_0) = \alpha _{\ell +} (r) C_{\ell +}^s(r) + \alpha _{\ell -} (r) C_{\ell -}^s(r), \end{equation} \tag{ 36 }$$

where C^s_ℓ± are given by (31). Here again to fix the osculating parameters, we add the condition

$$\begin{equation} \tilde{g}_{\ell }^{\prime } (r|r_0)= \alpha _{\ell +} (r) \big(C_{\ell +}^s\big)^{\prime }(r) + \alpha _{\ell -}(r) \big(C_{\ell -}^s\big)^{\prime } (r) \end{equation} \tag{ 37 }$$

from which it is found that

$$\begin{equation} \alpha _{\ell +} = \frac{i}{2} \big(\tilde{g}_{\ell } \big(C_{\ell -}^s\big)^{\prime } - \tilde{g}_{\ell }^{\prime } C_{\ell -}^s\big), \quad \alpha _{\ell -} = -\frac{i}{2} \big(\tilde{g}_{\ell } \big(C_{\ell +}^s\big)^{\prime } - \tilde{g}_{\ell }^{\prime } C_{\ell +}^s\big). \end{equation} \tag{ 38 }$$

Thus, for a given solution $\tilde{g}_{\ell }(r|r_0)$, the functions α_ℓ+(r) and α_ℓ−(r) can be computed. Note that $\tilde{g}_{\ell }(r|r_0)$] satisfies Equation (23) and since we are dealing with waves reflected in the interior, the solution must be regular at the center. With this condition alone, $\tilde{g}_{\ell }(r|r_0)$ is determined up to a constant factor. On the other hand, $\tilde{g}_{\ell } (r_0|r_0)$ can be taken as the one given by the asymptotic solution, namely $\tilde{g}_{\ell } (r_0)= \alpha _{\ell -}(r_0) C_{\ell -}^s(r_0) + \alpha _{\ell +}(r_0) C_{\ell +}^s(r_0)\,$, with α_ℓ±(r₀) given by Equation (30) and C^m_ℓ± by Equation (31). This fixes $\tilde{g}_{\ell }(r|r_0)$ completely and hence the inward and outward solutions, α_ℓ−(r)C^s_ℓ−(r) and α_ℓ+(r)C^s_ℓ+(r), respectively. Combining the outward solution with the prograde solution of the angular part, the outgoing wave is given by

$$\begin{eqnarray} G_{\ell +}^s (\vec{r}|\vec{r}_0) &=& -\frac{2l+1}{16\pi r}\csc ^{1/2} \theta \; \nonumber \\ &&\times \big\lbrace Q_{\ell }(\theta) B_{\ell -}^{\prime }(\theta) - Q_{\ell }^{\prime }(\theta) B_{\ell -}(\theta) \big\rbrace B_{\ell +}(\theta) \; \nonumber \\ &&\times \big\lbrace \tilde{g}_{\ell }(r) C_{\ell -}^{\prime }(r) - \tilde{g}_{\ell }^{\prime }(r) C_{\ell -}(r) \big\rbrace C_{\ell +} (r) \,. \end{eqnarray} \tag{ 39 }$$

The ingoing Green's function is the complex conjugate. G^s_ℓ− = G^s*_ℓ+ .

In the uppermost layers of the Sun, the asymptotic approximation fails. Moreover, non-adiabatic effects and the dynamics of convection that are not included in the wave equations are too large in these layers to be neglected. These effects were taking into account by Lindsey & Braun (2004) by using empirical control correlation functions. Here we use a completely different approach.

When computing the wavefunctions in the asymptotic approximation, these effects can be taken into account by introducing an additional phase shift at the upper turning point (e.g., Gough 1984; Deubner & Gough 1984). In the same way, Green's functions can be modified to include such effects, namely,

$$\begin{equation} \mathcal {G}^s_{\ell }(r,\theta |r_0,\theta _0) = G_{\ell }^s(r,\theta |r_0, \theta _0) \exp {(2i\tilde{s}\pi \delta \alpha),} \end{equation} \tag{ 40 }$$

where 2πδα − π/2 is the phase shift after one surface reflection. To a good approximation, δα is a function of frequency alone. If we were to include just the formal bounces at the surface, $\tilde{s}=s-1$. However, the pupil and, depending on the problem, the focus are so close to the upper turning points that each of them can account for an extra shift in the phase of, say, a half bounce. Thus, in practice, one should take $\tilde{s}=s-1/2$ if the focus is below the uppermost layers or $\tilde{s}=s$ if the focus is too close to the surface.

As shown by Christensen-Dalsgaard & Pérez Hernández (1992), a theoretical phase shift at every point, δα(r), can be computed by comparing the asymptotic approximation and the numerical solutions of the wave equations. This would allow us to use Equation (40) to correct Green's functions at any point in the interior, but, unfortunately, our tests have shown that this correction is not enough to achieve good seismic images. The reason is that they did not include any non-adiabatic or convective terms in the full wave equations, and this is not a simple task.

For this work, we have computed an r-independent phase shift based on the observed frequencies of the normal mode of oscillations. This way we include all the neglected aspects in the wave equations, but the focus must be assumed to be above or below these uppermost layers. Moreover, the method is limited to the asymptotic approximation for the radial factors of the solution. Specifically, δα(ω) can be estimated by inserting the observed frequencies into the asymptotic expression

$$\begin{equation} \int _{r_1}^{r_2} k_r {d}r = \pi \left(n- \frac{1}{2} - \delta \alpha \right) \,, \end{equation} \tag{ 41 }$$

where n is an integer, the radial order of the mode. In our case, the wavenumber is given by Equation (13) and hence

$$\begin{eqnarray} &&\int _{r_1}^{r_2} \left[ \frac{1}{c^2} \left(\omega ^2 -\omega _c^2 \right) - \left(1-\frac{N^2}{\omega ^2} \right) \frac{L^2}{r^2} \right]^{1/2} {d}r\nonumber\\ &&\quad = \pi \left(n- \frac{1}{2} - \delta \alpha \right). \end{eqnarray} \tag{ 42 }$$

To obtain δα(ω), eigenfrequencies from ν = 2500 μHz to ν = 4500 μHz and degrees ℓ = 20–80 were taken from MDI data. Modes of lower degrees are rejected because they do not satisfy the Cowling approximation as good as intermediate and high-degree modes. In addition, Model S from Christensen-Dalsgaard et al. (1996) was used. A fitted polynomial of degree 10 to Equation (42) was considered. The result is presented in Figure 1. Note that πδα changes by more than 1.5 rad in the frequency range of interest.

Zoom In Zoom Out Reset image size

To validate Green's functions calculated using this new approach, we apply them to create seismic maps of active regions at the far side of the Sun. In the case of far-side seismic holography, the focus is located at (or just below) the surface of the far hemisphere of the Sun, and the pupil is an annulus surrounding the antipodes of the focus. Although seismic maps of the full far-side hemisphere can be calculated using waves following one-, two-, and three-skip ray paths before arriving at the pupil (Braun & Lindsey 2001), for simplicity, we will concentrate here only on the central part of the maps. Traditionally, this is calculated considering waves following a two-skip ray path, as described by Lindsey & Braun (2000b).

6.1. Comparison of Ray Theory and Polar Expansion

As has been done in previous work, we normalize Green's functions at each frequency to its mean absolute value. Figure 2 shows these normalized Green's functions against the polar angle for different frequencies. Results for the polar expansion and the gravito-acoustic rays are compared. Equation (39) has been used for the polar expansion, but with the asymptotic approximation for the radial factors, while for the ray theory Equation (21) has been considered. In both cases the surface phase shift shown in Figure 1 has been added according to Equation (40). Finally, what we actually show are the derivatives of Green's functions that we will use to compute the ingression and egression fields following Equation (11). The focus is located at r₀ = 0.9995R_☉, θ₀ = 0, and the Green functions are shown at the upper turning point after two skips.

Zoom In Zoom Out Reset image size

In the ray theory, for a given number of skips, and a given frequency ω, each point in the pupil corresponds to a given value of the horizontal wavenumber at the focus, k_h(r₀), and hence of L = ℓ(ℓ + 1). Here ℓ is a real constant. Our pupil goes from θ_min = 2 to θ_max = 3 rad; this corresponds to an interval in "degrees" from l_min ≃ 15 to l_max ≃ 28 for ν = 2500 μHz and l_min ≃ 27 to l_max ≃ 50 for ν = 4500 μHz. In the polar expansion, there is not such relation; in fact in Figure 2 Green's functions for degrees up to l = 60 have been added up, the same for all the frequencies.

As can be seen in Figure 2, there is a systematic shift between Green's functions computed with the two different approaches, but the solutions are closer at low angles; in fact in the ray theory these angles correspond to higher l where this approximation is expected to work better. On the other hand, the solutions of the polar-expansion approximation look more irregular at low frequencies; this is related to the fact that at these frequencies the upper turning point r₂ is too close to the focus r₀.

We use both sets of Green's functions described in this paper and the traditional ones (acoustic ray approximation) to calculate maps of magnetic activity at the non-visible hemisphere of the Sun. Figure 3 presents far-side maps calculated using Global Oscillation Network Group (GONG; Harvey et al. 1996) data for 2005 September 1 and 2 for the three sets of Green's functions. The maps clearly show the signature of active region NOAA 10808 as a large negative phase shift. The Green functions used for the analysis include frequencies computed from ν = 2500 μHz up to 4500 μHz with a step of ∼8 μHz. Angles go from θ = 2 to 3 rad with a step of 0.0016 rad. Active region NOAA 10808 appeared at the East limb of the Sun 2005 September 7.

Zoom In Zoom Out Reset image size

The signature of the active region is very similar in the maps for all the cases. However, while the noise level is similar in all the three cases, the spatial distribution is only correlated between the maps calculated using the two types of Green's functions described in this paper. This may be to the different approaches used for the sub-surface phase shift correction, explained in Section 5. However, more research is necessary in order to understand the different contributors to the noise.

6.2. Green's Functions Combining Different Number of Bounces

As mentioned before, in the polar-expansion approximation, there is not a one to one relation between the angle in the pupil and the mode degree. Thus, even for the same pupil and the same angular degrees it is possible to combine Green's functions from different bounces. In Figure 4, we show Green's functions from the polar expansion against the polar angle. In this case, we have used degrees from ℓ = 15 to ℓ = 100. Green's functions are shown for s = 2 and 3 and the sum of both (with the normalization at each frequency done after the sum). We do not consider the s = 1 case since most of the signal for the first bounce in our pupil comes from very low degrees for which the Cowling approximation does not work properly. As can be seen, the maximum amplitude of these Green functions change both with s and frequency.

Zoom In Zoom Out Reset image size

We calculate a series of far-side maps using the polar-expansion Green's functions that combine two and three skips. In Figure 5, we compare these maps with those using Green's functions that include only s = 2. A sequence of three days spanning 2006 November 27–29 showing a weaker active region just south of the equator an past the central meridian is presented. This active region was at the same location in the front site the previous Carrington rotation (NOAA 10923) and reappeared at the East limb 2006 December 5 (NOAA 10930). We find that the signature associated with the magnetic area is enhanced when using the combination of two and three skips by approximately 20%–30% compared to the two-skips approach used for the center of the map.

Zoom In Zoom Out Reset image size

We have implemented a holographic, or migration technique, that includes a factorization of Green's functions in angular and radial parts. This technique is more general than the one provided by the ray theory. We show examples of how these new Green's functions can be used to improve far-side seismic images by combining Green's functions corresponding to travel waves with two and three skips or inner reflections.

Zhao (2007) and Ilonidis et al. (2009), using a time–distance technique, have demonstrated that adding more skips to the calculation of far-side maps improves the signalto noise. However, in their analysis, they combine seismic images obtained from different skip ray paths after they are calculated. Here, we show that by combining Green's functions associated with two and three skips before calculating the map in seismic holography we improve the signal to noise in a similar way. Since the spherical geometry of the observational pupil makes the calculation of the ingression and egression functions for far-side maps computationally expensive, the solution presented here is an ideal way of improving the signal to noise without the need for increasing the computation time. Our approach is also formally different, since we sum Green's functions from different bounces without a previous normalization.

There are some previous works where Green's functions beyond the ray approximation were used, for instance in Lindsey & Braun (2004), but there, a completely plane–parallel approximation was considered and hence the technique was limited to focus in the front side and close to the surface. On the other hand, for our Green's functions, although we have shown here the application to the particular case of far-side seismic holography, the technique is general enough to be applied in the front side and for focuses under the surface, even below the convection zone, at least formally. In that scenario, the signal to noise increases considerably, and more work would be required to test the inferences.

An interesting result of our work is the fact that the sub-surface phase shift that has been recognized as a fundamental correction to any Green's function in the past can be computed by using the observed frequencies of the normal modes of oscillations. This opens a link between local and global helioseismology and perhaps the issues of the local approach can be used to learn something about the global one and vice versa.

The authors thank C. Lindsey and D. Braun for their valuable comments on the paper and for making available large part of the far-side code used here. This research was supported by grant PNAyA2007-62650 from the Spanish National Research Plan and by NASA grant NAG 5-11703. SOHO is a project of international cooperation between ESA and NASA. This work utilizes data obtained by the GONG program, managed by the National Solar Observatory, which is operated by AURA, Inc. under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Interamerican Observatory.