Inversions for Deep Solar Meridional Flow Using Spherical Born Kernels
Vincent G. A. Böning Markus Roth Jason Jackiewicz Shukur Kholikov
astro-ph.SR













Kiepenheuer-Institut für Sonnenphysik, 79104 Freiburg, Germany


Kiepenheuer-Institut für Sonnenphysik, 79104 Freiburg, Germany

New Mexico State University, Las Cruces, NM 88001, USA

National Solar Observatory, Tucson, AZ 85719, USA









The solar meridional flow is a crucial ingredient in modern dynamo theory. Seismic estimates of this flow have, however, been contradictory in deeper layers below about 0.9 R_. 

Results from time-distance helioseismology have so far been obtained using the ray approximation. Here, we perform inversions using the Born approximation.


The initial result is similar to the result previously obtained by <cit.> using ray kernels while using the same set of GONG data and the SOLA inversion technique. However, we show that the assumption of uncorrelated measurements used in earlier studies may lead to inversion errors being underestimated by a factor of about two to four. 

In a second step, refined inversions are performed using the full covariance matrix and a regularization for cross-talk. 

As the results are found to depend on the threshold used in the singular value decomposition, they were obtained for a medium threshold (10^-7-10^-5, about 50% of the values used) and a threshold lower by a factor of 10 (about 70% of the values used).

The result obtained with the medium threshold is again similar to the original, with less latitudinal variation. However, using the lower threshold, the inverted flow in the southern hemisphere shows two or three cells stacked radially depending on the associated radial flows. 

Both the single-cell and the multi-cell profiles are consistent with the measured travel times. 

All our results confirm a shallow return flow at about 0.9 R_.






















§ INTRODUCTION
                                


Inferring the structure of the meridional flow in the deep solar interior has attracted considerable attention in recent years. While a conclusion on the flows in near-surface regions seems to have been reached <cit.>, the most recent measurements of the deep meridional flow <cit.> give seemingly contradictory results in deeper layers below about 0.9 R_, favoring a single-cell or multi-cell picture of the flow as summarized in the introduction of <cit.>.

Possible reasons for this discrepancy may include systematic effects like a center-to-limb effect in time-distance helioseismology <cit.>, perturbation of solar mode eigenfunctions by convection <cit.>, systematic effects introduced by magnetic fields <cit.>, B- or P-angle variations <cit.>, as well as differences in the instruments used and the time period considered.

In addition, there are several efforts underway to develop or validate new methods for inferring the meridional flow, in particular using local helioseismic techniques <cit.>.

As inversions for the meridional flow with time-distance helioseismology <cit.> have so far been modeled using the rather classical ray approximation <cit.>, the Born approximation has been brought forward very recently as an alternative <cit.>. Instead of assuming that travel times of acoustic waves can be perturbed by a flow along a ray path as in the ray approximation, the Born approximation <cit.> assumes that the whole wave field is scattered by the flow. Therefore, the travel time of a wave packet can be perturbed by a flow distant from the ray path. The Born approximation has been well tested and validated in Cartesian geometry for inferring small-scale flows <cit.>. As it is generally thought of as a more accurate model of the physics in the solar interior <cit.>, it is a method worth exploring for inferring large-scale flows such as the meridional flow.

Born approximation sensitivity functions (kernels) have very recently been validated in time-distance helioseismology of the meridional flow <cit.>. In this study, we will perform inversions for the deep meridional flow with these spherical Born kernels. As was shown by <cit.>, phase-speed filtered measurements seem particularly useful for this endeavor. We will therefore use the phase-speed filtered travel-time measurements obtained by <cit.> in this work. These measurements have been inverted for the meridional flow by <cit.> using the SOLA method <cit.>.

The first objective of this work is to answer the question whether the inversion results inferred by <cit.> using ray kernels can be confirmed using Born kernels or whether a different conclusion may be reached. Our second objective is to study in detail different sources of systematic errors in the inversion process, such as the error propagation and the cross-talk of the radial flow into the inversion for the horizontal flow component. For this work, we will employ an analytic formula for the covariance of travel-time measurements (, ). We will also use different strategies for analyzing the cross-talk and other systematics in SOLA inversions using Born kernels that were also employed by various authors in Cartesian inversions for small-scale flows <cit.>.

The paper is organized as follows. Section <ref> provides a short introduction to the data set and the SOLA inversion technique used. Section <ref> gives a summary of the computation of spherical Born kernels, which are used to model the travel-time measurements in this work. An analytic formula from the literature is used in Section <ref> to obtain a model for the covariance matrix of the measurements. A comparison of an initial inversion for the meridional flow using Born kernels to the result obtained by <cit.> is presented in Section <ref>. This inversion is performed under the assumption of uncorrelated measurements (i.e., a diagonal covariance matrix) and without including a regularization for cross-talk as it was done in previous inversions for the meridional flow using time-distance helioseismology <cit.>. In this section, we also provide a detailed study on the error propagation as well as the cross-talk between flow components. This information is used in Section <ref> to perform refined inversions for the meridional flow that include the full covariance matrix and a regularization for cross-talk. The results of this study are discussed in Section <ref> and conclusions are presented in Section <ref>.





§ DATA AND INVERSION TECHNIQUE



In this paper, we use the travel-time measurements obtained and described by <cit.>. The travel times were obtained using the GONG instrument and the data of 652 days between 2004 and 2012, where the duty cycle was high and the B-angle not too large. This period includes the declining phase of cycle 23 and the rising phase of cycle 24. An average sunspot number of 36 for the days used <cit.> indicates that the data are taken from times of moderate to low activity. North minus south (N – S) point-to-arc travel times were obtained for a total of 72 travel distances (Δ) and 384 latitude bins (λ). In addition, east minus west (E – W) travel times were measured in order to correct for a systematic center-to-limb effect detected by earlier studies (e.g., , ). The N – S travel times are corrected for this systematic effect by subtracting the E – W measurements.

The travel times were obtained using a Gabor fit <cit.>. <cit.> inverted this set of travel times for the horizontal component of the meridional flow using a standard SOLA technique <cit.>. This inversion technique will be the starting point of the work performed in this study. Our goal is to find the meridional flow with horizontal component v_θ(r,θ) and radial component v_r(r,θ) that satisfies for all travel-time measurements δτ_i

    δτ_i =    ∫_0^R_∫_0^π ∑_k=r,θ K_k^i(r,θ)  v_k(r,θ)   r sin(θ) 𝕀θ 𝕀 r  
        +  ϵ_i,

where K_k^i is the travel-time sensitivity function of the measurement i=(Δ,λ) to a flow in the direction k∈{r,θ}, and ϵ_i are measurement errors. Two-dimensional integrals are evaluated over the whole radial and latitudinal domains throughout this work. The meridional flow is sought at a series of target locations (r_T,θ_T) as a linear combination of the measured travel times, e.g., for the horizontal flow component,

    v_θ^inv(r_T,θ_T) = ∑_i w_i(r_T,θ_T)  δτ_i.

The aim is to construct averaging kernels _θ,

    _θ(r,θ;r_T,θ_T) = ∑_i K_θ^i(r,θ)w_i(r_T,θ_T),

which are sufficiently near a Gaussian target function T(r,θ;r_T,θ_T) with a certain radial and horizontal full width at half maximum (_r,_θ) centered at the target location. This is achieved by minimizing the misfit given by

    = ∬| _θ - T | ^2   r sin(θ) 𝕀θ 𝕀 r.

The minimization is performed with an additional regularization parameter μ that balances misfit and errors,

    ^2 =∑_i,j w_i(r_T,θ_T) Λ_i^j w_j(r_T,θ_T),

where the covariance matrix Λ is given by

    Λ_i^j = [ ϵ_i , ϵ_j].

The final cost function

    χ(w_i   (r_T,θ_T);μ)  =  + μ ^2

is then minimized at every target location subject to the constraint

    ∬_θ(r,θ;r_T,θ_T)   r sin(θ) 𝕀θ 𝕀 r = 1.

The resulting weights can be obtained using a matrix inversion <cit.>.




§ FORWARD-MODELING USING BORN KERNELS



The forward-modeling of the travel-time measurements in <cit.> was performed using ray kernels <cit.>. Instead, we will use the Born approximation <cit.> to model the effect of the flow on the travel times in this work.

While the ray approximation assumes that the acoustic waves in the solar interior are sensitive to a flow field only along a certain ray path, the Born approximation models the scattering of the full wave field due to advection in first order <cit.>. The method used in this work for computing Born kernels was developed by <cit.> and further refined by <cit.>. 

As the computation of the sensitivity kernels depends on an accurate match between model and data power spectra, the following free parameters in the model were adjusted <cit.>. 
For harmonic degree l<100, mode frequencies and damping rates from the GONG ftp site[<ftp://gong2.nso.edu/TSERIES/v1f/>] were used, and for l≥ 100, they were provided by Sylvain Korzennik <cit.>. 

For small harmonic degrees, the fitted mode widths are much smaller than the frequency resolution of a day-long time series. Therefore, small mode widths were increased to a minimum value 40% larger than the resolution of a day-long time series. This value was empirically found to produce a good match between data and kernel power spectra. As in <cit.>, the model power spectrum was further corrected by an l-dependent factor, which may be seen as an optical transfer function correcting instrumental effects; see <cit.> and <cit.>. In addition, the source correlation time, a free parameter that models the sources in the model, is adjusted to obtain a match between the mean frequency of the data and model power spectra. Finally, the phase-speed filters used in the data analysis procedure of <cit.> were applied to the model power spectrum.

As a result, zero-order model power spectra and cross-covariances agree well with those obtained from GONG data; see Figure <ref>. Only leakage effects visible in the GONG power spectrum at low harmonic degree (top-left panel) are not included in our model. However, no similar effect is visible in the cross-covariances (top middle panel), which are used to model the travel-time shift.

Inspection of kernel plots (right column in Figure <ref>) also reveals similar results to <cit.> with some changes introduced by differences in the filters and data power spectra.

For each of the 72 travel distances, one 3D kernel was computed. This kernel was then reprojected to the different latitudes and rotated to obtain the N – S point-to-arc geometry with 30 wide arcs used in the data analysis procedure. After integrating the kernels along the azimuthal domain, they were smoothed and rebinned with respect to distance in a similar way to that in the data analysis with a final number of 45 distances. Finally, both the travel times and the kernels were further rebinned by a factor of two to obtain 192 latitude bins, which are used in the inversions performed in this work.








§ FULL COVARIANCE MATRIX



In previous inversions for the deep meridional flow using time-distance helioseismology <cit.>, the travel-time measurements were assumed to be uncorrelated, which is equivalent to a diagonal covariance matrix. In order to study the effect of this assumption, we will use both a diagonal and the full covariance matrix in this work.

Analytic formulas for computing the covariance of travel-time measurements have been proposed by <cit.> and <cit.>. They were obtained on the basis of an empirical model for the noise in the power spectrum, which assumes uncorrelated noise for different frequency bins, and were shown to be in good agreement with data and Monte Carlo simulations.

For the purpose of this work, we compute the covariance matrix using formula (B.6) or equivalently (13) from <cit.>. As suggested in <cit.>, for longer time series (one day in our case), only the first term in the formula has to be taken into account. As an ingredient, only the zero-order or mean cross-covariances for all combinations of measurements are needed. In this work, we use the mean power spectrum of the data to obtain mean cross-covariances for every possible distance needed in the formula. For computational reasons, we compute the covariance for point-to-point measurements only. 

The resulting diagonal entries of the covariance matrix should in principle be equal to the squared errors of the travel-time measurements used in this study. For two reasons this is not the case. Firstly, the analytic formula for the diagonal entries only depends on the travel distance and not on the latitude, in contrast to the errors of the measurements, which increase toward the limb. Second, the analytic values have a different scale compared to the squared errors as the measurements have been further averaged, e.g. over a period of 652 days. This is taken into account by renormalizing the result from the analytic formula using

    Λ̃_i^j =  Λ_i^j σ_i  σ_j /√(Λ_i^i Λ_j^j),

where σ_i is the error of measurement i.

A cut through the covariance matrix obtained is displayed in Figure <ref>, which is in general agreement with the covariance matrix obtained by <cit.> from unfiltered data.






§ HORIZONTAL SOLA INVERSION: A FIRST COMPARISON




The aim of this section is to perform a first comparison of an inversion for the meridional flow using spherical Born kernels to an earlier inversion done with ray kernels. For this purpose, the inversion method is kept as close as possible to the one applied in <cit.>. Specifically, we only invert for the horizontal component of the flow, we assume the covariance matrix to be diagonal, and we do not take the cross-talk from the radial flow component into the horizontal flow component into account when obtaining the inversion weights.

Our intention is to first answer the question as to whether and how the inversion result is affected by simply exchanging ray kernels with Born kernels. Second, we will compare the magnitude of the errors propagated with a diagonal and with the full covariance matrix, and we will analyze the magnitude of the cross-talk.





 §.§ Choice of Regularization Parameters


The inversion procedure involves a number of free parameters, such as the full widths of the Gaussian target functions and the parameter μ, which may vary with target location. We therefore first performed a number of test computations in order to get an overview over the parameter space involved. These tests showed that a scaling of the full widths of the target functions with the sound speed, as proposed in <cit.>, is a reasonable choice. Near the surface (r_T≥0.95 R_), this scaling would imply very small target widths. We therefore set minimum target widths of _r,min = 0.03 R_ and _θ,min = 5. This limitation is justified by the minimum scales of the kernels involved in this study as the minimum travel distance is around three degrees and the kernels extend a little farther than that. For the scaling of the full widths, we find values at the bottom of the convection zone at r_T=0.7 R_ of _r = 0.09 R_ and _θ = 20. At first sight, these values may seem very large, especially when compared with the averaging kernels obtained by <cit.> using the ray approximation, which are much more localized. However, taking a closer look at the spatial scales involved in the kernels shown in Figure <ref>, which are much wider than ray kernels, the values obtained are plausible. Our test computations showed that, if the widths of the target functions are further decreased significantly, it is hard to obtain reasonable values for errors and misfit.


In addition, every inverse matrix in the inversion problem can be computed using different thresholds for the singular values (SVs). In practice, the inversion results do not show a large dependence on the choice of threshold in this inversion, where we use a diagonal covariance matrix.

In order to obtain an optimal choice for the remaining inversion parameter μ, which controls the relation between errors and misfit in the inversion, we first do a number of inversions for a series of values for μ on a coarse grid of target locations. Finally, we choose an optimal value for μ for every target depth, using maximum threshold values for the misfit and the errors of the inverted flows. In practice, we find that it is often hard to obtain a reasonable inversion error without compromising the fit of the averaging kernel to the target kernel. We thus chose a maximum error of 1, similar to <cit.>, and a maximum misfit of 0.2, when the misfit is normalized as 

    _norm  = /∬  T ^2   r sin(θ) 𝕀θ 𝕀 r.

We note here that it is possible to achieve much better values for the misfit, up to almost perfect agreement between the averaging and target kernels, but at an unacceptable expense in the errors. A maximum misfit of 0.2 is found to be about the largest possible value for obtaining an acceptable match between the averaging and target kernels. At the surface, this condition is relaxed by a factor of two, as we find that most of the misfit comes from locations just being a bit farther off, or a difficulty of the averaging kernel to achieve a Gaussian form. If no value for μ is found for the maximum error and misfit given, these conditions are relaxed step by step until possible inversion parameters are found. If several possible points are found, an optimal choice is made using the so-called L-curves <cit.>. In the optimizing procedure, we consider a small number of target latitudes within 20 from the equator.



Using the coarse grid of target locations, optimal inversion parameters are chosen for every target depth of the spatial grid. For the final inversion, a finer target grid is defined, and for each location, the inversion parameters are interpolated between the values from the preparatory inversion.






 §.§ Inversion Results


Inversion results are presented in Figure <ref>, where the inverted flow (left panel), the inversion errors (middle panel), and the misfit (right panel) are shown. We first note that the resulting meridional flow is similar to the inversion result obtained by <cit.> using ray kernels, at least qualitatively.


As this inversion result and the result obtained by <cit.> are obtained using a diagonal covariance matrix, it is questionable whether the error estimate is correct. Using the inversion weights obtained in this inversion and the full covariance matrix computed as in Section <ref>, it is possible to give more accurate estimates of the errors of the inverted flow and to estimate the impact of not taking the full covariance matrix into account.








In Figure <ref>, inversion errors using the diagonal covariance matrix (lines near the bottom) can be compared to errors obtained using the full covariance matrix (lines near the top) as a function of target depth for a series of target latitudes. Both errors are obtained using the same inversion weights. It can be seen that the errors from the diagonal covariance underestimate the errors from the full covariance by a factor of about two to four. For this comparison, the values on the diagonal of the two covariance matrices are identical prior to the final rebinning of the measurements by a factor of two in latitude; see Section <ref>. If they were to be set equal after this final rebinning, the errors from the diagonal covariance would increase on average by about 35 % and they would still underestimate the errors from the full covariance by a factor of 1.5 - 3.




Furthermore, although the cross-talk was not considered in the inversion procedure, it is possible to compute cross-talk averaging kernels, i.e., kernels for the influence of the radial flow on the inversion for horizontal flow,

    _r(r,θ;r_T,θ_T) = ∑_i K_r^i(r,θ)w_i(r_T,θ_T).

For a series of target depths, we show the averaging kernel for the horizontal flow (right subpanel) and the cross-talk kernel for the radial into the horizontal component (left subpanel) in Figure <ref>. We note that for each target depth plotted, both kernels are shown using the maximum values of the averaging kernel as a scale, with the maximum values of the cross-talk averaging kernels being 120, 100, 70, and 14 times larger for the given target depths r_T of 0.7, 0.8, 0.9, and 0.98 R_. These values are much larger than the ones obtained by <cit.> in an inversion for subsurface flows in Cartesian geometry.














We therefore study the nature and impact of the cross-talk in more detail, see Figure <ref>, where the normalized cross-talk,

    _norm  = ∬_r^2   r sin(θ) 𝕀θ 𝕀 r/∬  T ^2   r sin(θ) 𝕀θ 𝕀 r

is displayed as a function of depth, for a number of latitude bins. For example, a normalized cross-talk of 100 means that radial flows may contribute to the inversion for horizontal flow even if they are about √(100)=10 times smaller than the actual horizontal flow. In principle, even a small magnitude radial flow may thus leak into the inversion for the horizontal flow in such a case. In order to answer the question whether such values for the cross-talk introduce a large contribution of the radial flow to the inverted horizontal flow, we convolved the cross-talk averaging kernels with two choices of exemplary radial flow fields v_m^orig  (m=r),

    v_θ   ^conv(m)(r_T,θ_T)  
        =∬_m(r,θ;r_T,θ_T)   v_m^orig(r,θ)   r sinθ 𝕀θ 𝕀 r.

The resulting contribution of the radial flow signal to the inversion for v_θ is displayed in Figure <ref>. For the left panel, we chose v_r^orig to be the radial flow component of the single-cell meridional flow profile employed in the simulation of <cit.>, divided by a factor of 36 in order to mimic a realistic magnitude of a solar-type flow. For the right panel, we chose the radial flow component that was obtained by <cit.> from the inverted horizontal flow by applying mass conservation. In both cases, the contribution of the radial flow component to the inversion for v_θ has maximal values of around 2, with typical values of about 1 - 1.5. We may thus expect small contributions from the radial flow component, but it is likely that they do not alter the inversion result at a large scale. However, one should be aware of this contribution.













§ HORIZONTAL SOLA INVERSION INCLUDING FULL COVARIANCE AND CROSS-TALK



Refined inversions are now performed by including the full covariance matrix and a regularization term for cross-talk into the inversion problem. A similar problem has been studied for inversions for 3D  flows near the surface in Cartesian geometry by, e.g., <cit.> and <cit.>.




 §.§ The Full Inverse Problem


In the full SOLA inverse problem, we invert for a flow component k∈{r,θ},

    v_k^inv(r_T,θ_T) = ∑_i w_i^k(r_T,θ_T)  δτ_i,

by trying to match the averaging kernels,

    _m^k(r,θ;r_T,θ_T) = ∑_i K_m^i(r,θ)w_i^k(r_T,θ_T),

to a target kernel T_k^m = δ_m,k T, where T is the same target function from above. When inverting for a flow component k, we thus also intend to match the averaging kernel of the opposite component m≠ k to zero. Therefore, a regularization for the cross-talk of the flow component m≠ k into k,

    _k = ∬(_m≠ k^k)^2    r sin(θ) 𝕀θ 𝕀 r,

is added to the inversion problem outlined in Section <ref> (see, e.g, <cit.> for the equivalent formulation in Cartesian geometry). The cost function thus becomes

    χ_k(w_i^k;μ)    =  + μ ^2 + ν _k

which is to be minimized subject to the constraints

    ∬_m^k   r sin(θ) 𝕀θ 𝕀 r = δ_km.



As in the inversion without cross-talk, inversion parameters are chosen in a preparatory inversion using a coarser spatial grid for the target locations. The FWHMs of the target functions are the same.


In the following, we invert for the horizontal component (k=θ) and show the inversion results for parameters obtained using different strategies. In all cases, we first set an upper limit to the misfit of 0.2, which is the same as in the inversion without cross-talk in the previous section.



As the inversion results presented here are found to depend on the threshold used in the singular value decomposition (SVD), we show the SVs of a matrix used in the inversion at an example target depth of r_T=0.8 R_ in the left panel of Figure <ref>. As the distribution of SVs depends on target depth, the thresholds at each target depth are chosen relative to the point of highest negative curvature, which is found by fitting a high-order polynomial to the curve. Therefore, the thresholds and the fraction of SVs used in the matrix inversion also depend on target depth; see middle and right panel of Figure <ref>.









 §.§ Results: Controlling Misfit and Errors





In addition to the misfit, we aim to keep the errors under control in the first inversion. After choosing a threshold for the SVs, we also set a maximum threshold of 1.0 m s^-1 for the inversion errors as in Section <ref>. We then search the parameter space for the best cross-talk available. 

As a first case, we consider an SV threshold chosen to be just above the end of the rather flat plateau of SVs in the left panel of Figure <ref> (indicated in Figure <ref> by “medium SV”). Here, about half of the SVs are used in the matrix inversion.

In this case, we obtain an inverted flow profile that is very similar to the one obtained by <cit.>; see the left panel in Figure <ref>, subsequently termed “case 1”. Compared to the initial result shown in Figure <ref>, however, it shows less fluctuations as a function of latitude in each hemisphere. This difference is found to be due to the use of the full covariance matrix.






In a second case, we choose an SV threshold smaller by a factor of 10 just at the edge of the plateau of SVs in the left panel of Figure <ref> indicated by “low SV”. Here, about 70 % of SVs are used. In this case, the inverted flow profile shows some noteworthy features; see case 2 in Figure <ref>. 
Most notably, there is an additional extended poleward flow branch visible at about 0.8 R_ in the southern hemisphere. This flow component gives rise to a multi-cell structure of the flow stacked radially. Depending on the associated radial flow, the inverted meridional flow corresponds to a double-cell profile <cit.> or to three flow cells stacked radially <cit.>. At the same time, no additional flow cells are visible in the northern hemisphere. At the central location of the additional poleward flow structure in the southern hemisphere at about 0.8 R_ and 30 latitude, a patch with equatorward flow is also visible in the result for case 1, in the initial result in Figure <ref>, and in the result obtained by <cit.> using ray kernels.

Inversion errors using the full and a diagonal covariance matrix are shown in Figure <ref>. With the full covariance matrix, we now obtain errors that are generally smaller than when using a diagonal covariance matrix by about a factor of 1.5 in case 2, with a few exceptions. In the case of the medium SV threshold (case 1), however, the errors from the full covariance matrix are on average about twice as large as the ones from the diagonal covariance matrix, similar to the initial inversion result (see Figure <ref>). Note that the errors from the diagonal covariance matrix increase on average by about 35% if the diagonal values of the covariance matrices are assumed to be identical after rebinning in latitude as in Section <ref>.







Furthermore, one can see that the cross-talk decreased (see Figure <ref>, case 1 and 2), by about a factor of 3 especially at larger target depths compared to the cross-talk for the initial result; see Figure <ref>. As a consequence, the convolution of the two examples of radial flow profiles discussed in Section <ref> with the radial averaging kernels decreases by over 30% to below 1.4.

Concerning the cross-talk, it thus seems that the horizontal inversion without regularization for cross-talk already achieved results that are near optimum, or in other words, that the cross-talk is not a major concern for the horizontal inversion. A similar conclusion was reached by <cit.> in the case of a Cartesian inversion for the horizontal component of a 3D flow.

For misfit and errors, the distribution of values over the whole target grid is similar to that presented in Figure <ref> for all results presented in this section.










 §.§ Results: Controlling Misfit and Cross-talk





Still imposing the condition of a maximum misfit of 0.2 and the low SV threshold (same as in case 2 for a given matrix), we now take a slightly different strategy for choosing the trade-off parameters. Motivated by the analysis shown in Figure <ref>, we first choose a maximum value of 100 for the cross-talk (case 3 in the following). This means that the cross-talk averaging kernel _r^θ has a mean magnitude of about 10 times larger than the averaging kernel _θ^θ. If this condition cannot be met, both maximum misfit and cross-talk are increased step by step until a condition is found that can be met in the parameter space. We finally search for the regularization parameter that gives minimal errors.

The inversion results for this case are also displayed in Figure <ref>. In the whole region considered, the inverted flow is very similar to the one obtained in case 2 above. It is noteworthy, however, that the multi-cell structure in the southern hemisphere is even more pronounced in case 3 compared to case 2.

This result was obtained although the dependence of the errors and the cross-talk as a function of depth changed only slightly from case 2 to case 3; see Figures <ref> and <ref>. 
In the near-surface regions, we are able to obtain slightly smaller error bars that increase with depth, peaking at about 3 at the bottom of the convection zone. The cross-talk, on the other hand, decreased slightly in regions below 0.8 R_, and it increased slightly around 0.9 R_. As a result, the magnitude of the cross-talk averaging kernels changed accordingly, which is visible in the right panels of Figure <ref> with a smaller saturated region at r_T=0.9 R_⊙ in case 3 compared to case 1.


We note that to obtain the third result presented here, we allowed the misfit to increase near the bottom of the convection zone by about 20%. We observed that only a slight decrease of the misfit would have increased the errors or the cross-talk by a large amount there. This increase in the misfit is not expected to largely alter the quality of the inversion result compared with the one obtained in Section <ref>, as this change is barely visible in the averaging kernels (horizontal kernels in the left panels of Figure <ref>).





















§ DISCUSSION







The possible existence of a multi-cell meridional flow in the southern hemisphere is further underpinned by analyzing the convolution of the inverted flows v_θ^inv from Figure <ref> with the horizontal averaging kernels, that is, Equation (<ref>) with m=θ. The resulting convolved flow v_θ^conv(θ) is displayed in Figure <ref>. It gives an impression as to how the inverted flow would look like if there were no noise and the background flow would be equal to our inversion result v_θ^inv. As can be seen in Figure <ref>, the resulting flow profile would be washed out spatially given the large widths of the averaging kernels in all three cases. Locations with opposite signs near each other are especially no longer seen as clearly. Turning the argument around, it is possible that the locations with a sign change in the inversion result v_θ^inv could be even more pronounced in real solar flow.

In case 3, however, the convolved flow looks qualitatively very similar to the inverted flow profile, just with a somewhat lower amplitude. This may lead us to the conclusion that this flow profile is quite a robust result and that the original flow may be similar, just with a somewhat higher amplitude.











The match of the inverted flow to the data can be observed in Figures <ref> and <ref>, where forward travel times, i.e., the inverted flow convolved with the Born kernels, are compared to the measured travel times. The forward travel times are predominantly consistent with the measured ones within the measurement errors; see Figure <ref> for cases 1 and 3. We again notice for case 3 that the multi-cell circulation in fact does quite a nice job in reproducing the measured travel times in the southern hemisphere and that a slightly larger magnitude of the flow may explain the measured travel times even better. However, not all features visible in the measured travel times are recovered in the inverted flows, and the errors do not permit us to make a final decision in favor of a single or multiple cell profile of the meridional flow. On the other hand, given the above-mentioned relatively large values for misfit, cross-talk, and the width of the averaging kernels, one could say that the forward-modeled and measured travel times agree remarkably well.











Furthermore, we would like to put forward the somewhat more speculative idea that a multi-cell structure may also do a good job explaining the measured travel times in the northern hemisphere. A hint of an additional cell in the northern hemisphere may be seen in the initial inversion result shown in Figure <ref>. Admittedly, such a profile is not seen in the refined inversions of Section <ref> but we note that it may have been washed out by the large widths of the averaging kernels, especially if the magnitude of the flow is small. 

In addition, we note a few lessons learned in our test inversions, although to some extent of a more incomplete nature. Firstly, it is not clear whether the reversal of signs in the measured travel times at high latitudes – see travel distances above 9 and latitudes within about 10 of the edges in the right panel of Figure <ref> – is due to a signal from the flow or whether it is caused by a systematic effect. We therefore checked whether excluding this region from the data affects the inversion results for the three cases considered in Section <ref>. 
The inversion results from the reduced data set are qualitatively very similar to the results shown in Figure <ref>, especially regarding the general structure and direction of the flow. 
The magnitude of the flow at high latitudes and near the bottom of the convection zone, however, is reduced, which may be seen as physically more realistic regarding mass conservation.


Second, when increasing the threshold for the maximum errors in cases 1 and 2 step by step from 1 to 2.5, the multi-cell structure becomes less pronounced step by step and finally is nearly lost. We obtain a single-cell structure instead. At the same time, the match of the forward-modeled travel times to the measured ones becomes worse. This is not a surprising result, as in this case, the errors reach values similar to the magnitude of the flow in the considered regions. Furthermore, we note that further increasing the threshold for the SVD lets the inversion result further tend toward a smooth single-cell profile, although at some values hints of a second cell may be more pronounced than in case 1.



Finally, we checked that if we do not impose condition (<ref>) with m=r and k=θ, namely, that the integral of the cross-talk kernels be equal to zero, the inversion result does not change much, as the integrals of the cross-talk kernels stay near zero <cit.>. Similarly, some tests showed that further decreasing the cross-talk is only possible at the expense of a large increase in the errors and therefore is not an option considered in this work.








§ CONCLUSION




We have performed SOLA inversions for the deep solar meridional flow using Born kernels. We used GONG data obtained by <cit.> and inverted by <cit.> using ray kernels. In addition to performing a comparison to the ray kernel inversion, we performed several inversions in order to study systematic effects in the inversion process. 

A first comparison was made using the SOLA inversion technique employed by <cit.>, where the measurement errors were assumed to be uncorrelated (diagonal covariance matrix) and where no regularization for a cross-talk term was done. Our initial inversion result using Born kernels is qualitatively similar to the one obtained using ray kernels by <cit.> while displaying a little more fine structure.

Furthermore, we compared the errors of the inverted flows, which were propagated through the inversion using the diagonal covariance matrix, to the errors computed using the full covariance matrix. Previous results for the deep meridional flow using time-distance helioseismology (, ) were obtained under the assumption of uncorrelated measurements. In our inversion, we find that the errors may be underestimated by a factor of about two to four under such an assumption.


In a subsequent step, we performed refined inversions using the full covariance matrix, thus taking the correlation of the measurements into account, and we introduced a regularization for the cross-talk term as is standard practice in local helioseismology <cit.>. Three different cases for the choice of inversion parameters were studied. In all of these inversions, we find a shallow return flow at a depth of about 0.9 R_, in agreement with <cit.>.


In a first case, we obtained a single-cell meridional flow profile with a few opposite sign “bumps”, which is again similar to the result obtained by <cit.>, with less latitudinal variation compared to the result obtained with the diagonal covariance matrix. Here, we used a medium threshold for the SVD and about 50% of the SVs were kept.


When we lower the threshold by a factor of 10 and use about 70% of the SVs, the refined inversion results exhibit a multi-cell structure in the southern hemisphere and a single-cell meridional circulation structure in the northern hemisphere (cases 2 and 3). Depending on the associated radial flow, the inverted flow has two or three circulation cells stacked radially. At the locations of the additional cells, the inversion result for case 1 already showed a flow structure. This structure was, however, not as pronounced and no conclusion was made on this pattern earlier <cit.>.

In principle, a discussion of the inversion results (see Section <ref>) suggests that a multi-cell profile in the southern hemisphere is a suitable candidate as an inversion result, as it is reproduced when convolved with the averaging kernels. All other results are averaged out spatially by this procedure. Furthermore, we speculate that a multi-cell circulation structure may also be present in the northern hemisphere. Indications for this possibility can be seen in the initial result, although it is not present in the refined inversions. However, we point out that both the single-cell and the multi-cell profiles obtained in the refined inversions agree similarly well with the measured travel times within the measurement errors.


Furthermore, the cross-talk of the radial flow into the inversion result for the horizontal flow component was estimated to be below 1.4 . Therefore, it is only of minor importance, although the cross-talk averaging kernels have relatively large values in all inversions performed in this work.


Finally, a flow structure with multiple cells in latitude as obtained by <cit.> using global helioseismology is not seen in the refined inversion results of our study. However, some latitudinal variation can be observed in the travel times, in the initial inversion result, and in the original inversion result from <cit.> using ray kernels. Such a signal may have been averaged out by the relatively wide averaging kernels that have an FWHM of up to 20. We therefore cannot rule out the possibility of latitudinally stacked flow cells at this point.


Further insights in the inversion problem may be reached in future studies, e.g., by comparing to RLS/LSQR inversion techniques <cit.>, or to the Pinsker method newly introduced to helioseismology <cit.>. As pointed out by <cit.>, inversions for both radial and horizontal flows are most promising when a mass-conservation constraint is applied; see <cit.> for such an inversion. An inversion for radial flows without using mass conservation is not very promising; see <cit.>.

In addition, in order to reach an agreement on the nature of the deep meridional flow, understanding systematic effects such as the center-to-limb effect <cit.> and the effect of surface magnetic fields on the measurements <cit.> may play a key role in the future.







The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement n. 307117. This work was supported by the SOLARNET project (www.solarnet-east.eu), funded by the European Commission's FP7 Capacities Programme under the Grant Agreement 312495. J.J. acknowledges support from the National Science Foundation under Grant Number 1351311. S.K. was supported by NASA's Heliophysics Grand Challenges Research grant 13-GCR1-2-0036. The authors acknowledge fruitful discussions during the international team meeting on “Studies of the Deep Solar Meridional Flow” at ISSI (International Space Science Institute), Bern. This work utilizes data obtained by the Global Oscillation Network Group (GONG) Program, managed by the National Solar Observatory, which is operated by AURA, Inc., under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Interamerican Observatory. The authors thank Sylvain Korzennik for providing GONG mode frequencies and damping rates. V.B. thanks Damien Fournier for helpful discussions and Juan Manuel Borrero for reading an earlier version of the manuscript. The authors thank the referee for constructive comments that improved the manuscript.























[Baldner and Schou2012]Baldner2012
Baldner, C. S. and Schou, J.: 2012,
 760, L1,
  [http://adsabs.harvard.edu/abs/2012ApJ...760L...1BADS],
[http://dx.doi.org/10.1088/2041-8205/760/1/L1DOI]

[Birch and Felder2004]Birch2004
Birch, A. C. and Felder, G.: 2004,
 616, 1261,
  [http://adsabs.harvard.edu/abs/2004ApJ...616.1261BADS],
[http://dx.doi.org/10.1086/424961DOI]

[Birch and Gizon2007]BG2007
Birch, A. C. and Gizon, L.: 2007,
AN 328, 228,
  [http://adsabs.harvard.edu/abs/2007AN....328..228BADS],
[http://dx.doi.org/10.1002/asna.200610724DOI]

[Birch et al.2004]BirchKo2004
Birch, A. C., Kosovichev, A. G., and Duvall, Jr., T. L.: 2004,
 608, 580,
  [http://adsabs.harvard.edu/abs/2004ApJ...608..580BADS],
[http://dx.doi.org/10.1086/386361DOI]

[Birch et al.2001]Birch2001
Birch, A. C., Kosovichev, A. G., Price, G. H., and Schlottmann, R. B.:
  2001,
 561, L229,
  [http://adsabs.harvard.edu/abs/2001ApJ...561L.229BADS],
[http://dx.doi.org/10.1086/324766DOI]

[Böning et al.2017]Boening2017
Böning, V. G. A., Roth, M., Jackiewicz, J., and Kholikov, S.: 2017,
 838, 53,
  [http://adsabs.harvard.edu/abs/2017ApJ...838...53BADS],
[http://dx.doi.org/10.3847/1538-4357/aa6333DOI]

[Böning et al.2016]Boening2016
Böning, V. G. A., Roth, M., Zima, W., Birch, A. C., and Gizon,
  L.: 2016,
 824, 49,
  [http://adsabs.harvard.edu/abs/2016ApJ...824...49BADS],
[http://dx.doi.org/10.3847/0004-637X/824/1/49DOI]

[Couvidat et al.2006]Couvidat2006
Couvidat, S., Birch, A. C., and Kosovichev, A. G.: 2006,
 640, 516,
  [http://adsabs.harvard.edu/abs/2006ApJ...640..516CADS],
[http://dx.doi.org/10.1086/500103DOI]

[DeGrave and Jackiewicz2015]DeGrave2015
DeGrave, K. and Jackiewicz, J.: 2015,
 290, 1547,
  [http://adsabs.harvard.edu/abs/2015SoPh..290.1547DADS],
[http://dx.doi.org/10.1007/s11207-015-0693-0DOI]

[DeGrave et al.2014a]DeGrave2014Sunspot
DeGrave, K., Jackiewicz, J., and Rempel, M.: 2014a,
 794, 18,
  [http://adsabs.harvard.edu/abs/2014ApJ...794...18DADS],
[http://dx.doi.org/10.1088/0004-637X/794/1/18DOI]

[DeGrave et al.2014b]DeGrave2014QS
DeGrave, K., Jackiewicz, J., and Rempel, M.: 2014b,
 788, 127,
  [http://adsabs.harvard.edu/abs/2014ApJ...788..127DADS],
[http://dx.doi.org/10.1088/0004-637X/788/2/127DOI]

[Fournier et al.2014]Fournier2014
Fournier, D., Gizon, L., Hohage, T., and Birch, A. C.: 2014,
 567, A137,
  [http://adsabs.harvard.edu/abs/2014A%26A...567A.137FADS],
[http://dx.doi.org/10.1051/0004-6361/201423580DOI]

[Fournier et al.2016]Fournier2016
Fournier, D., Gizon, L., Holzke, M., and Hohage, T.: 2016,
Inverse Problems 32, 105002,
  [http://adsabs.harvard.edu/abs/2016InvPr..32j5002FADS],
[http://dx.doi.org/10.1088/0266-5611/32/10/105002DOI]

[Giles2000]GilesPhD
Giles, P. M.: 2000,
Ph.D. thesis, Stanford University,
[http://adsabs.harvard.edu/abs/2000PhDT.........9GADS]

[Gizon et al.2017]Gizon2017
Gizon, L., Barucq, H., Duruflé, M., Hanson, C. S., Leguèbe, M., Birch,
  A. C., Chabassier, J., Fournier, D., Hohage, T., and Papini, E.: 2017,
Astronomy and Astrophysics 600, A35,
  [http://adsabs.harvard.edu/abs/2017A%26A...600A..35GADS],
[http://dx.doi.org/10.1051/0004-6361/201629470DOI]

[Gizon and Birch2002]GB2002
Gizon, L. and Birch, A. C.: 2002,
 571, 966,
  [http://adsabs.harvard.edu/abs/2002ApJ...571..966GADS],
[http://dx.doi.org/10.1086/340015DOI]

[Gizon and Birch2004]GB2004
Gizon, L. and Birch, A. C.: 2004,
 614, 472,
  [http://adsabs.harvard.edu/abs/2004ApJ...614..472GADS],
[http://dx.doi.org/10.1086/423367DOI]

[Gizon and Birch2005]GB2005
Gizon, L. and Birch, A. C.: 2005,
Living Rev. Sol. Phys. 2, 6,
  [http://adsabs.harvard.edu/abs/2005LRSP....2....6GADS],
[http://dx.doi.org/10.12942/lrsp-2005-6DOI]

[Hartlep et al.2013]Hartlep2013
Hartlep, T., Zhao, J., Kosovichev, A. G., and Mansour, N. N.: 2013,
 762, 132,
  [http://adsabs.harvard.edu/abs/2013ApJ...762..132HADS],
[http://dx.doi.org/10.1088/0004-637X/762/2/132DOI]

[Hathaway2012]Hathaway2012
Hathaway, D. H.: 2012,
ApJ 760, 84,
  [http://adsabs.harvard.edu/abs/2012ApJ...760...84HADS],
[http://dx.doi.org/10.1088/0004-637X/760/1/84DOI]

[Hazra et al.2014]Hazra2014
Hazra, G., Karak, B. B., and Choudhuri, A. R.: 2014,
 782, 93,
  [http://adsabs.harvard.edu/abs/2014ApJ...782...93HADS],
[http://dx.doi.org/10.1088/0004-637X/782/2/93DOI]

[Jackiewicz et al.2012]Jackiewicz2012
Jackiewicz, J., Birch, A. C., Gizon, L., Hanasoge, S. M., Hohage, T.,
  Ruffio, J.-B., and Švanda, M.: 2012,
 276, 19,
  [http://adsabs.harvard.edu/abs/2012SoPh..276...19JADS],
[http://dx.doi.org/10.1007/s11207-011-9873-8DOI]

[Jackiewicz et al.2008]Jackiewicz2008
Jackiewicz, J., Gizon, L., and Birch, A. C.: 2008,
Solar Physics 251, 381,
  [http://adsabs.harvard.edu/abs/2008SoPh..251..381JADS],
[http://dx.doi.org/10.1007/s11207-008-9158-zDOI]

[Jackiewicz et al.2015]Jackiewicz2015
Jackiewicz, J., Serebryanskiy, A., and Kholikov, S.: 2015,
 805, 133,
  [http://adsabs.harvard.edu/abs/2015ApJ...805..133JADS],
[http://dx.doi.org/10.1088/0004-637X/805/2/133DOI]

[Kholikov et al.2014]Kholikov2014
Kholikov, S., Serebryanskiy, A., and Jackiewicz, J.: 2014,
 784, 145,
  [http://adsabs.harvard.edu/abs/2014ApJ...784..145KADS],
[http://dx.doi.org/10.1088/0004-637X/784/2/145DOI]

[Korzennik et al.2013]Korzennik2013ApJ
Korzennik, S. G., Rabello-Soares, M. C., Schou, J., and Larson, T. P.:
  2013,
 772, 87,
  [http://adsabs.harvard.edu/abs/2013ApJ...772...87KADS],
[http://dx.doi.org/10.1088/0004-637X/772/2/87DOI]

[Kosovichev and
  Duvall1997]KosovichevDuvall1997
Kosovichev, A. G. and Duvall, Jr., T. L.: 1997,
in F. P. Pijpers, J. Christensen-Dalsgaard, and C. S. Rosenthal
  (eds.), SCORe'96 : Solar Convection and Oscillations and their
  Relationship, Vol. 225 of Astrophysics and Space Science Library, pp
  241–260,
[http://adsabs.harvard.edu/abs/1997ASSL..225..241KADS],
[http://dx.doi.org/10.1007/978-94-011-5167-2_26DOI]

[Liang et al.2017]Liang2017
Liang, Z.-C., Birch, A. C., Duvall, T. L., Gizon, L., and Schou, J.: 2017,
Astronomy & Astrophysics 601, A46,
  [http://adsabs.harvard.edu/abs/2017arXiv170400475LADS],
[http://dx.doi.org/10.1051/0004-6361/201730416DOI]

[Liang and
  Chou2015a]Liang2015SurfaceMagnOnMeridFlow
Liang, Z.-C. and Chou, D.-Y.: 2015a,
 805, 165,
  [http://adsabs.harvard.edu/abs/2015ApJ...805..165LADS],
[http://dx.doi.org/10.1088/0004-637X/805/2/165DOI]

[Liang and
  Chou2015b]Liang2015ProbingMagneticFields
Liang, Z.-C. and Chou, D.-Y.: 2015b,
 809, 150,
  [http://adsabs.harvard.edu/abs/2015ApJ...809..150LADS],
[http://dx.doi.org/10.1088/0004-637X/809/2/150DOI]

[Miesch2005]Miesch2005LRSP
Miesch, M. S.: 2005,
Living Rev. Sol. Phys. 2, 1,
  [http://adsabs.harvard.edu/abs/2005LRSP....2....1MADS],
[http://dx.doi.org/10.12942/lrsp-2005-1DOI]

[Pijpers and Thompson1994]Pijpers1994
Pijpers, F. P. and Thompson, M. J.: 1994,
Astronomy and Astrophysics 281, 231,
  [http://adsabs.harvard.edu/abs/1994A%26A...281..231PADS]

[Rajaguru and Antia2015]Rajaguru2015
Rajaguru, S. P. and Antia, H. M.: 2015,
 813, 114,
  [http://adsabs.harvard.edu/abs/2015ApJ...813..114RADS],
[http://dx.doi.org/10.1088/0004-637X/813/2/114DOI]

[Roth et al.2016]Roth2016
Roth, M., Doerr, H.-P., and Hartlep, T.: 2016,
 592, A106,
  [http://adsabs.harvard.edu/abs/2016A%26A...592A.106RADS],
[http://dx.doi.org/10.1051/0004-6361/201526971DOI]

[Roth et al.2007]Roth2007
Roth, M., Gizon, L., and Beck, J. G.: 2007,
AN 328, 215,
  [http://adsabs.harvard.edu/abs/2007AN....328..215RADS],
[http://dx.doi.org/10.1002/asna.200610722DOI]

[Schad et al.2013]Schad2013
Schad, A., Timmer, J., and Roth, M.: 2013,
 778, L38,
  [http://adsabs.harvard.edu/abs/2013ApJ...778L..38SADS],
[http://dx.doi.org/10.1088/2041-8205/778/2/L38DOI]

[SILSO World Data
  Center2012]SSN_SILSO_2004_2012
SILSO World Data Center: 2004-2012,
International Sunspot Number Monthly Bulletin and online
  catalogue,
[http://www.sidc.be/silso/URL]

[Švanda2013]Svanda2013a
Švanda, M.: 2013,
 775, 7,
  [http://adsabs.harvard.edu/abs/2013ApJ...775....7SADS],
[http://dx.doi.org/10.1088/0004-637X/775/1/7DOI]

[Švanda et al.2011]Svanda2011
Švanda, M., Gizon, L., Hanasoge, S. M., and Ustyugov, S. D.: 2011,
 530, A148,
  [http://adsabs.harvard.edu/abs/2011A%26A...530A.148SADS],
[http://dx.doi.org/10.1051/0004-6361/201016426DOI]

[Švanda
  et al.2013a]Svanda2013Validation
Švanda, M., Roudier, T., Rieutord, M., Burston, R., and Gizon,
  L.: 2013a,
 771, 32,
  [http://adsabs.harvard.edu/abs/2013ApJ...771...32SADS],
[http://dx.doi.org/10.1088/0004-637X/771/1/32DOI]

[Švanda et al.2013b]Svanda2013c
Švanda, M., Schunker, H., and Burston, R.: 2013b,
Journal of Physics Conference Series 440(1), 012024,
  [http://adsabs.harvard.edu/abs/2013JPhCS.440a2024SADS],
[http://dx.doi.org/10.1088/1742-6596/440/1/012024DOI]

[Zhao et al.2013]Zhao2013
Zhao, J., Bogart, R. S., Kosovichev, A. G., Duvall, Jr., T. L., and
  Hartlep, T.: 2013,
 774, L29,
  [http://adsabs.harvard.edu/abs/2013ApJ...774L..29ZADS],
[http://dx.doi.org/10.1088/2041-8205/774/2/L29DOI]

[Zhao et al.2012a]Zhao2012
Zhao, J., Couvidat, S., Bogart, R. S., Parchevsky, K. V., Birch,
  A. C., Duvall, T. L., Beck, J. G., Kosovichev, A. G., and Scherrer,
  P. H.: 2012a,
 275, 375,
  [http://adsabs.harvard.edu/abs/2012SoPh..275..375ZADS],
[http://dx.doi.org/10.1007/s11207-011-9757-yDOI]

[Zhao et al.2012b]Zhao2012CtoL
Zhao, J., Nagashima, K., Bogart, R. S., Kosovichev, A. G., and
  Duvall, Jr., T. L.: 2012b,
 749, L5,
  [http://adsabs.harvard.edu/abs/2012ApJ...749L...5ZADS],
[http://dx.doi.org/10.1088/2041-8205/749/1/L5DOI]

[Zhao et al.2016]Zhao2016
Zhao, J., Stejko, A., and Chen, R.: 2016,
Solar Physics 291, 731,
  [http://adsabs.harvard.edu/abs/2016SoPh..291..731ZADS],
[http://dx.doi.org/10.1007/s11207-016-0864-7DOI]




