Combined helioseismic inversions for 3D vector flows and sound-speed perturbations

TL;DR

This paper presents an advanced helioseismic inversion method that simultaneously retrieves 3D vector flows and sound-speed perturbations in the solar interior, improving accuracy and reducing cross-talk effects.

Contribution

The authors develop a novel inversion approach using SOLA with Born approximation kernels to recover full 3D flows and sound-speed variations simultaneously, validated on synthetic data.

Findings

Successfully recover flows and sound-speed perturbations near the solar surface.

The methodology is sensitive to the entire vertical velocity component.

Cross-talk effects are negligible for horizontal flow inversions, but need minimization for vertical flows and sound-speed perturbations.

Abstract

Time-distance helioseismology is the method of the study of the propagation of waves through the solar interior via the travel times of those waves. The travel times of wave packets contain information about the conditions in the interior integrated along the propagation path of the wave. We introduce an improved methodology of the time-distance helioseismology which allows us to invert for a full 3D vector of flows and the sound-speed perturbations at once. Using this methodology one can also derive the mean value of the vertical component of flows and the cross-talk between the flows and the sound-speed perturbations. We used the SOLA method with a minimisation of the cross-talk as a tool for inverse modelling. In the forward model, we use Born approximation travel-time sensitivity kernels with the Model S as a background. The methodology was validated using forward-modelled travel times with both mean and difference point-to-annulus averaging geometries applied to a snapshot of fully self-consistent simulation of the convection. We tested the methodology on synthetic data. We demonstrate that we are able to recover flows and sound-speed perturbations in the near-surface layers. We have taken the advantage of the sensitivity of our methodology to entire vertical velocity, and not only to its variations as in other available methodologies. The cross-talk from both the vertical flow component and the sound-speed perturbation has only a negligible effect for inversions for the horizontal flow components. The inversions for the vertical component of the vector flows or for the sound-speed perturbations are affected by the cross-talk from the horizontal components, which needs to be minimised in order to provide valid results. It seems that there is a nearly constant cross-talk between the vertical component of the vector flows and the sound-speed perturbations.

Tables

Table 1. Table 1: Correlation coefficients and RMS of the inversions using the synthetic travel times. In the first two rows and the eighth and ninth rows there are the correlation coefficients of δ ​ q ^ α , α inv 𝛿 superscript subscript ^ 𝑞 𝛼 𝛼 inv \delta\hat{q}_{\alpha,\alpha}^{\mathrm{inv}} and the whole inversion with the ideal answer. In third, fourth, tenth and eleventh rows there are comparisons of the RMS value. In the fifth to seventh and twelfth to fourteenth rows there are RMS values of the individual components of the inversions.

GdG-type target function
cross-talk minimised
		$v_x$	$v_y$	$v_z$	$\delta c_s$
1.	$\mathrm{corr}(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}},\mathrm{id}.\mathrm{ans}.)$	0.97	0.97	0.87	0.90
2.	$\mathrm{corr}(\delta q_{\alpha }^{\mathrm{inv}},\mathrm{id}.\mathrm{ans}.)$	0.90	0.90	0.57	0.13
\hdashline3.	$\mathrm{RMS}\left(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}}\right)/\mathrm{RMS}(\mathrm{id}.\mathrm{ans}.)$	0.81	0.79	0.45	1.20
4.	$\mathrm{RMS}\left(\delta q_{\alpha }^{\mathrm{inv}}\right)/\mathrm{RMS}(\mathrm{id}.\mathrm{ans}.)$	0.86	0.85	0.44	1.50
\hdashline5.	$\mathrm{RMS}\left(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}}\right)$ [m s-1]	55.11	41.63	21.61	31.18
6.	$\mathrm{RMS}\left(\mathrm{cross}-\mathrm{talk}\right)$ [m s-1]	0.51	0.52	0.73	3.82
7.	$\mathrm{RMS}\left(\mathrm{noise}\right)$ [m s-1]	18.33	17.14	3.15	13.78
cross-talk ignored
8.	$\mathrm{corr}(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}},\mathrm{id}.\mathrm{ans}.)$	0.99	0.99	0.94	0.91
9.	$\mathrm{corr}(\delta q_{\alpha }^{\mathrm{inv}},\mathrm{id}.\mathrm{ans}.)$	0.93	0.94	0.36	0.42
\hdashline10.	$\mathrm{RMS}\left(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}}\right)/\mathrm{RMS}(\mathrm{id}.\mathrm{ans}.)$	0.92	0.93	0.47	1.25
11.	$\mathrm{RMS}\left(\delta q_{\alpha }^{\mathrm{inv}}\right)/\mathrm{RMS}(\mathrm{id}.\mathrm{ans}.)$	0.96	0.96	0.49	1.88
\hdashline12.	$\mathrm{RMS}\left(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}}\right)$ [m s-1]	62.76	48.86	22.13	31.39
13.	$\mathrm{RMS}\left(\mathrm{cross}-\mathrm{talk}\right)$ [m s-1]	2.55	2.30	8.92	16.19
14.	$\mathrm{RMS}\left(\mathrm{noise}\right)$ [m s-1]	18.01	15.98	5.25	13.17

Table 2. Table 2: Correlation coefficients and RMS of inversions using the synthetic travel times. The Az target function was used. See Table 1 for details.

Az-type target function
cross-talk minimised
		$\delta c_s$
1.	$\mathrm{corr}(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}},\mathrm{id}.\mathrm{ans}.)$	0.97
1.	$\mathrm{corr}(\delta q_{\alpha }^{\mathrm{inv}},\mathrm{id}.\mathrm{ans}.)$	0.44
\hdashline3.	$\mathrm{RMS}\left(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}}\right)/\mathrm{RMS}(\mathrm{id}.\mathrm{ans}.)$	0.85
4.	$\mathrm{RMS}\left(\delta q_{\alpha }^{\mathrm{inv}}\right)/\mathrm{RMS}(\mathrm{id}.\mathrm{ans}.)$	0.85
\hdashline5.	$\mathrm{RMS}\left(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}}\right)$ [m s-1]	48.49
6.	$\mathrm{RMS}\left(\mathrm{cross}-\mathrm{talk}\right)$ [m s-1]	1.27
7.	$\mathrm{RMS}\left(\mathrm{noise}\right)$ [m s-1]	5.22
cross-talk ignored
8.	$\mathrm{corr}(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}},\mathrm{id}.\mathrm{ans}.)$	0.98
9.	$\mathrm{corr}(\delta q_{\alpha }^{\mathrm{inv}},\mathrm{id}.\mathrm{ans}.)$	0.48
\hdashline10.	$\mathrm{RMS}\left(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}}\right)/\mathrm{RMS}(\mathrm{id}.\mathrm{ans}.)$	0.81
11.	$\mathrm{RMS}\left(\delta q_{\alpha }^{\mathrm{inv}}\right)/\mathrm{RMS}(\mathrm{id}.\mathrm{ans}.)$	0.71
\hdashline12.	$\mathrm{RMS}\left(\delta {\hat{q}}_{\alpha ,\alpha }^{\mathrm{inv}}\right)$ [m s-1]	46.28
13.	$\mathrm{RMS}\left(\mathrm{cross}-\mathrm{talk}\right)$ [m s-1]	7.18
14.	$\mathrm{RMS}\left(\mathrm{noise}\right)$ [m s-1]	4.77

Equations

$$\delta {\tau }^a\left(\mathbf{r}\right)=\underset{\odot }{\int }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z\sum _{\beta =1}^P{K}_{\beta }^a\left({\mathbf{r}}^{\prime }-\mathbf{r},z\right)\delta q_{\beta }\left({\mathbf{r}}^{\prime },z\right)+n^a\left(\mathbf{r}\right),$$

$$\delta q_{\alpha }\left({\mathbf{r}}_0;z_0\right)=\delta q_{\alpha }\left(\delta {\tau }^a,K_{\beta }^a\right)$$

$$\delta q_{\alpha }^{\mathrm{inv}}\left({\mathbf{r}}_0;z_0\right)=\sum _{i=1}^N\sum _{a=1}^M{w}_a^{\alpha }\left({\mathbf{r}}_i-{\mathbf{r}}_0;z_0\right)\delta {\tau }^a\left({\mathbf{r}}_i\right),$$

$${\chi }^2$$

$$+\mu \sum _{i,j,a,b}{w}_a^{\alpha }\left({\mathbf{r}}_i;z_0\right){\Lambda }^{ab}\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)w_b^{\alpha }\left({\mathbf{r}}_j;z_0\right)+$$

$$+\nu \sum _{\beta \ne \alpha }\underset{\odot }{\int }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{\left[{\mathcal{K}}_{\beta }^{\alpha }\left({\mathbf{r}}^{\prime },z;z_0\right)\right]}^2+ϵ\sum _{a,i}{\left[w_a^{\alpha }\left({\mathbf{r}}_i;z_0\right)\right]}^2+$$

$$+\sum _{\beta }{\lambda }^{\beta }\left[\underset{\odot }{\int }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{\mathcal{K}}_{\beta }^{\alpha }\left({\mathbf{r}}^{\prime },z;z_0\right)-{\delta }_{\beta }^{\alpha }\right].$$

$${\mathcal{K}}_{\beta }^{\alpha }\left({\mathbf{r}}^{\prime },z;z_0\right)=\sum _i\sum _a{w}_a^{\alpha }\left({\mathbf{r}}_i-{\mathbf{r}}_0;z_0\right)K_{\beta }^a\left({\mathbf{r}}^{\prime }-{\mathbf{r}}_i,z\right),$$

$$\delta q_{\alpha }^{\mathrm{inv}}\left({\mathbf{r}}_0;z_0\right)=\sum _{\beta }\underset{\odot }{\int }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathrm{d}z{\mathcal{K}}_{\beta }^{\alpha }\left({\mathbf{r}}^{\prime }-{\mathbf{r}}_0,z;z_0\right)\delta q_{\beta }\left({\mathbf{r}}^{\prime },z\right)+\mathrm{noise}.$$

$${\sigma }_{\alpha }^2=\sum _{i,j,a,b}{w}_a^{\alpha }\left({\mathbf{r}}_i;z_0\right){\Lambda }^{ab}\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)w_b^{\alpha }({\mathbf{r}}_j;z_0).$$

$$\delta c_s^2=c_s^2-c_{\mathrm{model}}^2=2c_{\mathrm{model}}\delta c_s+{\left(\delta c_s\right)}^2,$$

$$K_{\delta c_s^2/c_s^2}\frac{\delta c_s^2}{c_s^2}$$

$$=\frac{2K_{\delta c_s^2/c_s^2}}{c_{\mathrm{model}}}\delta c_s=K_{\delta c_s}\delta c_s,$$

$$K_{\delta c_s}\equiv \frac{2K_{\delta c_s^2/c_s^2}}{c_{\mathrm{model}}}.$$

$$\delta {\hat{q}}_{\alpha ,\beta }^{\mathrm{inv}}\left({\mathbf{r}}_0;z_0\right)$$

$$\mathrm{noise}\left({\mathbf{r}}_0;z_0\right)$$

$$\delta q_{\alpha }^{\mathrm{inv}}\left({\mathbf{r}}_0;z_0\right)$$

$$=\sum _{i,a}{w}_a^{\alpha }\left({\mathbf{r}}_i-{\mathbf{r}}_0;z_0\right)\delta {\tau }^a\left({\mathbf{r}}_i\right),$$

$$\mathrm{ideal}\mathrm{answer}\left({\mathbf{r}}_0;z_0\right)$$

Full text

*Aims. *We introduce an improved methodology of the time-distance helioseismology which allows us to invert for a full 3D vector of plasma flows and the sound-speed perturbations at once. Using this methodology one can also derive the mean value of the vertical component of plasma flows and the cross-talk between the plasma flows and the sound-speed perturbations.

*Methods. *We used the Subtractive Optimally Localised Averaging method with a minimisation of the cross-talk as a tool for inverse modelling. In the forward model, we use Born approximation travel-time sensitivity kernels with the Model S as a background. The methodology was validated using forward-modelled travel times with both mean and difference point-to-annulus averaging geometries applied to a snapshot of fully self-consistent simulation of the convection.

*Results. *We tested the methodology on synthetic data. We demonstrate that we are able to recover flows and sound-speed perturbations in the near-surface layers. We have taken the advantage of the sensitivity of our methodology to entire vertical velocity, and not only to its variations as in other available methodologies. The cross-talk from both the vertical flow component and the sound-speed perturbation has only a negligible effect for inversions for the horizontal flow components. Furthermore, this cross-talk can be minimised if needed. The inversions for the vertical component of the vector flows or for the sound-speed perturbations are affected by the cross-talk from the horizontal components, which needs to be minimised in order to provide valid results. It seems that there is a nearly constant cross-talk between the vertical component of the vector flows and the sound-speed perturbations.

††offprints: David Korda,

1 Introduction

The internal structure of the Sun and the dynamics of plasmas therein are key factors in understanding the solar activity, which in turn influences the interplanetary space (including effects on human infrastructures). Even though processes in the interior of the Sun may be modelled using numerical codes (e.g. Brun et al. 2004; Rempel et al. 2009), these models must be constrained by observations to ensure that they describe a physical reality to some extent.

Helioseismology is the only method that allows us to study subsurface layers of the Sun; these layers, at least to some extent, reflect mostly the properties of the convection zone. Knowledge of the structure and dynamics of plasmas in these regions allows us to put important constraints on theories of the solar dynamo and of the formation and evolution of magnetic fields, including the models of sunspots, and others. The models and inferences from helioseismology do not always agree. For instance, the estimate of convective velocity in the convection zone from helioseismology is two orders of magnitude less than that derived from theoretical studies (Hanasoge et al. 2012). These results, however, were not confirmed by an independent investigation by Greer et al. (2015) when using a different helioseismic analysis. Also, from theoretical considerations (Miesch et al. 2012) the results of Hanasoge et al. (2012) do not fit into the current view of convection in the Sun. On the other hand, these discrepancies have not been explained until now. Another surprising result is that there are helioseismic indications that the supergranules may have a different depth structure than predicted by theory (Duvall & Hanasoge 2013; Švanda 2012; Duvall et al. 2014). Even though both mentioned results are considered to be controversial, there is an obvious need for both theory and helioseismology to reconcile these existing discrepancies.

Helioseismology is used to study the solar interior via interpretation of the dispersion relations of acoustic and surface gravity waves. These waves are generated randomly via the vigorous convection in the convection zone. The waves manifest themselves in many different observables, such as intensity maps of the photosphere or Doppler shifts of photospheric spectral lines. The focus of global helioseismology is to infer the background structure of the solar interior. Local helioseismology focuses on the detection and quantification of small localised perturbations of plasma parameters with respect to the reference solar model.

In this study, we focused on the second approach, using the time–distance method (Duvall et al. 1993) to study the perturbations of travel times of waves caused by inhomogeneities of plasma parameters about the reference model. The time–distance method has been a standard method of solar research for several decades. A comprehensive review of methods of local helioseismology, including the time–distance method and results obtained by this method was presented for example in Gizon et al. (2010) or in Kosovichev & Zhao (2016). The time–distance helioseismology has been used in many applications: for inversions for the solar rotation (Schou 1997), meridional circulation (Zhao & Kosovichev 2004), perturbations of density, pressure, temperature, sound speed (Tong et al. 2003; Brüggen & Spruit 2000), and plasma flows (e.g. Švanda 2013).

The usual task of time–distance helioseismology is to invert for the horizontal components and/or spatial variations of the vertical component of the vector of plasma flows, or the sound-speed perturbations separately, using different sets of travel-time measurements. The separate inversions are justified under the assumption of independent action of various perturbers to the wave travel times. In the inverse modelling, however, it is not certain that this assumption holds. Also, in the realistic models of solar plasmas, the perturbers are correlated. For instance, in the mass-conserving convective flows, the upflows (in the vertical direction) in the cell centres are accompanied by outflows (in the horizontal direction). Thus, these two quantities are naturally correlated. This correlation between horizontal and vertical velocities was convincingly shown by Švanda et al. (2011); see their Fig. 8. A similar correlation plot between the vertical velocity and the sound-speed perturbation is given in Fig. 1, where a large correlation (correlation coefficient of 0.48 at 0.5 Mm depth) between these two quantities is seen. The correlation coefficient was calculated from a realistic solar convection simulation provided by Rempel (2014). This correlation has its origin in the fact that the large upflow velocities are usually located in the interior of the convective cells, where the temperature is greater compared to the surroundings. The greater temperature implies an increase in the sound speed, which naturally follows from the ideal-gas equation of state.

In time–distance helioseismology, the travel times are usually averaged over the annulus around the cell centre, which, together with the correlation described above, leads to a leakage – the cross-talk – of one perturber into the inversion for another one. Separate inversions do not allow one to quantify the possible cross-talk contributions. Information about the cross-talk is extremely important for proper interpretation of results, especially for the vertical velocity (Švanda et al. 2011). This is because the amplitudes of variations of the vertical velocity are less than 10 m s*-1*, but the amplitudes of the horizontal velocities are on the order of 100 m s*-1*. Even a small cross-talk can absolutely devalue the results.

Our aim is to combine inversion for the vector flows $\boldsymbol{v}=(v_{x},v_{y},v_{z})$, where $(x,y,z)$ are local Cartesian coordinates, and for the sound-speed perturbations $\delta c_{s}$. This approach will allow us to study, for the first time, the cross-talk between the flow perturbations and the sound-speed perturbations. The expected amplitude of the sound-speed perturbations known from the numeric simulations is approximately 50 m s*-1*, which is comparable with the amplitude of the vertical velocity, where the cross-talk has a very important effect.

Another remaining issue in the description of solar flows is that with the current methodology it is possible to invert only for the spatial variations of the vertical velocity. Due to the symmetries of the sensitivity kernels used in the flow inversions, the mean value of the vertical component vanishes. The improved methodology we propose in this paper, that is, merging inversions for the 3D flows and the sound-speed perturbations into one inversion, allows us to fully describe the vertical flow component.

2 Methodology

2.1 Forward problem

We measure the perturbed travel time $\delta\tau$ by cross-correlating the signal (the Doppler shift of a photospheric spectral line) in a point at a horizontal location $\boldsymbol{r}$ with the signal averaged over the surrounding annulus with the radius $\Delta$ ; we use the approach of Gizon & Birch (2004). For each $\Delta$, four different geometries are introduced by different weightings in the azimuthal angle (three difference types and one mean type). We use the outflow-inflow geometry denoted by o-i, which is the difference between the travel times of the waves travelling from the annulus centre to its rim and the waves travelling in the opposite direction. The east-west and north-south, denoted by e-w and n-s, are computed similarly, only the signal over the annulus is weighted by either cosine (e-w) or sine (n-s) of the polar angle about the central point. Thus, the measurement is sensitive to the waves travelling in one of these directions. The last one is the mean geometry, which averages travel times of the waves travelling from the centre to the rim of the annulus and of the waves travelling in the opposite direction. In the inversions seen in the literature, the difference geometries were used only to invert for flows, whereas the mean geometry was used to invert for the sound speed.

The perturbed travel time $\delta\tau$ can theoretically be computed from the model as

[TABLE]

where $\boldsymbol{r}$ and $\boldsymbol{r}^{\prime}$ are the horizontal positions, $z$ is the vertical position, and $K^{a}_{\beta}$ are the sensitivity kernels from the forward modelling (Gizon & Birch 2004; Birch & Gizon 2007). Here we use Model S (Christensen-Dalsgaard et al. 1996) as a reference background model. The $\delta q_{\beta}$ are the perturbations of the physical quantities (such as flows, sound-speed, density, etc.; indexed by $\beta$), $P$ is the number of these quantities, and $n^{a}$ is a realisation of the random noise (Jackiewicz et al. 2012). The superscript $a$ corresponds to the measurements of individual travel times and associates the selection of wave filters and averaging geometries as described above, including a selection of the radius of the annulus.

2.2 Inverse problem

The goal of an inverse problem is to invert Eq. (1) and compute

[TABLE]

at a position $\left(\boldsymbol{r}_{0};z_{0}\right)$, where $\alpha$ indicates the perturber (for instance the vertical flow velocity). For simplicity we use a subscript $\alpha$ to indicate the quantity in the direction of the inversion and a subscript $\beta$ as an index indicating a general physical quantity component. The above illustrated inversion is generally not possible because of high levels of noise. The time–distance helioseismology methods are used to derive an estimate of the given quantity, denoted as $\delta q_{\alpha}^{\mathrm{inv}}\left(\boldsymbol{r}_{0};z_{0}\right)$.

In our methodology we focus on a method called Optimally Localised Averaging (OLA), which corresponds to the method of an approximate inverse in mathematical terminology (Louis & Maass 1990). The OLA method (Backus & Gilbert 1968) was originally developed for Earth seismology but can also be used for helioseismology (Pijpers & Thompson 1992).

The SOLA method

For our study we use an improved version of the Subtractive OLA method (SOLA) described by Jackiewicz et al. (2012). The SOLA searches for $\delta q_{\alpha}^{\mathrm{inv}}$ at a given position $\left(\boldsymbol{r}_{0};z_{0}\right)$ in the form of a linear combination of the travel times and the unknown weight functions $w^{\alpha}$:

[TABLE]

where $N$ is the total number of the horizontal positions and $M$ is the number of the travel-time geometries. The functions $w^{\alpha}$ minimise the $\chi^{2}$ in the form

[TABLE]

We use $\mathcal{T}^{\alpha}_{\beta}=\mathcal{T}\delta_{\beta}^{\alpha}$, where $\mathcal{T}$ is a user-selected target function in the direction of the inversion and $\delta_{\beta}^{\alpha}$ is a Kronecker delta. The target function is localised around the area of investigation. $\Lambda^{ab}\left(\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\right)=\mathrm{cov}\left[n^{a}\left(\boldsymbol{r}_{i}\right),n^{b}\left(\boldsymbol{r}_{j}\right)\right]$ describes a noise covariance matrix between observations denoted by $a$ and $b$ in positions denoted by $r_{i}$ and $r_{j}$. Parameters $\mu$, $\nu$, and $\epsilon$ are user-selected trade-off parameters. The quantity $\mathcal{K}^{\alpha}_{\beta}\left(\boldsymbol{r}^{\prime},z;z_{0}\right)$ is defined as

[TABLE]

and is is referred to as the averaging kernel. It quantifies the level of smearing of the real quantity $\delta q_{\beta}$. By using the averaging kernels, the real quantity and the estimate from the inversion are linked by

[TABLE]

The first term of Eq. (4) is the misfit term between the user-selected localisation of the inverted quantity $\delta q_{\alpha}^{\mathrm{inv}}$ and the real localisation given by $\mathcal{K}^{\alpha}_{\beta}$. The second term is the regularisation of the propagation of the random noise. The third term is the regularisation of the cross-talk between the inverted quantity and the other quantities and its meaning is obvious from Eq. (6), where the quantity in the direction of the inversion and also other quantities contribute to the inverted estimate. The fourth term minimises the spatial spread of the weight functions and the last term is the normalisation of the averaging kernels added by the Lagrange multipliers $\lambda^{\beta}$ (for details see Švanda et al. 2011).

The minimisation of Eq. (4) with respect to $w^{\alpha}$ and $\lambda^{\beta}$ gives the solution for the weights. We assume the spatial homogeneity of the background model and therefore the problem is solved in the Fourier domain. There, it decouples and allows us to split a large realistic real-space problem into a series of small Fourier-space problems. This approach permits a realistic inverse problem to be solved even with a desktop computer. One can find the Fourier-space formulation in Jackiewicz et al. (2012) or Švanda et al. (2011).

Noise matrices

We use data-driven estimates for the noise matrices $\Lambda^{ab}$. They were computed by using a $1/T$ fitting introduced by Gizon & Birch (2004), where $T$ is the travel-time averaging time. We used a large set of travel-time maps measured for 31 days by the HMI/SDO instrument (Scherrer et al. 2012; Schou et al. 2012) in the quiet Sun period from August 16, 2010, to September 15, 2010. Following Fournier et al. (2014) we fitted covariances of the travel times $\delta\tau^{a}$ and $\delta\tau^{b}$ by using the function $\mathrm{cov}\left[\delta\tau^{a},\delta\tau^{b}\right]=c+d/T$ for $T=\{2,3,4,6,8,12,24\}$ hours, where $c$ and $d$ are the fitted parameters. The constant $d$ is then the element of $\Lambda^{ab}$. We also tried to fit the functions $d/T+e/T^{2}$ and $c+d/T+e/T^{2}$ with the fitting parameters $c$, $d$, and $e$ as suggested by Fournier et al. (2014), but the results were not as expected. For instance, the diagonal terms $\Lambda^{aa}$ (the RMS of the travel times) were overly smeared and did not correspond to the measured RMS values.

The full covariance matrix is important for getting a reasonable estimate of the variation $\sigma^{2}_{\alpha}$ of the inverted quantity $\delta q_{\alpha}^{\mathrm{inv}}$, which equals

[TABLE]

3 Improvements

The methodology described above is universal and has been used in various inversions. For instance, it was used to study the vector flows (Švanda et al. 2011, and following works) and was found useful for sound-speed inversions (Jackiewicz et al. 2012). These two types of inversions differ in the selection of the appropriate sensitivity kernels in the inverse problem and hence the corresponding travel-time measurements. For the flows the difference geometries were used (i.e. the o-i, e-w, and n-s geometries), whereas for the sound-speed inversions the mean geometry was used. This means that the subscript $\beta$ in Eq. (1) incorporated either the vector of the flow only with $\delta q_{\beta}=\left(v_{x},v_{y},v_{z}\right)$ or the fractional sound-speed perturbations only with $\delta q_{\beta}=\left(\delta c_{s}^{2}/c_{s}^{2}\right)$.

This is in agreement with recent findings by Burston et al. (2015) where the authors show in their Table 1 that the horizontal flow components $v_{x}$ and $v_{y}$ are largely sensitive to the difference types of the travel times (o-i, e-w, and n-s), whereas $\delta c_{s}$ is sensitive to the mean type. The vertical flow component $v_{z}$ is sensitive to both types however; to invert for the full $v_{z}$ and not only its spatial variations, the mean-type geometry needs to be included.

Inspired by such findings, we modified the working methodology so that the pipeline can be used to perform inversions for the vector flows and the sound-speed perturbations at once. In effect, the subscript $\beta$ in Eq. (1) now enumerates all components of the flow and the sound-speed perturbations at once, $\beta=\left(x,y,z,s\right)$. Thus, the vector $\delta q_{\beta}=\left(v_{x},v_{y},v_{z},\delta c_{s}\right)$.

Consequently, the travel-time geometries contain all four averaging geometries, the radii of the averaging annuli, and mode filtering – all of these choices are covered in the aggregate superscript $a$ in the equations.

This approach allows us to invert for the complete vector of the plasma flows and the sound-speed perturbations in one inversion, which has not been done before.

We note that variables were selected after consideration of extensive testing of four possible constructions of the vector $\delta q_{\beta}$:

The first one was $\left(\boldsymbol{v}/c_{s},\delta c_{s}/c_{s}\right)$. There is a problem with the interpretation of results, because in order to convert the results into physical units, one has to multiply them by the background sound-speed profile, smeared by the averaging kernel. The averaging kernels for the inverted horizontal velocities scaled by the sound speed had a deeper secondary lobe that could not be removed by varying the trade-off parameters, which made the interpretation of results difficult. Because of the deeper lobe the horizontal velocities were not localised around the target depth. The noise estimates targeted by the inversion were biased. 2. 2.

Another option was $\left(\boldsymbol{v},\delta c_{s}\right)$, which we took as the preferred one. This option has many advantages: The quantities are in physical units. The maps of the quantities, the cross-talk, and the localisation of $\delta q_{\alpha}^{\mathrm{inv}}$ are easy to interpret. The averaging kernels do not have the secondary lobe like in the previous case and are well-localised around the target depth. 3. 3,4.

The other options were $\left(\boldsymbol{v},\delta c^{2}_{s}/c_{s}\right)$ and $\left(\boldsymbol{v}/c_{s},\delta c^{2}_{s}/c^{2}_{s}\right)$. If we selected one of these options we would have to convert $\delta c^{2}_{s}/c_{s}$ or $\boldsymbol{v}/c_{s}$ and $\delta c^{2}_{s}/c^{2}_{s}$ to the physical units similarly to option #1.

The travel-time kernel codes111We use the Kc3 code provided by Aaron Birch. usually compute the kernels for the fractional squared perturbation of the sound speed, that is, for $\delta c^{2}_{s}/c^{2}_{s}$. Instead of modifying and extensively testing the kernel code, we rather compute the kernels for linear $\delta c_{s}$ from the fractional squared perturbations using the following approach.

We write the squared sound-speed perturbation as

[TABLE]

where $c_{\mathrm{model}}$ is the background sound-speed value. The second term has a negligible magnitude compared to the first term. We then compute $\delta c_{s}$ sensitivity kernels in the form

[TABLE]

where we neglected $\delta c^{2}_{s}$ against $c^{2}_{\mathrm{model}}$ (the uncertainty is less than 1 %). Finally, we obtain the relation for the sensitivity kernel for $\delta c_{s}$:

[TABLE]

4 Validation using synthetic data

4.1 Synthetic travel times

We validate our new methodology using synthetic data. A snapshot from a box222Available from http://download.hao.ucar.edu/pub/rempel/sunspot_models/Helioseismology/quiet_sun_98x98x18Mm_64x64x32km/. of realistic simulation of the convection zone (Rempel 2014; DeGrave et al. 2014) was used. It naturally contains both the complete vector of the plasma flows and the sound-speed perturbations.

The synthetic forward-modelled travel times were computed using Eq. (1). We further add a realistic representation of the random noise to the travel times. Having these two separate components, we were able to study various signal-to-noise-ratio situations.

The realisation of the random noise was computed in the Fourier space using a generator of multivariate normal random numbers using the covariance $\tilde{\Lambda}^{ab}(\boldsymbol{k})$ for each wave vector $\boldsymbol{k}$ separately, where $\tilde{\Lambda}^{ab}$ is the Fourier transform of $\Lambda^{ab}$. We forced the expected Fourier-space symmetries so that in the real space the noise realisation had only a real component as expected. By running a large number (10 000) of realisations and computing the noise covariance matrix again we verified that the random noise realisations have the same statistical distribution as the corresponding $\Lambda^{ab}$.

4.2 Comparison with the model

We applied our improved pipeline to the synthetic travel times and compared the results to the known inputs. By using the synthetic data, we have a complete understanding of the contributions of the individual components of the inversions, that is, all components of the averaging kernels, the noise, and also the expected ideal answer. Thanks to the use of the synthetic data, all these components may be validated separately for each quantity $\delta q_{\alpha}=\{v_{x},v_{y},v_{z},\delta c_{s}\}$. The inversion components to be validated are

[TABLE]

The ideal answer is what we expect to obtain by applying the inversion. The deviations between the ideal answer and our result are caused by the random noise and the other-than-expected localisation, described by the 4D averaging kernel.

For this particular study, we used only the $f$-mode ridge-filtered travel times. The reason for this is the strong localisation towards the surface of the Sun, which means that the results, if applied to the real data, can directly be confronted with the results obtained by other methods (similarly to e.g. Švanda et al. 2013). The disadvantage of using the $f$ mode is that it is not particularly sensitive to the sound-speed perturbations (Burston et al. 2015). According to the same paper, the acoustic $p$ modes should have an order-of-magnitude better sensitivity to the sound-speed perturbations. To confirm this we performed equivalent inversions using the $p_{1}$ mode. These results are discussed in Subsection 7.

4.3 The target function

In the SOLA inversions, the user is free to select the target function ${\cal T}$, which represents the first-order estimate of the averaging kernel. We describe the results of the comparison with a Gaussian-damped-Gaussian-type (GdG) target function, which consists of the 3D Gaussian with the horizontal full width at half maximum (FWHM) of 15 Mm and the vertical FWHM of 0.5 Mm. The centre of the target function is placed at a depth of 0.5 Mm. In the vertical direction, the target function is further multiplied by a smooth function so that ${\cal T}$ reaches zero at the solar surface, $z=0$, and does not extend to the domain of the solar atmosphere.

4.4 Trade-off parameters

The inverted $\delta q_{\alpha}^{\mathrm{inv}}$ depends on a selection of the trade-off parameters $\mu$, $\nu$, and $\epsilon$. We selected the trade-off parameters in the following manner: At first we set an upper limit of the level of the random noise, according to the inverted quantity (20 m s*-1* for $v_{x}$ and $v_{y}$, 15 m s*-1* for $\delta c_{s}$, 5 m s*-1* for $v_{z}$). This gave us a lower limit of the parameter $\mu$.

In the second step we plotted the averaging kernels and cuts through the weights for all combinations of the trade-off parameters which fulfilled the noise constraint. Our further aims were to minimise the cross-talk which is regularised by the parameter $\nu$. The parameter $\epsilon$ controls the spatial spread of weights to ensure that the inversion weights are spatially confined and to avoid ringing solutions. Larger values of the trade-off parameters are usually balanced by a worse fit of the target function by the averaging kernel. The optimal combination of trade-off parameters is chosen on the trial-and-error basis. We note that there is an L-curve (e.g. Hansen 1999) method which attempts to determine the values of the trade-off parameters in an exact way. Unfortunately, in our case this method usually leads to a large level of random noise and therefore we are not using it.

4.5 The method of evaluation

In the following sections we discuss the inversions for the horizontal and vertical flow and the sound-speed perturbations separately. When applicable, we compare the new results to the outputs of the Švanda et al. (2011) pipeline.

For all perturbations we always discuss the shape of the resulting averaging kernel to estimate the localisation in the Sun, the smearing, and the cross-talk contributions. In the plots we show the horizontal slices through the target depth and also the vertical slices perpendicular to the symmetry axis. We then describe the individual contributions to the inverted quantity and discuss the fulfilment of the inversion requirements, such as the spatial localisation of the weights.

5 Inversions for horizontal components of vector flows

Inversions for horizontal (that is $v_{x}$ and $v_{y}$) components of the vector flows have been typical tasks for the time-distance helioseismology. There is only a small amount of cross-talk with the vertical $v_{z}$ caused by natural correlations between the vertical upflow and the outflows in the mass-conserving flow. The averaging kernels are plotted in Fig. 2 (cross-talk minimised) and Fig. 3 (cross-talk not minimised). (We note that only plots for $v_{x}$ are given, plots for $v_{y}$ are similar except for the rotation by 90 degrees around the vertical axis.) By comparing the second, third, and fourth columns of the two figures, one can see that the cross-talk minimisation works as expected, when these contributions vanish in Fig. 2. The cross-talk minimisation is balanced by the worse fit of the component of the averaging kernel in the direction of the inversion to the target function. The averaging kernel ${\cal K}_{x}^{x}$ is more extended in the $x$ direction, reaching essentially an elliptical shape in the horizontal cut, whereas the target function is roundish. The use of the $f$-mode sensitivity kernels only also contributed to this issue because they do not represent the GdG function well.

In Fig. 4 one can see the individual contributions from the different quantities and the noise to both inversions (with and without the minimisation of the cross-talk). There is also a comparison with the ideal answer. Despite the appearance of the averaging kernel with cross-talk not being minimised, the total leakage from $v_{z}$ and $\delta c_{s}$ is negligible because the cross-talk components of the averaging kernel have a rather small extent. Moreover, the amplitudes of the leaking components are at least an order of magnitude smaller in the near-surface layers.

Due to the fact that the amplitudes of the horizontal flows are at least an order of magnitude larger than the amplitudes of the two remaining perturbers, there is no need to use the improved methodology for such inversions. However, it might be unnecessary to minimise the cross-talk at all. Such a claim will, however, lose its validity in the layers or regions, where the horizontal and vertical components reach a comparable magnitude.

6 Inversion for the vertical flow component

Until now, only variations of the vertical component $v_{z}$ were known. This is because only the difference travel-time geometries were used to invert for $v_{z}$ and the total horizontal integral of the sensitivity kernels vanished due to their symmetries (Burston et al. 2015). Therefore, such measurements are not sensitive to the horizontal average (the mean) $\langle v_{z}\rangle$ which is therefore identically zero in the inverted estimates. We only note that this property also requires a modification of the target function for such an inversion (Švanda et al. 2011) in order to fulfil the normalisation constraint in Eq. (4).

According to Burston et al. (2015) the mean travel-time geometry is sensitive to $\langle v_{z}\rangle$. Therefore, by adding the mean geometry we can invert for full $v_{z}$. From the previous works (e.g. Švanda et al. 2011) it is known that there is a large cross-talk with the horizontal flow components, which must be minimised in order to have credible results.

The inversion averaging kernels are plotted in Figs. 5 (cross-talk minimised) and 6 (cross-talk unconstrained). The cross-talk minimisation approach allows us to minimise the leakage from the horizontal flow components, however, it has only a limited effect on the leakage of the sound-speed perturbations. We found that this cross-talk with $\delta c_{s}$ has the origin in the contribution of the mean travel-time kernels. To estimate the severity of this cross-talk, we plotted the mean-geometry contribution of the four physical quantities to the averaging kernel. As one can see from Fig. 7, the mean contribution comes from the ${\cal K}^{z}_{z}$ kernel, with a small contribution from the ${\cal K}^{z}_{s}$ kernel. There is no contribution from the horizontal flow, because the mean travel-time sensitivity kernels for the horizontal components have vanishing horizontal integrals. From this plot we estimate that the leakage from the sound-speed perturbation to the inversion for the vertical flow is very small.

The effect of the minimisation of the cross-talk with the horizontal component is clearly demonstrated in Fig. 8, where individual inversion components are plotted. The cross-talk contributions seen in the bottom row have a significant impact on the sum and their value in the case of the horizontal-components leakage is higher than the typical variations of $\hat{v}_{z,z}^{\mathrm{inv}}$. For a better comparison between the components and the ideal answer, in Fig. 9 we plot $v_{z}^{\mathrm{inv}}$ with and without the minimisation of the cross-talk, $\hat{v}_{z,z}^{\mathrm{inv}}$, and the ideal answer. These plots show that the cross-talk minimisation indeed improves the inversion results. The correlation coefficients between the inverted value and the ideal answer is 0.36 when the cross-talk is ignored and 0.57 when it is minimised. By comparing our new results to the results of Švanda et al. (2011) we can say that even the supplement of the mean travel times in the vertical flow inversion helped to lower the cross-talk, because without them, the cross-talk unconstrained inversions ended with a correlation coefficient of close to $-1$ with the expected answer (also Zhao et al. 2007).

Obviously, the inversion scales the values by a factor of 0.44, which is due to a large misfit between the target function and the averaging kernel. Such an issue may be understood as being due to the fact that the inferred vertical flows are smoother than those expected from the setup. Such a result highlights the importance of having the information about the averaging kernel when interpreting the inversions.

The incorporation of the mean travel times allowed us to invert for the full vertical velocity and not only its horizontal variations. To demonstrate this, we additionally inverted for the vertical velocity using the old pipeline of Švanda et al. (2011) using only the difference travel-time kernels. It is necessary to modify the target function so that it has a vanishing horizontal mean, which is achieved by introducing a wide negative horizontal side lobe. The results shown in Fig. 10 demonstrate that the expected ideal answers obtained by using two different target functions are essentially indistinguishable, except for the systematic offset. By subtracting the two maps (seen in the third row of Fig. 10) we see that this offset actually has a structure with a very small variation compared to the actual variations of the ideal answer. The offset represents a large-scale structure of the vertical flow, which is not detectable using the Švanda et al. (2011) approach.

When the mean travel times are incorporated into the inversion, the overall structure of the inferred vertical flow does not change. The correlation coefficient of the improved vertical flow inversion and the inversion by the original pipeline is 0.91. However the large-scale offset is successfully measured by the improved methodology as one can see from the differences in the bottom row of Fig. 10. Except for the lowered amplitude (due to the imperfect averaging kernel), its structure is measured reasonably well.

7 Inversion for sound-speed perturbations

An example of the inversion for the sound-speed perturbations using principally the same methodology we use in this study was published by Jackiewicz et al. (2012). Their inversion used only the mean travel-time sensitivity kernels and thus did not take into account possible leakage from other quantities. We have the opportunity to do this now.

The sound-speed inversion averaging kernels are plotted in Figs. 11 (with the cross-talk minimisation) and 12 (without the cross-talk minimisation). The fits in the direction of the inversion (the right-most panels) are not very representative of the imposed GdG target function, which is the consequence of using the $f$ mode only. On the other hand, the cross-talk minimisation performs well by essentially setting the cross-talk components of the averaging kernels to zero. For that one has to expect a worse fit of the target function in the direction of the inversion, where even an extended small-amplitude surrounding annulus appears. We assume that the fit improves when more independent travel-time measurements are used in the inversion because of a larger number of independent observations. A wider averaging kernel causes a larger-than-expected smoothing of the results, which is demonstrated by the smaller amplitude of the inverted sound-speed perturbations compared to the ideal answer.

In Fig. 13 we plot individual contributions to $\delta c_{s}^{\mathrm{inv}}$. In both cases the inverted sound-speed-perturbations maps are dominated by the component in the direction of the inversion, which has the RMS value of 31 m s*-1*. The RMS of the total cross-talk is 16 m s*-1* when not minimised. The cross-talk reaches half of the magnitude of the desired signal, which we consider too large. After the minimisation, the cross-talk contribution lowers to only less than 4 m s*-1*. The cross-talk is positively correlated with the inverted value, that is, without its minimisation, we overestimate the real sound-speed perturbations by a factor of about 1.9.

In Fig. 14 one can compare $\delta c_{s}^{\mathrm{inv}}$ with and without the minimisation of the cross-talk, $\delta\hat{c}_{s,s}^{\mathrm{inv}}$, and the ideal answer. The correlation coefficient between $\delta\hat{c}_{s,s}^{\mathrm{inv}}$ and the ideal answer is 0.90. The difference between $\delta c_{s}^{\mathrm{inv}}$ (top left) and the ideal answer is mainly caused by the misfit because the sensitivity of the $f$ mode to $\delta c_{s}$ is weak.

Validation with the $p_{1}$ mode

According to Burston et al. (2015), the sensitivity of the acoustic modes is at least one order of magnitude higher than that of the $f$ mode. Therefore we ran the sound-speed inversion again, using only the kernels for the $p_{1}$ ridge. The first acoustic ridge samples similar depths to the $f$ mode, and therefore the inversions are directly comparable.

Furthermore, in the inversions presented above, we imposed a certain shape (the “GdG”) of the target function. We showed that in some cases the inverted results were far from optimal, which was due to the fact that the $f$-mode travel-time sensitivity kernels were not dependent enough on depth to fit the GdG target function. For this reason we recomputed the inversion with a target function that has a more natural profile in the vertical domain.

In our methodology, such a target function (called Arbitrary-z, Az) resembles the natural depth profile of the sensitivity kernels used in the inversion by representing their kinetic-energy profile. It is usually obtained by simply averaging the horizontal integrals as a function of depth of all sensitivity kernels used in the given inversion. Using the Az-class target function our inversion is almost equivalent to the 2D inversion and thus it is expected that the resulting averaging kernels will better fit this target function, while having a lower level of random noise at the same time. In the horizontal direction the target function is again a Gaussian with a FWHM of 15 Mm.

Figure 15 (to be compared with Fig. 11) demonstrates that the representation of the target function by the averaging kernel is indeed much better in the case of inversion using the Az-type target function and the $p_{1}$ mode. A similar conclusion may be drawn when comparing Figs. 16 and 12 for the inversions with an ignored cross-talk. The superiority of the $p_{1}$ inversion for $\delta c_{s}$ is even more visible when comparing the horizontally averaged averaging kernels plotted as a function of depth (Fig. 17).

8 Summary

The results shown in the figures presented here can also be represented by numbers. In Table 1 we give some statistical properties of the inversions for all physical quantities in question, that is, for three components of the flow velocity and the sound-speed perturbations.

For each quantity separately (by columns) we give the correlation coefficient of the inversion component in the direction of the inversion with the ideal (expected) answer, which is always positive and large. The correlation coefficient between the inverted estimate and the ideal answer is given, which is always positive, however, it is lower for the vertical flow and sound-speed perturbations. That is mainly due to the presence of the random noise. Finally, the RMS values of the component in the direction of the inversion, of the total cross-talk contribution, and of the random noise help to evaluate the signal-to-noise ratio of the inverted estimates and the significance of the cross-talk pollution. The comparison between the upper and lower tables demonstrates the expected performance of the cross-talk minimisation.

In Section 7 we performed an inversion for the sound-speed perturbations with a different target function and a different set of sensitivity kernels to improve the inversion. The improvements of such an approach are also apparent in Table 2 (that should be compared to the last column in Table 1). The correlation coefficient of the inverted estimate is much larger than before, the random-noise level is smaller, and also the cross-talk pollution is smaller even in the case when it is not minimised.

Moreover, in Fig. 18 we plot examples of the inversion weights for $\Delta\approx 19$ Mm. The weights are well localised around the central point.

9 Conclusions

We introduced the improved methodology of the SOLA time–distance inversions for local helioseismology. Our improvements consist of incorporating the mean travel-time geometry in the inversions and also of the use of the travel-time sensitivity kernels for all flow components and the sound-speed perturbations at once. This new methodology allows us to invert, for the first time, for the full value of the vertical velocity (not only for its variations as before) and to quantify the cross-talk between the vector flows and the sound-speed perturbations.

We validated our inversions by using the synthetic data based on the realistic numerical simulation of the solar convection provided by Matthias Rempel (Rempel 2014; DeGrave et al. 2014). The use of the synthetic data allowed us to properly quantify all the inversion components and to gain empirical experience of the use of the improved pipeline. To summarise the results:

• •

For $v_{x}$ and $v_{y}$ (the horizontal flow) it is not necessary to minimise the cross-talk (at least near the surface). Our methodology essentially does not constitute an improvement to what has been used recently by other authors. We only point out that the targeted RMS of the random noise should be 20 m s*-1* or less if one inverts the near-surface flows using travel times averaged over a day or so.

• •

For $v_{z}$ (the vertical flow) it is essential to minimise the cross-talk, otherwise the results are nullified by the leakage mainly from the horizontal components. The leakage from the sound-speed perturbations is acceptably small. Our findings confirm the results seen in the past. The RMS of the random noise should be lower than 2 m s*-1* for the near-surface flows (again, for one-day averaged travel times). An obvious advantage of our new methodology is that one can also quantify the large-scale (with scales larger than the typical horizontal extent of the averaging kernel) component of the vertical flow and not only its horizontal variations as before. The inversion is sensitive also to the large-scale offset because of the incorporation of the mean travel-time kernels. Also, it seems that the incorporation of the mean travel times helps to regularise the cross-talk contribution. By using the improved methodology, we achieved significant improvements in inversions for the vertical flow.

• •

At first sight it seems that the minimisation of the cross-talk in the case of the inversion for $\delta c_{s}$ does not improve the results. This is because the leakage component from the other quantities is positively correlated with the sound-speed perturbations and thus the leakage “helps” to “measure” $\delta c_{s}$. Such an approach leads to an overestimation of the magnitude of the sound-speed variations by almost a factor of two in the near-surface layers. Such an inversion, however, is not the inversion for $\delta c_{s}$. When the cross-talk minimisation is introduced, the agreement with the ideal answer decreases, which is in part due to the worse fit of the target function by the averaging kernel, and also by the lower signal-to-noise ratio achieved for the travel times averaged over 24 hours. The case we present here cannot be further improved because only the $f$ mode, which is not particularly sensitive to the sound-speed perturbations, was used. The reliability improves by a large extent when the first acoustic ridge is used in a similar setup.

Acknowledgements.

D.K. is supported by the Grant Agency of Charles University under grant No. 532217. M.Š. is supported by the project RVO:67985815 and the grant project 18-06319S awarded by Czech Science Foundation. We are grateful to Matthias Rempel for providing us with the snapshot of the numerical simulation. The sensitivity kernels were computed by the code Kc3 kindly provided by Aaron Birch. The authors would like to thank the anonymous referee for the valuable comments and suggestions, which greatly improved the quality of the paper. This research has made use of NASA Astrophysics Data System Bibliographic Services. Author contributions: M.Š. and D.K. designed the research. D.K. implemented and performed the research within his Ph.D projects under the supervision of M.Š. D.K. drafted the paper and both authors contributed to the final manuscript.