Higher-order quantum transformations of Hamiltonian dynamics

Efficiently simulating the dynamics of complex quantum systems is often stated as one of the main motivations of quantum computing. While such simulation is considered hard on classical computers, a range of efficient quantum algorithms have been developed for simulating Hamiltonian dynamics [1, 2, 3, 4, 5, 6, 7]. The core principle behind the standard Hamiltonian simulation algorithms is that the desired Hamiltonian dynamics can be well-approximated by a series of (arguably) simpler quantum operations. These algorithms rely on having a classical description of the desired Hamiltonian, which can often be used for obtaining a decomposition into a sum of easily implementable terms. This limits the way we can develop large-scale, complex quantum programs for dynamics simulation. Quantum algorithms which do not require detailed descriptions of quantum resources have a higher flexibility in quantum software development. This is related to the fundamental problem of understanding how much quantum algorithms need to rely on the classical description of their inputs in order to achieve quantum advantages in information processing.

In this work, we study Hamiltonian dynamics that can be implemented given a seed Hamiltonian $H$ without using a classical description of $H$. That is, we study transformations of black-box Hamiltonians. We present a quantum algorithm that simulates the dynamics of $f(H)$, where $f$ is any physically realizable linear function of $H$, given a description of $f$ and using a black-box Hamiltonian $H$ with a bounded energy range. This algorithm is an instance of a higher-order quantum transformation on the unitary operation realized by the seed Hamiltonian dynamics. The functions that the algorithm can implement include both the negative time-evolution and the time-reversal of an unknown Hamiltonian evolution by considering $f(H)=-H$ and $f(H)=H^{T}$ (transposition of $H$ in terms of the computational basis), respectively. Such general transformations have applications ranging from fundamental physics simulations to potential improvements in state-of-the-art algorithms, such as the Hamiltonian singular value transformation [8]. We also show an application of our algorithm for Hamiltonian learning [9], in particular, a task of efficiently estimating a parameter of a multi-parameter Hamiltonian using Hamiltonian dynamics by appropriately choosing $f(H)$.

Our work constitutes the first systematic study of higher-order quantum transformations in the context of Hamiltonian dynamics. Higher-order quantum transformations have attracted significant attention in recent years in the context of quantum circuit transformations, and are also known as superchannels, supermaps, quantum combs and process matrices [10, 11, 12, 13, 14, 15]. Higher-order algorithms for quantum computation can be seen as an analogue of functional programming in classical computing, where the possible inputs to an algorithm are quantum channels (for example, unitaries) specified “operationally” by their input-output description only (i.e. as black boxes).

Previous works on this topic have focused on the possible transformations that can be achieved when the input channels are taken to be a finite sequence of quantum gates [10, 16, 17, 18, 19, 20, 21, 22, 23, 15, 24]. Yet, the resources available in a given computation are not always best described by a finite sequence of gates, but rather by a continuously parameterized Hamiltonian evolution. In fact, it is known that certain functions such as controllization, which cannot be implemented on black box unitaries [25, 26, 27, 28], can in fact be implemented if access to the underlying Hamiltonian evolution is given [29, 17]. This is because it is possible to apply an arbitrary fractional power of an unknown Hamiltonian evolution by changing the evolution time, whereas applying a fractional power is not possible for black box unitaries.

Summary of algorithm.—We now present our algorithm in detail (see Algorithm 1). We represent Hilbert spaces of an $n$-qubit quantum system and a single-qubit auxiliary system by $\mathcal{H}$ and $\mathcal{H}_{c}$, respectively. We assume that we can invoke the Hamiltonian evolution $e^{-iH\tau}$ of a seed Hamiltonian $H\in\mathcal{L}(\mathcal{H})$ with an upper bound $\Delta_{H}$ of the difference between the maximum and the minimum energy eigenvalues is given, for any time $\tau>0$.

We assume that $f(I)\propto I$, which ensures that the resulting evolution $e^{-if(H)t}$ preserves the invariance under the global phase of $e^{-iH\tau}$. This class of $f$ covers all physically realizable linear transformations of $H$ as shown in Appendix C. In our setting, we are given the Pauli transfer matrices $\gamma$ [30] as in Eq. (1) of a hermitian-preserving linear map $f:\mathcal{L}(\mathcal{H})\to\mathcal{L}(\mathcal{H})$. Our algorithm simulates the Hamiltonian evolution $e^{-if(H)t}$ for any $t>0$ representing the time for the transformed Hamiltonian dynamics up to an error $\epsilon>0$ and variance $4\epsilon$ (see proof in Appendix A, which relies on more general results proven in Appendix B. A similar analysis of variance is obtained in probabilistic state synthesis [31]).

The runtime of our algorithm is upper-bounded as $O(\beta^{2}t^{2}\Delta_{H}^{2}n/\epsilon)$ in terms of $\beta:=2\sum_{\vec{w},\vec{u}}|\gamma_{\vec{w},\vec{u}}|$, which is a function of elements of the Pauli transfer matrix. The total evolution time of the input dynamics $e^{-iH\tau}$ is $\beta t$ which can be shown from step 3 and step 6 of the Algorithm 1.

In Algorithm 1, the gate sequence $V_{f,j}$ is constructed only from controlled-Pauli gates, which are Clifford gates. The only element which may be non-Clifford is the black box dynamics $e^{-iH\tau}$. Dependence on the transformation $f$ is specified only through the probability distribution $p_{\vec{u},\vec{w}}^{(2)}$ in choosing $(\vec{u},\vec{w})$ in Step 4 and through the gate $X^{s_{f}}$ in Step 5. The total runtime $O(\beta^{2}t^{2}\Delta_{H}^{2}n/\epsilon)$ is calculated by multiplying the number of iterations $N$ with the runtime $O(n)$ for implementing the controlled-Pauli gates in $V_{f,j}$ using CNOT gates and single-qubit Clifford gates. Note that $N$ is independent of $n$, even though the set of parameters $j\in(\{0,1,2,3\}^{n})^{4}$ has exponentially many terms. The procedure of Algorithm 1 is summarized in Figure 1.

To understand how the gate sequence $V_{f,j}$ transforms the Hamiltonian at each iteration, Figure 2 shows the explicit evolution of an arbitrary seed Hamiltonian $H=\sum_{\vec{v}}c_{\vec{v}}\sigma_{\vec{v}}$ after pre- and post-processing with each successive gate in the (random) sequence $V_{f,j}$ averaged over $\vec{v}$ and $\vec{v}^{\prime}$, namely, $\frac{1}{16^{n}}\sum_{\vec{v},\vec{v}^{\prime}}V_{f,j}(I\otimes e^{-iHt\beta/N})V_{f,j}^{\dagger}$. For simplicity, $H$ is assumed to be traceless (any trace-full part is proportional to the identity and is therefore invariant under the overall transformation $f$, by construction). The gate sequence of $\frac{1}{16^{n}}\sum_{\vec{v},\vec{v}^{\prime}}V_{f,j}(I\otimes e^{-iHt\beta/N})V_{f,j}^{\dagger}$ is constructed in a functional programming approach, namely, by concatenations of a series of higher-order transformations, here called Processes \scriptsize1⃝ to \scriptsize7⃝. Each of these processes is designed to implement a Hamiltonian dynamics whose Hamiltonian is given by

Applying the first controlled-$\sigma_{\vec{v}}$ gate before and after the seed Hamiltonian evolution $e^{-iHt\beta/N}$ with $\vec{v}$ chosen independently from the uniform distribution in each iteration but perfectly correlated between the pre- and post-processing within each iteration (Process \scriptsize1⃝) implements Hamiltonian controllization [17]. That is, the effective evolution $({\tt ctrl}\sigma_{\vec{v}})e^{-i(I\otimes H)t\beta/N}({\tt ctrl}\sigma_{\vec{v}})$ averaged over $\vec{v}$ simulates a Hamiltonian of the form $H\oplus 0$.

Process \scriptsize4⃝ is based on the identity

where $H_{00},\ H_{01},\ H_{11}\in\mathcal{L}(\mathcal{H})$ are arbitrary operators, noting that for all $\vec{u}$, $\mathrm{tr}(H)=\mathrm{tr}(\sigma_{\vec{u}}H\sigma_{\vec{u}})=0$ and $\mathrm{tr}(\sigma_{\vec{u}}H)=\mathrm{tr}(H\sigma_{\vec{u}})=2^{n}c_{\vec{u}}$. The two gates $\sigma_{\vec{v}^{\prime}}$ are chosen independently from the uniform distribution at each iteration.

Algorithm 1 is universal in the sense that it transforms the dynamics of any seed Hamiltonian $H$ to that of the Hamiltonian $f(H)$ for any choice of a physically realizable linear transformation $f$, even if $H$ is only given as a black box. Therefore the algorithm is an instance of higher-order quantum transformations of Hamiltonian dynamics. The algorithm makes use of a general approximation technique for simulating Hamiltonians of the form $g(H)=\sum_{j}h_{j}U_{j}HU_{j}^{\dagger}$, where $\{U_{j}\}_{j}$ is a set of unitaries, $\{h_{j}\}_{j}$ is a set of positive numbers and $H$ is a seed Hamiltonian. This approximation is represented by the following circuit

Applications of the algorithm.—We describe three applications of our algorithm: the negative time-evolution of Hamiltonian dynamics $e^{-iH\tau}\mapsto e^{iHt}\ (\tau,t>0)$, the time-reversal of Hamiltonian dynamics $e^{-iH\tau}\mapsto e^{-iH^{T}t}\ (\tau,t>0)$, and a Hamiltonian single-parameter learning task of estimating one of the parameters represented by a Pauli coefficient $c_{\vec{v}}$ $(|c_{\vec{v}}|\leq 1\ ,\vec{v}\in\{0,1,2,3\}^{n})$ of a Hamiltonian $H=\sum_{\vec{u}}c_{\vec{u}}\sigma_{\vec{u}}$ with Heisenberg-limited precision scaling using its dynamics $e^{-iH\tau}$ $(\tau>0)$.

In general, all three applications can be performed even if the dynamics $e^{-iH\tau}$ is given as a black box, apart from knowledge of $\Delta_{H}$. However, given the knowledge that $H$ belongs to a subspace of $\mathcal{L}(\mathcal{H})$ spanned by the set $\{\sigma_{\vec{v}}\}_{\vec{v}\in J}$ for some $J\subset{\{0,1,2,3\}^{n}}$, negative time-evolution and time-reversal can be performed in a runtime of $O[\mathrm{poly}(|J|)]$. This property is useful when the Hamiltonian is known to be $k$-local for some constant $k$, in which case $J=\left\{\vec{w}\ :||\vec{w}||_{0}\leq k\right\}$ satisfies $|J|\sim O(n^{k})$, so that the overall runtime is polynomial in the system size $n$, based on the fact that $\Delta_{H}$ is also ${\rm poly}(n)$.

In quantum algorithms that make direct use of Hamiltonian dynamics, both the positive and negative time-evolution are often assumed to be readily accessible. For example, this is required in the recent Hamiltonian singular value transformation [8]. However, in practice, a Hamiltonian evolution being native to a given hardware does not automatically guarantee that the same is true for the corresponding negative time-evolution. Therefore, the ability to efficiently simulate the negative time-evolution of any Hamiltonian given as a black box can decrease the resources required for such algorithms. On the more foundational side, given access to a black box Hamiltonian evolution, one might be interested in simulating the corresponding time-reversed evolution. For example, the evolution of an antiparticle can be described by the time-reversal of the corresponding particle evolution [33, 34].

The simulations of both negative time-evolution and time-reversal are performed by choosing the function $f$ as $f^{\rm neg}(H):=-H$ and $f^{\rm rev}(H):=H^{T}$, respectively, which are specified by

where $s_{(w_{1},\ldots,w_{n})}:=|\{j\in\{1,\ldots,n\}\mid w_{j}=2\}|$. In the definition of $\gamma^{\rm rev}_{\vec{w},\vec{u}}$, the fact that $I^{T}=I$, $X^{T}=X$, $Y^{T}=-Y$, and $Z^{T}=Z$ are used.

In both of these cases, $\beta=2\sum_{\vec{w},\vec{u}}|\gamma_{\vec{w},\vec{u}}|=2(4^{n}-1)$, thus the runtime $O(\beta^{2}t^{2}\Delta_{H}^{2}n/\epsilon)$ is exponential in $n$ in general. However, when $H$ is in a subspace of $\mathcal{L}(\mathcal{H})$ spanned by the set $\{\sigma_{\vec{v}}\}_{\vec{v}\in J}$ , we can define

since $f(H)$ does not depend on values of $\gamma_{\vec{w},\vec{u}}$ for $\vec{u}\notin J$. In this case, $\beta=2|J|$ so the runtime scales as $O(|J|^{2}t^{2}\Delta_{H}^{2}n/\epsilon)$, which is $O[\mathrm{poly}(n)]$ for a realistic Hamiltonian whose number of terms $|J|$ is polynomial in the system size $n$. For a general Hamiltonian linear transformation $f$, if both the seed Hamiltonian and the transformed Hamiltonian have a polynomial number of terms in $n$, then the non-zero elements of $f$ can be truncated so that the runtime $O(\beta^{2}t^{2}\Delta_{H}^{2}n/\epsilon)$ has a polynomial dependence on $n$.

We note that the runtime scales as $t^{2}$, meaning that in order to perform the time-reversal or negative time-evolution by this algorithm, the dynamics is slowed down quadratically. An application of simulating the negative time-evolution to Hamiltonian block encoding [8] is described in Appendix D.

Finally, we consider an application of our algorithm to Hamiltonian single-parameter learning. Estimation techniques of parameters of unknown Hamiltonians for Hamiltonian learning have many applications in quantum sensing [35], analyzing properties of quantum many-body physics [36], and quantum device calibration [37]. Recently, an estimation technique achieving the Heisenberg limit for the precision scaling in the estimation of parameters of a low-interaction Hamiltonian utilizing transformations of Hamiltonian dynamics has been proposed [32]. Our algorithm can be used to extend similar techniques to a more general class of $n$-qubit Hamiltonians.

Our estimation algorithm consists of two steps. The first step simulates $e^{-if_{\vec{v}}(H_{0})t}$ $(t>0)$ using the Hamiltonian dynamics $e^{-iH\tau}$ $(\tau>0)$ where $\vec{v}$ specifies $c_{\vec{v}}$ that we want to estimate and $f_{\vec{v}}$ is a hermitian-preserving linear map chosen as $f_{\vec{v}}(H)=c_{\vec{v}}Y\otimes I\otimes\cdots\otimes I$. This function $f_{\vec{v}}$ “filters” to keep only the coefficient $c_{\vec{v}}$ and changes all other coefficients to be zero, and then sends the coefficient $c_{\vec{v}}$ to the coefficient of $Y\otimes I\otimes\cdots\otimes I$, which is chosen for the convenience of the second step. The corresponding $\gamma$ is given by $\gamma_{\vec{w}\vec{u}}:=\delta_{\vec{w},(2,0,\ldots,0)}\delta_{\vec{u},\vec{v}}$. The second step performs robust phase estimation [38] using $e^{-if_{\vec{v}}(H)t}$ similarly to the technique in [32] to obtain an estimate for $c_{\vec{v}}$, by measuring only the first qubit in our algorithm. The total evolution time is $O((\log\delta)/\epsilon)$ where $\epsilon$ is precision and $\delta$ is the failure probability, which achieves the Heisenberg-limited precision scaling. The detailed procedure and analysis of the total evolution time are given in Appendix E.

For parameter estimation of low-interaction Hamiltonians, the method of [32] can perform the full-parameter estimation in a single run with total evolution time $O((\log\delta)/\epsilon)$, while our method requires polynomially longer total evolution time as we need to repeat the single parameter estimation for every parameter to perform the same task. However, the method of [32] requires exponential total evolution time for estimating a high-interaction coefficient (a coefficient of a $k$-local Hamiltonian term with $k=O(n)$), while our algorithm requires the same total evolution time for any coefficient. Therefore, our algorithm is suitable for estimating a single parameter of non-local Hamiltonians.

Summary and outlook.— We presented a universal algorithm that can simulate any linear physically realizable hermitian-preserving transformation of any Hamiltonian dynamics given as a black box. Our algorithm requires only a finite number of calls to the black box Hamiltonian dynamics and random pairs of correlated controlled-Pauli gates. We showed how our algorithm can simulate both the time-reversal and negative time-evolution of any unknown Hamiltonian dynamics, as well as an application to Hamiltonian single-parameter learning, efficiently estimating a single parameter of a multi-parameter Hamiltonian.

In our algorithm, the probability distribution for choosing multiple gates at different time steps are correlated in the sense that the gate $V_{f,j}$ is always used together with its adjoint $V_{f,j}^{\dagger}$, and the probabilities for picking its component controlled-Pauli gates are correlated via a joint probability distribution. This algorithm demonstrates how multiply correlated randomness can be leveraged to construct unitary operators without introducing decoherence. Our algorithm is a starting point for the emerging field of black box Hamiltonian simulation. One possible future direction is to extend higher-order quantum transformations of Hamiltonian dynamics to Hamiltonian transformations beyond hermitian-preserving linear transformations.

We would like to thank Seiseki Akibue, Mile Gu, Toshinori Itoko, Antonio Mezzacapo, Kunal Sharma, Philip Taranto, and Satoshi Yoshida for fruitful discussions. This work was supported by the MEXT Quantum Leap Flagship Program (MEXT QLEAP) Grant No. JPMXS0118069605, the Japan Society for the Promotion of Science (JSPS) KAKENHI Grants No. 21H03394 and No. 18K13467, and partly by IBM-UTokyo lab. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Colleges and Universities.