Finite-frequency sensitivity kernels in spherical geometry for time-distance helioseismology

TL;DR

This paper develops and validates finite-frequency sensitivity kernels in spherical geometry for time-distance helioseismology, enabling more accurate inference of the Sun's internal properties from surface wave travel times.

Contribution

It introduces a method to compute finite-frequency sensitivity kernels in spherical geometry using parallelized Green's function calculations, improving speed and accuracy.

Findings

Kernels produce travel time estimates within 0.47% of true values for uniform flows.

Travel-time errors are less than 1 ms for flows similar to meridional circulation.

Validated sensitivity kernels for sound-speed perturbations enhance helioseismic inferences.

Abstract

The inference of internal properties of the Sun from surface measurements of wave travel times is the goal of time-distance helioseismology. A critical step toward the accurate interpretation of travel-time shifts is the computation of sensitivity functions linking seismic measurements to internal structure. Here we calculate finite-frequency sensitivity kernels in spherical geometry for two-point travel-time measurements. We numerically build Green's function by solving for it at each frequency and spherical-harmonic degree and summing over all these pieces. These computations are performed in parallel ("embarrassingly"), thereby achieving significant speedup in wall-clock time. Kernels are calculated by invoking the first-order Born approximation connecting deviations in the wavefield to perturbations in the operator. Validated flow kernels are shown to produce travel times within 0.47% of the true value for uniform flows up to 750 m/s. We find that travel-time can be obtained with errors of 1 millisecond or less for flows having magnitudes similar to meridional circulation. Alongside flows, we also compute and validate sensitivity kernel for sound-speed perturbations. These accurate sensitivity kernels might improve the current inferences of sub-surface flows significantly.

Tables

Table 1. Table 1: Values of the coefficients in the polynomial expansion of damping scale γ 𝛾 \gamma (Equation ( 16 )).

$a_0$ $\left(\mu \text{Hz}\right)$	$a_1$ ${\left(\mu \text{Hz}\right)}^0$	$a_2$ ${\left(\mu \text{Hz}\right)}^{-1}$	$a_3$ ${\left(\mu \text{Hz}\right)}^{-2}$	$a_4$ ${\left(\mu \text{Hz}\right)}^{-3}$	$a_5$ ${\left(\mu \text{Hz}\right)}^{-4}$	$a_6$ ${\left(\mu \text{Hz}\right)}^{-5}$
$1.33\times {10}^{-6}$	$-3.20\times {10}^{-3}$	$3.12$	$-1.57\times {10}^3$	$4.30\times {10}^5$	$-6.19\times {10}^7$	$3.69\times {10}^9$

Equations

$$-{\rho }_0(\mathbf{r}){\left(\omega +i\gamma \right)}^2\mathbf{\xi }(\mathbf{r},\omega )=-\nabla p^{\prime }(\mathbf{r},\omega )+{\rho }^{\prime }(\mathbf{r},\omega ){\mathbf{g}}_0(\mathbf{r})+\mathbf{F}(\mathbf{r},\omega ),$$

$${\partial }_r{p}^{\prime }$$

$${\nabla }_h{p}^{\prime }$$

$$\rho^{\prime}=-\frac{1}{r^{2}}\partial_{r}\left(r^{2}\rho_{0}\xi_{r}\right)+\rho_{0}\mbox{$\bm{\nabla}$}_{h}\cdot\bm{\xi}_{h}.$$

$${\rho }^{\prime }=\frac{p^{\prime }}{c^2}+\frac{{\rho }_0}{g_0}{N}^2{\xi }_r.$$

$${\xi }_i(\mathbf{r},\omega )=\int G_{ij}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )F_j({\mathbf{r}}^{\prime },\omega )d{\mathbf{r}}^{\prime },$$

$$F_i(\mathbf{r},\omega )=\delta (\mathbf{r}-{\mathbf{r}}_s){\delta }_{ir}.$$

$${\xi }_r(\mathbf{r},\omega )$$

$${\xi }_r(\mathbf{r};\omega )$$

$$p^{\prime }(\mathbf{r};\omega )$$

$$F_r(\mathbf{r},\omega )$$

$$\mathrm{M}\left(\begin{array}[]{c}\alpha_{\ell\omega}\left(r\right)\\ \beta_{\ell\omega}\left(r\right)\end{array}\right)=\left(\begin{array}[]{c}0\\ \delta(r-r_{s})\end{array}\right),$$

$$\mathrm{M}$$

$$G_{rr}(\mathbf{r},{\mathbf{r}}_s,\omega )$$

$${\mathbf{G}}_{hr}(\mathbf{r},{\mathbf{r}}_s,\omega )$$

$$\gamma =a_0+a_1\omega +a_2{\omega }^2+a_3{\omega }^3+a_4{\omega }^4+a_5{\omega }^5+a_6{\omega }^6,$$

$${\xi }_r(r=0.2R_{\odot },\theta ,\varphi ;\omega )=0.$$

$$p^{\prime }(r=r_{\text{out}},\theta ,\varphi ;\omega )=0.$$

$${\alpha }_{\ell \omega }(r=r_{\text{in}})$$

$${\beta }_{\ell \omega }(r=r_{\text{out}})$$

$$\delta (r-r_s)\approx \frac{\text{exp}[-(r-r_s{)}^2/(2{\Delta }^2)]}{\sqrt{2\pi }{r}^2\Delta },$$

$${\xi }_i(\mathbf{r},t)={\int }_{-\infty }^{\infty }dt{G}_{ir}(\mathbf{r},{\mathbf{r}}_s,\omega )F_r({\mathbf{r}}_s,\omega )e^{i\omega t}.$$

$$F({\mathbf{r}}_s,\omega )=\text{exp}\left(-\frac{(\omega -{\omega }_0{)}^2}{2{\sigma }^2}\right),$$

$${\xi }_{\ell m}\left(r_{\text{obs}};\omega \right)=\int {\xi }_r\left(r_{\text{obs}},\theta ,\varphi ;\omega \right)Y_{\ell m}(\theta ,\varphi )d\Omega ,$$

$$P_{\ell }\left(r_{\text{obs}};\omega \right)=\frac{1}{2\ell +1}\sum _{m=-\ell }^{\ell }{\left|{\xi }_{\ell m}\left(r_{\text{obs}};\omega \right)\right|}^2.$$

$$\delta \tau ({\mathbf{r}}_r,{\mathbf{r}}_s)={\int }_{-\infty }^{\infty }dth(t)({\xi }_r({\mathbf{r}}_r,{\mathbf{r}}_s,t)-{\xi }_r^0({\mathbf{r}}_r,{\mathbf{r}}_s,t)),$$

$$h(t)=\frac{-W(t){\dot{\xi }}_{r,0}({\mathbf{r}}_r,{\mathbf{r}}_s,t)}{{\int }_{-\infty }^{\infty }d{t}^{\prime }W(t){\left[{\dot{\xi }}_{r,0}({\mathbf{r}}_r,{\mathbf{r}}_s,t)\right]}^2},$$

$$\delta \tau ={\int }_{\odot }{K}_q(\mathbf{r})\delta q\left(\mathbf{r}\right)d\mathbf{r},$$

$$K_q(\mathbf{r})\delta q(\mathbf{r})=\int d\omega G_{rj}({\mathbf{r}}_r,\mathbf{r},\omega ){\left[\delta \mathcal{L}G(\mathbf{r},{\mathbf{r}}_s,\omega )\right]}_{jr}{h}^{\ast }(\omega )F({\mathbf{r}}_s,\omega ),$$

$$c\left(\mathbf{r}\right)=c_0\left(r\right)+\delta c\left(\mathbf{r}\right),$$

$$\delta \mathcal{L}{\mathbf{G}}_r=-2\nabla ({\rho }_{0}c\delta c\nabla \cdot {\mathbf{G}}_r),$$

$$K_c(\mathbf{r})={\int }_{-\infty }^{{}^{\infty }}d\omega 2{\rho }_{0}c\nabla \cdot {\mathbf{G}}_r(\mathbf{r},{\mathbf{r}}_r,\omega )\nabla \cdot {\mathbf{G}}_r(\mathbf{r},{\mathbf{r}}_s,\omega )h^{\ast }(\omega )F({\mathbf{r}}_{s,}\omega ).$$

$$G_{ij}({\mathbf{r}}_1,{\mathbf{r}}_2,\omega )=G_{ji}({\mathbf{r}}_2,{\mathbf{r}}_1,\omega ),$$

Full text

1 Introduction

Seismic waves are observed on the solar surface by studying Doppler shifts of specific spectral lines produced in the photosphere. These waves are produced by vigorous turbulence near the solar surface, and they travel through the solar interior before resurfacing. Measuring the wave velocity field on the surface therefore opens up a window into the solar subsurface that is otherwise opaque to electromagnetic observations. Seismic waves are sensitive to subsurface features that either change the wave speed; in turn, information gleaned from surface observations of these waves can be inverted to image the interior that the wave has traversed. Local helioseismology can be used to infer, among other things, flows of various length scales inside the Sun, magnetic fields and active regions, and thermal anomalies leading to deviations in sound speed from that in the stratified hydrostatic background.

There are various approaches of relating seismic observations to subsurface features (for reviews see e.g. Gizon & Birch (2005); Gizon et al. (2010); Hanasoge et al. (2016)), one among them being time-distance helioseismology (Duvall et al., 1993) where we relate travel-time maps on the solar surface — obtained from wave cross-correlations — to interior features . Wave travel-times, as measured on the surface, will change if the wave encounters sound-speed perturbations or flows as it passes through the solar interior. Among several other approaches to relate change in travel times with perturbations in the background medium, the formalism proposed by Birch & Kosovichev (2000); Gizon & Birch (2002) using first-order Born approximations has been widely adopted. Key to the relationship between travel time shifts and perturbations in the medium is travel-time sensitivity kernel which describes how sensitive travel times are to changes in model parameters. Several authors, e.g. Jackiewicz et al. (2007); Birch & Gizon (2007); Burston et al. (2015) have used the formalism of Gizon & Birch (2002) to compute sensitivity kernels for sound speed and flows in Cartesian geometry. Cartesian formulations of the inverse problem are limited to spatial scales much smaller than the solar radius e.g. for the studies of sunspots, supergranulation etc.

Since the Sun is spherical, it is important to extend this formalism to spherical geometry to reliably image large-scale structures e.g. meridional flows (Duvall, 1979; Giles et al., 1997), differential rotation, tachocline etc. Computing these kernels is expensive and due to this limitation, several authors e.g. Zhao et al. (2013); Jackiewicz et al. (2015); Rajaguru & Antia (2015) have applied the ray approximation in place of the first-order Born approximation to compute flow-sensitivity kernels. Ray theory is an infinite frequency limit in which the travel time is sensitive only to perturbations along the ray path. Results from ray theory are reliable only if the length scale of the perturbation is significantly greater than the wavelength (Birch et al., 2001; Birch & Felder, 2004). Since the length scales over which perturbations vary are not known a priori in these inverse problems, it is important to perform inversions using the best-possible kernels. Recently, Böning et al. (2016); Gizon et al. (2016) have computed sensitivity kernels in spherical geometry. Böning et al. (2016) use a normal-mode expansion to compute Green’s function. This approach converges slowly and is therefore computationally expensive (personal communication, A. C. Birch, Gizon et al. (2016)). Gizon et al. (2016) reduce a gravity-free wave equation to a scalar equation and solve it using a finite element analysis method in an axisymmetric background.

In this work, we propose a different approach. We follow the measurement process described in Gizon & Birch (2002) to derive expressions for sensitivity kernels for sound-speed, flows and stream function in terms of Green’s function and its derivative. We numerically solve for Green’s functions in a spherically symmetric background using a finite-difference based scheme and compute kernels with high accuracy. We also show kernels for an azimuthal stream function which takes into account continuity and therefore appropriate for meridional-flow inversions. Since kernels are computed about a spherically symmetric background, so the inversions have to be linear, but we show that linearity is a good assumption for flows having magnitudes similar to meridional circulations.

2 Computing Green’s function

We consider a temporally stationary, spherically symmetric, non-rotating, non-magnetic solar model at hydrostatic equilibrium parametrized through material composition and thermodynamic properties at each point. Assuming spherical symmetry, material properties such as density and acceleration due to gravity, and thermal properties such as pressure and sound-speed depend only on the radial distance $r$ from the center of the Sun. We choose Model S (Christensen-Dalsgaard et al., 1996) as our background solar model. In further analysis, we use the symbol $\rho_{0}\left(r\right)$ to denote the radial density profile, $p_{0}\left(r\right)$ to denote the radial pressure profile, $\mathbf{g}_{0}\left(r\right)$ to denote acceleration due to gravity and $c\left(r\right)$ to denote the sound-speed. Seismic waves result in small deviations of these parameters about their equilibrium values, we denote these deviations using unsubscripted and primed variables. $\bm{\xi}(\mathbf{r},\omega)$ which is displacement vector of seismic waves follows the wave equation, where $\omega$ is temporal frequency,

[TABLE]

where $\mathbf{F}\left(\mathbf{r},\omega\right)$ denotes sources excitation, $\gamma$ is attenuation. $p^{\prime}$ and $\rho^{\prime}$ pressure and density perturbation respectively. Splitting (1) into tangential and radial components, we obtain

[TABLE]

where $\bm{\nabla}_{h}$ represents the lateral component of the gradient $\bm{\nabla}$, and $N$ is the Brunt-Väisälä frequency. The perturbed parameters are also constrained by the continuity equation,

[TABLE]

We also assume that the perturbations are adiabatic in nature, so the pressure perturbation $p^{\prime}$ and the density perturbation $\rho^{\prime}$ and radial displacement $\xi_{r}$ are related through

[TABLE]

This set of equations forms a well-determined system that we solve for quantities $\bm{\xi}$, $\rho^{\prime}$ and $p^{\prime}$. We simplify the system by eliminating $\rho^{\prime}$ and the tangential components of $\bm{\xi}$, therefore reducing the system to two equations in two unknowns: $\xi_{r}$ and $p^{\prime}$.

The displacement vector, $\bm{\xi}(\mathbf{r},\omega)$ is related to Green’s function through

[TABLE]

where indices $i$, $j$ denote $r,\,\theta,\,\phi$. We use Einstein’s summation convention here. $G_{ij}(\mathbf{r},\mathbf{r}^{\prime})$ is the seismic response of the $j$ th component of the point source, located at $\mathbf{r}^{\prime}$, measured in the $i$ th component of the displacement vector, at position $\mathbf{r}$. In order to obtain Green’s function, a radially directed point source, placed at $\mathbf{r}_{s}$ is considered as a source function in the wave equation

[TABLE]

Applying Equation (7) to Equation (6), we obtain

[TABLE]

which means that the radial and horizontal components of the displacement vector $\bm{\xi}$ for a radially directed delta function point source describe Green’s function $G_{rr}$ and $\mathbf{G}_{hr}$ respectively, where $\mathbf{G}_{hr}=\left(G_{\theta r},G_{\phi r}\right)$. We expand $\xi_{r}$, $p^{\prime}$ and source $F_{r}$ in the spherical-harmonic basis

[TABLE]

where $Y_{\ell m}(\theta,\phi)$ is the spherical harmonic of degree $\ell$ and azimuthal order $m$. Substituting Equation (9) into Equation (2) and (3), we obtain a coupled system of ordinary differential equations

[TABLE]

where

[TABLE]

Equation (10) has to be solved numerically as a function of radius for each temporal frequency $\omega$ and harmonic degree $\ell$, yielding the pair $\left(\alpha_{\ell\omega}\left(r\right),\beta_{\ell\omega}\left(r\right)\right)$. Using Equation (8), we construct components of Green’s function from $\alpha_{\ell\omega}$ and $\beta_{\ell\omega}$,

[TABLE]

2.1 Model for wave damping

Waves in the Sun have finite lifetimes, and are attenuated over a period of a few days. The decay of modes results from dynamical origins such as coupling with turbulent pressure and leakage into the atmosphere, as well as thermal ones such as radiative losses and interaction of waves with turbulent heat flux (see Houdek et al., 1999; Bhattacharya et al., 2015, and references therein). Damping of wave modes is usually modeled by adding a small imaginary component to the mode frequency, that is by setting $\omega_{n\ell}=\omega_{n\ell}^{0}+i\gamma$, where $\omega_{n\ell}^{0}$ represents the frequency of the ideal adiabatic undamped wave. Observational studies (Schou, 1999) show that the damping parameter $\gamma$ is primarily dependent on mode eigenfrequency $\omega_{n\ell}$, and to a lesser extent on the harmonic degree $\ell$ of the mode. We have plotted the measured damping parameter as a function of frequency in Fig. 1. Ignoring the $\ell$-dependence of $\gamma$, we find that it can be approximately represented as a sixth-order polynomial of frequency as

[TABLE]

where the value of the coefficients are noted in Table (1).

2.2 Boundary conditions

The system in Equation (10) has to be augmented with appropriate boundary condition to obtain solutions. We are interested in trapped modes, that is waves with frequency lying in the range $2\,\text{mHz}$ to $5.5\,\text{mHz}$; these waves are reflected back into the solar interior at the surface. The inwards reflection takes place because of a sharp increase in the acoustic cutoff frequency close to the surface. While propagating into the interior, seismic waves are refracted away from the center because of increasing sound speed, and at a specific depth — referred to as the turning point — these modes are totally internally reflected back towards the solar surface. The depth at which total internal reflection occurs, depends on the frequency and horizontal wavenumber $k_{\text{h}}=\sqrt{\ell(\ell+1)}/R_{\odot}$. This picture of waves being totally reflected back, however, is inherently ray-theoretic in nature; waves of a finite frequency are exponentially damped beyond the turning point and have a finite non-zero — albeit decaying — amplitude deeper in the interior.

We choose $0.2R_{\odot}$ as the inner boundary and we do not consider modes whose turning points are below $0.2R_{\odot}$. With no loss of generality, we can push the lower boundary closer to the core. We assume that waves corresponding to harmonic degrees greater than $20$ have turning points above $0.2R_{\odot}$ and choose $20$ as the lower cutoff of harmonic degrees in our analysis. We set the radial component of wave displacement to zero at the lower boundary, that is

[TABLE]

Beyond the outer surface, the waves with frequencies below the acoustic cutoff are exponentially damped. The pressure perturbation corresponding to the wave would rapidly decay to zero with height, which is why we peg its value to zero at the upper boundary of our domain, that is at $r=r_{\text{out}}$ the pressure perturbation $p^{\prime}$ satisfies

[TABLE]

The Equations (17) and (18) hold for all $(\theta,\,\phi)$ and from Equations (15) and (14) that is possible only if

[TABLE]

We use boundary condition (19) to solve Equation (10) for $\alpha_{l\omega}$ and $\beta_{l\omega}$.

2.3 Numerical technique

Evaluating Green’s function requires us to solve Equation (10) for each (discretized) frequency $\omega$ and harmonic degree $\ell$ that encompass the spectrum of solar seismic eigenmodes. We choose a frequency range from $2\,\text{mHz}$ to $4.5\,\text{mHz}$, split into $1250$ bins. We choose harmonic degree $\ell$ lying in a range $\left[20,\ell_{\text{max}}\right]$. The choice of the upper cutoff $\ell_{\text{max}}$ is primarily governed by the convergence of the final sensitivity kernel, since increasing the cutoff $\ell_{\text{max}}$ would also necessitate increasing the resolution of the discretized angular $\left(\theta,\phi\right)$ grid to avoid aliasing while computing wave travel-times using the kernel. Evaluating Green’s function using Equation (10) involves discretizing the radius $r$ and generating the matrix on the left-hand side; each $\left(\omega-\ell\right)$ pair leads to one matrix, leading to one set of solutions $\left(\alpha_{\ell}\left(r,\omega\right),\beta_{\ell}\left(r,\omega\right)\right)$. We use Model S to evaluate matrix elements. We choose $1596$ radial points distributed evenly in acoustic distance, leading to matrices of size $3192\times 3192$. Spherical symmetry and linearity dictates that the solutions for different $\left(\omega,\ell\right)$ pairs are independent, a fact that we utilize to compute the different solutions in parallel on a computer cluster. We solve Equation (10) using the linalg module implemented in numpy, and subsequently evaluate various components of Green’s function listed in Equation (15) in $\omega-\ell$ space. We construct the matrix in Equation (10) by discretizing derivatives on the radial grid using various stencils.

We apply a second-order backward finite difference scheme to evaluate the first derivative in Equation (10) at the boundary points. Close to the boundary except for boundary points, we use second-order central differences. Farther away from the boundary, we increase the accuracy of the central-difference scheme up to sixth order. We approximate the delta function by the following Gaussian:

[TABLE]

where $\Delta$ is the width of the function. We have chosen $\Delta=8\,\mathrm{km}$ and we place our source at $75$ km below the surface. The reason for this particular choice of $\Delta$ is to use $30$ points to resolve the Gaussian. We have considered $\ell_{\text{max}}=300$ for all the plots of sensitivity kernels in the following sections.

3 Validation of Green’s function

3.1 Time-distance diagram

The primary observation in seismology is the line-of-sight projected velocity at each point on the solar disk. Waves in the Sun are stochastically excited by turbulent convection near the surface, and the sources that excite waves are distributed randomly over the solar disk. In our analysis, we place a point source and study waves emanating from it. We record the waves as they pass through specific points on the surface that we label as “receivers”. Each source-receiver pair yields information about the sub-surface medium that the wave travels through. Waves recorded at each receiver over the entire period of observation is referred to as a time-distance diagram (for a description of time-distance diagrams and how they are obtained from observations of solar disk, see Duvall et al. (1993)).

The time-distance diagram acts as a validation test for Green’s function as we may compare it with the diagram obtained separately in the ray theory limit. In our case, we study the wave displacement instead of velocity, the former being a time-integral of the latter. The wave displacement is given by

[TABLE]

We assume a Gaussian frequency dependence of the source, that is

[TABLE]

where $\omega_{0}=2\pi\times 3.2$ mHz, and $\sigma=2\pi\times 0.4$ mHz. We use the same parameters for the computation of kernel. We compare the time-distance diagram from our simulation with ray-theory (e.g. Giles (2000)) in Fig. 2.

3.2 Power spectrum

We compute the power spectrum of the waveform in temporal and spatial frequency space. Time series of velocity amplitudes of seismic waves, recorded by the Michelson Doppler Imager (MDI, Scherrer et al., 1995) onboard the Solar and Heliospheric Observatory (SOHO, Domingo et al., 1995), have been used to generate high-resolution seismic power spectra (Rhodes et al., 1997, 1998; Schou, 1999). This provides us with a ready test for Green’s functions, in that the resonant ridges in the numerically computed function should match those observed in the Sun.

The first step in computing the power spectrum is to carry out a spherical harmonic transform of the wave displacement to obtain

[TABLE]

where $d\Omega$ is the spherical solid angle and $r_{\text{obs}}$ is the radial coordinate of the height at which observations are carried out. For simplicities, we set $r_{\text{obs}}=r_{s}$ here. Since the background model is spherically symmetric, the spectrum does not depend on the azimuthal degree $m$, therefore we average over it to obtain power at each angular mode $\ell$. The $m-$averaged power spectrum of the wave displacement is given by

[TABLE]

In Fig. 3, we compare the numerical spectrum computed from our analysis with that obtained from $72$ days MDI mode-parameter measurements by Schou (1999). We notice small mismatch between simulated and measured frequencies in Fig. 3. This may be attributed to inaccuracies in our choice of surface boundary conditions as compared to the Sun (Rhodes et al., 2001) and the imperfect modeling of surface layers in model S (Rosenthal et al., 1999).

4 Sensitivity kernels

A change in the background that the wave propagates through results in a variation in seismic waves measured at the surface. This in turn may change in the wave travel times between the source and receivers. Gizon & Birch (2002) developed a technique to compute travel times from wave cross-correlations by minimizing the misfit between the observed and model cross correlations. Their formulation, however, is not specific to cross-correlations and can be extended to other wave measurables. The use of cross-correlations is necessary for solar observations since the measured wave velocity is inherently a stochastic quantity. This is because the location of sources and excitation of waves is random. In our analysis, however, we assume that the location and excitation of the wave source is entirely deterministic. Under this assumption, we relate the wave displacement directly to travel-time shifts. Denoting the radial component of wave displacement in spherically symmetric Model S by $\xi_{r}^{0}$ and that in a different background — possibly with reduced symmetry — by $\xi_{r}$, the difference in source-receiver travel times for these two wavefields can be expressed as

[TABLE]

where $\mathbf{r}_{r}$ and $\mathbf{r}_{s}$ are the receiver and source location respectively, and the function $h(t)$ is defined as

[TABLE]

where $W(t)$ is the window function that, in our case, selects only the first arrival of the waves at the receiver point $\mathbf{r}_{r}$.

The background model can change because of various reasons, for example a local bump in the thermal properties resulting in an altered sound speed, or there being small or large scale flows that the waves propagate through and are advected by. These perturbations will leave their imprint on wave travel times. Key to seismic inference is a linear relation between wave travel-times and the model perturbation. Given a generic three-dimensional local perturbation $\delta q\left(\mathbf{r}\right)$ in the solar model, the impact it has on the travel time can be quantified as

[TABLE]

where $K_{q}$ is referred to as the sensitivity kernel. This kernel encodes information about the local impact of a perturbation on measured travel times. Viewed from the vantage of an inverse problem, the kernel also represents the gradient of travel-times in the parameter-space of the perturbation $\delta q$. In the first-order Born approximation, the sensitivity kernel $K_{q}$ obeys

[TABLE]

where $\delta\mathcal{L}$ is the change in wave operator $\mathcal{L}$ due to the change in parameter $q$ and $h^{*}(\omega)$ is the complex conjugate of the Fourier transform of the function $h(t)$. In this work, we propose an efficient way to evaluate sensitivity kernels in spherical geometry.

4.1 Sensitivity kernel for sound speed

We assume that the wave propagates through a background that has a sound speed given by

[TABLE]

where $\delta c\left(\mathbf{r}\right)$ is a small three-dimensional perturbation to the spherically symmetric sound speed $c_{0}\left(r\right)$ in Model S. The corresponding change $\delta\mathcal{L}$ in the wave operator $\mathcal{L}$ takes the form

[TABLE]

where $\mathbf{G}_{r}=(G_{rr},G_{\theta r},G_{\phi r})$. Substituting Equation (30) in Equation (28), we obtain the expression for the sound-speed kernel

[TABLE]

We have used the reciprocity relation derived from the adjoint nature of the operator (Hanasoge et al., 2011)

[TABLE]

to arrive at the Equation (31). The expression for kernel $K_{c}$ is symmetric on the interchange of the source and receiver location, $\mathbf{r}_{s}$ and $\mathbf{r}_{r}$ and this symmetry can be seen in Fig. 4. The value of the kernel is small near the ray path — as seen in “banana-doughnut” kernels in geophysics literature (Marquering et al., 1999) — and peaks near the source and receiver locations. Fresnel zones surrounding the ray path oscillate between positive and negative values.

4.2 Validation of sound-speed kernel

In order to validate the sound-speed kernel, we consider the simple scenario where the perturbation in sound speed $\delta c$ is only function of the radius $r$, and the background remains spherically symmetric. In that case, we can solve for Green’s function numerically in a manner similar to that described in Section 2, the only change being $c_{0}\rightarrow c_{0}+\delta c$. After obtaining the Green’s function for the perturbed model, we can compute the $\xi_{r}(\mathbf{r}_{r},\mathbf{r}_{s},t)$ from Equation (21). We also obtain the wave displacement $\xi_{r}^{0}(\mathbf{r}_{r},\mathbf{r}_{s},t)$ for Model S through a similar computation. Once we have both $\xi_{r}(\mathbf{r}_{r},\mathbf{r}_{s},t)$ and $\xi_{r}^{0}(\mathbf{r}_{r},\mathbf{r}_{s},t)$, we estimate the change in travel time from Equation (25) and compare it with that obtained from sound-speed kernel (27). In Fig.  5, we plot the results for several different distances between source and receiver for a particular case in which the sound speed of the model is perturbed by $10^{-3}\%$. We find the two estimates of $\delta\tau$ to be in good agreement, demonstrating that the sensitivity kernel has been computed accurately. We also compare the accuracy of the sound-speed kernel by varying the perturbation in sound speed in Fig. (6).

4.3 Sensitivity kernel for flow

In presence of a temporally stationary flow with a velocity field $\mathbf{v}\left(\mathbf{r}\right)$, there will be an advection term in the wave equation given by

[TABLE]

If the velocity field $\mathbf{v}$ is small compared to the sound-speed $c$, the change in travel time $\delta\tau$ is linearly related to $\mathbf{v}$,

[TABLE]

where $\mathbf{K}_{\mathbf{v}}$ is the sensitivity kernel for velocity. The expression for $\mathbf{K}_{\mathbf{v}}$ — in the first-order Born approximation — is

[TABLE]

where index $j$ is summed over. We compute the $\theta$ and $\phi$ component of the velocity kernel $\mathbf{K}_{\mathbf{v}}$. The basic features of the flow kernel are same as the sound-speed kernel. The expression of kernel $K_{\mathbf{v}}$ is not symmetric in the source and receiver locations, $\mathbf{r}_{s}$ and $\mathbf{r}_{r}$ and this asymmetry is reflected in Fig. (7). The flow kernel also has a small value along the ray path.

Realistic inversions for flows in the Sun should ensure mass conservation. In temporally stationary backgrounds the condition for mass conservation can be expressed as \bm{\nabla}\cdot\text{\left(\rho\mathbf{v}\right)}=0. The constraint can be enforced automatically if we derive the velocity field from a stream function $\bm{\chi}$. As we are interested in meridional circulation, we follow the approach of Rajaguru & Antia (2015) and consider an azimuthal stream function $\bm{\chi}=\chi\left(r,\theta\right)\hat{\bm{\phi}}$. The corresponding velocity field is

[TABLE]

We assume that $\chi=0$ at the solar surface. Substituting Equation (36) in Equation (34), we obtain

[TABLE]

where $K_{\chi}=\hat{\bm{\phi}}\cdot\bm{\nabla}\times\left(\frac{1}{\rho}\mathbf{K}_{\mathbf{v}}\right)$ is the sensitivity kernel for the stream function. We have computed $K_{\chi}$ and it is shown is Fig. (8). The values of this kernel increases rapidly close to the surface, therefore we multiply it with density before plotting to highlight the functional variation with depth. The kernel, $K_{\chi}$ is shown in Fig (8). The grainy pattern near the surface is reminiscent of those observed by Böning et al. (2016) and Gizon et al. (2016). It appears due to the finite cutoff in $\ell_{\text{max}}$ chosen to compute the Green’s function. Increasing $\ell_{\text{max}}$ appears to further localize the pattern to shallower layers.

4.4 Validation of flow kernel

To test the accuracy of the kernel, we consider a flow field equivalent to a solid-body rotation, thereby retaining spherical symmetry. In this case the velocity field will have the form

[TABLE]

where $\Omega$ is the angular velocity of the rotation. The perturbed wave field $\bm{\xi}$ will be related to unperturbed wave field $\bm{\xi}_{0}$ through a change in reference frame

[TABLE]

where $\Delta$ is the angular distance between source and receiver. We place both the source and receiver on the equator. Using Equation (39), we compute $\delta\tau$ from Equation (25) and compare it with that obtained from Equation (34). We plot the dependence $\delta\tau$ on the strength of flow velocity in Fig. 9. In Fig. 10, we show the dependence of $\delta\tau$ with source-receiver distance when the flow speed is $20\,\text{m/s}$ at the surface, and compare the difference between the values computed through the two techniques mentioned above.

5 Summary and discussion

In this paper, we have developed a technique to compute seismic sensitivity kernels in spherical geometry using the first-order Born approximation. Computation of spherical sensitivity kernels are typically expensive. We have shown that assuming a spherically symmetric background, Green’s function decouples in frequency and harmonic degree and therefore computation of each frequency and harmonic degree can be done efficiently in parallel on a computer cluster. It takes around $16$ seconds to compute the displacement vector and pressure perturbation of equation (10) for each $(\ell,\,\omega)$ pair on a single processor. For the parameters chosen in this work, the entire Green’s function takes around six hours to compute when evaluated in parallel using $300$ processors on a computer cluster. It takes a further hour to compute the sensitivity kernel from the Green’s function.

We have studied in this work how weak flow has to be in order for linear relationship between travel-time delay and flow to hold. We have found that travel times can be obtained within $0.47\%$ accuracy using the flow kernel computed through our approach for uniform flows up to $750$ m/s. Since the observed velocity of meridional circulation on the solar surface is around $20$ m/s, we expect that linearity might be an appropriate assumption for the study of meridional circulation.

We have considered a single deterministic source in our work.** **In the case of uniformly distributed sources and for certain types of wave damping, it can be shown that (Snieder, 2004, 2007) the positive and negative branches of the cross-correlation measurement may be interpreted as waves originating from one measurement pixel to the other and vice versa. This equivalence between cross correlations and Green’s function, while possibly not very accurate in the Sun owing to line-of-sight projection and a complicated damping mechanism (among other effects), represents a useful starting point. Indeed, travel-time inversions of meridional circulation are typically performed using kernels computed in the ray approximation (e.g. Giles et al. (1997); Zhao et al. (2013); Rajaguru & Antia (2015)). Ray theory assumes that the wave frequency is infinite, relies on a single-source picture and does not take into account line-of-sight projection. In contrast, the Born approximation can account for line-of-sight projection and because it is a finite-frequency model, is more accurate than ray theory. Therefore kernels based on the Born approximation, computed in the single-source picture, though not the best, are still better to use for inversions than kernels computed using ray theory. A more complete theory would aim to model the cross-correlation measurement and take into account line-of-sight projection effects, which will be a part of our future work.

KM, JB & SMH acknowledge the financial support provided by the Department of Atomic Energy, India. SMH also acknowledges support from Ramanujan fellowship SB/S2/RJN-73/2013, the Max-Planck partner group program and thanks the Center for Space Science, New York University at Abu Dhabi.