An algorithm for solving unconstrained unitary quantum brachistochrone problems

TL;DR

This paper presents an iterative algorithm to find time-optimal Hamiltonians for driving quantum systems between states, with rapid convergence and applications in quantum state preparation and gate design.

Contribution

The paper introduces a novel iterative method for solving quantum brachistochrone problems that converges quickly and saturates minimal evolution time bounds.

Findings

Method converges rapidly for large systems

Solutions reach the minimal evolution time bound

Provides a geometric interpretation linked to quantum phases

Abstract

We introduce an iterative method to search for time-optimal Hamiltonians that drive a quantum system between two arbitrary, and in general mixed, quantum states. The method is based on the idea of progressively improving the efficiency of an initial, randomly chosen, Hamiltonian, by reducing its components that do not actively contribute to driving the system. We show that our method converges rapidly even for large dimensional systems, and that its solutions saturate any attainable bound for the minimal time of evolution. We provide a rigorous geometric interpretation of the iterative method by exploiting an isomorphism between geometric phases acquired by the system along a path and the Hamiltonian that generates it, and discuss resulting similarities with Grover's quantum search algorithm. Our method is directly applicable as a powerful tool for state preparation and gate design problems.

Figures

Equations

$$O(\mathbf{\phi })=\sum _k{e}^{i{\phi }_k}|s_k⟩⟨r_k|,$$

$$\displaystyle\mathcal{M}_{\rho}[H]=D\big{(}M\circ D^{\dagger}HD)\big{)}D^{\dagger},\;\mbox{with}\;M_{ij}=\delta_{\lambda_{i}\lambda_{j}},$$

$$O^{(j+1)}=O^{(j)}{e}^{-i{\mathcal{M}}_{\rho }[H^{(j)}]},\text{mbox}with{H}^{(j)}=i\log O^{(j)}.$$

$$T_{QSL}\le {\tau }_{\text{min}}\le {\tau }^{(n)},$$

$${\partial }_t{U}_t=-i{O}_t{H}_{\perp }{O}_t^{†}{U}_t.$$

$$\begin{array}{rl}{\partial }_t{\tilde{U}}_t& =-i{H}_{\perp }{\tilde{U}}_t+({\partial }_t{O}_t^{†})U_t\\ & =-i(H_{\perp }-H){\tilde{U}}_t\\ & =i{H}_{\Vert }{\tilde{U}}_t.\end{array}$$

$${\eta }^{\star }(H,\rho )=\frac{\sqrt{\text{mbox}tr[{\rho }^2{H}^2]-\text{mbox}tr[(\rho H{)}^2]}}{\sqrt{\text{mbox}tr[{\rho }^2{H}^2]-(\text{mbox}tr[(\rho H{)}^2]-\text{mbox}tr[\rho H{]}^2)}},$$

$$\rho \to {\rho }^{\prime }=(1-\delta )\rho +\delta \chi ,\text{and}$$

$$\rho \to {\rho }^{\prime }=e^{iV\delta }\rho e^{-iV\delta },$$

Full text

Introduction — Inspired by the problem of finding the curve of fastest descent between two points, quantum brachistochrone problems (QBPs) aim to find the Hamiltonian that generates the time-optimal evolution between two given quantum states. Such problems have been considered, in various forms, to obtain accurate minimum-time protocols for the control of quantum systems Wang et al. (2011); Geng et al. (2016), clarify the role of entanglement and quantum correlations in time-optimal evolution Borras et al. (2008); Basilewitsch et al. (2017), study the speed of Hermitian and non-Hermitian Hamiltonians Bender et al. (2007); Assis and Fring (2008); Mostafazadeh (2007), and improve bounds on the minimal time of evolution, known as quantum speed limits Russell and Stepney (2017); Brody et al. (2015); Campaioli et al. (2017a). QBPs have a special role in quantum information theory and technology, where they are of fundamental importance to accurately perform tasks such as preparing a desired state of a system, or implementing a specific gate, while satisfying the strict physical and fault tolerance requirements imposed by the locality of interactions and short decoherence times Caneva et al. (2009); Wang et al. (2015). For this reason, their solutions have found applications in information processing Werschnik and Gross (2007); Reich et al. (2012); Jäger et al. (2014); Dür et al. (2014); Huang and Goan (2014); Schulte-Herbrüggen et al. (2005), quantum state preparation Vandersypen and Chuang (2005); van Frank et al. (2016); Islam et al. (2011); Senko et al. (2015); Jurcevic et al. (2014); Zhou (2016), cooling Wang et al. (2011); Machnes et al. (2012), metrology Kessler et al. (2014); Giovannetti et al. (2004); Burkard et al. (1999), and quantum thermodynamics Rezakhani et al. (2009); del Campo (2013); Binder et al. (2015); Campaioli et al. (2017b). Moreover, QBPs give a physical interpretation to the complexity of quantum algorithms, which emerges from the minimisation of the time required to obtain the desired unitary transformation Carlini et al. (2006).

Solving QBPs is generally hard, and accurate analytic and numerical solutions are only known for some special cases, such as the unconstrained unitary evolution of pure states, and control problems of well structured quantum systems Carlini et al. (2006); Levitin and Toffoli (2009); Carlini et al. (2011); Carlini and Koike (2012, 2013); Brody et al. (2015); Wang et al. (2015). There are methods that can be used to address a relatively large class of QBPs, but obtaining precise solutions becomes increasingly challenging as the dimension of the system grows, and the constraints on the control Hamiltonian become more complex. For instance, the quantum brachistochrone equations proposed in Ref. Carlini et al. (2006) involve the solution of ordinary differential equations with boundary conditions, for which there are no efficient numerical methods when high-dimensional systems are considered Krotov (1993); Carlini et al. (2007); Rezakhani et al. (2009); Carlini et al. (2008); Wang et al. (2015). One crucial open problem is the case of unconstrained unitary evolution between two mixed states. As opposed to the case of pure states, the solution to this QBP is not known, except for special cases with highly degenerate spectra. This is akin to the problem of finding tight bounds on the minimal time of evolution, known as quantum speed limits (QSLs). These well-known bounds are attainable for pure states, but are often loose for mixed states due to the complicated structure of the space of density operators Deffner and Lutz (2013); del Campo et al. (2013); Sun et al. (2015), though we have recently proposed geometric methods to derive simple, efficient, and tight bounds for the time of evolution of mixed states Campaioli et al. (2017a, 2018).

In this Letter we take a similar approach to solving the complementary problem, that of finding the optimal unconstrained unitary evolution between two mixed states. That is, we look for the generator $H$ of a unitary operator $U=\exp(-iHt)$ that takes a mixed state $\rho$ to $\sigma=U\ \rho\ U^{\dagger}$, such that the transition time $t$ is minimised. We introduce an iterative method to search for the optimal time-independent Hamiltonian, while respecting an energetic constraint. The method progressively improves the efficiency of the Hamiltonian Uzdin et al. (2012) that generates the evolution, by applying an algorithm that rapidly converges even for large dimensional systems.

We investigate the efficacy of our method using bounds on the minimal time of evolution, i.e., QSLs. By comparing the achievable upper bound, provided by the iterative method, with the inviolable lower bound offered by QSLs, we demonstrate the strong synergy of the two results: When the two are found to be close, we can be sure that both results are close to optimal. We study the performance of our method with respect to the size of the system, and its dependence on parameters such as the convergence threshold, before discussing direct applications of the algorithm, its potential combination with time-optimal gate design and Monte-Carlo methods, and its geometric interpretation, juxtaposing it to that of Grover’s quantum search algorithm.

Time-optimal evolution and Hamiltonian efficiency — Let us begin by considering the QBP for a two-level system, defined by the initial state $\rho=(\mathbb{1}+p\,\Lambda_{x})/2$, and the target state $\sigma=(\mathbb{1}+p\,\Lambda_{y})/2$. Here $p\in(0,1]$, and $\bm{\Lambda}=(\Lambda_{x},\Lambda_{y},\Lambda_{z})$ are Pauli matrices. The Bloch vectors for the two states are given by $\bm{r}=(p,0,0)$ and $\bm{s}=(0,p,0)$, respectively. Any unitary $O(\bm{\varphi})=e^{i\varphi_{1}}|{s_{1}}\rangle\!\langle{r_{1}}|+e^{i\varphi_{2}}|{s_{2}}\rangle\!\langle{r_{2}}|$ can be used to map $\rho\to\sigma$, where $|{r_{k}}\rangle$ and $|{s_{k}}\rangle$ are the eigenvectors of $\rho$ and $\sigma$, respectively, and where $\bm{\varphi}=(\varphi_{1},\varphi_{2})$ represent the different geometric phases that the state gathers from $O(\bm{\varphi})$. The Hamiltonians $H(\bm{\varphi})=i\log O(\bm{\varphi})$ defines different evolutions depending on the choice of $\bm{\varphi}$.

Let us focus on two possible choices for these phases, $\bm{\varphi}_{z}=(\pi/4,\pi/4)$, and $\bm{\varphi}_{xy}=(\pi/4,-3\pi/4)$. The associated Hamiltonians are given by $H(\bm{\varphi}_{z})=\Lambda_{z}\pi/4$, and $H(\bm{\varphi}_{xy})=(\Lambda_{x}+\Lambda_{y})\sqrt{2}\pi/4$. These two Hamiltonians generate different paths on the Bloch sphere. Under energetic constraints, such as fixing the standard deviation $\Delta H_{\rho}=\sqrt{\mbox{tr}{\rho H^{2}}-\mbox{tr}{\rho H}^{2}}$ with respect to the initial state, the lengths of these paths meaningfully correspond to evolution times. The path generated by $H(\bm{\varphi}_{z})$ connects $\rho$ to $\sigma$ with an arc of great circle, i.e., the geodesic with respect to the Fubini-Study (FS) metric 111The FS metric is the unique, unitarily invariant metric on the space of pure states Bengtsson and Zyczkowski (2008), thus it constitutes a solution to the considered QBP for any homogeneous energy constraint Russell and Stepney (2017), whereas $H(\bm{\varphi}_{xy})$ draws a longer path, which is not time optimal.

A heuristic explanation for the variable performance of different Hamiltonians is that the Bloch vectors $\bm{h}$ associated with them $H=\bm{h}\cdot\bm{\Lambda}$ have a different orientation with respect to $\bm{r}$ and $\bm{s}$. In particular, when $\bm{h}$ is orthogonal to $\bm{r}$, it generates a rotation on a plane that passes through the origin of the Bloch sphere, while when $\bm{h}$ is not perpendicular to $\bm{r}$, it generates a rotation on a plane that does not. The slower evolution can be seen as due to a less efficient use of the resources encoded in constraints on the energy, which are wasted on parts of the Hamiltonian that not actively drive the system Uzdin et al. (2012). The notion of efficiency can be formalised for pure states as in Ref. Uzdin et al. (2015), and generalised to mixed states, as discussed in the Supplemental Material.

This intuitive geometric argument can be generalised for dimension $d>2$ by replacing the notion of orthogonality between vectors with commutation relations between states and Hamiltonians. Given an operator $\rho$, the space of Hamiltonians (the Lie algebra $\mathfrak{su}(d)$) splits into a ‘parallel’ subspace, closed under multiplication, that commutes with $\rho$, and an orthogonal ‘perpendicular’ subspace that does not. This allows us to decompose the Hamiltonian $H$ into components $H_{\|}$ and $H_{\perp}$ which are elements, respectively, of these two subspaces, such that $[H_{\|},\rho]=0$ and $[H,\rho]=[H_{\perp},\rho]$ 222These correspond to elements of the vertical and horizontal tangent bundles in a fibre bundle representation where the state space (with a given spectrum) is the base manifold and the phases $\mathbf{\varphi}$ are the fibres.. This observation leads to the idea at the core of the iterative method that we will introduce in the next section, where efficient Hamiltonians for a QBP are achieved by requiring their parallel component to be vanishing at all times during the evolution.

Iterative method for efficient Hamiltonians — We can now consider the more general QBP defined by an arbitrary isospectral pair of initial and final states $\rho=\sum_{k}\lambda_{k}|{r_{k}}\rangle\!\langle{r_{k}}|$ and $\sigma=\sum_{k}\lambda_{k}|{s_{k}}\rangle\!\langle{s_{k}}|$ of finite dimension $d$. When non-degenerate states are considered, the operator

[TABLE]

represents the most general unitary that connects initial and final states $\rho$ and $\sigma$, while $\varphi_{k}$ represent the geometric phases gained evolving along the path generated by $i\log O(\bm{\varphi})$. When degenerate states are considered, the phases $\varphi_{k}$ are replaced by unitary operators $U_{k}\in SU(m)$ on the subspace associated with eigenvalues $\lambda_{k}$ with multiplicity $m$. Since any unitary $O(\bm{\varphi})$ maps $\rho\to\sigma$, we can choose an arbitrary initial phase $\bm{\varphi}^{(0)}$ to obtain the initial unitary $O^{(0)}=O(\bm{\varphi}^{(0)})$. The Hamiltonian $H^{(0)}=i\log O^{(0)}$, canonically associated with $O(\bm{\varphi}^{(0)})$, is then split into a parallel and perpendicular component $\mathcal{M}_{\rho}[H^{(0)}]=H^{(0)}_{\|}$ and $H^{(0)}_{\perp}=H^{(0)}-H^{(0)}_{\|}$, such that $[H^{(0)}_{\|},\rho]=0$, via a map $\mathcal{M}_{\rho}$ that depends on $\rho$, and which will be referred to as mask. A possible choice for such a mask is given by

[TABLE]

where $D$ diagonalises $\rho$, i.e., its columns are given by the eigenstates $\{|{r_{k}}\rangle\}$ of $\rho$, and where $\circ$ is the Hadamard, i.e., element-wise, product Horn and Johnson (2012). Note that, when degenerate states are considered, the eigenvalues’ multiplicity defines the structure of the mask via $\delta_{\lambda_{i}\lambda_{j}}$.

We now notice that the unitary $U^{(0)}=\exp[-iH_{\|}^{(0)}]$ can be composed with $O^{(0)}$ to obtain another unitary $O^{(1)}=O^{(0)}U^{(0)}$, formally equivalent to evolving with the time-dependent Hamiltonian $e^{-iH^{(0)}t}H^{(0)}_{\perp}e^{iH^{(0)}t}$ (as we detail in the Supplemental Material), which also drives $\rho\to\sigma$ since $[\rho,H_{\|}]=0$. In general, the unitary $O^{(1)}$ is associated with a new geometric phase, and the Hamiltonian $H^{(1)}=i\log O^{(1)}$ draws a different path, which might be shorter (or longer) than that generated by $H^{(0)}$. In a second iteration, we apply the mask to the new Hamiltonian to obtain $H^{(1)}_{\|}=\mathcal{M}_{\rho}[H^{(1)}]$, and define another gate $O^{(2)}$ analogously. By iterating this method we obtain a sequence of Hamiltonians $\{H^{(j)}\}_{j=0}^{n}$ that drive $\rho\to\sigma$ via $\exp[-iH^{(j)}]$. The unitary at each step is related to the previous one by

[TABLE]

A necessary and sufficient condition for the routine to reach a fixed point $H^{(n)}$, such that $H^{(n+1)}=H^{(n)}$, is given by $\mathcal{M}_{\rho}[H^{(n)}]=0$, i.e., when the parallel part of the $n$-th Hamiltonian vanishes under the action of the mask $\mathcal{M}_{\rho}$, which follows trivially from the fact that in this case $e^{-i\mathcal{M}_{\rho}[H^{(n)}]}=\mathbb{1}$.

When implementing the method numerically one has to fix a convergence threshold to stop the routine as soon as the parallel component $H_{\|}^{(n)}$ becomes small enough with respect to the full Hamiltonian $H^{(n)}$. We have chosen to quantify this threshold with a positive number $\varepsilon\ll 1$, such that the method stops when $\lVert H_{\|}^{(n)}\rVert_{HS}\leq\varepsilon\lVert H^{(n)}\rVert_{HS}$. Accordingly, the number of iterations $n$ required for the method to converge implicitly depends on $\varepsilon$. Numerical evidence suggests that the sequence always converges towards a Hamiltonian $H^{(n)}$ that is fully perpendicular with respect to $\rho$ along the whole evolution, in the sense that the parallel components eventually vanish within the precision defined by $\varepsilon$. This iterative method is summarised in its simplest form by Algorithm 1, and can be interpreted as an optimisation of energy cost associated with the different geometric phases $\bm{\varphi}$, accomplished via the recursive suppression of ineffective components of the Hamiltonians.

There is also a geometric interpretation of our method analogous to that of Grover’s famous quantum search algorithm Grover (1997). The latter aims to find the unique input, encoded in a quantum state, of a function that produces a particular output value. Its action on an initialisation state $|{\psi^{(0)}}\rangle$ can be interpreted as a sequence of rotations that converges with high probability to the desired state. Similarly, the iterative method introduced here can be seen as a sequence of rotations of some initialisation unit vector $\hat{\bm{h}}^{(0)}$, associated with the Hamiltonian via $H^{(0)}=\bm{h}^{(0)}\cdot\bm{\Lambda}$, where $\bm{\Lambda}=(\Lambda_{1},\dots,\Lambda_{d^{2}-1})$ forms a Lie Algebra for $SU(d)$. However, unlike Grover’s algorithm, the vectors $\{\bm{h}^{(n)}\}$ do not span a two-dimensional real plane, but a high-dimensional subspace of $\mathbb{R}^{d^{2}-1}$.

Performance of the method — As anticipated in the introduction, our iterative method can be directly applied to solve QBPs defined by the unconstrained time-optimal unitary evolution of any isospectral pair of density operators. A sub-class of these problems with known solutions are those of unconstrained unitary evolution between pure states or between mixed states whose eigenvalues are all degenerate except one. To demonstrate the performance of our method, we have tested it on Bures-random pairs of pure states $|{\psi}\rangle\!\langle{\psi}|$, $|{\phi}\rangle\!\langle{\phi}|$ of dimension $d=2,\dots,100$, successfully obtaining Hamiltonians that are time-optimal and fully efficient with respect to the notion of efficiency introduced in Ref. Uzdin et al. (2015). We can be confident that the solutions we obtain are globally optimal by exploiting the fact that the quantum speed limit for unitary evolution of pure states is attainable Levitin and Toffoli (2009). Comparing the evolution time $\tau^{(n)}$ of the optimized Hamiltonian $H^{(n)}$ obtained from Algorithm 1 with the minimal evolution time given by the Mandelstam-Tamm bound $T_{QSL}(|{\psi}\rangle,|{\phi}\rangle)=d_{FS}(|{\psi}\rangle,|{\phi}\rangle)/\Delta H^{(n)}$ Deffner and Campbell (2017), where $d_{FS}(|{\psi}\rangle,|{\phi}\rangle)=\cos^{-1}(|\langle{\psi}||{\phi}\rangle|)$ is the FS distance, we find $\tau^{(n)}\approx T_{QSL}$, within the precision imposed by $\varepsilon$, for all cases, as shown in Fig. 1.

A more challenging test was run on random pairs of mixed states of dimension $d=3,\dots,100$ (for initial states with both degenerate and non-degenerate spectra), for which a general solution to the unconstrained unitary QBP is not known. Since the quantum speed limit for the unitary evolution of mixed states Campaioli et al. (2017a) is in general not tight, it is harder to benchmark the quality of the solutions provided by our method in this case. On the other hand, the natural synergy between this iterative method and the QSLs introduced in Ref. Campaioli et al. (2017a) can be used to gauge the performance of the former and the tightness of the latter. The actual minimal time of evolution $\tau_{\text{min}}$ for a given choice of $\rho$ and $\sigma$ is bounded as

[TABLE]

where $T_{QSL}$ corresponds to the inviolable lower bound offered by QSLs of Ref. Campaioli et al. (2017a), and $\tau^{(n)}$ to an achievable upper bound, provided by the solutions obtained with our iterative method. The first inequality can be saturated for pure states and for states with $d-1$ degenerate eigenvalues. For such case, the second inequality can be saturated up to the precision imposed by $\varepsilon$, due to the fact that the optimized time $\tau^{(n)}$ inherently depends on the convergence threshold $\varepsilon$, which can be seen as an implicit trade-off between $n$ and $\varepsilon$. By looking at the difference between $\tau^{(n)}$ and $T_{QSL}$ for the given solution $H^{(n)}$, we can assess the quality of both results, which we find to coincide within the precision set by $\varepsilon$ for some choices of $\rho$ and $\sigma$, even when their degeneracy structure differs from that of pure states.

Moreover, the performance of the solutions $H^{(n)}$ obtained with this iterative method is stable under small perturbations of the initial state $\rho$. These could be introduced through convex mixing $\rho\to\rho^{\prime}=(1-\delta)\rho+\delta\>\chi$, where $\chi$ is a random state of the same dimension as $\rho$, or unitary rotation $\rho\to\rho^{\prime}=e^{iV\delta}\rho\;e^{-iV\delta}$, where $V$ is a random Hamiltonian of unit Hilbert-Schmidt norm. As shown in the Supplemental Material, the difference between the solution to the original problem and that for the perturbed version does not grow quickly as a function of the perturbation strength $\delta$. We now discuss the rate at which the algorithm converges.

Convergence of the method — The number $n$ of iterations required for convergence generally grows with the dimension $d$ of the system, and the strictness of the precision set by $\varepsilon$. Recalling that the number of elements of the matrices associated with non-commuting density operators $\rho$ and $\sigma$ grows quadratically with $d$, one might expect conservatively that the average number of iterations $\bar{n}$ would grow in the same way. However, $\bar{n}$ grows logarithmically with $d$, as shown in Fig. 2, with a slower growth for the case of pure states and highly degenerate mixed states. Equivalently, for composite systems $\bar{n}$ grows linearly with the number of constituent subsystems.

The choice of initial phase vector $\bm{\varphi}^{(0)}$ also affects the number of iterations required to converge, without noticeably affecting the performance of the end point $H^{(n)}$; we have never observed the case where two different choices for $\bm{\varphi}^{(0)}$ lead to converged solutions with inconsistent evolution times. Moreover, while $\bar{n}$ grows with $d$ and $\varepsilon^{-1}$, a good choice of initial geometric phase $\bm{\varphi}^{(0)}$ can still lead to quick convergence, returning runs that can take less then 10 iterations even for $d=100$ and $\varepsilon=10^{-4}$. For these reasons, there is the potential to combine Algorithm 1 with Monte-Carlo sampling methods Hastings (1970); Vanderbilt and Louie (1984); Kalos and Whitlock (2008), in order to deploy many parallel runs of the iterative method, for the same pair of initial and final states $\rho$ and $\sigma$, with each run initialised with a different, randomly chosen geometric phase $\bm{\varphi}$. Since the whole computation is stopped as soon as one of these runs converges to a solution, the number $n$ of iterations is guaranteed to be smaller than the average one. We give evidence in the Supplemental Material that this approach is feasible. There we show that, for fixed $\rho$ and $\sigma$, the fraction of initial phases that converges quickly is still almost always significant, even with large $d$ and $\varepsilon^{-1}$. The use of Monte-Carlo methods comes at the expense of requiring a much larger number of simultaneous tasks, which form an embarrassingly parallel workload Herlihy and Shavit (2012).

Discussion — In this Letter, we have introduced an iterative algorithm to obtain efficient time-independent Hamiltonians that generate fast evolution between two $d$-dimensional isospectral states $\rho$ and $\sigma$. This constitutes a solution to the quantum brachistochrone problems for unitary evolution of mixed states, with direct applications to state preparation and quantum optimal control problems. We have shown that a necessary and sufficient condition for our method to converge is to reach an efficient Hamiltonian, and provided extensive numerical evidence that our method rapidly reaches a solution even for large Hilbert space dimension $d\sim 100$. In particular, we have shown that the average number of iterations required for convergence grows logarithmically with the dimension of the system and that this number can be further reduced by running multiple instances of the algorithm in parallel. For this reason our method provides a viable and powerful approach to the design of time-optimal Hamiltonians of high-dimensional systems.

A possible application of the iterative method introduced here could be found in its combination with quantum optimal control methods. Algorithm 1 can be used to identify a unitary $\widetilde{O}$, associated with an efficient Hamiltonian $\widetilde{H}$, to be implemented by means of time-optimal gate design methods, such as those experimentally realised in Ref. Geng et al. (2016). As opposed to QBPs defined by the specific choice of a pair of state $\rho$ and $\sigma$, time-optimal gate design problems aim to realize a certain unitary $U$ (independently of the input states) by means of a time-dependent Hamiltonian, achieved by optimizing some control pulses. More generally, our method has broad application in further elucidating the geometric structure of quantum control, since it has been shown that quantum brachistochrone problems can be recast as those of finding geodesics in the space of unitary operators Wang et al. (2015).

Our method could also be extended to more general control problems. For instance, a similar approach could be taken to the open evolution quantum systems, where the Hamiltonian is replaced by a Liouvillian generating a family of CPTP maps. Moreover, whether open or not, constraints on the form of the generators could be included. Such constraints are naturally imposed by physical restrictions on the order and range of interactions, and can dramatically change the bounds on minimal evolution time, as recently shown in Ref. Campaioli et al. (2017b). Even for unconstrained unitary evolution, there may be alternative formulations of our algorithm which perform better under certain circumstances. In the Supplemental Material we briefly outline a version (which does not qualitatively differ in its performance) where the part of the Hamiltonian that commutes with the final state $\sigma$ is additionally removed in each iteration. It is possible that there are further modifications along these lines that could improve the rate of convergence on the solution of the QBP. Finally, the study of the stability of these methods on experimental errors, such as the precision on the control parameters and on time measurements, would be of paramount practical importance.

I Derivation of Algorithm 1

With reference to the problem defined in the second section of this Letter, let $O=\sum_{k}|{s_{k}}\rangle\!\langle{r_{k}}|$, and $i\log O=H=H_{\|}+H_{\perp}$, where $H_{\|}$ commutes with $\rho$. Consider the unitary $U_{t}$ generated by the time dependent Hamiltonian $O_{t}H_{\perp}O^{\dagger}_{t}$, where $O_{t}=\exp[-iHt]$, which satisfies the equation

[TABLE]

To leading order in time, this will rotate the perpendicular part of the Hamiltonian to follow the system’s evolution. Now consider the unitary $\widetilde{U}_{t}=O^{\dagger}_{t}U_{t}$, whose equation of motion is given by

[TABLE]

Eq. (6) implies that $\widetilde{U}_{t}=\exp[iH_{\|}t]$, thus that $\widetilde{U}_{t}$ commutes with $\rho$. Accordingly, $U=U_{1}=O_{1}\widetilde{U}_{1}$ transforms $\rho\to\sigma$, just as well as $O$ does. The procedure can be iterated until it reaches a fixed point, i.e., until $H=H_{\perp}$ and $H_{\|}=0$.

The iterative method introduced here involves removing the part of the Hamiltonian operator that commutes with the initial state $\rho$. A similar approach could be taken, where the part of the Hamiltonian that commutes with the final state $\sigma$ is removed, and the unitary joining the two states is acted on from the left in each iteration, instead of from the right. It also possible to remove both components simultaneously, replacing Eq. (3) with $O^{(j+1)}=e^{i\mathcal{M}_{\sigma}[H^{(j)}]}O^{(j)}e^{-i\mathcal{M}_{\rho}[H^{(j)}]}$. While we find this approach generally changes the convergence rate of the algorithm, it does not seem to perform significantly better than the algorithm outlined in the main text.

II Hamiltonian Efficiency for Density Operators

When pure states are considered, we refer to the the notion of Hamiltonian efficiency introduced in Ref. Uzdin et al. (2015), given by $\eta[H,\rho]:=\Delta H_{\rho}/\lVert H\rVert_{op}$, with $\lVert\cdot\rVert_{op}$ being the operator norm. It is possible to consider different energetic constraints, such as $\lVert H\rVert_{HS}=\omega$, known as finite energy bandwidth, or $\lVert H\rVert_{op}=\omega$, which corresponds to bounding the largest eigenvalue of the Hamiltonians. Different notions of efficiency can be obtained for any choice of speed and total energy measure. The efficiency measures how much of the energy available, quantified by $\lVert H\rVert_{op}$, is transferred to the system to drive the evolution, and thus converted into speed $\Delta H_{\rho}$. If we assess the performance of the two Hamiltonians $H(\bm{\varphi}_{z})$ and $H(\bm{\varphi}_{xy})$ considered in the Letter using $\eta$, we obtain $\eta(H(\bm{\varphi}_{z}),\rho)=p$, and $\eta(H(\bm{\varphi}_{xy}),\rho)=p/2$. The dependence on $p$ (and thus on the purity of the states) implies that this notion of efficiency cannot not be saturated in general.

A possible notion of Hamiltonian efficiency that can be saturated for the case of density operators is given by

[TABLE]

which is a positive function smaller than 1, that is saturated for $H_{\|}=0$. The numerator of $\eta^{\star}$ is proportional to the speed $\lVert\dot{\rho}\rVert_{HS}$ of the generator $\dot{\rho}=-i[H,\rho]$ del Campo et al. (2013); Deffner and Lutz (2013); Campaioli et al. (2017a), and reduces to $\Delta H_{\rho}$ when pure states are considered, for which $\eta^{\star}=1\leftrightarrow\eta=1$. Moreover, like $\Delta H_{\rho}$, $\lVert\dot{\rho}\rVert_{HS}$ is also invariant under the addition of a parallel component $H_{\|}$ with respect to $\rho$. Accordingly, the Hamiltonians generated by the iterative method are of unit efficiency for all the considered QBPs, which reflect the ability to converge to a fully-perpendicular Hamiltonian. However, the convergence of this iterative method is, in general, non-monotonic with respect to this notion of efficiency, which can decrease before reaching $\eta^{\star}=1$ after several iterations, as shown in Fig. 3.

III Performance and convergence of Algorithm 1

The two parameters that we considered to study the performance of the iterative method introduced in this Letter are $n$ and $\tau^{(n)}/T_{QSL}$, respectively given by the number of iterations required to converge, and the ratio between the optimised time $\tau^{(n)}$ and the QSL for the problem defined by the chosen initial and final state and $H^{(n)}$. Eq. 4 implies that $\tau^{(n)}/T_{QSL}\geq 1$, but since the $T_{QSL}$ is not tight in general, $\tau^{(n)}/T_{QSL}>1$ does not imply that $H^{(n)}$ is not a time-optimal Hamiltonian.

To study how the performance of the iterative method varies with $d$, we sampled a large number of pairs $\rho,\sigma$ from the Bures-ensemble of pure and mixed states of dimension up to 100. We have considered a sample size of $10^{4}$ pairs of states $\rho$ and $\sigma$ for dimension $2,3,4,5,6,7,8,9,10,20,30,40,50$, and a sample of $10^{2}$ pairs for $d=60,70,80,90,100$. Numerical evidence suggests that the average number of iterations required for convergence $\bar{n}$ grows linearly with $\varepsilon^{-1}$, and, as shown in Fig. 2, logarithmically with $d$.

The choice of initial geometric phase $\bm{\varphi}^{(0)}$ strongly affects $n$ without affecting $\tau^{(n)}$ (within the precision set by $\varepsilon$). We studied the dependence of $n$ on the choice of initial geometric phases by running several instances of the method for the same pair of states $\rho$ and $\sigma$, uniformly sampling each phase $\varphi^{(0)}_{j}$ in the interval $[0,2\pi]$. As shown in Fig. 4, the slow growth of the 20th percentile of the number of iterations required for converges suggests that Monte-Carlo sampling methods could combined with this method to speed up the computation.

IV Stability under perturbations

As mentioned in the Letter, the iterative method defined by Algorithm 1 is stable under small perturbations of the initial and final states $\rho$ and $\sigma$. We considered convex and unitary perturbations, respectively given by

[TABLE]

where $\chi$ is a random state of the same dimension as $\rho$, and $V$ is a random Hamiltonian of unit Hilbert-Schmidt norm. For each pair $\rho,\sigma$, the perturbed final state $\sigma^{\prime}$ is then obtained by applying any unitary $O$ that maps $\{|{r_{k}}\rangle\}_{k=1}^{d}\to\{|{s_{k}}\rangle\}_{k=1}^{d}$ on the perturbed initial state $\rho^{\prime}$. We look for the effect of perturbations calculating the relative deviations from the unperturbed solutions, given by $\lVert H^{(n)}-{H^{\prime}}^{(m)}\ \rVert_{HS}/\lVert H^{(n)}\rVert_{HS}$, where ${H^{\prime}}^{(m)}$ is the solutions to the perturbed problem. Numerical evidence suggest that the relative deviations grow slowly as a function of the perturbation strength $\delta$, and are negligible for $\delta\ll\varepsilon$, as shown in Fig. 5.