Quantum mechanics and speed limit of ultrafast local control in spin chains

TL;DR

This paper investigates how to optimize the speed and fidelity of ultrafast local control in spin chains, revealing system-independent features and quantum speed limits that can enhance quantum system manipulation.

Contribution

It introduces a method to significantly reduce infidelity in ultrafast spin chain transformations by designing optimal control profiles, supported by analytical predictions and numerical validation.

Findings

Optimal control profiles drastically decrease infidelity.

Qualitative features of control are system-independent.

Restrictions on control shape persist up to adiabatic limits.

Abstract

We study optimization of fidelity for ultrafast transformation of a spin chain via external control of a local exchange coupling. We show that infidelity of such a process can be dramatically decreased by choosing a proper control profile in nonadiabatic time domain, predict main features of this profile analytically, corroborate them numerically with a gradient search algorithm, and discuss the corresponding quantum speed limit. For ultrafast transformations, the qualitative features of the obtained optimal control are system-independent. Moreover, the main restrictions on its shape do not depend on the transformation time and remain valid up to the adiabatic limit. Our results can be applied to control a broad variety of quantum systems.

Tables

Table 1. Table 1: Numerical evidence of linearity of  γ T ​ ( 𝐚 opt ) subscript 𝛾 𝑇 subscript 𝐚 opt \gamma_{T}(\mathbf{a}_{\rm opt}) .

$T$	0.005	0.02	0.05	0.1
$a_{2opt}$	648.3	162.4	65.2	32.84
${\gamma }_T\left({𝐚}_{\mathrm{opt}}\right)$	0.00322	0.0129	0.0322	0.0642
${\gamma }_T\left({𝐚}_{\mathrm{opt}}\right)/T$	0.644	0.644	0.643	0.642

Equations

$$H(g(t))=J\sum _{n=1}^{N-1}{X}_n{X}_{n+1}+B\sum _{n=1}^N{Z}_n+g(t)J{X}_1{X}_N,$$

$$U(t)=\mathcal{T}\exp \left(-i{\int }_0^{t}H(g(s))ds\right),$$

$$f_T\equiv |⟨{\phi }_f⟩\psi (T)|,R_T\equiv 1-f_T.$$

$$R_T\to 0,\mathrm{for}T\to \infty ,$$

$$g(\mathbf{a},t)=\frac{t}{T}+a_1\sin \left(\frac{\pi t}{T}\right)+a_2\sin \left(\frac{2\pi t}{T}\right),$$

$$U(T)$$

$$G(t)$$

$$R_0-R_T\approx |f_0|\alpha (\mathbf{a})T^2,$$

$$U(T)\approx e^{-iVG(T)}\left(\mathbb{I}-i{\int }_0^T{e}^{iVG(t)}{H}_0{e}^{-iVG(t)}dt\right),$$

$$e^{\pm iG(t)X_1{X}_N}=\mathbb{I}\cos G(t)\pm i{X}_1{X}_N\sin G(t).$$

$$\begin{gathered} {\int }_0^T{e}^{iVG(t)}{H}_0{e}^{-iVG(t)}dt={\int }_0^T[H_0-B\left(Z_1+Z_N\right)+ \\ B{e}^{iVG(t)}(Z_1+Z_N)e^{-iVG(t)}]dt=\left(H_0-B(Z_1+Z_N)\right)T+ \\ B(Z_1+Z_N){\beta }_T(\mathbf{a})+B(Y_1{X}_N+X_1{Y}_N){\gamma }_T(\mathbf{a}), \end{gathered}$$

$${\beta }_T(\mathbf{a})\equiv {\int }_0^T\cos (2G(t))dt,{\gamma }_T(\mathbf{a})\equiv {\int }_0^T\sin (2G(t))dt.$$

$${\int }_0^{T/2}\sin (2g_{1}t)dt=\frac{T}{2}\cdot (g_{1}T{)}^{-1}(1-\cos (g_{1}T)),$$

$$\underset{T\to 0}{\lim }|⟨{\phi }_f⟩U_{\mathrm{opt}}(T){\phi }_i|\ge |⟨{\phi }_f⟩{\phi }_i|.$$

$$\underset{T\to 0}{\lim }{a}_{1\mathrm{opt}}=C_1,\underset{T\to 0}{\lim }{a}_{2\mathrm{opt}}=\frac{C_2}{T},$$

$$R_0-R_T\approx |f_0|B\left[\left({\beta }_T(\mathbf{a})-T\right)F_1+{\gamma }_T(\mathbf{a})F_2\right],$$

$$F_1=\mathrm{Im}\frac{f_Z{f}_0^{\ast }}{|f_0{|}^2},F_2=\mathrm{Im}\frac{f_{XY}{f}_0^{\ast }}{|f_0{|}^2},$$

$$f_Z=⟨{\phi }_f⟩Z_1+Z_N{\phi }_i,f_{XY}=⟨{\phi }_f⟩X_1{Y}_N+Y_1{X}_N{\phi }_i.$$

$$\begin{gathered} ⟨{\phi }_f⟩U(T){\phi }_i=⟨{\phi }_f|[\mathbb{I}-iB(Z_1+Z_N)({\beta }_T(\mathbf{a})-T)- \\ iB(Y_1{X}_N+X_1{Y}_N){\gamma }_T(\mathbf{a})]e^{-i{H}_{0}T}|{\phi }_i⟩+\mathcal{O}(T^2). \end{gathered}$$

$$R_0-R_T=|f_0|B{K}_{\gamma }{F}_{2}T.$$

$$\frac{T_{\mathrm{QSL}}}{T}=\frac{\arccos \left|⟨{\phi }_i⟩U_{\mathrm{opt}}(T){\phi }_i\right|}{{\int }_0^T\left|⟨{\phi }_i⟩H(g({\mathbf{a}}_{\mathrm{opt}},t){\phi }_i\right|dt},$$

$$\dot{g}(t\to 0)>0,\dot{g}(t\to T)>0,$$

$$-\left(1/\pi +2a_2\right)<a_1<1/\pi +2a_2.$$

Full text

I Introduction

Recent progress in experimental research on quantum systems described by moderate-size Hilbert spaces, such as ensembles of qubits, posed fascinating problems of optimal quantum control Peirce et al. (1988); Brif et al. (2010); Borzi et al. (2017); Caneva et al. (2011a, 2014) of these systems. The quantum control aims at achieving desired quantum states or certain quantum operations with maximum possible fidelity using limited resources such as time or energy. The dynamics of quantum systems under external control can be unitary or non-unitary. The unitary dynamics is driven by a time-dependent controllable Hamiltonian $H\left({\mathbf{g}}(t)\right)$, where ${\mathbf{g}}(t)$ is a multicomponent control function. The controllable non-unitary dynamics is achievable by system measurements Li et al. (2011); Luchnikov and Filippov (2017); Wu et al. (2004); Pyshkin et al. (2016); Torres et al. (2017), via a controllable interaction with a non-Markovian environment Verstraete et al. (2009); Jing and Yu (2010); Luo et al. (2015) or via control of the unitary part of the evolution of open system Schmidt et al. (2011).

We consider driving a quantum system from a ground state of initial Hamiltonian $H_{i}$ to achieve at time $T$ ground state of a final Hamiltonian $H_{f}$ with $H(t)=H_{i}+g(t)(H_{f}-H_{i}),$ where $g(0)=0$ and $g(T)=1.$ Although a high fidelity can be obtained by an adiabatic process Born and Fock (1928) driven by a slowly varying $H(t)$ with, e.g., $g(t)=t/T$, this method requires a long evolution while optimized $g(t)$ can permit achieving a demanded quantum state for a relatively short $T.$

A possible approach to the quantum control, where the transitions occur between the ground states of $H(t)$, is based on the shortcut to adiabaticity  Torrontegui et al. (2013); Guéry-Odelin et al. (2019); Demirplak and Rice (2003); Berry (2009); Ren et al. (2017). However, this technique requires a control of the all parts of a complex quantum system. Implementation of a such a shortcut can be a part of quantum computation in arrays of quantum dots Loss and DiVincenzo (1998); Burkard et al. (1999), or in quantum annealing Das and Chakrabarti (2008), such as applied in D-Wave computer McGeoch (2014). Here by focusing on high fidelity ultrafast processes, we analytically obtain properties of optimal local control in the ultra-short time domain for a particular many-body system and corroborate our reasoning by a direct numerical optimization. We show that several properties of the finite time quantum control (even for the ultra-short time) can be explained by requiring a smooth passage to the adiabatic protocols, thus, connecting these two limits. Although the reported results are obtained for spin chains, the proposed heuristic reasoning and numerical approach can be extended to a much broader class of quantum systems.

II Ultrafast local control: the problem setting

We concentrate on a local control acting only on a small part of a complex system, being a natural tool for cutting or stitching links between its parts, thus, modifying its size and/or topology. We consider an Ising chain with $N$ spins, as shown in Fig. 1, described by the Hamiltonian

[TABLE]

where $X_{k},Z_{k}$ are corresponding Pauli matrices of the  $k$-th spin and $B$ is a magnetic field. It is useful to rewrite Hamiltonian (1) in a short form: $H(t)=H_{0}+g(t)V,$ where $V\equiv JX_{1}X_{N}.$ We assume antiferromagnetic interaction and set $J=1.$ The last term in (1) connects the first and last spins in the chain. By assuming $g(0)=0$ and $g(T)=1,$ we perform a transformation from an open to a ring-shaped chain via “stitching” a single link between the spins (see Fig. 1). Now we define initial and final Hamiltonians: $H_{i}\equiv H(g(0))$, and $H_{f}\equiv H(g(T))$. The corresponding ground states of these Hamiltonians are $\ket{\varphi_{i}}$ and $\ket{\varphi_{f}}$ : $H_{i}\ket{\varphi_{i}}=\varepsilon_{i}\ket{\varphi_{i}}$, and $H_{f}\ket{\varphi_{f}}=\varepsilon_{f}\ket{\varphi_{f}}$. The state of the system during the evolution is $\ket{\psi(t)}=U(t)\ket{\psi(0)},$ where

[TABLE]

$\mathcal{T}$ is the time-ordering operator, and we set  $\hbar\equiv 1$, and the time unit as $1/J$. We assume that the initial state $\ket{\psi(0)}=\ket{\varphi_{i}}$ and study the controlled state-transition process with the following target fidelity $f_{T}$ and infidelity $R_{T}$:

[TABLE]

Adiabatic theorem allows us to have an ideal state-transition protocol:

[TABLE]

when $g(t)=t/T$. Note that the protocol (4) is valid only in the absence of level crossing for an arbitrary $g(t)\in(0,1)$. In this work we assume that the evolution time belongs to one of three domains: ultrashort one ($T\ll 1$), short time ($T\sim 1$), and adiabatic one ($T\gg 1$), where the adiabatic theorem is valid. We will concentrate mainly on the physics of the optimal control in the ultrashort domain and analyze the general features being common for all three domains. Such ultrafast control can be achieved, e.g., by electrical manipulation of the bonds connecting quantum dots on the time scale much shorter than global change in the magnetic field strength Burkard et al. (1999).

If ground state of $H_{i}$ or $H_{f}$ is degenerate we select $\ket{\varphi_{i}}$ or/and $\ket{\varphi_{f}}$ from some subspace. In such a case we assume that $\ket{\varphi_{i}}$ ($\ket{\varphi_{f}}$) is a non-degenerate ground state of Hamiltonian $H(\delta g)$ ($H(1-\delta g)$) for $\delta g\rightarrow+0$.

Our task is to find the optimal $g(t)$ to minimize the target infidelity functional $R_{T}[g(t)]$ for a finite time $T$. Note that our system (1) doesn’t have complete controllability Schirmer et al. (2002); Ramakrishna et al. (1995) because of locality of our control. Although the local control can, in general, be complete (see Ref. Burgarth et al. (2009)), the Hamiltonian (1) doesn’t satisfy assumptions made in Ref. Burgarth et al. (2009) since $iX_{1}X_{N}$ does not generate a Lie algebra in the subspace of 1st and $N$th spins. Therefore, the result of our optimization is the minimal nonzero $R_{T}.$

In order to deal with a function instead of a functional we parameterize as follows,

[TABLE]

with $\mathbf{a}\equiv(a_{1},a_{2})$. Now the target infidelity is a function of two parameters $R_{T}=R_{T}(a_{1},a_{2})$. Parametrization (5) can be considered as a simple particular case of chopped-random-basis optimization Caneva et al. (2011b); Rach et al. (2015). Our task is to find optimal $\mathbf{a}=\mathbf{a}_{\rm opt}$, with $R_{T}(\mathbf{a}_{\rm opt})=\min_{\mathbf{a}}\{R_{T}(\mathbf{a})\}$. Also we can rewrite statement (4) as $R_{T}\rightarrow 0$ for $T\rightarrow\infty$ and $\mathbf{a}\rightarrow\mathbf{0},$ and therefore we expect that $\mathbf{a}=\mathbf{0}$ is a good starting point for gradient numerical search.

III Optimal control for ultrafast evolution

III.1 Analytical results for singular control parameters

It is possible to use the first two terms of Dyson series for approximation of $U(T)$ when $T\rightarrow 0$. One can write

[TABLE]

and $\mathbb{I}$ is the identity operator. By using $\braket{\varphi_{f}}{V}{\varphi_{i}}=\braket{\varphi_{f}}{H_{f}-H_{i}}{\varphi_{i}}=(\varepsilon_{f}-\varepsilon_{i})f_{0}$, where $f_{0}=\braket{\varphi_{f}}{\varphi_{i}}$, we obtain that this approximation leads to a quadratic $T-$dependence of the target fidelity

[TABLE]

where $R_{0}\equiv 1-|f_{0}|$, and $\alpha$ is a coefficient, e.g., $\alpha(\mathbf{0})=(\varepsilon_{i}-\varepsilon_{f})^{2}/8.$ The quadratic behavior with $dR_{T}/dT|_{T=0}=0$ is an understandable feature of the sudden approximation Sakurai and Napolitano (2010) where the initial state remains almost intact after fast change in the Hamiltonian. However, (7) being valid only for  $|\mathbf{a}|\ll T^{-1},$ provides an inefficient optimization and, therefore, one needs to go beyond this condition.

A more convenient way to go beyond the simple sudden approximation is to apply the following interaction picture:

[TABLE]

with the validity of the integral expression (8) being not explicitly related to the magnitude of $g(t).$ For $V\equiv X_{1}X_{N}$ we simplify matrix exponents in (8) as:

[TABLE]

Since only two terms, that is $BZ_{1}$ and $BZ_{N},$ in $H_{0}$ do not commute with $\exp(\pm iX_{1}X_{N}G(t)),$ by using algebra of Pauli matrices with $(X_{1}X_{N})^{2}=\mathbb{I}$, $X_{k}Z_{k}X_{k}=-Z_{k}$ and $Z_{k}X_{k}=-X_{k}Z_{k}=iY_{k},$ we simplify integral in (8) as:

[TABLE]

where

[TABLE]

Now we analyze the integral terms in (10). If we assume the control is very weak when for any given $t,$ $G(t)\rightarrow 0$, we arrive at Eq. (7). The opposite case of a strong control should be considered in more detail. To get insight into the evolution, we assume for the moment a pulsed control: $g(t)=g_{1}$, for $0<t<T/2$ and $g(t)=g_{2}$, for $T/2<t<T$. In this case we have $G(t)=g_{1}t$, $t<T/2$. Part of the last integral in (10) can be written as

[TABLE]

where we picked out linear proportionality on $T$ as in the other terms in (10). Now we see that the contribution of this integral goes to zero in two limits: 1) $g_{1}T\rightarrow 0$ and 2) $g_{1}T\rightarrow\infty$. Thus, we have reached an important conclusion that for effective control one must have $g_{1}T={\rm const}$ for $T\rightarrow 0$. The same conclusion can be made for $g_{2}$ by analysis of the $T/2<t<T$ interval.

As the next step in our reasoning we require that the optimal control enhances the fidelity, that is:

[TABLE]

In order to satisfy (13) we can require $\lim_{T\rightarrow 0}G(T)=0$ for optimal control (see Eq.(8)). This means that for two pulses one must have $g_{2{\rm opt}}=-g_{1{\rm opt}}$. All these conclusions now can be applied for smooth optimal control function (5) in the following way:

[TABLE]

where constants $C_{1,2}$ to be obtained by numerical calculations based on Eq. (10).

Using result (10) we write the target infidelity as:

[TABLE]

where

[TABLE]

Since without loss of generality, we can assume that the states $\ket{\varphi_{i}}$ and $\ket{\varphi_{f}}$ are real, we obtain $F_{1}=0$, ${\rm Re}F_{2}=0$ and thus our result does not depend on $\beta_{T}(\mathbf{a}).$

Expression (15) was derived using

[TABLE]

Here we also assume $G(T)=0$ for all values of $T$. This assumption means that we chose $a_{1}=-\pi/4$ as optimal value for any small nonzero $T$ 2019-1 (201). It is important that optimized function $\gamma_{T}(\mathbf{a}_{\rm opt})=K_{\gamma}T$ is linear for small $T$, where $\gamma_{T}(\mathbf{a}_{\rm opt})=\max_{\mathbf{a}}\{\gamma_{T}(\mathbf{a})\}$, $K_{\gamma}$ is a system-dependent coefficient, and $\mathbf{a_{{\rm opt}}}$ satisfies (14). Thus, the linear approximation to optimal infidelity is

[TABLE]

Linearly decreasing behavior of infidelity (19) under optimal control gives a big advantage in comparison with quadratic (7) for short time $T$. The spatial symmetry of Hamiltonian (1) assures that $F_{2}$ is an odd function of $B,$ corresponding to the fidelity independent of the direction of the magnetic field.

III.2 Numerical examples

To illustrate the above arguments, we study a chain with $N=6$ spins and $B=0.9$, and relate the results of direct numerical simulations to expressions (7), (15), and (19). Although the Ising chain in a transverse field is exactly solvable Pfeuty (1970), we obtain the states $\ket{\varphi_{i}}$ and $\ket{\varphi_{f}}$ by direct numerical diagonalization of the corresponding Hamiltonians. Next, we use gradient Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm Fletcher (1988) in direct numerical search of $\mathbf{a}_{\rm opt}(T)$. Exact numerical diagonalization was used in order to calculate $\mathbf{a}-$dependent propagators (2). In Fig. 2 we show non-optimized ($g(t)=t/T$) and optimized infidelities obtained by direct numerical simulations. To quantify linear approximation for the optimized fidelity, we first numerically obtain $f_{0}=0.9525$, $f_{Z}=-1.4090$, $f_{XY}=0.389i$, resulting in  $F_{1}=0,$ as expected, and $F_{2}=0.408$. It is easy to numerically find a maximum of a function $\gamma_{T}(\mathbf{a})$ for fixed $T$ and $a_{1}=-\pi/4$; the example of dependence $\gamma_{T}(a_{2})$ for $T=0.005$ is depicted in Fig. 3. From the last line of Table 1 we find $K_{\gamma}=0.644$ and obtain $BK_{\gamma}F_{2}=0.237$. Additional numerical checks show that adding next harmonic in (5) doesn’t considerably change the optimal fidelity value. Note, that expression (19) being linear in $T,$ is not linear in $B$ since $F_{2}$ is $B-$dependent.

In Fig. 4 (main panel) we show numerically optimized values of $a_{2{\rm opt}}$, and approximation $a_{2}(T)=3.24/T,$ which follows from expression (14). The coefficient $3.24$ can be obtained by taking product $Ta_{2{\rm opt}}$ from the first column of Table 1. We see a good agreement between analytical and exact numerical optimization up to $T\approx 0.1$. In the inset of Fig. 4 we show numerically optimized values of $a_{1{\rm opt}}(T)$ and see that $a_{1}$ remains finite in the limit $T\rightarrow 0$ in agreement with (14), and thus our simplification $G(T)=0$ for $T\rightarrow 0$ in (15) was rational. Note, in Fig. 4 we do not have $\lim_{T\rightarrow 0}a_{1{\rm opt}}=-\pi/4$ as we used in our analysis. However, the approximation we made is valid because $G(T)=T(1/2+2a_{1{\rm opt}}/\pi)\ll 1$ for any finite $|a_{1{\rm opt}}|$ and $T\rightarrow 0$.

Important input of (19) is that the optimal control parameters $\mathbf{a}_{\rm opt}$ can be evaluated from analysis of function $\gamma_{T}(\mathbf{a})$ (11), which, in turn, does not contain information about $N$ and $B$. Therefore we expect that the optimal control $g(t)$ for different quantum Ising chains is only very weakly system-dependent in the ultrashort time limit. This universality is confirmed by Fig. 5 (a) where we present the BFGS optimization results for chains with different parameters.

III.3 Relation to quantum speed limit

The optimal control (14) produces a strong perturbation, where the characteristic energy given by the time-energy uncertainty Mandelstam and Tamm (1991) is proportional to $T^{-1}.$ Therefore, it is instructive to consider the quantum speed limit (QSL) time $T_{\rm QSL},$ that is the minimal possible time required to transform the initial $|\varphi_{i}\rangle$ into the final $U_{\rm opt}(T)|\varphi_{i}\rangle$ state Margolus and Levitin (1998). The $T_{\rm QSL}$ time is computed as Deffner and Lutz (2013):

[TABLE]

and the ratio $T/T_{\rm QSL}$ can be considered as the efficiency of the quantum control (see, e.g., Murphy et al. (2010); Caneva et al. (2009)). In Fig. 5(b), where we show $T-$dependence of $T/T_{\rm QSL}$, one can see the ratio $T/T_{\rm QSL}\gtrsim 5$, and this relatively large value can be related to the locality of the control.

IV From ultrafast to adiabatic evolution

Here we briefly discuss features of optimal control when the evolution time $T$ runs from ultrashort values to the adiabatic domain. In Fig. 6 we show non-optimized and optimized infidelity for an extended interval of $T.$ As can be seen, the non-optimized infidelity goes to zero as a consequence of adiabaticity (4). Also, the maximal difference between optimized and non-optimized infidelity appears at short time $T\lesssim 1$.

As we have shown in Refs. Pyshkin et al. (2018, 2019) the optimal shape of control function $g(t)$ is restricted by two requirements: (1) continuous transition from non-adiabatic to adiabatic time domain, and (2) nonzero time derivative $\dot{g}(t)$ at $t=0$ and $t=T.$ These two assumptions lead to the following conditions:

[TABLE]

written in parametrization (5) as:

[TABLE]

Remarkably, the conditions (22) are $T-$independent. It turns out that $\mathbf{a}_{\rm opt}$ from the numerical calculations satisfies (22) for ultrashort, short, and adiabatic processes, regardless of the time domain (this statement is as well $T-$independent when we take more than two harmonics in (5)). Note that first two harmonics in (5) provide the simplest efficient parametrization, which satisfies (21), for a strong energy pumping, with $|g(t)|\gg 1$ at a finite time interval.

In Fig. 7 we show the landscape of output fidelity as a function of $\mathbf{a}$ for $T=0.2$. One can see that the high-fidelity “islands” form horizontal equidistant lines. Appearance of these lines is related to possibility of the satisfaction of (13) by letting  $U(T\rightarrow 0)=\mathbb{I}$ with $a_{1}=\pi^{2}l/2T$ $(\approx 24.7l,\,l=\pm 1,\pm 2,\dots)$ (see Eqs.(8) and (9)). In comparison with the fidelity of different local maximums we see that initial point $\mathbf{a}=\mathbf{0}$ is the valid choice for numerical BFGS search in order to avoid traps  Rabitz et al. (2004); Ho and Rabitz (2006); Pechen and Tannor (2011); Zhdanov and Seideman (2015); Zhdanov (2018), and our assumption $\lim_{T\rightarrow 0}G(T)=0$ is corroborated. Moreover, the universal initial point $(0,0)$ (which satisfies (22)) for numerical search connects together adiabatic and non-adiabatic time domains.

V Conclusion

We demonstrated that properly-designed incomplete local control can greatly decrease infidelity of unitary evolution in the non-adiabatic time domain, even for ultrafast transition processes. We presented an approximate analytical solution for finding the optimal control parameters in the ultrashort $T$ domain and showed that optimization can lead to a linear in  $T$ decrease in the infidelity. Rather than achieving zero infidelity, this linearity is the main benefit of using the unrestricted energy resource in the case of incomplete local control.

The main features of the optimal control found by heuristic reasoning and analytical derivations have been confirmed by direct numerical simulations. Our results show that optimal control parameters for short $T$, being system-independent, are somehow universal. Surprisingly, in our approach one needs to only analyze one of extrema of a single-variable analytical function to find the optimal control parameters instead of the conventional numerical algorithm for computing propagators. We hope that our findings and approaches will be useful for further improvements of efficiency in realistic quantum control in broad variety of systems.

Acknowledgements

We gratefully acknowledge National Research, Development and Innovation Office of Hungary (Project Nos. K124351 and 2017-1.2.1-NKP-2017-00001), the Basque Country Government (Grant No. IT472-10), the Spanish Ministry of Economy, Industry, and Competitiveness (MINECO) and the European Regional Development Fund FEDER Grant No. FIS2015-67161-P (MINECO/FEDER, UE).