© ESO 2022

1. Introduction

Helioseismology is a commonly used tool to infer large-scale subsurface features in the Sun, such as zonal and meridional flows. This works by relating surface measurements of seismic observables – such as wave travel-times – to subsurface features such as flows, through a sensitivity kernel that describes the response of the observable to a change in the background through which the waves propagate. The sensitivity kernel depends both on the details of the background model such as the thermal structure and the presence of inhomogeneous features such as flows, as well as on specifics of the measurement procedure. Recently, considerable effort has been made to develop frameworks to evaluate sensitivity kernels that describe wave propagation in the first Born approximation while accounting for the spherical geometry of the Sun (Böning et al. 2016; Mandal et al. 2017; Gizon et al. 2017; Fournier et al. 2018). The present work builds on the results of Bhattacharya et al. (2020) and Bhattacharya (2020), in which the authors described the procedure to evaluate sensitivity kernels by accounting for line-of-sight projections that are inherent in Doppler measurements of wave velocity, as well as the differences in line-formation height between the solar disk center and the limb.

In this work, we describe the numerical scheme in computing the sensitivity kernel in detail. Secondly, we improve upon the algorithm presented in the previous work to evaluate the angular dependence of the sensitivity kernels, which makes it possible to evaluate the kernels for waves with arbitrarily high spherical harmonic degrees, which was not possible in the approach presented by Bhattacharya (2020). We also investigated the impact of line-of-sight projection on the kernels, and demonstrate that the kernels for line-of-sight projected measurements may differ significantly from that computed for the radial components. Furthermore, we illustrate how the assumed spherical symmetry in the background solar model makes the computation of kernels numerically more efficient even if line-of-sight projections are included, as pairs of points that are related through a rotation have the corresponding kernels related through a unitary transformation.

We provide a Julia implementation of the algorithm described here and in Bhattacharya (2020) that may be used to evaluate helioseismic sensitivity kernels in spherical geometry. We chose Julia (Bezanson et al. 2017) as the language because it merges superior performance with the productivity of a high-level language. The code is modular and features various subsections that might be useful in other related analyses of functions in spherical geometry.

2. Theory

We denote a point in the Sun by x, which, in spherical polar coordinates, may be represented in terms of the radial coordinate r, the co-latitude θ and the azimuthal angle ϕ. The sensitivity kernel for flows in the Sun, K(x,x₁,x₂) is a vector function that depends on two observation points x₁ and x₂ and a third point x where the flow velocities are to be inferred. We follow Bhattacharya (2020) and expand the sensitivity kernels in a basis of vector spherical harmonics in x. The coefficients of expansion are functions of the two angular coordinates x₁ and x₂, which we further expand in a basis of bipolar spherical harmonics. In the following sections, we expand upon the approach to evaluate the angular functions involved in the analysis, while we use the radial functions as described by Bhattacharya (2020).

2.1. Angular functions

2.1.1. Vector spherical harmonics

Vector spherical harmonics (VSH) form a complete basis that may be used to decompose vector fields on a sphere. Analogous to scalar spherical harmonics, the vector spherical harmonics are eigenfunctions of the angular momentum operator and may be uniquely labeled by the angular degree j, the azimuthal order m, and an index α that indicates the direction of the vector. The choice of the vector indices is not unique, as linear combinations of one basis may also form another complete basis. This provides us with the freedom to choose the most convenient set to work with. The wave equation may be written most conveniently in terms of the Hansen vector spherical harmonics (Hansen 1935):

$$\begin{array}{cc}\hfill {\mathit{H}}_{\mathit{jm}}^{(-1)}\left(\widehat{\mathit{n}}\right)& ={\mathit{e}}_r{Y}_{\mathit{jm}}\left(\widehat{\mathit{n}}\right),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\mathit{H}}_{\mathit{jm}}^{\left(0\right)}\left(\widehat{\mathit{n}}\right)& =\frac{-i}{\sqrt{j(j+1)}}({\mathit{e}}_r\times {\mathbf{\nabla }}_{\mathrm{\Omega }})Y_{\mathit{jm}}\left(\widehat{\mathit{n}}\right),\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill {\mathit{H}}_{\mathit{jm}}^{(+1)}\left(\widehat{\mathit{n}}\right)& =\frac{1}{\sqrt{j(j+1)}}{\mathbf{\nabla }}_{\mathrm{\Omega }}{Y}_{\mathit{jm}}\left(\widehat{\mathit{n}}\right),\hfill \end{array}$$(3)

where e_r is the radial unit vector, ∇_Ω is the angular gradient operator, $\widehat{\mathit{n}}=(\theta ,\varphi )$ represents a point on the unit sphere, and $Y_{\mathit{jm}}\left(\widehat{\mathit{n}}\right)$ is the scalar spherical harmonic. This is because conservative restoring forces that appear in the helioseismic wave equation may be expressed in terms of the harmonics ${\mathit{H}}_{\mathit{jm}}^{(-1)}$ and ${\mathit{H}}_{\mathit{jm}}^{(+1)}$ with no component directed along ${\mathit{H}}_{\mathit{jm}}^{\left(0\right)}$.

The evaluation of the sensitivity kernels is easier in a linear combination of the Hansen harmonics that we refer to as the Phinney-Burridge harmonics, following their use by Phinney & Burridge (1973):

$$\begin{array}{cc}\hfill {\mathit{P}}_{\mathit{jm}}^{+1}\left(\widehat{\mathit{n}}\right)& =\frac{1}{\sqrt{2}}({\mathit{H}}_{\mathit{jm}}^{(+1)}\left(\widehat{\mathit{n}}\right)-{\mathit{H}}_{\mathit{jm}}^{\left(0\right)}\left(\widehat{\mathit{n}}\right)),\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill {\mathit{P}}_{\mathit{jm}}^0\left(\widehat{\mathit{n}}\right)& ={\mathit{H}}_{\mathit{jm}}^{(-1)}\left(\widehat{\mathit{n}}\right),\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill {\mathit{P}}_{\mathit{jm}}^{-1}\left(\widehat{\mathit{n}}\right)& =\frac{1}{\sqrt{2}}({\mathit{H}}_{\mathit{jm}}^{(+1)}\left(\widehat{\mathit{n}}\right)+{\mathit{H}}_{\mathit{jm}}^{\left(0\right)}\left(\widehat{\mathit{n}}\right)).\hfill \end{array}$$(6)

The advantage of the PB harmonics over the Hansen ones is that these are diagonal at each point in the covariant helicity basis, defined as

$$\begin{array}{cc}\hfill {\mathit{e}}_{+1}& =-\frac{1}{\sqrt{2}}({\mathit{e}}_{\theta }+i{\mathit{e}}_{\varphi }),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill {\mathit{e}}_0& ={\mathit{e}}_r,\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill {\mathit{e}}_{-1}& =\frac{1}{\sqrt{2}}({\mathit{e}}_{\theta }-i{\mathit{e}}_{\varphi }),\hfill \end{array}$$(9)

with the diagonal components being generalized spherical harmonics ${\mathit{e}}_{\beta }·{\mathit{P}}_{\mathit{jm}}^{\alpha }\left(\widehat{\mathit{n}}\right)={\delta }_{\alpha \beta }{Y}_{\mathit{jm}}^{\alpha }\left(\widehat{\mathit{n}}\right)$, where δ_αβ is the Kronecker delta function, and the generalized spherical harmonics $Y_{\mathit{jm}}^{\alpha }\left(\widehat{\mathit{n}}\right)$ are normalized elements of the Wigner D-matrix (Dahlen & Tromp 1998). This implies that on this basis, we only require three components instead of nine to describe the vector harmonics, which provides a significant computational advantage.

2.1.2. Bipolar spherical harmonics

Analogous to using a basis of scalar or vector spherical harmonics to decompose one-point functions on a sphere, a two-point function may be decomposed on a basis of bipolar spherical harmonics. These may be evaluated through angular-momentum coupling of monopolar PB harmonics as

$$\begin{array}{c}\hfill {\mathit{B}}_{\ell m}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2)={\displaystyle \sum _{m_1{m}_2}{C}_{j_1{m}_1{j}_2{m}_2}^{\ell m}{\mathit{P}}_{j_1{m}_1}^{{\alpha }_1}\left({\widehat{\mathit{n}}}_1\right){\mathit{P}}_{j_2{m}_2}^{{\alpha }_2}\left({\widehat{\mathit{n}}}_2\right),}\end{array}$$(10)

where $C_{j_1{m}_1{j}_2{m}_2}^{\ell m}$ is a Clebsch-Gordan coefficient and ${\mathit{P}}_{\mathit{jm}}^{\alpha }\left({\widehat{\mathit{n}}}_i\right)$ represents a monopolar PB vector harmonic at the point ${\widehat{\mathit{n}}}_i$. These bipolar harmonics are rank-2 tensors, which are eigenfunctions of the angular momentum operator J_z and satisfy similar properties such as orthogonality and completeness. We may use these bases in evaluating cross-covariances that are two-point functions on the solar surface. In our analysis, we sought the projection of the bipolar harmonics along the line of sight at each observation point. Assuming that the line of sight at each point on the Sun is directed along e_x, we obtain the following projected harmonic:

$$\begin{array}{c}\hfill B_{\ell m,xx}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2)={\mathit{e}}_x{\mathit{e}}_x:{\mathit{B}}_{\ell m}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2),\end{array}$$(11)

where the double contraction operator “:” is defined using the Einstein summation convention as A : B = A_ijB_ij.

2.1.3. Rotation of harmonics

One advantage of expressing two-point functions on the basis of bipolar spherical harmonics is that their behavior under a rotation of coordinates may be described as a unitary transformation. We assume that the points ${\widehat{\mathit{n}}}_1$ and ${\widehat{\mathit{n}}}_2$ have coordinates (θ₁, ϕ₁) and (θ₂, ϕ₂) in a frame S, and $({\theta }_1^{\prime },{\varphi }_1^{\prime })$ and $({\theta }_2^{\prime },{\varphi }_2^{\prime })$ in a frame S′, where the frames S and S′ are related through S′=RS for some rotation R. We refer to the points ${\widehat{\mathit{n}}}_1$ and ${\widehat{\mathit{n}}}_2$ as ${\widehat{\mathit{n}}}_1^{\prime }$ and ${\widehat{\mathit{n}}}_2^{\prime }$ in the frame S′, with the understanding that the points with and without the primes are identical, and the primes refer to the difference in coordinates that arise from a change in basis between the frames. Under such a rotation, the components of bipolar spherical harmonics satisfy

$$\begin{array}{cc}\hfill B_{\ell m,xx}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathbf{n}}}_1^{\prime },{\widehat{\mathbf{n}}}_2^{\prime })& =\left(R^{-1}{\mathbf{e}}_x\right)\left(R^{-1}{\mathbf{e}}_x\right):\hfill \\ \hfill & {\displaystyle \sum _{m^{\prime }}{D}_{m^{\prime }m}^{\ell }\left(R\right){\mathit{B}}_{\ell m^{\prime }}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2),}\hfill \end{array}$$(12)

where $D_{m^{\prime }m}^{\mathcal{l}}\left(R\right)$ is the Wigner D-matrix corresponding to the rotation. In our analysis, we evaluate the dot products in the helicity basis, in which case the rotation matrices become $R^{\prime -1}=A{(\widehat{\mathit{n}})}^{\ast }{R}^{-1}A{(\widehat{\mathit{n}})}^T$, where the matrix A transforms between the Cartesian and helicity bases and has elements $A_{ij}={\mathit{e}}_j\cdot {\mathit{e}}_j^{\prime }$, and the superscript T represents a transpose, while an asterisk as a superscript denotes complex conjugation.

2.2. Sensitivity kernels

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 x by v(x, t). We assume that the e_x direction coincides with the direction of the line of sight, so the component of velocity that is measured is e_x ⋅ v(x, t) The frequency-domain components of the time evolution of the velocity may be denoted by v_ω(x). Propagation of seismic waves are best studied in terms of the cross-covariance of wave velocity at two points x₁ and x₂, which is given by

$$\begin{array}{cc}\hfill {\mathit{C}}_{\omega }({\mathit{x}}_1,{\mathit{x}}_2)& =\langle {\mathit{v}}_{\omega }^{\ast }\left({\mathit{x}}_1\right){\mathit{v}}_{\omega }\left({\mathit{x}}_2\right)\rangle ,\hfill \end{array}$$(13)

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

$$\begin{array}{cc}\hfill C_{\omega }({\mathit{x}}_1,{\mathit{x}}_2)& ={\mathit{e}}_x{\mathit{e}}_x:{\mathit{C}}_{\omega }({\mathit{x}}_1,{\mathit{x}}_2).\hfill \end{array}$$(14)

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

$$\begin{array}{cc}\hfill C_{\omega }({\mathit{x}}_1^{\prime },{\mathit{x}}_2^{\prime })& =\left(R^{-1}{\mathit{e}}_x\right)\left(R^{-1}{\mathit{e}}_x\right):{\mathit{C}}_{\omega }({\mathit{x}}_1,{\mathit{x}}_2).\hfill \end{array}$$(15)

The knowledge of the covariance tensor C_ω(x₁,x₂) therefore suffices to evaluate the line-of-sight-projected cross-covariance $C_{\omega }({\mathit{x}}_1\prime ,{\mathit{x}}_2\prime )$ between any two pairs of points related by a rotation.

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

$$\begin{array}{cc}\hfill \mathit{u}\left(\mathit{x}\right)& ={\displaystyle \sum _{\ell m\gamma }{u}_{\ell m}^{\gamma }\left(r\right){\mathit{P}}_{\ell m}^{\gamma }\left(\widehat{\mathit{n}}\right).}\hfill \end{array}$$(16)

Seismic waves interact with this velocity field, and their local propagation speeds are altered. This shows up in the measured cross-covariances C_ω(x₁,x₂). 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 x₁ and x₂ 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

$$\begin{array}{c}\hfill \delta \tau ({\mathit{x}}_1,{\mathit{x}}_2)={\displaystyle \int }\mathrm{d}\mathit{x}\mathit{K}(\mathit{x},{\mathit{x}}_1,{\mathit{x}}_2)·\mathit{u}\left(\mathit{x}\right).\end{array}$$(17)

We may use Eq. (16) and expand the sensitivity kernel K(x,x₁,x₂) analogously on a basis of the PB VSH to rewrite Equation (17) as

$$\begin{array}{c}\hfill \delta \tau ({\mathit{x}}_1,{\mathit{x}}_2)={\displaystyle \sum _{\ell m\gamma }{\int }_0^{R_{\odot }}{r}^2\mathrm{d}r{K}_{\ell m\gamma }(r,{\mathit{x}}_1,{\mathit{x}}_2)u_{\ell m}^{\gamma }\left(r\right).}\end{array}$$(18)

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

$$\begin{array}{cc}\hfill K_{\gamma ,\ell m}(r,{\mathit{x}}_1,{\mathit{x}}_2)=& {\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\omega }{2\pi }\sum _{j_1{j}_2}\sum _{{\alpha }_1{\alpha }_2}{K}_{\ell j_1{j}_2\omega ;{\alpha }_1{\alpha }_2}^{\gamma }(r,{\mathit{x}}_1,{\mathit{x}}_2)\times }\hfill \\ \hfill & {\mathit{e}}_x{\mathit{e}}_x:{\mathit{B}}_{\ell m}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2),\hfill \end{array}$$(19)

where j₁ and j₂ are wave numbers associated with the waves, and α₁ and α₂ 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 }(r,{\mathit{x}}_1,{\mathit{x}}_2)$ 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 }(r,r_1,r_2)$, and a window function h_ω(r₁,r₂) that selects an appropriate wave arrival, as

$$\begin{array}{c}\hfill K_{\ell j_1{j}_2\omega ;{\alpha }_1{\alpha }_2}^{\gamma }(r,{\mathit{x}}_1,{\mathit{x}}_2)=2\mathfrak{R}\left[h_{\omega }^{\ast }({\mathbf{x}}_1,{\mathit{x}}_2)C_{\ell j_1{j}_2\omega ;{\alpha }_1{\alpha }_2}^{\gamma }(r,r_1,r_2)\right].\end{array}$$(20)

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_γ, ℓm(r,x₁,x₂) and $K_{\gamma ,\ell m}(r,{\mathit{x}}_1^{\prime },{\mathit{x}}_2^{\prime })$ 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 ${\mathit{B}}_{\ell m}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2)$ and ${\mathit{B}}_{\ell m}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1^{\prime },{\widehat{\mathit{n}}}_2^{\prime })$ transform linearly between each other following Eq. (12). The window function may be evaluated as a function of frequency while the remaining sums over j₁, j₂, α₁ and α₂ may be evaluated separately, and the two may be multiplied at the final stage and integrated over frequency to obtain the kernel components.

3. Method

We evaluated the vector harmonics by noting that these may be generated by coupling unit vectors and scalar spherical harmonics through the mechanism of angular momentum addition. Specifically, we constructed the following linear combination of the Cartesian basis vectors:

$$\begin{array}{cc}\hfill {\mathit{\chi }}_{-1}& =-\frac{1}{\sqrt{2}}({\mathit{e}}_x+i{\mathit{e}}_y),\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill {\mathit{\chi }}_0& ={\mathit{e}}_z,\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill {\mathit{\chi }}_{+1}& =\frac{1}{\sqrt{2}}({\mathit{e}}_x-i{\mathit{e}}_y).\hfill \end{array}$$(23)

We followed Varshalovich et al. (1988) and refer to this basis as the covariant spherical basis. We coupled this basis with scalar spherical harmonics to construct the vector spherical harmonics:

$$\begin{array}{cc}\hfill {\mathit{Y}}_{\mathit{jm}}^{\ell }\left(\widehat{\mathit{n}}\right)& ={\displaystyle \sum _{\gamma =-1}^1\sum _{n=-\ell }^{\ell }{C}_{\ell n1\mu }^{\mathit{jm}}{Y}_{\ell n}\left(\widehat{\mathit{n}}\right){\mathit{\chi }}_{\mu }.}\hfill \end{array}$$(24)

The Clebsch-Gordan relations involved in this sum may be evaluated explicitly, and we list these in Appendix A. For a detailed discussion on these harmonics, we invite the reader to consult James (1976).

We may show (Varshalovich et al. 1988) that these harmonics are related to the Hansen ones through

$$\begin{array}{cc}\hfill {\mathit{H}}_{\mathit{jm}}^{\left(1\right)}\left(\widehat{\mathit{n}}\right)& =\sqrt{\frac{j+1}{2j+1}}{\mathit{Y}}_{\mathit{jm}}^{j-1}\left(\widehat{\mathit{n}}\right)+\sqrt{\frac{j}{2j+1}}{\mathit{Y}}_{\mathit{jm}}^{j+1}\left(\widehat{\mathit{n}}\right),\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill {\mathit{H}}_{\mathit{jm}}^{\left(0\right)}\left(\widehat{\mathit{n}}\right)& ={\mathit{Y}}_{\mathit{jm}}^j\left(\widehat{\mathit{n}}\right),\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill {\mathit{H}}_{\mathit{jm}}^{(-1)}\left(\widehat{\mathit{n}}\right)& =\sqrt{\frac{j}{2j+1}}{\mathit{Y}}_{\mathit{jm}}^{j-1}\left(\widehat{\mathit{n}}\right)-\sqrt{\frac{j+1}{2j+1}}{\mathit{Y}}_{\mathit{jm}}^{j+1}\left(\widehat{\mathit{n}}\right).\hfill \end{array}$$(27)

This provides us with a route to get from scalar spherical harmonics to the vector ones that we are interested in. For numerical efficiency, we used the conjugation relation

$$\begin{array}{cc}\hfill {\mathit{Y}}_{\mathit{jm}}^{\ell \ast }\left(\widehat{\mathit{n}}\right)& ={(-1)}^{j+\ell +m+1}{\mathit{Y}}_{j-m}^{\ell }\left(\widehat{\mathit{n}}\right)\hfill \end{array}$$(28)

to evaluate the harmonics for negative m.

We developed a Julia implementation of this algorithm to evaluate vector spherical harmonics, which is available freely under the MIT license as the package VectorSphericalHarmonics.jl¹. In this, we used the spherical harmonics code by Limpanuparb & Milthorpe (2014), which can generate accurate values of the harmonics till angular degrees of around 1000, with an absolute and relative accuracy of 10⁻¹⁰. Since the components of the vector harmonics are linear combinations of normalized spherical harmonics, the accuracy of these functions are similar. This approach of computing vector spherical harmonics differs from the previous approach by Bhattacharya (2020), where the authors had chosen to evaluate them through the computation of Wigner D-matrices, where the accuracy was limited to degrees of around 100. The present approach offers accuracy to a much higher degree, as well as being computationally more efficient.

The monopolar scalar or vector harmonics may be coupled to evaluate bipolar harmonics. We recognize that the expressions for the sensitivity kernels may be cast in a way that involves the harmonics $B_{\ell m,xx}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2)$ as well as $B_{\ell m,xx}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_2,{\widehat{\mathit{n}}}_1)$. We may use the symmetries of the Clebsch-Gordan coefficients to obtain

$$\begin{array}{cc}\hfill B_{\ell m,xx}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2)& ={(-1)}^{j_1+j_2-\ell }{B}_{\ell m,xx}^{\left(j_2{\alpha }_2\right)\left(j_1{\alpha }_1\right)}({\widehat{\mathit{n}}}_2,{\widehat{\mathit{n}}}_1).\hfill \end{array}$$(29)

We used this relation to avoid the need to evaluate certain bipolar harmonics, and we used pre-computed values instead.

We developed a Julia code to evaluate bipolar spherical harmonics that is freely available under the MIT license as the package BipolarSphericalHarmonics.jl². We evaluated the Clebsch-Gordan coefficients using the freely available FORTRAN package SHTOOLS (Wieczorek & Meschede 2018) for low angular degrees, and the Julia implementation of the public code wigxjpf (Johansson & Forssén 2016) for high angular degrees. The former is computationally more efficient, using a double precision implementation of the algorithm presented by Luscombe & Luban (1998) and accurate up to angular degrees of around 160, whereas the latter is accurate at higher angular degrees as it carries out computation using arbitrary-precision arithmetic. The Julia implementation of the wigxjpf algorithm is available as the package WignerSymbols.jl. The ability to evaluate the Clebsch-Gordan coefficients for arbitrarily high degrees is an improvement over previous work by Bhattacharya (2020), as the kernels may now be computed using a wider range of wave modes.

We followed Feng et al. (2015) to evaluate Wigner D-matrices that feature in Eq. (12). A Julia implementation of this code is available freely under the MIT license as the package WignerD.jl³. This uses an exact diagonalization of the angular momentum operator J_y to evaluate the matrix elements, an approach that is tested to be accurate in absolute values up to degrees of around 100, although relative errors may be dominated by machine precision for values close to zero. The latter does not pose a significant challenge in our analysis, given that the harmonics are normalized, and the rotated harmonics depend overwhelmingly on the D-matrix elements with absolute values on the order of 1. Further tests might be necessary to use the rotation relation for bipolar harmonics at higher angular degrees.

4. Results

4.1. Sensitivity kernels

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

Fig. 1. Radial profiles of spherical-harmonic components of the sensitivity kernel for the radial component (γ = 0) of flow velocity. The gray line is computed using radial components of wave velocity, whereas the black doted line is computed using the line-of-sight projected components. The title indicates the mode indices (ℓ,m).

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.

4.2. Rotation of coordinates

Given two sets of points $({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2)$ and $({\widehat{\mathit{n}}}_1^{\prime },{\widehat{\mathit{n}}}_2^{\prime })$ that are related through a rotation, the bipolar spherical harmonics evaluated at these points are related through Eq. (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 $({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2)$ to $({\widehat{\mathit{n}}}_1^{\prime },{\widehat{\mathit{n}}}_2^{\prime })$ may be expressed as the product of two rotations, where the first – referred to as R₁ – maps ${\widehat{\mathit{n}}}_1$ to ${\widehat{\mathit{n}}}_1^{\prime }$, and the second maps $R_1{\widehat{\mathit{n}}}_2$ to ${\widehat{\mathit{n}}}_2^{\prime }$. The operation of the first rotation may be represented as

$$\begin{array}{cc}\hfill R_1{\widehat{\mathit{n}}}_1& ={\widehat{\mathit{n}}}_1^{\prime },\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill R_1{\widehat{\mathit{n}}}_2& ={\widehat{\mathit{n}}}_2^{″},\hfill \end{array}$$(31)

whereas the second one – referred to as R₂ – performs the transformation

$$\begin{array}{c}\hfill R_2{\widehat{\mathit{n}}}_2^{″}={\widehat{\mathit{n}}}_2^{\prime }.\end{array}$$(32)

The point ${\widehat{\mathit{n}}}_2^{″}$ is an intermediate position that depends on the choice of R₁. As an aside, we note that the choice of R₁ is not unique. To illustrate this, we denote a rotation in the axis-angle formulation for an axis $\widehat{\mathit{n}}$ and angle γ as $R(\widehat{\mathit{n}},\gamma )$. We see that the rotation $R_{\prime }^1=R({\widehat{\mathit{n}}}_1^{\prime },\psi )R_{1}R({\widehat{\mathit{n}}}_1,\varphi )$ for arbitrary ψ, and ϕ also map ${\widehat{\mathit{n}}}_1$ to ${\widehat{\mathit{n}}}_1^{\prime }$. One choice for the rotation R₁ is $R_z({\varphi }_1\prime )R_y({\theta }_1\prime -{\theta }_1)R_z(-{\varphi }_1)$.

The rotation R₂ may be denoted in the axis-angle formulation as $R_2=R({\widehat{\mathit{n}}}_1^{\prime },\omega )$, where ω satisfies

$$\begin{array}{cc}\hfill {|{\widehat{\mathit{n}}}_1^{\prime }\times {\widehat{\mathit{n}}}_2^{″}|}^{2}sin\omega & ={\widehat{\mathit{n}}}_2^{\prime }·({\widehat{\mathit{n}}}_1^{\prime }\times {\widehat{\mathit{n}}}_2^{″}),\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill {|{\mathit{n}}_2^{″}-({\widehat{\mathit{n}}}_1^{\prime }·{\widehat{\mathit{n}}}_2^{″}){\widehat{\mathit{n}}}_1^{\prime }|}^{2}cos\omega & =[{\mathit{n}}_2^{\prime }-({\widehat{\mathit{n}}}_1^{\prime }·{\widehat{\mathit{n}}}_2^{\prime }){\widehat{\mathit{n}}}_1^{\prime }]·{\widehat{\mathit{n}}}_2^{″}.\hfill \end{array}$$(34)

Collectively, the rotation R₂R₁ maps the set of points $({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2)$ to $({\widehat{\mathit{n}}}_1^{\prime },{\widehat{\mathit{n}}}_2^{\prime })$. Given this rotation, we may use Eq. (12) to evaluate the bipolar harmonics ${\mathit{B}}_{\ell m}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1^{\prime },{\widehat{\mathit{n}}}_2^{\prime })$ by transforming ${\mathit{B}}_{\ell m}^{\left(j_1{\alpha }_1\right)\left(j_2{\alpha }_2\right)}({\widehat{\mathit{n}}}_1,{\widehat{\mathit{n}}}_2)$ instead of an explicit evaluation.

We demonstrate the transformation of sensitivity kernels in Fig. 2 by choosing the pairs of points ${\widehat{\mathit{n}}}_1=(0,0)$, ${\widehat{\mathit{n}}}_2=(\pi /4,0)$, ${\widehat{\mathit{n}}}_1^{\prime }=(\pi /2,0)$ and ${\widehat{\mathit{n}}}_2^{\prime }=(\pi /2,\pi /4)$. The solid line represents the kernels computed using Eq. (20) for the points ${\widehat{\mathit{n}}}_1^{\prime }$ and ${\widehat{\mathit{n}}}_2^{\prime }$, while the dots represent that computed by rotating the kernel components for ${\widehat{\mathit{n}}}_1$ and ${\widehat{\mathit{n}}}_2$ using Eq. (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.

Fig. 2. Radial profiles of the imaginary part of the kernel for (ℓ = 1, m = 0) and vector index γ = 1, computed directly for ${\widehat{\mathit{n}}}_1=(\pi /2,0)$ and ${\widehat{\mathit{n}}}_2=(\pi /2,\pi /4)$ using the equation (solid line) and by rotating the kernels for ${\widehat{\mathit{n}}}_1=(0,0)$ and ${\widehat{\mathit{n}}}_2=(\pi /4,0)$ (dots). The kernels are computed for line-of-sight projected measurements. Bottom panel: represents the absolute difference between the two results.

4.3. Evaluation time

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 ℓ_max = 30. The absolute difference using ℓ_max = 30 is 17.5 h of computation time. The timings were measured on 224 Intel Xeon CPU E5-2680 v4 cores running at 2.40 GHz.

Fig. 3. Computational time required to evaluate kernels Kℓm for all modes for 0 ≤ ℓ ≤ ℓmax and 0 ≤ m ≤ ℓmmax, where ℓmax is the angular degree cut-off. Top panel: the solid line indicates the time required using the radial component of the displacement, while the dashed line indicates the time required using the line-of-sight projected component of the wave velocity.

5. Discussion and conclusions

In this work, we improved upon the previous work by Bhattacharya (2020) to evaluate helioseismic sensitivity kernels for travel-time shifts computed using line-of-sight projected seismic wave velocities. The ability to transform between kernels for pairs of points that are related by a rotation may be helpful if one is interested in inverting for averaged travel times to improve the signal-to-noise ratio. The Julia implementation of the approach presented in this work is available freely under the MIT license⁴. The modular structure of the code ensures that parts of the entire project, such as vector spherical harmonics, bipolar spherical harmonics, or Wigner D-matrices, may be reused by other projects as necessary. The code developed is embarrassingly parallel in wave modes and temporal frequency and uses the message passing interface (MPI) for communication between processes, so it may be directly run on computing clusters to evaluate sensitivity kernels.

In this work, we did not focus on the impact of the differences in line-formation heights on the sensitivity kernels, although it is expected to play a part in observed travel-time differences as well (Zhao & Chen 2020). While this may also be incorporated into our approach, it requires a more careful modeling of the observation heights and the atmospheric boundary condition experienced by seismic waves. Further refinements in this direction might be possible. We also assumed that the wave excitation occurs at one specific radius, while in reality the excitation is perhaps spread over a radial range. A different model of source covariance following Gizon et al. (2017) might be used to obtain a better match to observations, although this would significantly increase the computational burden. Furthermore, the temporal spectrum of the source covariance was chosen arbitrarily in this work to be a Gaussian centered at 3.2 mHz, although a better profile may be obtained from the spectral envelope of the observed wave covariances.

Acknowledgments

This work was supported by NYUAD Institute Grant G1502 “NYUAD Center for Space Science”. This research was carried out on the High-Performance Computing resources at New York University Abu Dhabi.

References

Appendix A: Vector spherical harmonics

We numerically evaluated the vector spherical harmonics following Equation (24). We list the components of the VSH in the spherical basis, as detailed by Varshalovich et al. (1988):

$$\begin{array}{cc}\hfill {\left[{\mathbf{Y}}_{\mathit{jm}}^{j+1}\left(\widehat{\mathbf{n}}\right)\right]}^{+1}& ={\left[\frac{(j-m+1)(j-m+2)}{2(j+1)(2j+3)}\right]}^{\frac{1}{2}}{Y}_{j+1m-1}\left(\widehat{\mathbf{n}}\right),\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill {\left[{\mathbf{Y}}_{\mathit{jm}}^{j+1}\left(\widehat{\mathbf{n}}\right)\right]}^0& =-{\left[\frac{(j-m+1)(j+m+1)}{(j+1)(2j+3)}\right]}^{\frac{1}{2}}{Y}_{j+1m}\left(\widehat{\mathbf{n}}\right),\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill {\left[{\mathbf{Y}}_{\mathit{jm}}^{j+1}\left(\widehat{\mathbf{n}}\right)\right]}^{-1}& ={\left[\frac{(j+m+1)(j+m+2)}{2(j+1)(2j+3)}\right]}^{\frac{1}{2}}{Y}_{j+1m+1}\left(\widehat{\mathbf{n}}\right),\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill {\left[{\mathbf{Y}}_{\mathit{jm}}^j\left(\widehat{\mathbf{n}}\right)\right]}^{+1}& =-{\left[\frac{(j+m)(j-m+1)}{2j(j+1)}\right]}^{\frac{1}{2}}{Y}_{jm-1}\left(\widehat{\mathbf{n}}\right),\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill {\left[{\mathbf{Y}}_{\mathit{jm}}^j\left(\widehat{\mathbf{n}}\right)\right]}^0& =\frac{m}{\sqrt{j(j+1)}}{Y}_{\mathit{jm}}\left(\widehat{\mathbf{n}}\right),\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill {\left[{\mathbf{Y}}_{\mathit{jm}}^j\left(\widehat{\mathbf{n}}\right)\right]}^{-1}& ={\left[\frac{(j-m)(j+m+1)}{2j(j+1)}\right]}^{\frac{1}{2}}{Y}_{jm+1}\left(\widehat{\mathbf{n}}\right),\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill {\left[{\mathbf{Y}}_{\mathit{jm}}^{j-1}\left(\widehat{\mathbf{n}}\right)\right]}^{+1}& ={\left[\frac{(j+m)(j+m-1)}{2j(2j-1)}\right]}^{\frac{1}{2}}{Y}_{j-1m-1}\left(\widehat{\mathbf{n}}\right),\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill {\left[{\mathbf{Y}}_{\mathit{jm}}^{j-1}\left(\widehat{\mathbf{n}}\right)\right]}^0& ={\left[\frac{(j-m)(j+m)}{j(2j-1)}\right]}^{\frac{1}{2}}{Y}_{j-1m}\left(\widehat{\mathbf{n}}\right),\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill {\left[{\mathbf{Y}}_{\mathit{jm}}^{j-1}\left(\widehat{\mathbf{n}}\right)\right]}^{-1}& ={\left[\frac{(j-m)(j-m-1)}{2j(2j-1)}\right]}^{\frac{1}{2}}{Y}_{j-1m+1}\left(\widehat{\mathbf{n}}\right).\hfill \end{array}$$(A.8)

The spherical harmonics $Y_{\ell m}(\widehat{\mathbf{n}})$ are zero for indices (ℓ,m) that do not satisfy |m| ≤ ℓ.

All Figures

	Fig. 1. Radial profiles of spherical-harmonic components of the sensitivity kernel for the radial component (γ = 0) of flow velocity. The gray line is computed using radial components of wave velocity, whereas the black doted line is computed using the line-of-sight projected components. The title indicates the mode indices (ℓ,m).
In the text

Fig. 2. Radial profiles of the imaginary part of the kernel for (ℓ = 1, m = 0) and vector index γ = 1, computed directly for ${\widehat{\mathit{n}}}_1=(\pi /2,0)$ and ${\widehat{\mathit{n}}}_2=(\pi /2,\pi /4)$ using the equation (solid line) and by rotating the kernels for ${\widehat{\mathit{n}}}_1=(0,0)$ and ${\widehat{\mathit{n}}}_2=(\pi /4,0)$ (dots). The kernels are computed for line-of-sight projected measurements. Bottom panel: represents the absolute difference between the two results.

In the text

	Fig. 3. Computational time required to evaluate kernels Kℓm for all modes for 0 ≤ ℓ ≤ ℓmax and 0 ≤ m ≤ ℓmmax, where ℓmax is the angular degree cut-off. Top panel: the solid line indicates the time required using the radial component of the displacement, while the dashed line indicates the time required using the line-of-sight projected component of the wave velocity.
In the text