© The Authors 2024

1. Introduction

The photospheric layers of the Sun are optically thick, which prevents us from directly studying the interior of our mother star. In the past seven decades, some of the information hiding below the impenetrable solar skin may have serendipitously been revealed indirectly with the helioseismic method.

Helioseismology (Gizon & Birch 2005; Basu 2016) in general encompasses a set of methods that enable us to examine the inner composition and dynamics of the Sun through the analysis of solar oscillations. Many aspects of the solar interior structure have been successfully revealed through helioseismology, such as the depth of the convection zone, the character of the rotation of the convective envelope, the global meridional mass circulation, or indications of a pre-emergence of active regions.

Helioseismology is based on the observation and interpretation of the manifestations of seismic waves that propagate through the interior. These waves are excited near the solar surface, likely by the vigorous convection within solar granules. The travel path (in the first approximation) is determined by the 1D structure of the Sun. Perturbations in this background structure cause the waves to deviate in propagation and change speed (Christensen-Dalsgaard 2008). The waves manifest themselves in some surface observables, such as the continuum intensity, or in Doppler shifts and the line depths of photospheric spectral lines and others. They are best seen in the above mentioned Doppler shifts. Hence, current helioseismic instruments record spatially resolved Doppler shifts of a particular photospheric spectral line, the Dopplergrams. The sequence of the Dopplergrams that records dynamical changes at a particular spatial position is then processed in order to learn about the behaviour of the solar seismic waves. These Dopplergram sequences are available for instance from the Michelson Doppler Imager (MDI; Scherrer et al. 1995) or from the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012; Schou et al. 2012 instruments.

A specific approach called local helioseismology (see a review by Gizon et al. 2010) focuses on observation and analysis of the waves in a spatially limited region. The wave field within the chosen region is influenced by the properties of that particular area at the surface and beneath it. Consequently, by analysing the wave patterns in this local region, it becomes possible to deduce local 3D structures and dynamics beneath the surface (Christensen-Dalsgaard 2002). The standard procedure used in all methods of local helioseismology involves selecting a small region on the solar surface and tracking it within a co-rotation frame, that is, the region does not move in time relative to the surface of the Sun.

A specific method of local helioseismology, the time-distance helioseismology (Duvall 1993), studies the travel times of the waves that are tracked when they travel between two places on the Sun. In its simplicity, the travel time τ is measured between two measurement points (pixels) in the surface measurements, r₁ and r₂. We indicate a helioseismic observable (e.g. a sequence of Dopplergrams that might be even filtered) at these points as Ψ(r₁, t) and Ψ(r₂, t). The key mathematical quantity to be evaluated is then the cross-covariance, which is evaluated as

$$\begin{array}{c}\hfill C({\mathbf{r}}_1,{\mathbf{r}}_2,t)=\frac{h_t}{T-|t|}\sum _{t^{\prime }}\mathrm{\Psi }({\mathbf{r}}_1,t^{\prime })\mathrm{\Psi }({\mathbf{r}}_2,t^{\prime }+t)\end{array}$$(1)

in its discrete form (Gizon & Birch 2002). Here, h_t is the temporal sampling, and T is the duration of the observation.

In realistic applications, the cross-covariance function that is measured for two points contains much realisation noise that is due to the stochastic nature of solar oscillations. Hence, it is almost impossible to measure the travel time of the waves that travel between two individual pixels in the helioseismic observables. Some averaging is usually required. Most often, this averaging is applied both in time and space. In time, cross-covariances are usually averaged over a time period T longer than about 6 hours. To further improve the signal-to-noise ratio, it was first suggested by Duvall (1993) to average the cross-covariance C(r₁, r₂, t) over points r₂ that belong to the annulus or quadrants centred at r₁ with a distance Δ = |r₂ − r₁|. To study the waves travelling from the central point towards the surrounding annulus and in the opposite direction, averaging over a full annulus is suitable. To investigate the waves that travel in the direction parallel to the equator (east-west), an averaging over 90-degree quadrants of the annulus in the corresponding direction is more suitable. In the method used in this study, quadrants were replaced by sections of the annulus multiplied by the cosine (east-west direction) or sine (south-north direction) of the horizontal polar angle. These geometries constitute a continuous transition from quadrant geometries.

For a given set of r₁ and Δ, the cross-covariance function oscillates around a characteristic travel time in both the positive and negative part of the time axis. It is a goal of the time-distance methods to search for the travel time that maximises the cross-covariance (1). One of the approaches to measuring the travel time is to fit a Gaussian wavelet to the cross-covariance function (Kosovichev & Duvall 1997),

$$\begin{array}{c}\hfill C_{+}(t;{\mathbf{r}}_1,\mathrm{\Delta })=Aexp[-{\gamma }^2{(t-t_{\mathrm{g}})}^2]cos[{\omega }_0(t-{\tau }_{+})],\end{array}$$(2)

where A is the amplitude, t_g represents the group travel time, and τ₊ is the positive phase travel time. Parameters A and γ describe the amplitude and decay rate of the wavelet envelope. An analogous equation may be written for a negative travel time. The time-distance analysis relies on the measurement and interpretation of the phase travel times τ₊ and τ₋.

The fitting of the Gaussian wavelet is demanding and often fails due to the presence of noise. Therefore, alternative definitions of travel times were developed that are more robust in the presence of noise. They stem from the fact that the cross-covariance of the waves may be computed from the solar model by solving the equations of the wave propagation. These cross-covariances represent some sort of a reference for a spherically symmetrical 1D quiet Sun. We can safely assume that the changes to the cross-covariances caused by the perturbations of the model introduce a modification of the background-model cross-covariances. Then, instead of fitting a five-parameter Gaussian wavelet, this reference cross-covariance can be used as a template, and only a one-parameter fitting in time is required. This approach is also used in geoseismology (Zhao & Jordan 1998), and for the need of time-distance helioseismology, it was devised by Gizon & Birch (2002). Following this paper, we refer to these travel times as to GB02 travel times henceforth. In this approach, not the full travel time τ is measured, but the deviation δτ of the travel time from the reference template.

Gizon & Birch (2004) later simplified this definition even further for the case of small perturbations. This robust definition allows the measurement of travel times averaged over a short time (e.g. around T = 2 hours) and with a high spatial resolution. This linearised approach is called GB04 henceforth.

We note that the GB02 and GB04 definitions of the travel times both result in a travel time that is a mixture of the group and phase travel times. It is not straightforward to decouple the two. A proper choice of the reference cross-covariance C⁰ is essential to minimise the contribution of the group travel time to the results. In some cases, a reference cross-covariance from the forward modelling may not be accurate enough. In these cases, the construction of the reference cross-covariance as a spatial average of different realisations over the field of view in the quiet-Sun regions may be a better choice.

We considered that the solar model is fully described by a multidimensional vector of quantities Q^α(x), where the index α lists the considered physical quantities. These quantities can, for instance, be vector flows, the sound speed, the density, or the adiabatic exponent. In our study, we primarily focused on the flows v. The quantities Q^α are functions of 3D spatial coordinates x and, in principle, also of time. However, we only considered stationary models in terms of the averaging over the observation interval T.

Similarly to the discussion of the travel times, the values of the quantities can be split into the unperturbed background Q₀^α(x), which in most cases is represented by the background solar model, and the perturbation q^α(x).

To compute the linear adiabatic oscillations of the Sun, methods have been available for a very long time (Lynden-Bell & Ostriker 1967). This computation is slow and expensive even for the computational power of current computers. This approach is not viable for realistic problems. When we focus on linear perturbations of the model alone, we can relate travel-time deviations δτ with perturbations of the solar model q^α by the equation

$$\begin{array}{c}\hfill \delta {\tau }^a(\mathbf{r})=\sum _{\alpha }{\int }_{\odot }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{\mathbf{K}}_{\alpha }^a({\mathbf{r}}^{\prime }-\mathbf{r},z)·q^{\alpha }({\mathbf{r}}^{\prime },z)+n^a(\mathbf{r}).\end{array}$$(3)

The index a represents the travel-time geometry and combines the choice of the k − ω filter, point-to-annulus, or point-to-quadrant averaging, and the radius of the annulus Δ. The measured travel times are subject to random realisation noise n^a. The quantity K_α^a is called a sensitivity kernel and represents a function that translates changes in the solar model into travel-time deviations. The assumption of small perturbations allows us to split the total travel-time deviation into individual contributions by individual perturbers α.

The sensitivity kernel functions are calculated numerically using the background solar model. Various approaches have been developed by different authors. The simplest approximation ignores the finite-wavelength effects and assumes that the waves propagate along their optimal rays (see Kosovichev & Duvall 1997, and references therein). Finite-wavelength effects are considered by approximating the wave propagation in terms of scattering. Single-source sensitivity kernels were proposed and computed, for example, by Birch & Kosovichev (2000). Sensitivity kernels considering randomly distributed sources were proposed and developed by Gizon & Birch (2002). Later, a linearised approach consistent with the definition of linearised travel times was proposed by Gizon & Birch (2004). Gizon and Birch emphasised the importance of incorporating data processing steps (not only the k − ω filtering, but also an accurate estimate of the telescope point spread function) into the calculation of the kernels and ensuring that they align with the measured travel times. For details of the kernel calculation, we refer to Burston et al. (2015).

It needs to be mentioned that there are other more sophisticated methods that allow the computation of sensitivity kernels, namely the adjoint method (Hanasoge et al. 2011). The method relies on the computation of the wave field that is driven by two sources. One source is the excitation spectrum resulting in the forward wavefield, and the other source is the difference between the predicted wavefield and observations (the adjoint wavefield). The method thus minimises the misfit between the modelled and observed wavefield, whereas the sensitivity kernels are produced naturally within the process.

2. Motivation for this work

In real applications, the problem is usually posed in the opposite direction: We wish to learn from the measured travel-time maps about the properties of the solar interior, and therefore, we usually wish to invert (3) to infer q^α. Not to mention issues with the chosen inversion schemes used by the community, the procedure critically depends on the accuracy of the sensitivity kernels. The kernels represent the background model of the Sun, around which the problem is linearised and solved. If the kernels were problematic, the recovered tomographic maps of the solar interior would also be affected by these issues.

Recently, one of us indicated inconsistencies in the special type of inversions for flows, where the technique of ensemble averaging was applied (Švanda 2015). This technique uses the assumption that when many realisation of the same phenomenon are averaged (it was the supergranular cells in that case), the averaging decreases the noise term in Eq. (3). In the inversion procedure, the term containing the sensitivity kernels is then of the most importance, and the results of the statistical inversion are constrained by these kernels alone. The study showed that various choices of sensitivity kernels in the inversion lead to very different models of the average supergranule. It was not possible to tell from the results and travel times alone which of the models was better, that is, closer to reality.

A similar issue was subsequently confirmed by DeGrave & Jackiewicz (2015), who showed that the observed travel-time maps are inconsistent with the forward-modelled travel times using various supergranular models. Their study directly tested Eq. (3) and pointed out possible issues with the kernels.

On the other hand, Duvall (2006) directly measured sensitivity kernels for flows for the surface gravity modes using small magnetic points as scatterers and obtained a reasonable match with the models.

A direct validation of the sensitivity kernels with observations is not possible. It is possible to test the kernels with respect to the Sun-like magnetohydrodynamic simulations, but most of the available simulations are not exactly Sun-like in terms of the properties of the surface flows and spectrum of solar oscillations. This motivated us to test at least some of the properties of the kernels with observations.

3. Testing procedure

The aim of this study is to test the validity of velocity (flow) perturbation sensitivity kernel integrals with a model-independent method. As a tool, we injected an artificial constant longitudinal flow (velocity perturbation) v₀ into data from observations. This was achieved by tracking the selected region on the solar disc with a (known) constant velocity different from the solar rotation. This effectively means that to all perturbations q^α mentioned in Eq. (3), an additional perturbation v₀ representing the constant injected velocity was added, transforming the equation into

$$\begin{array}{c}\hfill \delta {\tau }^a(\mathbf{r})=\sum _{\alpha }{\int }_{\odot }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{\mathbf{K}}_{\alpha }^a({\mathbf{r}}^{\prime }-\mathbf{r},z)·[q^{\alpha }({\mathbf{r}}^{\prime },z)-{\mathbf{v}}_0]+n^a(\mathbf{r}).\end{array}$$(4)

The minus sign in the term [q^α(r′,z)−v₀] arises because we injected the artificial flow by moving the coordinate frame. However, moving the frame in one direction corresponds to an apparent flow of the same velocity, but in the opposite direction, similarly to watching a landscape passing by through the window of a moving train. The term with v₀ in the integrand is bound with a corresponding vector sensitivity kernel where K_v^a. Since v₀ is constant and independent of position, it is not affected by the volume integral and can be taken out of it,

$$\begin{array}{cc}\hfill \delta {\tau }^a(\mathbf{r})=& \sum _{\alpha }{\int }_{\odot }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{\mathbf{K}}_{\alpha }^a({\mathbf{r}}^{\prime }-\mathbf{r},z)·q^{\alpha }({\mathbf{r}}^{\prime },z)\hfill \\ \hfill & -{\mathbf{v}}_0·{\int }_{\odot }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{\mathbf{K}}_v^a({\mathbf{r}}^{\prime },z)+n^a(\mathbf{r}).\hfill \end{array}$$(5)

Strictly speaking, the additional velocity v₀ should not be mistaken for the change in the rotational rate of the reference frame because it is constant with depth. A constant angular velocity corresponding to a different rate of solar rotation would manifest itself with a decreasing tangential velocity with depth. At the surface, the two approaches are equal. However, most of the time-distance kernels we have at our disposal are sensitive only in a depth of a few Mm below the surface, hence the inaccuracies from exchanging these two points of view in an interpretation of v₀ in Eq. (5) are negligible. What is not negligible at all is the contribution of the random-noise realisation n^a. It also is not known a priori. One of the ways to suppress this term is by means of spatial averaging. At the position of neighbouring pixels, the random realisation is correlated (e.g. Gizon & Birch 2004), but the correlation length-scale, on the other hand, is only a few pixels. By averaging over many pixels (100×100 in our case), the contribution of the random noise is therefore significantly lowered.

Equation (5) then becomes

$$\begin{array}{cc}\hfill \langle \delta {\tau }^a(\mathbf{r})\rangle =& \langle \sum _{\alpha }{\int }_{\odot }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{\mathbf{K}}_{\alpha }^a({\mathbf{r}}^{\prime }-\mathbf{r},z)·q^{\alpha }({\mathbf{r}}^{\prime },z)\rangle \hfill \\ \hfill & -{\mathbf{v}}_0·{\int }_{\odot }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{\mathbf{K}}_v^a({\mathbf{r}}^{\prime },z),\hfill \end{array}$$(6)

where angle brackets indicate the spatial averaging. The first term on the right side of Eq. (6) can be viewed as a certain average background travel-time perturbation that is unaffected by the implanted velocity. The second term is independent of position, and therefore, the averaging does not affect it and it remains unchanged. The average travel-time perturbation in a selected central area of the tracked region is clearly directly proportional to the volume integral of the velocity perturbation kernel, with the injected velocity acting as a constant of proportionality. The injected velocity has only a longitudinal (east-west, usually called v_x) component. East-west quadrant-like difference travel times are most sensitive to this type of perturbation (Burston et al. 2015). Then, finally,

$$\begin{array}{c}\hfill \langle \delta {\tau }^a(\mathbf{r})\rangle =\langle \delta {\tau }_{\mathrm{back}}^a(\mathbf{r})\rangle -v_0{\int }_{\odot }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{K}_{v_x}^a({\mathbf{r}}^{\prime },z).\end{array}$$(7)

By applying a set of implanted velocities v₀, we can independently measure the spatial integral of the sensitivity kernel K_vx^a using a simple linear fit stemming from Eq. (7).

4. Data and methods

4.1. Data

The travel times were measured from a series of Dopplergram datacubes. The Dopplergrams were observed by the HMI instrument and stored within the Joint Science Operations Center (JSOC¹). HMI scans the Fe^I 617.3 nm photospheric absorption line at six positions and measures the intensity across different polarisation states. The Dopplergrams were computed using a fast MDI-like algorithm (Couvidat et al. 2012) that employs the use of the discrete approximations of the Fourier coefficients describing the spectral line and correction by factors from pre-computed look-up tables. The final data products are stored within the series hmi.V_45s at JSOC, from which they are read by our pipelines using a module in the SUNPY package.

4.2. Tracking pipeline

The tracking pipeline is responsible for preparing datacubes that are then used for measuring the travel times. We wrote a versatile pipeline from scratch using routines from the SUNPY (SunPy Community et al. 2020) package so that it exactly matched our needs. In addition to the traditional steps such as selecting the local area, applying the azimuthal equidistant (Postel) projection, and creating the datacube in FITS format with a correct header, this pipeline offers an additional feature specific to this work: it implants artificial velocities by apparently moving the selected area at a predefined speed. This is a crucial feature used to determine the model-independent velocity kernel integrals as explained in previous sections. Since this is not part of standard pipelines for datacube preparation, our own solution was necessary.

The pipeline is highly versatile and prepared to be executed both on a single computer and in a cluster environment. Individual consecutive steps were written as separate modules, so that the user can verify the outcome of each step before proceeding further. At a high level, the tracking pipeline performs the steps listed below.

1. Prepare the folder structure together with all necessary configuration files (for the datacube creation step and for the subsequent travel-time pipeline) based on the input configuration (JSOC Data Record Management System (DRMS) queries, data paths, projection origins, and implanted velocities).

2. (OPTIONAL) Download data (Dopplergrams) from JSOC. This step is unnecessary when the required data are already downloaded.

3. Create a datacube for each combination of DRMS query and velocity and save it in the form of a FITS file into its designated folder prepared in step 1.

The configuration of the pipeline is given in the form of a text file. Among others, it contains the DRMS request string specifying the time range to be downloaded and processed, the positions of the centre of the tracked region in Carrington coordinates, lower and upper limits for implanted velocities, and the velocity sample count. The set of implanted velocities is then randomly sampled from within the given velocity bounds. Subsequently, a loop is executed over the realisation of the implanted velocities. Within this loop, the full-disc Dopplergrams are downloaded (they are downloaded only once for the given DRMS query), tracked and mapped following the requests, and stacked to a datacube.

Before a final saving into an FITS file with headers compatible with the following processing, invalid values (NaNs or blank pixels) are removed from the datacube frames and replaced by a median of non-empty values. Then, from each frame, a quadratic smooth surface is removed to subtract large-scale trends in Dopplergrams. A large-scale trend that smoothly varies across the field of view is naturally produced by the solar rotation, for instance.

One potential question is whether the injected longitudinal velocity truly results in an artificial flow within the data. To address this, ring diagrams were constructed from the datacube power spectrum that show that the rings are shifted in the datacubes with an implanted velocity with respect to the situation in the datacube that was tracked with the solar rotation (see Fig. 1). It qualitatively demonstrates the presence of a horizontal flow in the east-west direction, as evidenced by the horizontal distortion of the rings.

Fig. 1. Cuts through the power spectra at a frequency of 4 mHz in the datacube that was tracked with solar rotation (upper panel) and the same datacube tracked 500 m/s faster than solar rotation. The distortion of the rings due to the artificial velocity is clearly visible.

A direct comparison of the horizontal velocity with the injected velocity is also possible. In the Dopplergrams, large-scale convective cells of supergranulation (Rincon & Rieutord 2018) are clearly visible. The supergranular structures can therefore be used to determine the horizontal flow field on the solar surface (e.g. Švanda et al. 2006, and the follow-up papers). This goal may be achieved by applying the method of local correlation tracking (LCT; November 1986) to Dopplergram datacubes.

The local correlation tracking algorithm searches for the optimal displacement that minimises the differences of the floating spatial windows that capture slowly evolving features in the time sequence. Knowing the sampling of the cross-correlated frames, we can convert the detected displacement into the horizontal velocity vector. In the test performed here, the horizontal velocity vector is dominated by the longitudinal component, which is the only one we compare with the implanted velocity.

The results of our pipeline test are shown in Fig. 2. The measured surface horizontal velocities clearly correspond very well to the implanted velocity, and thereby validate the performance of the tracking pipeline.

Fig. 2. Comparison of the implanted velocity and the longitudinal component of the horizontal velocity vector determined from the manipulated datacubes using the LCT algorithm. The fit has a slope of −1, which is a consequence of the fact that the implanted velocity was injected by artificially changing the reference frame, which naturally introduced an apparent flow in the opposite direction.

4.3. Travel-time pipeline

In general, measurements of helioseismic travel times consist of three consecutive steps. The travel-time pipeline operated on the data-processing computer cluster at the Astronomical Institute of the Czech Academy of Sciences in Ondejov is no exception. First, the input datacube, that is, the output of the tracking pipeline, is detrended in time, and a mean frame (average of the frames over time) is subtracted. This step removes existing large-scale trends in the observations. In Dopplergrams, these large-scale trends may be produced by a residual signal from the differential rotation, for instance, and also by long-living velocity features, such as most of the signal from solar supergranules. The detrended datacube is apodised in both space and time and transformed into Fourier space for further processing.

The Fourier-transformed datacube is filtered to keep only certain types of waves. We use two types of filtering that are traditional in time-distance helioseismology, that is, the ridge filters (isolating only the surface gravity f modes or acoustic p-modes to the fourth order) and the phase-speed time-distance filters (indicated as td1 to td11, following Duvall 1997).

These filtered datacubes are inputs for the travel-time measurements procedure. Following the usual recipes, cross-covariance maps are computed numerically (see Equation (1)) for spatially smoothed datacubes using the generalised point-to-annulus geometry. Then, the travel times that minimise the cross-covariance of the signals at the given point are measured using a user-selected definition. In our case, we used both the GB02 and the GB04 definition.

This workflow was implemented in a pipeline coded in the MATLAB language. Each of the three steps is represented in a separate module, and the pipeline is executed in the cluster environment.

4.4. Sensitivity kernels

The sensitivity kernels that were tested were computed in the KC3 code by Aaron Birch (Birch & Gizon 2007). The code uses a mode summation approach in Cartesian geometry and a Born approximation with randomly distributed sources. It uses eigenfunctions and eigenfrequencies of the oscillations from the solar model S (Christensen-Dalsgaard et al. 1996). KC3 is a reference kernel code that produced sensitivity kernels for flow or sound-speed perturbation that were used in many studies before.

5. Results

To achieve the goals of this study, we analysed 716 separate datacubes in total. An overview is given in Table 1. The processed data volume consisted of 705 GB in the Dopplergram datacubes and 8.5 TB in the resulting travel-time maps and necessary intermediate data products. The processing took about two months using the Radegast cluster at the Solar Department of the Astronomical Institute of the Czech Academy of Sciences in Ondřejov, which has 102 CPU cores in total.

For each ridge filter (f and p₁ to p₄ modes), we evaluated the travel times for 16 different distances Δ, starting at 5 px and ending at 20 px, with a pixel size of 1.46 Mm. The reason for the degraded MDI-like resolution is that the sensitivity kernels we have at our disposal and that we regularly use for helioseismic inversions within our group (e.g. Švanda et al. 2013; Švanda 2013; Švanda et al. 2014; Korda & Švanda 2019, 2021) have exactly this pixel size. As for the phase-speed filters td1 to td11, we used five distances Δ for each, following the usual definition of these filters by Couvidat & Birch (2006) (also Table 1 in Švanda 2013). For each mode (filter) and distance, we calculated travel-time maps for all three point-to-annulus and point-to-quadrant geometries. In the following, we only used the east-west averaging geometry, because it is most sensitive to the zonal flows (Burston et al. 2015) that we introduced artificially. Hence, we verified kernel integrals for 135 different sensitivity kernels in total.

The procedure is demonstrated for example in Fig. 3. There, for a chosen sensitivity kernel (index a), we plot the dependence given by Eq. (7), that is, the dependence of the average travel-time deviation in the central region ⟨δτ^a(r)⟩ of the datacube as a function of the implanted zonal velocity v₀. The dependence is monotonous (this is the case for all considered as) and very close to linear, as expected. For higher values of v₀ above around 300 m/s the average travel times start to deviate slightly from the linear fit, which is also over-plotted. This non-linearity is expected for high velocities and is consistent with previous works (see e.g. Jackiewicz et al. 2007). Therefore, we evaluated the fit of the linear function only in the range of ( − 300, +300) m/s. For higher modes (e.g. p₄) and higher phase speeds, the non-linearity is slightly more pronounced. Following Eq. (7), the slope of the linear fit corresponds to the total spatial integral of the corresponding sensitivity kernel.

Fig. 3. Measured average travel time around the disc centre as a function of the implanted velocity. Top: Plot for mode td4 at a distance of 9.9 px. The linear fit and its properties are also indicated. The slope of the linear fit corresponds to the total integral of the appropriate sensitivity kernel. Bottom: Similar plot for the p3 mode and Δ = 15 px when using several dates. Different dates are grouped around virtual lines that are shifted vertically, as indicated by the time stamps.

For the different dates, a different realisation of the background travel-time term ⟨δτ_back^a(r)⟩ in (Eq. 7) occurs (see Fig. 3 bottom panel). For instance, an uncorrected differential rotation contributes to this term. This term changes the intercept of the linear dependence. We are only interested in the slopes. In order to learn about the characteristic slope value, we fitted the linear dependence individually for each date and then obtained the representative slope value s as a weighted average,

$$\begin{array}{c}\hfill s=\frac{\sum _i{w}_i{s}_i}{\sum _i{w}_i},\mathrm{with}{w}_i=\frac{1}{{\sigma }_{s,i}^2},\end{array}$$(8)

where s_i are the individual slopes with their error bars σ_s, i.

The comparison of the resulting slope that represents the model-independent value of the total spatial integral of the sensitivity kernel with the value from the forward modelling for all modes and all distances is plotted in Fig. 4. The colours represent different modes, and the total integral values usually increase with the distance Δ. An expected line with a slope of −1 is over-plotted for reference.

Fig. 4. Volume integrals of the velocity perturbation kernels obtained from linear regression slopes (model independent) vs. volume integrals of the velocity perturbation kernels calculated from forward modelling. The data were obtained from 578 different datacubes for various ridge and phase-speed filters. The dashed line indicates the expected line with a slope of −1. The values on both axes are in units of s/(km/s).

The figure shows that the measured and modelled travel-time sensitivity kernel integrals correspond to each other for most of the ridge filters. This is also demonstrated further in Table 2, where the slopes for the f and p₁ to p₄ modes are close to the expected value of −1. The fits show a small spread, low intercept values, and a large (anti)correlation, and they are statistically significant. The slope deviates from −1 for the p₄ mode, which penetrates to deeper layers, where the assumption of implanting the surface velocity that is applicable to the whole sub-surface region no longer holds.

The situation is different for the phase-speed filters (see e.g. Fig. 5). Even when the correlation coefficient is close to −1 and the fits are significant (all cases except for td1 and perhaps td4, when we take the level of statistical significance 0.01), the slopes differ significantly from the expected value of −1. Moreover, the values of the intercepts are mostly high. The most notable differences occur for low phase speeds, for filters td1 to td5, for instance. Large deviations from the expectations are registered also for the highest phase speeds (filters td9 to td11), where the slope is much larger than the expected value. Some of this disagreement is probably due to the reason described above when discussing the p₄ mode. On the other hand, p-modes with low phase speeds only penetrate shallow near-surface layers, and our testing approach should be valid for them.

Fig. 5. Same as Fig. 4 with details of phase-speed filtered kernels. Moreover, this plot is obtained for a particular date, not for a combination of several dates. The dashed line indicates the expected line with a slope of −1. The values on both axes are in units of s/(km/s).

We note that the linear fits by mode are rather indicative. A complete set of kernel integrals from the forward modelling and their comparison to the model-independent determination is given in Table A.1.

6. Discussion

Our results indicate an unresolved issue with the travel-time sensitivity kernels using phase-speed filtering. Some of the disagreement might be a consequence of the chosen testing procedure, where the implanted velocity is under control only at the surface. Deeper down, its proper interpretation becomes more difficult. The waves with high phase speeds and higher modes are more sensitive to the state of deeper modes. This property is prominent especially for the phase-speed-filtered waves (see e.g. Fig. 2 in Švanda 2013). Therefore, it is difficult to properly interpret the disagreement between the modelled and observed travel-time sensitivity kernel integrals for these modes. On the other hand, the disagreement is also present for very shallow waves (short distances in ridge-filtered data and lowest phase-speed filtered data with modes td1 to td4), where high-degree modes prevail. The disagreement of amplitudes and line widths of the models for the high-degree modes with the observations is known from the literature (see e.g. Korzennik et al. 2013), which affects the Green functions and cross-covariance functions used in the kernel models, among others.

It would seem at the first glance that the different conclusions drawn for ridge- and phase-speed-filtered kernels might indicate different approaches taken in processing the ridge- and phase-speed-filtered data. We wish to stress that this is not the case. In our testing procedure, an identical kernel code and an identical travel-time pipeline was used to process different filtering approaches. The corresponding k − ω filters are supplied to the code as external files with the same format. Both codes and pipelines are written in the programming language MATLAB and share as many common sub-routines as necessary. We would like to stress that our travel times and sensitivity kernels are mutually consistent in terms of the considered data-processing steps. Filtering, spatial resolutions, and so on are taken into account consistently in both the computation of the sensitivity kernels and measuring the travel times.

There is a handfull of choices in running the travel-time pipeline that might possibly influence the results of the presented comparison.

Choice of the travel-time fitting method. Our travel-time pipeline offers three definitions in total of the travel time to be used: Gabor-wavelet fitting, GB02, and GB04. We only used GB02 and GB04 in this study because our travel-time sensitivity kernels are compatible with these definitions. We tested that both approaches provided essentially indistinguishable results for implanted velocities in the range of about ( − 200, 200) m/s, then non-linearities started to appear for GB04 definition for some travel-time geometries. This is expected because the GB04 definition is only valid for small perturbations. For the GB04 definition, the points start to deviate from the linear fit sooner than shown in Fig. 3, for example, which was computed with the GB02 definition. For this reason, we preferred the GB02 definition for our final results even though it is significantly more demanding computationally. It allowed a wider range of the implanted velocities to be examined.

Averaging time. Averaging of the cross-covariances over a longer time span leads to a lower contribution of the random noise; the signal-to-noise ratio scales approximately as $\sqrt{T}$ with T being an averaging time. We computed several testing runs using T = 24 hours and T = 6 hours and compared the results. Averaging of the travel times in the spatial window around the centre of the tracked patch (we used 100 × 100 px window) is much more effective in suppressing the random noise than the choice of the longer averaging time. The averaging time therefore has no effect on the results of this study.

Choice of the reference cross-covariance. The GB02 and GB04 definitions of the travel times rely heavily on the reference cross-covariance that is used as a template to determine the offset δτ of the measured travel time. As pointed out already in the introduction, the reference cross-covariance should be close to the measured one. In the case of low implanted velocities, a usual quiet-Sun reference cross-covariance may be used, for instance the one from the forward modelling. With increasing v₀, the deviations from the template became significant and break the assumption of the perturbation, which is a small change. When the surface is tracked at a rate that differs significantly from the background solar rotation, even the position of the emergent point changes, which again influences the measured cross-covariance when applying the non-modified pipeline. For a wave travelling about 20 minutes from the central point to the surrounding annulus (which is a characteristic time of the propagation in our choice of the waves), for example, the emergent point shifts by about 600 km with a tracking speed of 500 m/s with respect to the background solar rotation. This is about 0.4 pixels, which may already play a role in the measurement procedure.

We found that the use of the pre-computed cross-covariance from the solar model is insufficient to represent the changes introduced by the tracking. The use of the reference cross-covariance computed as an average over all representations in the field of view solves the issue.

Heliographic latitude. The results presented above were obtained with datacubes whose centres were positioned on the solar equator. We also ran the same procedure for two datacubes that were located at heliographic latitudes of +20 and −20 degrees. The results were very similar and only differed by the value of the intercept on the implanted velocity versus averaged travel time plot (see Fig. 3). The shift of the intercepts is a consequence of the differential rotation, which is not compensated for, as the basic rotation rate was the Carrington rotation rate. The comparisons between the measured sensitivity kernel integrals and the models were comparable to those presented in Table 2.

Other travel-time geometries. We focused primarily on the east-west geometry from the point-to-quadrant choices because it is naturally sensitive to the implanted zonal velocity (Burston et al. 2015). Other geometries, such as the point-to-annulus (also known as out-minus-in) or south-north geometry, can be used. The volume integrals of these geometries are zero for the zonal (v_x) velocity. Therefore, the averaged travel times around the disc centre have a low value and vary quasi-randomly with implanted velocity (see Fig. 6).

Fig. 6. Figure similar to Fig. 3, but for the f mode, point-to-annulus geometry with Δ = 11 px.

7. Concluding remarks

We tested the volume integrals of the time-distance sensitivity kernels for flows and pointed out a possible issue in the sensitivity kernels that use the phase-speed filtering scheme. We need to stress that our study does not necessarily have to invalidate results and interpretations of the published helioseismic inversions in the past for the reasons listed below.

• This test represents a form of a necessary condition sensitivity kernels should fulfil, but not a sufficient condition, as the procedure only tested the volume integrals of the kernels.

• The disagreement of the total volume integral may lie in details that might in fact be marginal for a proper analysis. For instance, if there is a slight disagreement of the kernel value at pixels that are far from the central point and if this disagreement has a large area even when the difference is small, it will affect the total integral by a large amount due to the 2D integration.

• Similarly, a shifted mean horizontal value will behave in the same way. A small shift but over a large area will significantly affect the total integral of the sensitivity kernel.

• Švanda (2015) speculated that the reason of the inconsistency of the flow models derived from the inversion using different combinations of travel-time geometries was that during the inversion, the contributions of some kernels were cancelled by contributions from the others. The subtraction to zero critically depends on precise values, and when the values have small errors, the results are largely unconstrained. When the combination leads to the summation of the kernels, these results should be more robust in the presence of inconsistencies. If this hypothesis were to be true, the inversion schemes call for a method for determining the optimal combination of the sensitivity kernels that are to be used in the scheme that would avoid large-scale subtractions.

• Comparisons between different codes for computing the sensitivity kernels were made in the past (e.g. Parchevsky et al. 2010; Böning et al. 2016), and they resulted in a generally good match between different techniques. The disagreement between various codes was to within a few percent. Even the comparison of the total integrals (Parchevsky et al. 2010) of kernels using very different wave-field approximations (ray approximation vs. Born approximation) showed an acceptable agreement. It should be possible with the numbers in Table A.1 to verify the kernel integrals also using other codes, provided that the data processing involved in the kernel computation is consistent with the travel-time measurements presented in this study.

• We would also like to point out that except for the phase-speed filters td1–td4, the measured kernel integrals are somewhat proportional to the values obtained from the forward models (see Fig. 5). This is merely an observation, and no conclusion is drawn from it at the moment.

Many of the helioseismic inferences were achieved using methods involving sensitivity kernels. Our test indicates that there might be an issue with some of them, and therefore, more work on verification and validation is needed. Some helioseismic methods do not require sensitivity kernels and use the full waveform information about the solar oscillations (see e.g. Woodard 2007; Hanson et al. 2021; Mani et al. 2022; Hanson et al. 2024).

Acknowledgments

This paper presents the results obtained within completing the MSc. studies of D.Ch. under the supervision of M.Š. at the Faculty of Mathematics and Physics, Charles University. This work was partly done with a support of the institutional support ASU:67985815 of the Czech Academy of Sciences. M.Š. further acknowledges the support by Czech-German common grant, funded by the Czech Science Foundation under the project 23-07633K and by the Deutsche Forschungsgemeinschaft under the project BE 5771/3-1 (eBer23-13412). We are grateful to Aaron Birch for making his kernel code free to use.

References

Appendix A: All investigated kernels

All Tables

All Figures

	Fig. 1. Cuts through the power spectra at a frequency of 4 mHz in the datacube that was tracked with solar rotation (upper panel) and the same datacube tracked 500 m/s faster than solar rotation. The distortion of the rings due to the artificial velocity is clearly visible.
In the text

	Fig. 2. Comparison of the implanted velocity and the longitudinal component of the horizontal velocity vector determined from the manipulated datacubes using the LCT algorithm. The fit has a slope of −1, which is a consequence of the fact that the implanted velocity was injected by artificially changing the reference frame, which naturally introduced an apparent flow in the opposite direction.
In the text

Fig. 3. Measured average travel time around the disc centre as a function of the implanted velocity. Top: Plot for mode td4 at a distance of 9.9 px. The linear fit and its properties are also indicated. The slope of the linear fit corresponds to the total integral of the appropriate sensitivity kernel. Bottom: Similar plot for the p3 mode and Δ = 15 px when using several dates. Different dates are grouped around virtual lines that are shifted vertically, as indicated by the time stamps.

In the text

Fig. 4. Volume integrals of the velocity perturbation kernels obtained from linear regression slopes (model independent) vs. volume integrals of the velocity perturbation kernels calculated from forward modelling. The data were obtained from 578 different datacubes for various ridge and phase-speed filters. The dashed line indicates the expected line with a slope of −1. The values on both axes are in units of s/(km/s).

In the text

	Fig. 5. Same as Fig. 4 with details of phase-speed filtered kernels. Moreover, this plot is obtained for a particular date, not for a combination of several dates. The dashed line indicates the expected line with a slope of −1. The values on both axes are in units of s/(km/s).
In the text

	Fig. 6. Figure similar to Fig. 3, but for the f mode, point-to-annulus geometry with Δ = 11 px.
In the text