Numerical evaluation of time-distance helioseismic sensitivity kernels in spherical geometry

The fundamental measurement in time-distance helioseismology is line-of-sight projected Doppler velocity measurements of oscillations at the solar surface. We denote the time-evolution of the velocity of oscillation at the point $\mathbf{x}$ by $\mathbf{v}(\mathbf{x},t)$. We assume that the $\mathbf{e}_{x}$ direction coincides with the direction of the line of sight, so the component of velocity that is measured is $\mathbf{e}_{x}\cdot\mathbf{v}(\mathbf{x},t)$ The frequency-domain components of the time evolution of the velocity may be denoted by $\mathbf{v}_{\omega}(\mathbf{x})$. Propagation of seismic waves are best studied in terms of the cross-covariance of wave velocity at two points $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$, which is given by

where the angular brackets indicate an ensemble average. The measured cross-covariance accounting for projection may be expressed as

While the cross-covariances of the radial components of the wave velocity depends solely on the angle between the two observation points $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$, the line-of-sight projection breaks this spherical symmetry. We chose a second pair of points $\mathbf{x}^{\prime}_{1}$ and $\mathbf{x}^{\prime}_{2}$ that are related to the pair $(\mathbf{x}_{1},\mathbf{x}_{2})$ through a rotation $R$. The cross-covariances measured at these pairs of points are related through

The knowledge of the covariance tensor $\mathbf{C}_{\omega}\left(\mathbf{x}_{1},\mathbf{x}_{2}\right)$ therefore suffices to evaluate the line-of-sight-projected cross-covariance $C_{\omega}\left(\mathbf{x}_{1}^{\prime},\mathbf{x}_{2}^{\prime}\right)$ between any two pairs of points related by a rotation.

Following Bhattacharya (2020), we express the flow velocity within the Sun as

Seismic waves interact with this velocity field, and their local propagation speeds are altered. This shows up in the measured cross-covariances $\mathbf{C}_{\omega}\left(\mathbf{x}_{1},\mathbf{x}_{2}\right)$. It is common in time-distance helioseismology to apply a filter and follow the propagation of a subsection of waves; for example, those that arrive in a short time interval. This is achieved by choosing an appropriate window that filters out wave arrivals at a particular observation point. We may obtain the difference in the time taken by seismic waves to travel between the observation points $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ by comparing the measured wave covariance with that in a model that does not account for the flow. Gizon & Birch (2002) demonstrated a way to relate the shift in travel times linearly to the flow velocity as

We may use Equation (16) and expand the sensitivity kernel $\mathbf{K}\left(\mathbf{x},\mathbf{x}_{1},\mathbf{x}_{2}\right)$ analogously on a basis of the PB VSH to rewrite Equation (17) as

The radial profiles of the spherical harmonic coefficients of the sensitivity kernel, following Bhattacharya (2020), may be expressed as

where $j_{1}$ and $j_{2}$ are wave numbers associated with the waves, and $\alpha_{1}$ and $\alpha_{2}$ indicate vector components of the wave velocity in the helicity basis. In practice, we may limit the frequency integral to the range corresponding to waves that are bound within the Sun, although we may include higher frequencies corresponding to evanescent waves by choosing appropriate boundary conditions (Gizon et al., 2017). The function $K_{\ell j_{1}j_{2}\omega;\alpha_{1}\alpha_{2}}^{\gamma}\left(r,\mathbf{x}_{1},\mathbf{x}_{2}\right)$ may be expressed as a product of a radial function of Green’s function, which we refer to as $C_{\ell j_{1}j_{2}\omega;\alpha_{1}\alpha_{2}}^{\gamma}\left(r,r_{1},r_{2}\right)$, and a window function $h_{\omega}\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)$ that selects an appropriate wave arrival, as

The detailed expressions for these radial functions may be found in Bhattacharya (2020). In this work, we focused on the angular dependence of the kernel.

The window function (Gizon & Birch, 2002) is a nonlinear function of the modeled cross-covariance, which, in turn, implies that the kernel components $K_{\gamma,\ell m}\left(r,\mathbf{x}_{1},\mathbf{x}_{2}\right)$ and $K_{\gamma,\ell m}\left(r,\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2}\right)$ are related to each other nonlinearly through the rotation operator $R$. We may, however, split the evaluation of the kernel into two parts by noting that the components of the bipolar harmonics $\mathbf{B}_{\ell m}^{\left(j_{1}\alpha_{1}\right)\left(j_{2}\alpha_{2}\right)}\left(\hat{\mathbf{n}}_{1},\hat{\mathbf{n}}_{2}\right)$ and $\mathbf{B}_{\ell m}^{\left(j_{1}\alpha_{1}\right)\left(j_{2}\alpha_{2}\right)}\left(\hat{\mathbf{n}}^{\prime}_{1},\hat{\mathbf{n}}^{\prime}_{2}\right)$ transform linearly between each other following Equation (12). The window function may be evaluated as a function of frequency while the remaining sums over $j_{1}$, $j_{2}$, $\alpha_{1}$ and $\alpha_{2}$ may be evaluated separately, and the two may be multiplied at the final stage and integrated over frequency to obtain the kernel components.

We compared the sensitivity kernel for the radial component of the flow velocity computed using line-of-sight wave velocities and compared it with that computed using a radial component of wave velocities, choosing the two observation points to lie on the equator at a height of $150$ km above the photosphere, at azimuths of $30$ and $90$ degrees. We plot the result in Fig. 1. The kernels were computed using Green’s functions with angular degrees from $1$ to $100$, and frequencies from $2$ to $4.5$ mHz. The limitation on angular degrees arises from numerical artifacts in evaluating Green’s functions using the finite-difference scheme that was presented by Bhattacharya et al. (2020) that arise at higher degrees; however, this is not a fundamental limitation of this approach as Green’s functions computed using alternate numerical schemes may be used in its place to increase the accuracy (e.g., the finite-element scheme described by Gizon et al. (2017)).

We see that the components of the kernels may be significantly different if line-of-sight projections are not considered. While a simple relation between the differences in the kernels and the angular distance from the disk center is difficult to obtain, our tests indicate that the difference is significant for observation points closer to the solar limb, where the line-of-sight projection effects are increasingly important to account for.

Given two sets of points $(\hat{\mathbf{n}}_{1},\hat{\mathbf{n}}_{2})$ and $(\hat{\mathbf{n}}^{\prime}_{1},\hat{\mathbf{n}}^{\prime}_{2})$ that are related through a rotation, the bipolar spherical harmonics evaluated at these points are related through Equation (12). We made use of this in reverse by treating it as an active rotation of the points instead of as a passive rotation of coordinate frames. The Euler angles may be computed by expressing the rotation in terms of the points. We note that for the rotation $R,$ maps $(\hat{\mathbf{n}}_{1},\hat{\mathbf{n}}_{2})$ to $(\hat{\mathbf{n}}^{\prime}_{1},\hat{\mathbf{n}}^{\prime}_{2})$ may be expressed as the product of two rotations, where the first — referred to as $R_{1}$ — maps $\hat{\mathbf{n}}_{1}$ to $\hat{\mathbf{n}}^{\prime}_{1}$, and the second maps $R_{1}\hat{\mathbf{n}}_{2}$ to $\hat{\mathbf{n}}^{\prime}_{2}$. The operation of the first rotation may be represented as

whereas the second one — referred to as $R_{2}$ — performs the transformation

The point $\hat{\mathbf{n}}^{\prime\prime}_{2}$ is an intermediate position that depends on the choice of $R_{1}$. As an aside, we note that the choice of $R_{1}$ is not unique. To illustrate this, we denote a rotation in the axis-angle formulation for an axis $\hat{\mathbf{n}}$ and angle $\gamma$ as $R(\hat{\mathbf{n}},\gamma)$. We see that the rotation $R^{\prime}_{1}=R\left(\hat{\mathbf{n}}^{\prime}_{1},\psi\right)R_{1}R\left(\hat{\mathbf{n}}_{1},\phi\right)$ for arbitrary $\psi,$ and $\phi$ also map $\hat{\mathbf{n}}_{1}$ to $\hat{\mathbf{n}}^{\prime}_{1}$. One choice for the rotation $R_{1}$ is $R_{z}\left(\phi_{1}^{\prime}\right)R_{y}\left(\theta_{1}^{\prime}-\theta_{1}\right)R_{z}\left(-\phi_{1}\right)$.

The rotation $R_{2}$ may be denoted in the axis-angle formulation as $R_{2}=R\left(\hat{\mathbf{n}}^{\prime}_{1},\omega\right)$, where $\omega$ satisfies

Collectively, the rotation $R_{2}R_{1}$ maps the set of points $(\hat{\mathbf{n}}_{1},\hat{\mathbf{n}}_{2})$ to $(\hat{\mathbf{n}}^{\prime}_{1},\hat{\mathbf{n}}^{\prime}_{2})$. Given this rotation, we may use Equation (12) to evaluate the bipolar harmonics $\mathbf{B}_{\ell m}^{\left(j_{1}\alpha_{1}\right)\left(j_{2}\alpha_{2}\right)}\left(\hat{\mathbf{n}}^{\prime}_{1},\hat{\mathbf{n}}^{\prime}_{2}\right)$ by transforming $\mathbf{B}_{\ell m}^{\left(j_{1}\alpha_{1}\right)\left(j_{2}\alpha_{2}\right)}\left(\hat{\mathbf{n}}_{1},\hat{\mathbf{n}}_{2}\right)$ instead of an explicit evaluation.

We demonstrate the transformation of sensitivity kernels in Fig. 2 by choosing the pairs of points $\hat{\mathbf{n}}_{1}=(0,0)$, $\hat{\mathbf{n}}_{2}=(\pi/4,0)$, $\hat{\mathbf{n}}^{\prime}_{1}=(\pi/2,0)$ and $\hat{\mathbf{n}}^{\prime}_{2}=(\pi/2,\pi/4)$. The solid line represents the kernels computed using Equation (2.2) for the points $\hat{\mathbf{n}}^{\prime}_{1}$ and $\hat{\mathbf{n}}^{\prime}_{2}$, while the dots represent that computed by rotating the kernel components for $\hat{\mathbf{n}}_{1}$ and $\hat{\mathbf{n}}_{2}$ using Equation (12). The bottom panel of Fig. 2 depicts the absolute difference between the kernels computed using the two approaches, and the magnitude of the difference indicates that the error is almost purely numerical. The close match shows that the components of sensitivity kernels for line-of-sight projected wave velocities for points that are related through a rotation may be transformed among each other, and kernels for several pairs of points may be computed simultaneously. In the presence of center-to-limb factors that break spherical symmetry (e.g., differences in line-formation heights), such a transformation may only be carried out among pairs of points for which symmetry-breaking factors remain unchanged on rotation.

The numerical implementation of the algorithm described in Bhattacharya (2020) has three parts — the first where the angular functions are computed, the second where the radial functions are evaluated, and the third where the products of the two are summed over. In this work, we described an approach to perform the first part more efficiently than before. We also improved upon the second and third parts by harnessing Julia’s loop vectorization capabilities. We demonstrate one comparison of the evaluation times in Fig. 3, where we plot the time required to evaluate the kernels using the line-of-sight projections of wave velocity and that using the radial component of wave velocity with the maximum angular degree up to which which the kernels are evaluated. We used Green’s functions corresponding to wave numbers in the $1$ to $100$ range and frequency in the $2$ mHz to $4.5$ mHz range uniformly divided over $4000$ bins. The relative difference between the evaluation times using line-of-sight projected velocities and that using radial velocities decreases with an increase in the cutoff in angular degree, reaching $15\%$ at $\ell_{\text{max}}=30$. The absolute difference using $\ell_{\text{max}}=30$ is $17.5$ hours of computation time. The timings were measured on $224$ Intel Xeon CPU E5-2680 v4 cores running at $2.40$GHz.