Helioseismic Inversion to Infer the Depth Profile of Solar Meridional Flow Using Spherical Born Kernels

Observing solar meridional circulation is challenging because of its small magnitude compared to other flows at the solar surface. Improvements in observational techniques have made it feasible to reliably infer the profile of meridional circulation at the surface and near-surface layers through a variety of techniques, e.g., Doppler-shift measurement (Duvall 1979; Hathaway et al. 1996; Ulrich 2010), tracking of surface features like small-scale magnetic remnants (Komm et al. 1993; Hathaway 2012), ring-diagram analysis (Schou & Bogart 1998; Basu et al. 1999; Haber et al. 2002; González Hernández et al. 2008; Basu & Antia 2010), and time–distance helioseismology (Giles et al. 1997). All these studies show that the magnitude of meridional circulation at mid-latitudes in the near-surface region is 10–20 m s⁻¹ at mid-latitudes, depending on the phase of the solar cycle, and the magnitude is directed poleward in both hemispheres. Mass conservation requires there to be an equatorward return flow below the solar surface. The location of the return flow plays a crucial role in some models of the solar dynamo. There have been several attempts to understand the internal structure of meridional circulation using helioseismic techniques, e.g., global helioseismology (Schad et al. 2013; Woodard et al. 2013), time–distance helioseismology (Zhao et al. 2013; Jackiewicz et al. 2015; Rajaguru & Antia 2015; Böning et al. 2017; Chen & Zhao 2017), etc. Such studies are very important for constraining dynamo models (Choudhuri et al. 1995; Dikpati et al. 2010) because meridional circulation is believed to carry magnetic fields as well as transport energy and angular momentum in the convection zone.

These studies have arrived at a variety of conclusions. From the analysis of two years of continuous data taken by the Solar Dynamics Observatory (SDO)/Helioseismic and Magnetic Imager (HMI), Zhao et al. (2013) reported multiple cells in depth with a return flow starting at about 0.9R_⊙ and the second cell below 0.82 R_⊙. Using 2 years of GONG data, Jackiewicz et al. (2015) obtained results that agreed with Zhao et al. (2013) about the shallow return flow but could not find the second cell in the deeper convection zone. Recently, Chen & Zhao (2017), using 7 years of SDO/HMI data covering 2010 May 1 to 2017 April 30, found a qualitatively similar profile as Zhao et al. (2013). These studies did not apply the mass conservation constraint. Using 4 years of HMI data, Rajaguru & Antia (2015) found a single-cell profile with a return flow below 0.77 R_⊙ using stream functions, which automatically ensure mass conservation. On the other hand, Schad et al. (2013) found multiple cells in both latitude and depth, which is very different from those found using time–distance helioseismology.

Systematics and noise significantly affect results obtained below 0.9 R_⊙ using time–distance helioseismology. All studies have considered data obtained from different instruments covering different periods of time. Sensitivity kernels computed using ray theory were used by Zhao et al. (2013), Jackiewicz et al. (2015), Rajaguru & Antia (2015), and Chen & Zhao (2017) to invert for flows. Ray kernels are sensitive only to perturbations along a ray path. If the length scale of the perturbation is smaller than the acoustic wavelength, results using ray kernels may not be reliable (Birch et al. 2001; Birch & Felder 2004).

The first-Born approximation, an alternative to ray theory, does not suffer from the above limitation. Recently, Mandal et al. (2017) computed travel-time sensitivity kernels in the Born limit in spherical geometry (two other alternative approaches were proposed independently by Böning et al. 2016; Gizon et al. 2017). Using Born kernels, Böning et al. (2017) inverted for meridional circulation and found a return flow at 0.9 R_⊙ using a SOLA (Pijpers & Thompson 1994) inversion technique. They found that both single and multiple-cell profiles are consistent with their measured travel times.

In this work, we use a stream-function approach, which automatically takes into account mass conservation, and compute relevant sensitivity kernels. The other advantage of using stream functions is that both radial and latitudinal components of the meridional flow may be simultaneously derived from it.

Helioseismic travel times capturing the signals due to meridional flows are estimated the same way as described in detail in Rajaguru & Antia (2015). The basic data used are the full-disk Doppler observations made by the HMI on board the SDO, and we have added two more years of data to that of Rajaguru & Antia (2015), covering a total length of six years (2010 May 1 through 2016 April 30). Each estimate of travel time is made from a Gabor wavelet (Kosovichev & Duvall 1997) fit to the monthly average of time–distance correlations computed in point-to-point deep-focus arc geometry employing 60 travel distances. Travel distances ranging between 216 and 4464 in steps of 072, covering depths from near the surface down to about 0.7 R_⊙. A major uncertainty in travel-time estimates is the center-to-limb systematics—CLS (Baldner & Schou 2012; Zhao et al. 2012, 2013; Rajaguru & Antia 2015), whose magnitude is several times the signal, due to meridional flows. We correct the north to south (N–S; or the meridional direction) travel times for this systematics in the same way as originally proposed by Zhao et al. (2012) and also as done by Rajaguru & Antia (2015): CLS is estimated from the antisymmetric component of west to east (W–E) travel times (Figure 1; the symmetric part corresponds to the solar rotation), which are then subtracted from the N–S travel times. These W–E travel times are estimated from the central W–E equatorial belt spanning 15° on either side of the equator. Recently, Chen & Zhao (2017) have proposed a new empirical method to remove the systematics. They measure travel-time shifts between two points located on the solar disk for many different azimuthal angles and skip distances, then solve a system of linear equations containing both center-to-limb systematics and travel-time shifts due to meridional flow. In addition, Chen & Zhao (2017) removed data pixels containing fields stronger than a threshold of 10 G, following a method suggested by Liang & Chou (2015). In our travel-time measurement process, we have not removed the oscillation signals over the strongly magnetized surface regions. It should be noted that surface magnetic fields have been shown to corrupt the meridional flow measurements (Liang & Chou 2015); although such influence of surface regions is not expected to majorly affect the overall inferences on the deep structure of meridional circulation inferred from 6 years of observations, a detailed analysis with and without the surface magnetic regions' removal is necessary to quantify the surface effects. However, detailed analyses of the character of center-to-limb systematics and its dependence on the frequency of acoustic waves (S. P. Rajaguru & H. M. Antia 2018, in preparation, and Chen & Zhao 2018), indicate that further careful study is needed to fully account for them in the current estimates of meridional flow travel times. Measured travel times induced by solar meridional flow from our analyses and corresponding errors have been shown in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An efficient approach to computing sensitivity kernels in spherical geometry was proposed by Mandal et al. (2017). They showed that Green's function can be expressed as (Equations (12) and (13) from Mandal et al. 2017)

$$\begin{eqnarray}&&{G}_{{rr}}({\boldsymbol{r}},{{\boldsymbol{r}}}_{s},\omega )\\ &&=\displaystyle \sum _{{\ell }}\displaystyle \frac{(2{\ell }+1)}{4\pi }{\alpha }_{{\ell }\omega }(r){P}_{{\ell }}(\cos (\hat{{\boldsymbol{r}}}\cdot {\hat{{\boldsymbol{r}}}}_{s})),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\boldsymbol{G}}}_{{hr}}({\boldsymbol{r}},{{\boldsymbol{r}}}_{s},\omega )\\ &&=\displaystyle \sum _{{\ell }}\displaystyle \frac{(2{\ell }+1)}{4\pi {\omega }^{2}{\rho }_{0}}{\beta }_{{\ell }\omega }(r){{\boldsymbol{\nabla }}}_{h}{P}_{{\ell }}(\cos (\hat{{\boldsymbol{r}}}\cdot {\hat{{\boldsymbol{r}}}}_{s})),\end{eqnarray} \tag{ 2 }$$

where G_rr(r, r_s, ω) and ${{\boldsymbol{G}}}_{{hr}}({\boldsymbol{r}},{{\boldsymbol{r}}}_{s},\omega )=({G}_{\theta r}({\boldsymbol{r}},{{\boldsymbol{r}}}_{s},\omega ),{G}_{\phi r}({\boldsymbol{r}},{{\boldsymbol{r}}}_{s},\omega ))$ are radial and tangential components of the displacement of a wave with temporal frequency ω, measured at r due to a point source placed at r_s, P_ℓ the Legendre polynomial of degree ℓ. The terms α_ℓω, β_ℓω are obtained by solving a coupled differential equation (Equation (10) in Mandal et al. 2017) using a finite-difference based scheme for each harmonic degree ℓ and frequency ω. We compute the pair $({\alpha }_{{\ell }\omega },{\beta }_{{\ell }\omega })$ once and use them to evaluate Green's function, which is used subsequently for obtaining kernels. The expression of sensitivity kernels for a stream function, K_ψ, can be written in terms of Green's function and its derivatives. The kernel computation time depends on the size of the (r, θ, ϕ) grid used in the analysis, the maximum harmonic degree for the computation of Green's function and resolution in frequency. In this work, we consider 1398 × 200 × 400 grid points spanning 0.7 R_⊙–1.0 R_⊙ in the radial direction and the entire horizontal domain. We choose the harmonic degree in the range 20 ≤ ℓ ≤ 100 to compute Green's function using Equations (1) and (2). We choose the frequency range 2.0–4.5 mHz, split uniformly into 1250 bins. It takes approximately 1.5 hr to compute a stream-function sensitivity kernel for one source-receiver distance using 200 processors on a computer cluster. It is indeed very expensive to compute sensitivity kernels for all 60 travel distances centered around 329 latitude points. Since we do not consider line-of-sight effects in this work, we need to compute sensitivity kernels only for 60 travel distances. We then translate the kernels at different latitude points because these kernels depend only on the source-receiver distance and not on their exact location. We have validated these sensitivity kernels by choosing simple flow profiles.

Flows differently alter the travel times of upstream and downstream propagating waves. The sensitivity kernels for flows relate them to the travel-time difference, thus

$$\begin{eqnarray}&&\delta \tau ={\int }_{\odot }{K}_{{v}_{r}}{v}_{r}d{\boldsymbol{r}}+{\int }_{\odot }{K}_{{v}_{\theta }}{v}_{\theta }d{\boldsymbol{r}}.\end{eqnarray} \tag{ 3 }$$

The contribution of the second integral is much larger than that of the first. The rising and falling part of the v_r component cancels out, whereas the contribution of two branches of the v_θ component add to each other. Furthermore, the magnitude of v_r is also smaller. Because of the small contribution from the v_r component, it is almost impossible to directly determine it from inversions. We use the stream function, ψ, instead of velocity for inversion, which automatically takes into account mass conservation and both components of velocity, v_r and v_θ, can be determined from this single function using the following relations:

$$\begin{eqnarray}&&\rho {v}_{r}=\displaystyle \frac{1}{r}\displaystyle \frac{\partial \psi }{\partial \theta }+\displaystyle \frac{\cos \theta }{r\sin \theta }\psi ,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\rho {v}_{\theta }=-\displaystyle \frac{\partial \psi }{\partial r}-\displaystyle \frac{\psi }{r},\end{eqnarray} \tag{ 5 }$$

where ρ is the density of the medium. Substituting Equations (4) and (5) into Equation (3), we obtain

$$\begin{eqnarray}&&\delta \tau ={\int }_{\odot }{K}_{\psi }\psi d{\boldsymbol{r}}.\end{eqnarray} \tag{ 6 }$$

To compute sensitivity kernels for stream function ψ, we use model S Christensen-Dalsgaard et al. (1996) as the background solar model (Mandal et al. 2017). These kernels are highly sensitive to surface layers and are especially sensitive to the base of the convection zone, making it difficult to determine the profile of meridional circulation at depth. Furthermore, the density varies by 6 orders of magnitude over the convection zone, while the velocity of meridional flow does not vary to that extent. To account for this variation we invert for the quantity ψ' = ψ/ρ.

The size of the problem will become very large and ill-posed if we perform inversions for ψ at all spatial points. B-spline basis functions for both radial and latitude coordinates were used by e.g., Rajaguru & Antia (2015) to represent the stream function. Recently, Bhattacharya et al. (2017) used B-spline basis functions successfully for synthetic inversions of supergranules. This motivates us to use B-spline functions to represent variations in the radial direction. In the latitudinal direction, we use a derivative of the Legendre polynomial for the inversion. The total number of parameters is reduced drastically in this approach:

$$\begin{eqnarray}&&\psi (r,\theta )=\displaystyle \sum _{i,{\ell }}{a}_{i{\ell }}{B}_{i}(r)\displaystyle \frac{{{dP}}_{{\ell }}(\cos \theta )}{d\theta },\end{eqnarray} \tag{ 7 }$$

where a_iℓ are the coefficients of expansion, which we need to determine. The reason for choosing Legendre polynomials instead of B-splines in the latitudinal direction is that fewer coefficients are needed. In addition, we assume hemispheric symmetry of the meridional circulation and consider only the derivative of even Legendre polynomials, which reduces the number of parameters in the inversion and presents the problem more clearly. In all inversions, we aim to minimize the misfit function, defined here as

$$\begin{eqnarray}&&{\chi }_{\mathrm{MF}}\\ &&=\,\displaystyle \sum _{i}{\left(\displaystyle \frac{\displaystyle \int {K}_{i}\psi (r,\theta ){rdrd}\theta -\delta {\tau }_{i}}{{\sigma }_{i}}\right)}^{2}+\mathrm{Regularization},\end{eqnarray} \tag{ 8 }$$

where i indicates a pair of observation points. The terms K_i and δτ_i denote the corresponding kernel and travel-time measurement. Despite the reduction in parameter space, we still need to apply regularization because travel-time measurements are subject to systematic and realization noise, which may be amplified in the inversion. We use second-derivative-smoothing both in r and θ:

$$\begin{eqnarray}\mathrm{Regularization} & = & \displaystyle \int \left[{\lambda }_{r}{\left(\displaystyle \frac{{\partial }^{2}\psi (r,\theta )}{\partial {r}^{2}}\right)}^{2}\right.\\ & & \left.+\ {\lambda }_{\theta }{\left(\displaystyle \frac{{\partial }^{2}\psi (r,\theta )}{\partial {\theta }^{2}}\right)}^{2}\right]{r}^{2}\sin \theta d\theta {dr},\end{eqnarray} \tag{ 9 }$$

where λ_r and λ_θ are regularization parameters that can be tuned to obtain a smooth solution. We apply the constraint ψ = 0 at the upper boundary. We use two different methods based on the regularized least squares technique to obtain unknown coefficients in the Equation (7) by minimizing the misfit function.

4.1. First Method

After substituting Equation (7) into Equation (8) and taking derivatives with respect to unknown coefficient a_iℓ, we obtain a system of equations that is written in matrix form,

$$\begin{eqnarray}&&({M}^{T}{\rm{\Lambda }}M+{M}_{\mathrm{Reg}})a={M}^{T}{\rm{\Lambda }}\tau ,\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{M}_{{ij}}=\int {K}_{i}(r,\theta ){B}_{k}(r)\displaystyle \frac{{{dP}}_{{\ell }}(\cos \theta )}{d\theta }{rdrd}\theta ,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}{\rm{\Lambda }}=\left(\begin{array}{cccc}\displaystyle \frac{1}{{\sigma }_{1}^{2}} & 0 & 0 & ..\\ 0 & \displaystyle \frac{1}{{\sigma }_{2}^{2}} & 0 & ..\\ .. & .. & .. & ..\\ .. & .. & 0 & \displaystyle \frac{1}{{\sigma }_{n}^{2}}\end{array}\right),\end{eqnarray} \tag{ 12 }$$

and M_Reg is the Regularization matrix, obtained from Equation (9). Each pair of (k, ℓ) determines a unique index j in Equation (11). The column vector a is composed of expansion coefficients {a_iℓ} of Equation (7) and is obtained by solving Equation (10).

$$\begin{eqnarray}&&a={({M}^{T}{\rm{\Lambda }}M+{M}_{\mathrm{Reg}})}^{-1}{M}^{T}{\rm{\Lambda }}\tau .\end{eqnarray} \tag{ 13 }$$

The inverse of the matrix in Equation (13), which may be obtained using many methods, is computed here using singular value decomposition.

4.2. Second Method

A regularization misfit from Equation (8) can be satisfied if we solve the following system of equations to obtain the unknown coefficients in Equation (7):

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\sigma }_{i}}{\int }_{\odot }{K}_{i}(r,\theta )\psi (r,\theta )r\,{drd}\theta =\displaystyle \frac{1}{{\sigma }_{i}}\delta {\tau }_{i},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\lambda }_{r}\displaystyle \frac{{\partial }^{2}\psi ({r}_{k},{\theta }_{m})}{\partial {r}^{2}}=0,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\lambda }_{\theta }\displaystyle \frac{{\partial }^{2}\psi ({r}_{k},{\theta }_{m})}{\partial {\theta }^{2}}=0.\end{eqnarray} \tag{ 16 }$$

$({r}_{k},{\theta }_{m})$ is one of many points in the grid onto which we place constraints (15) and (16) in order to obtain a smooth solution. The condition number of the matrix for the second method is better than that of the first method.

Systematic errors can result from using finite numbers of Legendre polynomials and B-spline knots to expand the stream function. The chance of such an error decreases with increasing numbers of knots for a B-spline and degree of Legendre polynomials. However, the tradeoff is that the condition number of the problem becomes correspondingly larger. We define the misfit function for this case as

$$\begin{eqnarray*}{\chi }_{\mathrm{MF}} & = & \displaystyle \sum _{i}{\left(\displaystyle \int {K}_{i}\psi (r,\theta ){rdrd}\theta -\delta {\tau }_{i}\right)}^{2}\\ & & +\ \mathrm{Regularization}.\end{eqnarray*}$$

In our work, we consider 40 knots in the radial direction and all even values for ℓ with the highest degree 16 for Legendre polynomials in Equation (7). In order to validate our inversion algorithm, we first consider a few profiles, e.g., single cell or double-cell, and estimate the travel time corresponding to that profile by forward-modeling. We perform inversions with these travel-time measurements and compare the retrieved profile with the original one.

After ensuring that systematic errors do not affect the inversion, we add realization noise into the travel-time measurements (obtained in the previous section). We assume that noise is Gaussian and randomly perturb the travel time in proportion with errors found in observations. We then perform inversions with these travel-time data and compare the retrieved profile with the original one. In this case, we consider the misfit function (8). We are only going to present inversion results for noisy synthetic data in this paper. As expected, inverted profiles obtained using noise-free measurements agree better than those obtained using noisy travel times.

In order to get reliable noise estimates, we follow the approach of Rajaguru & Antia (2015). We randomly perturb the travel-time according to the error in the observation. We estimate the velocity and the associated standard deviation of these values, thereby propagating measurement to inferential errors.

We test our inversion algorithm with artificial single-cell and double-cell profiles to verify that we can retrieve both profiles, even in the presence of the observational error. The velocity amplitudes of both the cells are 20 m s⁻¹ at the surface, close to what was observationally found. The single-cell profile extends down to 0.75 R_⊙ before it changes signs. For the double-cell profile, the first cell extends down to 0.85 R_⊙ and the second one ends at 0.70 R_⊙. We generate artificial travel-time data with these two cells and add random errors (see Figure 3 for example). We then choose smoothing parameters in order to obtain a smooth solution from noisy travel-time data. The results for a single-cell profile are shown in Figure 4. We have shown both input and inverted profiles together for comparison. It is seen that our inversion is able to recover the input profile fairly well. The result for the double-cell profile is shown in Figure 5. Again, we find good agreement between inverted and input profiles. In both cases, we are also able to recover the radial component of the velocity, v_r, well. The depth dependence of the velocity profiles, v_θ and v_r, averaged over latitude is shown in Figure 6 for single-cell and double-cell cases. We can see that our inversion accurately recovers the depth profile of v_θ. In Figure 3, we compare the input travel time with that obtained from the inverted profile by forward-modeling. All these results give us the confidence to proceed further and perform inversions with observed travel-time measurement data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use travel-time measurements of meridional circulation obtained from analyses of 6 years of observational data by SDO/HMI in the periods between 2010 and 2016. As mentioned above, we perform inversions to determine the stream function ψ and evaluate components of flow, v_r and v_θ, using Equations (4) and (5). We find results similar to the results from using the two methods described in Section 4. Inversion results are shown in Figures 7 and 8. The results are qualitatively similar to those obtained by Rajaguru & Antia (2015) using ray theory. We plot depth variations of horizontal and radial flows, v_θ and v_r, averaged over the latitude range described in the Figure 8. We find that the return flow is below 0.78 R_⊙ at all latitudes (less than 20^◦; see Figure 7). Similar to Rajaguru & Antia (2015), we also see a sign change in v_θ at low latitudes below 0.9 R_⊙. However, this is not significant when considering the size of the error bar (see Figure 8). We are unable to completely rule out the picture of multiple-cell meridional circulation because of this. At higher latitudes, it clearly represents a single-cell profile. The inverted profile for radial velocity v_r (right panel of Figure 7) suggests that flow is directed outward at lower latitudes and inward at higher latitudes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because meridional circulation is important for understanding solar dynamo processes, with different techniques there have been several attempts to infer meridional circulation. Due to systematic errors and associated small magnitudes, the inferred models of this circulation tend to vary widely from one study to another (Zhao et al. 2013; Jackiewicz et al. 2015; Rajaguru & Antia 2015; Böning et al. 2017; Chen & Zhao 2017). In this work, we introduce a combination of techniques to improve the accuracy of these inferences and the condition number of the inverse problem: (1) to ensure mass conservation we use a stream function, (2) we compute sensitivity kernels for stream functions in the first-order Born approximation, (3) we project the solution on a cubic B-spline, Legendre-polynomial derivative basis to reduce the number of free parameters, and (4) we assume hemispheric symmetry and only consider even Legendre polynomials. In order to validate this method, we test it on single-cell and double-cell models of the meridional circulation. The surface-velocity amplitudes of these profiles are chosen to be close to the observational value. We compute synthetic travel times by integrating Born kernels against the stream function and add randomly generated noise in proportion to observational errors.

Upon successful validation, we apply the method to infer meridional circulation from travel-time measurements obtained from 6 years of SDO/HMI observational data. We find qualitatively similar results to those of Rajaguru & Antia (2015). Our inversion results suggest a single-cell profile covering ±60° in latitude and up to 0.7 R_⊙ in depth in the radial direction. We find a sign change in the horizontal component of the velocity, v_θ beneath 0.9 R_⊙, within ±20° in latitude, but it is inconclusive considering the size of the error bar. We find an equatorward return flow below 0.78 R_⊙.

Our analysis period almost exactly overlaps with the analysis period of Chen & Zhao (2017). Differences in the two inversion results might be due to the following factors. We have not removed surface magnetic regions, in order to obtain the travel time induced by solar meridional flow. Liang & Chou (2015) demonstrated that not removing surface magnetic regions affects the travel-time measurements. We have accounted for local mass conservation and also employed Born kernels, whereas Chen & Zhao (2017) have not considered mass conservation and used ray kernels.

Braun & Birch (2008) have shown that in order to infer the depth profile of a weak signal like the solar meridional flow in the entire convection zone, many years of data are needed. However, during this period, the flow might change as a function of the phase of the solar activity cycle. In that case, we can only obtain an average meridional flow profile for the time period of the analysis. Improvements such as including a full covariance matrix, as done by Böning et al. (2017), have been suggested and can account for line-of-sight effects (the latter is only possible in a wave-theoretic approach); their viability is reserved for our future studies.

K.M. acknowledges the financial support provided by the Department of Atomic Energy, India. S.M.H. acknowledges support from Ramanujan Fellowship SB/S2/RJN-73 and the Max-Planck Partner group program.