The Kitaev–Feynman clock for open quantum systems

The notion that quantum computers can efficiently simulate the time-evolution of other quantum systems was originally pioneered by Feynman [1–3] and has since spurred a plethora of experimental and theoretical work in the field of quantum computation and simulation. Feynman and later Kitaev [4] envisioned a quantum simulator, in which the simulated system is entangled with a clock particle and the entire history of the simulation is encoded as the ground state of a Hamiltonian. This construction, which we refer to as the 'Kitaev–Feynman clock' is particularly appealing, since it enables simulation of time-dependent quantum mechanics on a quantum computer using a time-independent setup [5, 6, 10]. The Kitaev–Feynman clock has also been an important tool in proving theorems, such as the equivalence of the adiabatic and gate models of quantum computation [6].

The Kitaev–Feynman clock in itʼs original formulation is restricted to isolated quantum systems evolving unitarily. However, many quantum systems of interest in chemistry and physics are not isolated, but undergo energy exchange and decoherence due to interaction with a thermal environment or measuring apparatus. Examples include energy transfer in photosynthetic and excitonic complexes [11], condensed phase spectroscopy and cavity quantum electrodynamics to name a few. In all these systems, there still exists sufficient quantum coherence that one expects quantum mechanics to be important, but interactions with the environment are certainly not negligible. Several important studies have explored the possibility of using an open quantum system to simulate the evolution of another open quantum system [7–9].

In the present manuscript, we will present a construction analogous to the Kitaev–Feynman clock, but which unlike the original construction is applicable to open quantum systems. This construction is useful for a number of reasons. First, by mapping the open-system dynamics onto a ground state problem, it becomes time-independent and variational. This can be used to develop computational methods, which simulate open quantum systems on classical computers. This use of the Kitaev–Feynman clock for closed quantum systems was presented in [5]. Second, it is a useful tool for proving theorems about quantum computation, when the dynamics are no longer assumed to be unitary [8]. In this way, many of the proofs that use the Kitaev–Feynman clock for unitary evolution might be extendable to open quantum systems. Third, as large-scale quantum computing devices are experimentally realized, it may be possible to build a Kitaev–Feynman clock as a quantum simulator. In fact, this is what Feynman originally had imagined. For many realistic applications, a Kitaev–Feynman clock would need to be able to not only simulate a quantum system, but also its interaction with an environment.

The manuscript is organized as follows. In section 2, we review the Kitaev–Feynman clock for unitary evolution and also the stochastic unraveling of the Lindblad master equation. Section 3 presents the formal theory behind the Kitaev–Feynman clock for the Lindblad master equation. In section 4, the formal theory is demonstrated with a numerical study of a two-level atom undergoing spontaneous emission. Section 5 provides a conclusion by discussing experimental implementations and extensions of the theory to non-Markovian systems. We have set $\hbar =1$ throughout, unless specified otherwise.

Our goal in section 3 will be to construct an ensemble of Kitaev–Feynman clock Hamiltonians, which will have as their ground states the stochastic trajectories that unravel a Lindblad master equation. As a prelude, in this section we will separately review the Kitaev–Feynman clock for unitary evolution and the stochastic Schrodinger equation (SSE) method of evolving the Lindblad equation.

2.1. The Kitaev–Feynman clock

Feynmanʼs original clock construction assumes an ideal quantum simulator, which simulates another system described by a wave function evolving under the time-dependent Schrödinger equation, ${\rm i}\frac{\partial }{\partial t}|\psi (t)\rangle =\hat{H}|\psi (t)\rangle$, whose solution is $|\psi (t)\rangle ={{{\rm e}}^{-{\rm i}\hat{H}t}}|\psi (0)\rangle$, given an initial state $|\psi (0)\rangle$. The entire simulation is then encoded as a superposition in a 'history state' $|\eta \rangle =\sqrt{\frac{\delta t}{T+\delta t}}\sum _{t=0}^{T}|\psi (t)\rangle \otimes |t\rangle$, where $|t\rangle$ denotes the state of an auxiliary quantum system known as the 'clock,' used to keep track of the evolution over the time of the simulation from t = 0 to t = T. The clock can be any quantum degree of freedom, discrete or continuous, that can be coupled with the system we wish to simulate. By performing a projective measurement of the clock at a specific time of interest t, the history state collapses to the wave function at that time, $|\psi (t)\rangle$.

In the construction of Kitaev [4], $|\eta \rangle$ can be encoded as the ground state of the Hamiltonian

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{H} & = & \mathop{\sum }\limits_{t=0}^{T-\delta t}-U\otimes |t+\delta t\rangle \langle t|-{{U}^{\dagger }}\otimes |t\rangle \langle t+\delta t|+|t\rangle \langle t| \\ {} & {} & +|t+\delta t\rangle \langle t+\delta t|+(1-|\psi (0)\rangle \langle \psi (0)|)\otimes |0\rangle \langle 0|, \\ \end{array}\end{eqnarray} \tag{ 1 }$$

where $U={{{\rm e}}^{-{\rm i}\hat{H}\delta t}}$ and $\delta t$ is the time-step which represents the distance between sites of the clock, assumed to be discrete. The first four terms in equation (1) ensure the history state encodes the correct time evolution, while the last term enforces the correct initial state. It can be readily verified that $\mathcal{H}|\eta \rangle =0$ and because $\mathcal{H}$ is positive semidefinite, $|\eta \rangle$ is the unique ground state of $\mathcal{H}$ with eigenvalue 0. In general, while the ground state encodes the history state, excited states of the Kitaev–Feynman clock do not have an obvious physical interpretation.

2.2. Stochastic unraveling of the Lindblad master equation

We wish to use a Kitaev–Feynman clock construction to simulate an open quantum system, described by a density matrix evolving under the widely used Markovian Lindblad master equation [15]

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm dt}}{{\hat{\rho }}_{s}}={\rm i}[{{\hat{\rho }}_{s}},{{\hat{H}}_{s}}]-\frac{1}{2}\mathop{\sum }\limits_{m}({{\hat{C}}_{m}}\hat{C}_{m}^{\dagger }{{\hat{\rho }}_{s}}+{{\hat{\rho }}_{s}}{{\hat{C}}_{m}}\hat{C}_{m}^{\dagger })+\mathop{\sum }\limits_{m}{{\hat{C}}_{m}}{{\hat{\rho }}_{s}}\hat{C}_{m}^{\dagger }.\end{eqnarray} \tag{ 2 }$$

Here, ${{\hat{\rho }}_{s}}$ is the density matrix of the system, ${{\hat{H}}_{s}}$ is the system Hamiltonian, ${{\hat{C}}_{m}}$ and $\hat{C}_{m}^{\dagger }$ describe interactions of the system with its environment. Because equation (2) describes the evolution of a mixed-state density matrix and not a pure-state wave function, a straightforward generalization of the Kitaev–Feynman clock (equation (1)) is not possible. However, the SSE procedure deals directly with wave functions and therefore serves as an ideal starting point for extending the Kitaev–Feynman clock to open quantum systems [12–14].

The SSE simulates a set of m individual realizations of the open quantum system, whose wave functions $\{|{{\psi }^{i}}(t)\rangle \}$ yield the ensemble averaged density matrix in equation (2) by averaging over stochastic trajectories according to

$$\begin{eqnarray}&&{{\hat{\rho }}_{s}}(t)=\frac{1}{m}\mathop{\sum }\limits_{i=1}^{m}|{{\psi }^{i}}(t)\rangle \langle {{\psi }^{i}}(t)|.\end{eqnarray} \tag{ 3 }$$

It can be shown that in the limit $m\to \infty$, equation (3) converges to the exact density matrix [14]. Each realization can be thought of as simulating a single experiment on an individual member of the ensemble, such as a single-molecule fluorescence measurement or single-atom–photon detection experiment. The actual number of SSE trajectories required to yield a good approximation in equation (3) depends on the particular system and the type of operators whose expectation values one wishes to calculate. For local operators, which depend for instance on only a single atom in a large sample, one must have $m\gg n$, where n is the size of the atomʼs Hilbert space. For global operators such as the average kinetic energy of an entire sample of atoms, one must have $m\gg 1$. Therefore, the SSE becomes optimal for calculating global properties of many-body systems [12]. For a detailed discussion of the convergence properties of the SSE with rigorous bounds, we refer the reader to section 7 of [12].

In the SSE, the set of wave functions $\{|{{\psi }^{i}}(t)\rangle \}$ are simulated as follows. At time t, one evolves with a non-Hermitian Hamiltonian to $t+\delta t$ according to

$$\begin{eqnarray}&&|{{\psi }^{i}}(t+\delta t)\rangle =\left[ 1-{\rm i}{{{\hat{H}}}_{s}}\delta t-\frac{1}{2}\mathop{\sum }\limits_{m}\hat{C}_{m}^{\dagger }{{{\hat{C}}}_{m}} \right]\frac{|{{\psi }^{i}}(t)\rangle }{\sqrt{1-\delta p(t)}},\end{eqnarray} \tag{ 4 }$$

with probability $1-\delta p(t)$, where $\delta p(t)={{\sum }_{m}}\delta {{p}_{m}}(t)$ and $\delta {{p}_{m}}(t)=\delta t\langle {{\psi }^{i}}(t)|\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}|{{\psi }^{i}}(t)\rangle$. On the other hand, with probability $\delta p(t)$, one instead collapses the wave function to the state

$$\begin{eqnarray}&&|{{\psi }^{i}}(t+\delta t)\rangle =\frac{{{{\hat{C}}}_{m}}|{{\psi }^{i}}(t)\rangle }{\sqrt{\delta {{p}_{m}}(t)/\delta t}},\end{eqnarray} \tag{ 5 }$$

chosen from among the bath operators ${{\hat{C}}_{m}}$ with probability $\delta {{p}_{m}}(t)$. $\delta p(t)$ is typically small, so the majority of the time the system evolves freely without interaction with the environment. The Hamiltonian is non-Hermitian, because by learning that the system has not interacted with the environment, we have gained information, i.e. indirectly measured the system. Occasionally, with a small probability $\delta p(t)$, the system interacts with its environment causing wave function collapse.

In this section we develop a formalism where history states are constructed, which correspond to the non-deterministic evolution of the stochastic trajectories discussed in section 2.2. Each history state is encoded as the ground state of a non-Hermitian Hamiltonian, chosen probabilistically to enforce the correct stochastic jump probabilities. By ensemble averaging these history states, one can recover the entire history of the density matrix in equation (2). Before including stochastic jumps, however, we will begin with a deterministic, but non-Hermitian clock describing the free evolution in equation (4).

3.1. The non-Hermitian Kitaev–Feynman clock

Our goal is to construct a history state, $|\eta \rangle =\sqrt{\frac{\delta t}{T+\delta t}}\sum _{t=0}^{T}|\psi (t)\rangle \otimes |t\rangle$, where $|\psi (t)\rangle$ is identical to the evolution described by equation (4). Such a history state would encode the non-Hermitian evolution of an open system on the interval [0, T], when no stochastic jumps are generated. For an atom–photon experiment for example, $|\eta \rangle$ encodes the history of an atom, which is observed to evolve freely with no photons emitted. This history state is explicitly given by

$$\begin{eqnarray}&&|\eta \rangle =\sqrt{\frac{\delta t}{T+\delta t}}\mathop{\sum }\limits_{t=0}^{T}\frac{{{(\hat{R})}^{t}}}{\sqrt{1-\delta p(t)}}|\psi (0)\rangle \otimes |t\rangle ,\end{eqnarray} \tag{ 6 }$$

where $\hat{R}=1-{\rm i}{{\hat{H}}_{s}}\delta t-\frac{1}{2}\delta t{{\sum }_{m}}\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}$, up to corrections of order δt². This is equivalent to keeping only the first-order term in a Suzuki–Trotter expansion [16, 17]. In order to construct a Kitaev–Feynman clock, we need to write this history state as the ground state of a Hamiltonian. It can be readily verified that the history state in equation (6) satisfies $\mathcal{H}|\eta \rangle =0$ where

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{H} & = & \mathop{\sum }\limits_{t=0}^{T-\delta t}-\hat{R}|t+\delta t\rangle \langle t|-{{{\hat{R}}}^{-1}}|t\rangle \langle t+\delta t|+|t\rangle \langle t| \\ {} & {} & +|t+\delta t\rangle \langle t+\delta t|+(1-|\psi (0)\rangle \langle \psi (0)|)\otimes |0\rangle \langle 0|, \\ \end{array}\end{eqnarray} \tag{ 7 }$$

and ${{\hat{R}}^{-1}}=1+{\rm i}{{\hat{H}}_{s}}\delta t+\frac{1}{2}\delta t{{\sum }_{m}}\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}$. The Hamiltonian in equation (7) is non-Hermitian, so it is not immediately obvious that $|\eta \rangle$ is the ground state or if this even has any meaning if the spectrum were complex. However, in the appendix we show that $\mathcal{H}$ has a complete, non-degenerate and real spectrum of positive eigenvalues. Therefore, $|\eta \rangle$ is in fact the ground state and we have succeeded in constructing a Kitaev–Feynman clock, which encodes the free evolution of an open system when no jumps occur.

3.2. The stochastic Kitaev–Feynman clock

So far we have constructed a non-Hermitian Kitaev–Feynman clock describing the free SSE evolution, without stochastic jumps. We now show that it is possible to generate an ensemble of history states $\{|{{\eta }^{i}}\rangle \}$ where

$$\begin{eqnarray}&&|{{\eta }^{i}}\rangle =\sqrt{\frac{\delta t}{T+\delta t}}\mathop{\sum }\limits_{t=0}^{T}|{{\psi }^{i}}(t)\rangle \otimes |t\rangle ,\end{eqnarray} \tag{ 8 }$$

is the history state of the ith stochastic trajectory and the set $\{|{{\psi }^{i}}(t)\rangle \}$ are identical to those obtained from the SSE procedure in equations (4) and (5). This is done by choosing an ensemble of stochastic Kitaev–Feynman clock Hamiltonians $\{{{\mathcal{H}}^{i}}\}$ according to the following procedure. For each i, write ${{\mathcal{H}}^{i}}$ as a sum of local terms according to,

$$\begin{eqnarray}&&{{\mathcal{H}}^{i}}=\mathop{\sum }\limits_{t=0}^{T-\delta t}{{h}^{i}}(t+\delta t)+(1-|\psi (0)\rangle \langle \psi (0)|)\otimes |0\rangle \langle 0|.\end{eqnarray} \tag{ 9 }$$

Then for each value of t one chooses these terms to be

$$\begin{eqnarray}\begin{array}{rcl} {{h}^{i}}(t+\delta t) & = & -\sqrt{\frac{1-\delta {{p}^{i}}(t)}{1-\delta {{p}^{i}}(t+\delta t)}}\hat{R}(\delta t)|t+\delta t\rangle \langle t| \\ {} & {} & -\sqrt{\frac{1-\delta {{p}^{i}}(t+\delta t)}{1-\delta {{p}^{i}}(t)}}{{{\hat{R}}}^{-1}}(\delta t)|t\rangle \langle t+\delta t| \\ {} & {} & +|t\rangle \langle t|+|t+\delta t\rangle \langle t+\delta t|, \\ \end{array}\end{eqnarray} \tag{ 10 }$$

with probability $1-{{p}^{i}}(t)$, where $\delta {{p}^{i}}(t)=\delta t{{\sum }_{m}}\langle {{\psi }^{i}}(t)|\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}|{{\psi }^{i}}(t)\rangle$. On the other hand, with probability ${{p}^{i}}(t)$ one instead chooses

$$\begin{eqnarray}&&{{h}^{i}}(t+\delta t)=(1-\frac{\delta t}{\delta p_{m}^{i}}{{\hat{C}}_{m}}|{{\psi }^{i}}(t)\rangle \langle {{\psi }^{i}}(t)|\hat{C}_{m}^{\dagger })\otimes |t+\delta t\rangle \langle t+\delta t|,\end{eqnarray} \tag{ 11 }$$

from among the various bath operators $\hat{C}_{m}^{\dagger }$ with probability $\delta p_{m}^{i}=\delta t\langle {{\psi }^{i}}(t)|\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}|{{\psi }^{i}}(t)\rangle$. The terms in equation (10) force the ground state to have free evolution at the specified times, while the terms in equation (11) force the ground state to have the appropriate collapsed wave function corresponding to a jump. For a given Hamiltonian ${{\mathcal{H}}^{i}}$ generated in this way, it can be verified by direct substitution that the state $|{{\eta }^{i}}\rangle =\sqrt{\frac{\delta t}{T+\delta t}}\sum _{t=0}^{T}|{{\psi }^{i}}(t)\rangle \otimes |t\rangle$ satisfies ${{\mathcal{H}}^{i}}|{{\eta }^{i}}\rangle =0,$ provided the states $|{{\psi }^{i}}(t)\rangle$ are generated from the SSE with the same realization of jumps. The set of history states, $\{|{{\eta }^{i}}\rangle \}$, encode the entire evolution of individual trajectories and can be used to compute ensemble averages. The density matrix of equation (2) can be obtained by making projective measurements of the clock and then ensemble averaging these measurements according to

$$\begin{eqnarray}&&{{\hat{\rho }}_{s}}(t)=\frac{1}{m}\mathop{\sum }\limits_{i=1}^{m}Tr[|{{\eta }^{i}}\rangle \langle {{\eta }^{i}}||t\rangle \langle t|]=\frac{1}{m}\mathop{\sum }\limits_{i=1}^{m}|{{\psi }^{i}}(t)\rangle \langle {{\psi }^{i}}(t)|.\end{eqnarray} \tag{ 12 }$$

We see that the Hamiltonian in equation (9) is a nonlinear functional of the state $|{{\eta }^{i}}\rangle$. This nonlinearity enters in two ways. First, there is an implicit nonlinearity, because the choice of terms in the Hamiltonian is determined probabilistically from the state $|{{\eta }^{i}}\rangle$. Second, there is an explicit dependence on $|{{\eta }^{i}}\rangle$ appearing in the terms in equations (10) and (11). In the appendix, we show that despite this nonlinearity, the spectrum of each ${{\mathcal{H}}^{i}}$ is strictly real and there exists a corresponding $|{{\eta }^{i}}\rangle$ that is the ground state with eigenvalue 0. Because finding the spectrum of the stochastic Kitaev–Feynman clock is a nonlinear eigenvalue problem, it needs to be solved self-consistently to obtain an exact solution. However, in the next section we consider a perturbative expansion, which yields a valid approximation when the system–bath interaction is weak.

3.3. Perturbative expansion of the stochastic Kitaev–Feynman clock

As discussed in the previous section, the exact stochastic history states are solutions to a nonlinear eigenvalue problem, which must be solved self-consistently. However, it is often the case that the environment interacts only weakly with the system. In these situations, one expects that a majority of terms in the Hamiltonian will be of the form in equation (10) describing free evolution. Only occasionally does a jump occur, with a term of the form in equation (11) appearing. We can therefore linearize the stochastic Kitaev–Feynman clock and develop a perturbative expansion about the free evolution.

One first solves the linear eigenvalue problem for the non-Hermitian Hamiltonian in equation (7), and obtains the history state in equation (6), which describes free evolution when no jumps occur. We denote this state $|{{\eta }_{0}}\rangle$, which serves as the zeroth-order history state in our perturbative expansion. From $|{{\eta }_{0}}\rangle$, the set of jump probabilities $\{\delta {{p}_{m}}(0)\},\{\delta {{p}_{m}}(\delta t)\},\{\delta {{p}_{m}}(2\delta t)\},\ldots \{\delta {{p}_{m}}(T-\delta t)\}$ for all bath operators m at each time-step are generated as follows. First, the initial set of jump probabilities $\{\delta {{p}_{m}}(0)\}$ are obtained with a projective measurement of the clock particle at t = 0 and simultaneous measurement of the bath operator $\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}$ according to $\delta {{p}_{m}}(0)=\delta t\langle {{\eta }_{0}}|(\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}\otimes |0\rangle \langle 0|)|{{\eta }_{0}}\rangle =\delta t\langle \psi (0)|\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}|\psi (0)\rangle$. The remaining jump probabilities are obtained recursively, since the jump probabilities $\delta {{p}_{m}}(t+\delta t)$ at time $t+\delta t$, can be obtained from the history state $|{{\eta }_{0}}\rangle$ and jump probabilities at earlier times $\{\delta {{p}_{m}}(0)\},\{\delta {{p}_{m}}(\delta t)\},\ldots ,\{\delta {{p}_{m}}(t)\}$ through the relation

$$\begin{eqnarray}&&\delta {{p}_{m}}(t+\delta t)=\frac{\delta t}{1-\mathop{\sum }\limits_{t^{\prime} =0}^{t}\delta p(t^{\prime} )}\langle {{\eta }_{0}}|(\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}\otimes |t+\delta t\rangle \langle t+\delta t|)|{{\eta }_{0}}\rangle ,\end{eqnarray} \tag{ 13 }$$

where $t^{\prime} \lt t$. This entire procedure necessitates storage of $M\frac{T}{\delta t}$ copies of the state $|{{\eta }^{0}}\rangle$, where M is the number of bath operators, i.e. the algorithm is polynomial in the size of the systemʼs Hilbert space and the run time of the simulation.

From the set of jump probabilities just obtained and ground state $|{{\eta }^{0}}\rangle$, we now generate an ensemble $\{{{\mathcal{H}}^{i}}\}$ of 'single-jump' stochastic clock Hamiltonians. This is done by writing the ith Hamiltonian in the ensemble as a sum of terms acting locally as in section 3.2, ${{\mathcal{H}}^{i}}\equiv \sum _{t=0}^{T-\delta t}{{h}^{i}}(t+\delta t)$. Using the probabilities obtained in equation (13), we generate the ensemble by choosing

$$\begin{eqnarray}\begin{array}{rcl} {{h}^{i}}(t+\delta t) & = & -\sqrt{\frac{1-\delta p(t)}{1-\delta p(t+\delta t)}}\hat{R}(\delta t)|t+\delta t\rangle \langle t| \\ {} & {} & -\sqrt{\frac{1-\delta p(t+\delta t)}{1-\delta p(t)}}{{{\hat{R}}}^{-1}}(\delta t)|t\rangle \langle t+\delta t| \\ {} & {} & +|t\rangle \langle t|+|t+\delta t\rangle \langle t+\delta t|, \\ \end{array}\end{eqnarray} \tag{ 14 }$$

with probability $\delta p(t)=1-{{\sum }_{m}}\delta {{p}_{m}}(t)$. On the other hand, one chooses

$$\begin{eqnarray}&&{{h}^{i}}(t+\delta t)=(1-\frac{\delta t}{\delta {{p}_{m}}(t)}){{\hat{C}}_{m}}|{{\psi }_{0}}(t)\rangle \langle {{\psi }_{0}}(t)|\hat{C}_{m}^{\dagger }\otimes |t+\delta t\rangle \langle t+\delta t|\end{eqnarray} \tag{ 15 }$$

with probability $\delta {{p}_{m}}(t)$, where $|{{\psi }_{0}}(t)\rangle \equiv \langle t|{{\eta }_{0}}\rangle$ . Once a jump has occurred, terms at later times are chosen to have the form in equation (10), which enforces the one-jump approximation. The terms in equation (14) force the ground state of ${{\mathcal{H}}^{i}}$ to have the correct free evolution before and after the jump, while the terms in equation (15) generate energy penalties that enforce the jumps with the correct probabilities. The Hamiltonian ${{\mathcal{H}}^{i}}$ generated this way is block diagonal, with each block corresponding to the free evolution before and after the jump has occurred. The ground state of ${{\mathcal{H}}^{i}}$ will have eigenvalue zero and be two-fold degenerate, with one eigenstate corresponding to evolution before the jump and the other eigenstate to evolution after the jump. The physical history state, given as an equal superposition of these two degenerate states yields a single stochastic trajectory, $|\eta _{1}^{i}\rangle =\sqrt{\frac{\delta t}{T+\delta t}}\sum _{t=0}^{T}|\psi _{1}^{i}(t)\rangle \otimes |t\rangle$. The set of history states $\{|\eta _{1}^{i}\rangle \}$ are first-order in the system–bath interaction. The set of wave functions $\{|\psi _{1}^{i}(t)\rangle \}$ are precisely the wave functions from a SSE evolution, when only one jump has occurred.

It is clear that the single-jump stochastic clock Hamiltonians depend only on $|{{\eta }_{0}}\rangle$, so we need only solve a set of linear eigenvalue equations to to obtain the first-order history states $\{|\eta _{1}^{i}\rangle \}$. Similarly, the set of first-order history states can be used to generate linear eigenvalue equations for a set of second-order history states $\{|\eta _{2}^{i}\rangle \}$, which describe stochastic trajectories in which two jumps occur. The expansion can be continued up to order $N=\frac{T}{\delta t}$, which is equivalent to the solution of the nonlinear eigenvalue equations in section 3.2. In the appendix, we bound errors involved in truncating the expansion. Presently, we demonstrate the formal theory with an exactly solvable model system, in which the perturbation expansion to first-order is exact.

As a simple numerical demonstration, we construct the stochastic Kitaev–Feynman clock for a two-level atom undergoing spontaneous emission, as determined by its interaction with a photon detector [12, 13]. If no photon is detected at time t, spontaneous emission has not occurred and the atom continues to evolve freely to time $t+\delta t$. If on the other hand a photon is detected, spontaneous emission has occurred and the measurement causes the atomic wave function to jump to the ground state. In this situation the one-jump approximation is exact, because measurement causes the wave function to collapse to an eigenstate of the atomic Hamiltonian, leaving no possibility for a second jump to occur. The atomic wave function can be expanded in terms of the ground state $|0\rangle$ and excited state $|1\rangle$ as $|\psi (t)\rangle =\alpha (t)|0\rangle +\beta (t)|1\rangle$. The system Hamiltonian is then given by ${{\hat{H}}_{s}}=\omega |1\rangle \langle 1|$, where ω is the excitation energy. There is a single jump operator resulting from the measurement process given by $\hat{C}=\sqrt{\Gamma }|0\rangle \langle 0|$, where Γ is the emission rate.

An ensemble of two-level atoms in this setup is described by the Lindblad master equation

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}{{\hat{\rho }}_{s}}={\rm i}[{{\hat{\rho }}_{s}},{{\hat{H}}_{s}}]-\frac{\Gamma }{2}({{\hat{\sigma }}^{+}}{{\hat{\sigma }}^{-}}{{\hat{\rho }}_{s}}+{{\hat{\rho }}_{s}}{{\hat{\sigma }}^{+}}{{\hat{\sigma }}^{-}})+\Gamma {{\hat{\sigma }}^{-}}{{\hat{\rho }}_{s}}{{\hat{\sigma }}^{+}},\end{eqnarray} \tag{ 16 }$$

where ${{\hat{\sigma }}^{+}}=|1\rangle \langle 0|$ and ${{\hat{\sigma }}^{-}}=|0\rangle \langle 1|$. The above density matrix is obtained by averaging over stochastic trajectories of the SSE as described in section 2.2.

For the simulation, we choose the initial state to be an equal superposition of the atomic ground state and first excited state, i.e. $\alpha (0)=\beta (0)=\frac{1}{\sqrt{2}}$. The deterministic non-Hermitian clock Hamiltonian in equation (6) is constructed and the ground state $|{{\eta }_{0}}\rangle =\sqrt{\frac{\delta t}{T+\delta t}}\sum _{t=0}^{T}|{{\psi }_{0}}(t)\rangle \otimes |t\rangle$ is obtained via an exact numerical diagonalization of the clock Hamiltonian, ${{\mathcal{H}}_{0}}$. The populations and coherence obtained from the coefficients $|{{\psi }_{0}}(t)\rangle =\langle t|{{\eta }_{0}}\rangle$ are shown in figure 1 as a function of $|t\rangle$. As expected, the free non-Hermitian evolution causes the excited state population and coherence to decay, while the ground state population of the atom remains unchanged. Although we have solved for the ground state of a time-independent Hamiltonian, we get exactly what we would expect from a propagation of the SSE with no jumps. Figure 2 shows the spectrum of ${{\mathcal{H}}_{0}}$, which has a single non-degenerate zero eigenvalue corresponding to $|{{\eta }_{0}}\rangle$. In general, the excited states of ${{\mathcal{H}}_{0}}$ are unphysical and only the ground state encodes a valid evolution, which is true of the unitary Kitaev–Feynman clock as well [5]. Figure 3 shows populations and coherences of a single stochastic Kitaev–Feynman clock, in which a jump penalty (equation (15)) corresponding to a photon detection has been imposed at $t=1.0(\hbar \omega )$. The ground state is two-fold degenerate, with one ground state corresponding to free evolution before the jump and the other after the jump has caused the wave function to collapse to the atomic ground state. In figure 4, one sees that in general, the spectrum is divided into pairs of degenerate states for every energy level. For each pair, one state is localized before the jump occurs, while the second is localized afterward. This happens because the clock Hamiltonian assumes a block diagonal form, with one block acting only on clock states $\{|t\rangle \}$ before the jump and the other block acting only on states after the jump.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5(a) shows the elements of the density matrix obtained from averaging the ground states of 20 stochastic Kitaev–Feynman clock Hamiltonians, for a runtime of $T=1.0(\hbar \omega )$. The density matrix is seen to be nearly identical to that obtained from a conventional propagation of the SSE using the same realization of the random number generator, shown in figure 5(b). This serves to demonstrate that the stochastic Kitaev–Feynman clock exactly reproduces the unraveling of the Lindblad equation with the SSE (equation (16)).

Zoom In Zoom Out Reset image size

In a realistic experimental setup one expects imperfections in the stochastic clock Hamiltonian to cause the history state to deviate from the open-system dynamics simulated in real-time by a SSE. To study the effects of such imperfections, we include diagonal static disorder by choosing the stochastic clock Hamiltonians to have the form ${{\mathcal{H}}^{i}}+\delta$, where δ is a positive-semidefinite diagonal random matrix. In figure 6, we see that for small runtimes the static disorder has little effect and the density matrix produced by ensemble averaging stochastic clock history states still agrees faithfully with that from the SSE evolution. However, as the runtime increases, the density matrix from the stochastic clock history states are seen to deviate appreciably from the SSE evolution. This occurs due to contamination from excited states of the stochastic clocks, which do not correspond to physical evolution. Specifically, the gap between the ground and first excited history state of the clock is proportional to $\frac{1}{{{T}^{2}}}$ [5], so the effects of noise become appreciable when ${{\delta }_{{\rm max} }}{{T}^{2}}\sim 1$, where ${{\delta }_{{\rm max} }}$ is the maximum eigenvalue of the matrix δ. However, the error grows only polynomially in the runtime and in general static disorder can be minimized more easily than dynamic disorder. From a quantum simulation standpoint, this might offer an advantage of the stochastic clock construction over conventional simulation in real-time, in which decoherence results during evolution.

Zoom In Zoom Out Reset image size

Experimental implementation of the stochastic Kitaev–Feynman clock should be possible with a variety of highly tunable quantum systems, such as ultracold atoms [18], superconducting qubits [19] and quantum dots [20]. For instance, with ultracold atoms in an optical lattice, the optical lattice can be tuned to a very large local potential at a randomly chosen site to impose stochastic jumps. In contrast, a perfectly periodic optical lattice in the superfluid regime will have delocalized eigenstates corresponding to a history state with free evolution. Non-Hermitian evolution can be generated by coupling sites representing the history state to a large number of other sites in the lattice, effectively generating a continuum. Furthermore, many stochastic history states can be generated in parallel in a single experiment and stored in disconnected sites of the optical lattice. The coefficients of the history states can be read and manipulated as required by our procedure using a high resolution quantum gas microscope as presented in [18].

From an experimental standpoint, the Kitaev–Feynman clock formulation of quantum simulation may offer an advantage over conventional simulation in real-time, for certain systems. By recasting the simulation as a time-independent problem, all quantum gates can be constructed as ground-state interaction terms and no ultrafast, real-time manipulations are needed. Such real-time manipulations can introduce unwanted decoherence, therefore the time-independent setup used in a Kitaev–Feynman clock may offer an advantage. However, as demonstrated in section 4, imperfections in the couplings introduce mixing of unphysical excited states, which become increasingly problematic as the gap between the ground and first-excited state grows with the simulation time. Therefore, one trades the problem of decoherence in conventional real-time quantum simulation, for contamination from excited states in the Kitaev–Feynman clock method. Which method is more beneficial depends on the specific system and experimental tunability of the couplings.

Many theoretical results in quantum computation rely on the unitary Kitaev–Feynman clock to prove theorems. We anticipate that the stochastic Kitaev–Feynman clock can be a useful tool for extending many of these results to open quantum systems. Examples include the equivalence of the adiabatic and circuit models of quantum computing [6, 21] and the complexity of k-local Hamiltonians [22]. We have focused specifically on Markovian environments, but it is also possible to construct non-Markovian stochastic Kitaev–Feynman clocks starting from a non-Markovian quantum jump model [23]. Also, we have chosen to work directly with wave functions by using the SSE as our starting point, rather than the master equation for the density matrix. This is the most straightforward generalization of the unitary Kitaev–Feynman clock, since one still works within a Hilbert space. However, this approach has the obvious disadvantage that the stochastic Kitaev–Feynman clock Hamiltonians are non-linear functionals acting in the Hilbert space. In future work, we will explore the possibility of formulating the Kitaev–Feynman clock in Liouville space, which allows one to derive linear equations for the density matrix [24].

We are grateful to Jarrod McClean, Joel Yuen-Zhou and Gian Giacomo Guerreschi for valuable discussions. We acknowledge NSF CHE-1152291 for financial support.

In section 3.1 we introduced the non-Hermitian clock $\mathcal{H}$ and found that the history state $|\eta \rangle =\sqrt{\frac{\delta t}{T+\delta t}}\sum _{t=0}^{T}|\psi (t)\rangle \otimes |t\rangle$ satisfies $\mathcal{H}|\eta \rangle =0$. Our goal is to show that the spectrum of $\mathcal{H}$ is real and positive so that $|\eta \rangle$ is in fact the ground state.

We first show that the operator $\mathcal{H}$ is normal and therefore admits a spectral decomposition [25]. We can write $\mathcal{H}=\mathcal{L}+\mathcal{T}$ as a sum of a Hermitian part $\mathcal{L}$ and anti-Hermitian part, $\mathcal{T}$. Explicitly,

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{L} & = & \mathop{\sum }\limits_{t=0}^{T-\delta t}-(1-{\rm i}\delta t\hat{H})|t+\delta t\rangle \langle t|-(1+{\rm i}\delta t\hat{H})|t\rangle \langle t+\delta t|+|t\rangle \langle t| \\ {} & {} & +|t+\delta t\rangle \langle t+\delta t|+(1-|\psi (0)\rangle \langle \psi (0)|)\otimes |0\rangle \langle 0| \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

and

$$\begin{eqnarray}&&\mathcal{T}=\delta t\mathop{\sum }\limits_{t=0}^{T-\delta t}[\hat{D}|t+\delta t\rangle \langle t|-\hat{D}|t\rangle \langle t+\delta t|],\end{eqnarray} \tag{ A.2 }$$

where $\hat{D}=\frac{1}{2}\delta t{{\sum }_{m}}\hat{C}_{m}^{\dagger }{{\hat{C}}_{m}}$. Working out the commutator of $\mathcal{H}$ with its adjoint ${{\mathcal{H}}^{\dagger }}$ and keeping only terms linear in $\delta t$ one finds

$$\begin{eqnarray}&&[{{\mathcal{H}}^{0}},{{({{\mathcal{H}}^{0}})}^{\dagger }}]\approx \delta t[\hat{D}(1-|{{\psi }_{0}}\rangle \langle {{\psi }_{0}}|)+(1-|{{\psi }_{0}}\rangle \langle {{\psi }_{0}}|)\hat{D}]\otimes |0\rangle \langle 0|.\end{eqnarray} \tag{ A.3 }$$

We notice two things. First, while the operator $\mathcal{H}$ is not normal strictly speaking, it becomes normal in the limit $\delta t\to 0$. Second, the commutator is non-zero only if the bath operator $\hat{D}$ does not commute with the projector ${{P}_{\psi 0}}\equiv |{{\psi }_{0}}\rangle \langle {{\psi }_{0}}|\otimes I$. Consequently

$$\begin{eqnarray}&&{{P}_{\psi 0}}[{{\mathcal{H}}^{0}},{{({{\mathcal{H}}^{0}})}^{\dagger }}]{{P}_{\psi 0}}=0,\end{eqnarray} \tag{ A.4 }$$

and $\mathcal{H}$ is normal in the subspace of history states satisfying the initial condition. Therefore, $\mathcal{H}$ has a complete spectrum of eigenstates $\{|{{\eta }_{k}}\rangle \}$ satisfying

$$\begin{eqnarray}&&{{\mathcal{H}}^{0}}|{{\eta }_{k}}\rangle ={{\lambda }_{k}}|{{\eta }_{k}}\rangle ,\end{eqnarray} \tag{ A.5 }$$

either in the limit $\delta t\to 0$ or to all orders in $\delta t$ within the projected space of ${{P}_{\psi 0}}$.

We now prove that all of the eigenvalues ${{\lambda }_{k}}$ are real and positive. To this end, we introduce an operator $\hat{O}$, which acts on the clock states as $\hat{O}|t\rangle =(1+t\delta t\hat{D})|t\rangle$. It has an inverse whose action is ${{\hat{O}}^{-1}}|t\rangle =(1-t\delta t\hat{D})|t\rangle$. We may then re-write equation (A.5) as

$$\begin{eqnarray}&&\hat{O}\mathcal{H}{{\hat{O}}^{-1}}\hat{O}|{{\eta }_{k}}\rangle ={{\lambda }_{k}}\hat{O}|{{\eta }_{k}}\rangle ,\end{eqnarray} \tag{ A.6 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} \hat{O}\mathcal{H}{{{\hat{O}}}^{-1}} & = & \mathop{\sum }\limits_{t=0}^{T-}\left[ -(1-{\rm i}\hat{H})|\delta t+\delta t\rangle \langle t| \right. \\ {} & {} & -\left. (1+{\rm i}\hat{H})|t\rangle \langle t+\delta t|+|t\rangle \langle t|+|t+\delta t\rangle \langle t+\delta t| \right] \\ {} & {} & +(1-|\psi (0)\rangle \langle \psi (0)|)\otimes |0\rangle \langle 0|. \\ \end{array}\end{eqnarray} \tag{ A.7 }$$

This transformed Hamiltonian is Hermitian, positive semidefinite and has eigenvalues between 0 and 4. Since it has an identical spectrum to the original Hamiltonian in equation (A.5), this is true of the Hamiltonian $\mathcal{H}$ as well.

We now show that despite being nonlinear functionals of their respective ground states, the stochastic Kitaev–Feynman clock Hamiltonians $\{{{\mathcal{H}}^{i}}\}$ have a real and positive spectrum. Consequently, for each realization the state $|{{\eta }^{i}}\rangle =\sqrt{\frac{\delta t}{T+\delta t}}\sum _{t=0}^{T}|{{\psi }^{i}}(t)\rangle \otimes |t\rangle$, which satisfies ${{\mathcal{H}}^{i}}|{{\eta }^{i}}\rangle =0$ is the ground state. To show this, consider a particular realization

$$\begin{eqnarray}&&{{\mathcal{H}}^{i}}=\mathop{\sum }\limits_{t=0}^{T-\delta t}{{h}^{i}}(t+\delta t)+(1-|\psi (0)\rangle \langle \psi (0)|)\otimes |0\rangle \langle 0|,\end{eqnarray} \tag{ B.1 }$$

which has local jump Hamiltonians of the form

$$\begin{eqnarray}&&h({{s}_{j}})=(1-\frac{\delta t}{\delta {{p}_{m}}}{{\hat{C}}_{m}}|\psi ({{s}_{j}}-\delta t)\rangle \langle \psi ({{s}_{j}}-\delta t)|\hat{C}_{m}^{\dagger })\otimes |{{s}_{j}}\rangle \langle {{s}_{j}}|,\end{eqnarray} \tag{ B.2 }$$

at a set of n times $\{{{s}_{j}}\}$ and free non-Hermitian evolution at other times. The Hamiltonian in equation (B.1) can be written as

$$\begin{eqnarray}&&\mathcal{H}=\mathop{\sum }\limits_{j=0}^{n}\mathcal{H}({{s}_{j}}),\end{eqnarray} \tag{ B.3 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{H}({{s}_{j}}) & = & \mathop{\sum }\limits_{{{s}_{j+1}}}^{t={{s}_{j}}}\left[ -\sqrt{\frac{1-\delta {{p}^{i}}(t)}{1-\delta {{p}^{i}}(t+\delta t)}}\hat{R}(\delta t)|t+\delta t\rangle \langle t| \right. \\ {} & {} & -\left. \sqrt{\frac{1-\delta {{p}^{i}}(t+\delta t)}{1-\delta {{p}^{i}}(t)}}{{{\hat{R}}}^{-1}}(\delta t)|t\rangle \langle t+\delta t| \right] \\ {} & {} & +|t\rangle \langle t|+|t+\delta t\rangle \langle t+\delta t|+1 \\ {} & {} & -\frac{\delta t}{\delta {{p}_{m}}}{{{\hat{C}}}_{m}}|\psi ({{s}_{j}}-\delta t)\rangle \langle \psi ({{s}_{j}}-\delta t)|\left. \hat{C}_{m}^{\dagger } \right)\otimes |{{s}_{j}}\rangle \langle {{s}_{j}}|. \\ \end{array}\end{eqnarray} \tag{ B.4 }$$

We see from the above expression, that each $\mathcal{H}({{s}_{j}})$ has the form of a free non-Hermitian clock Hamiltonian, starting from the initial state $\sqrt{\frac{\delta t}{\delta {{p}_{m}}}}{{\hat{C}}_{m}}|\psi ({{s}_{j}}-\delta t)\rangle$ instead of $|\psi (0)\rangle$. Furthermore, the various $\mathcal{H}({{s}_{j}})$ commute with one another. Therefore, $\mathcal{H}$ has a block diagonal structure, and the eigenvalue equation $\mathcal{H}|{{\eta }_{k}}\rangle ={{\lambda }_{k}}|{{\eta }_{k}}\rangle$ separates into separate eigenvalue equations in each block

$$\begin{eqnarray}&&\mathcal{H}({{s}_{j}})|{{\eta }_{k}}({{s}_{j}})\rangle ={{\lambda }_{k}}({{s}_{j}})|{{\eta }_{k}}({{s}_{j}})\rangle .\end{eqnarray} \tag{ B.5 }$$

Applying the results of (appendix A), each ${{\lambda }_{k}}({{s}_{j}})$ is real and positive. Therefore, ${{\lambda }_{k}}={{\sum }_{j}}{{\lambda }_{k}}({{s}_{j}})$ is real and positive as well. The state

$$\begin{eqnarray}\begin{array}{rcl} |\eta \rangle & = & \sqrt{\frac{\delta t}{T+\delta t}}\mathop{\sum }\limits_{{{s}_{0}}-\delta t}^{t=0}\frac{{{{\hat{R}}}^{\frac{t}{\delta t}}}}{\sqrt{1-p(t)}}|\psi (0)\rangle \otimes |t\rangle \\ {} & {} & +\mathop{\sum }\limits_{n}^{j=0}\mathop{\sum }\limits_{{{s}_{j+1}}-\delta t}^{t={{s}_{j}}}\frac{{{{\hat{R}}}^{\frac{t-{{s}_{j}}}{\delta t}}}}{\sqrt{1-p(t)}}{{{\hat{C}}}_{m}}|\psi ({{s}_{j}}-\delta t)\rangle \otimes |t\rangle \\ {} & {} & +\mathop{\sum }\limits_{T}^{t={{s}_{n+1}}}\frac{{{{\hat{R}}}^{\frac{t-{{s}_{n+1}}}{\delta t}}}}{\sqrt{1-p(t)}}{{{\hat{C}}}_{m}}|\psi ({{s}_{n+1}}-\delta t)\rangle \otimes |t\rangle , \\ \end{array}\end{eqnarray} \tag{ B.6 }$$

satisfies $\mathcal{H}|\eta \rangle =0$, so it is the ground state. From inspection, we see that this state has exactly the form $|{{\eta }^{i}}\rangle =\sqrt{\frac{\delta t}{T+\delta t}}\sum _{t=0}^{T}|{{\psi }^{i}}(t)\rangle \otimes |t\rangle$, with $|{{\psi }^{i}}(t)\rangle$ being the wave functions generated by an SSE propagation for the same realization of jumps. This proves the desired result, that the ground states of the stochastic Kitaev–Feynman clock Hamiltonians are history states of the stochastic trajectories.

Thus far, we have written the non-Hermitian evolution operator $\hat{R}$ appearing in equation (7) in a completely general form, which in principle can involve interactions involving any number of variables. In order for the stochastic clock Hamiltonian to be experimentally implemented on a quantum computer or simulator, it must be possible to construct $\hat{R}$ from terms that involve only a few qubits at a time. In this section, we write the stochastic clock Hamiltonian in an explicitly local form and provide examples.

Let the system consist of N variables represented by qubits and the dimension of the Hilbert space be ${{2}^{N}}$. Following Lloyd [26], we first decompose the system Hamiltonian as ${{\hat{H}}_{s}}=\sum _{i=1}^{l}{{\hat{H}}_{i}}$, where each ${{\hat{H}}_{i}}$ acts on a space of dimension ${{m}_{i}}\lt {{2}^{N}}$. Similarly, we decompose the bath operator as $\hat{D}=\sum _{i=1}^{p}{{\hat{D}}_{i}}$, where each ${{\hat{D}}_{i}}$ acts on a space of dimension ${{m}_{i}}\lt {{2}^{N}}$. Therefore, $\hat{R}=1-{\rm i}{{\hat{H}}_{s}}\delta t-\hat{D}\delta t$ can be written as $\hat{R}=1-{\rm i}\sum _{i=1}^{l}{{\hat{H}}_{i}}\delta t-\sum _{i=1}^{p}{{\hat{D}}_{i}}\delta t$, which is a sum of $l+p$ operators acting on a space of dimension ${{m}_{i}}\lt {{2}^{N}}$ and containing at most k of the N variables. Using this decomposition, the non-Hermitian clock Hamiltonian in equation (7) can be written as

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{H} & = & {{I}_{s}}\otimes {{\left[ \mathop{\sum }\limits_{t=0}^{T-\delta t}-|t+\delta t\rangle \langle t|-|t\rangle \langle t+\delta t|+|t\rangle \langle t|+|t+\delta t\rangle \langle t+\delta t| \right]}_{c}} \\ {} & {} & +\delta t{{\left[ i\mathop{\sum }\limits_{i=1}^{l}{{{\hat{H}}}_{i}}+\mathop{\sum }\limits_{i=1}^{p}{{{\hat{D}}}_{i}} \right]}_{s}}\otimes {{\left[ \mathop{\sum }\limits_{t=0}^{T-\delta t}|t+\delta t\rangle \langle t|-|t\rangle \langle t+\delta t| \right]}_{c}} \\ {} & {} & +(1-|\psi (0)\rangle {{\langle \psi (0)|)}_{s}}\otimes |0\rangle \langle 0{{|}_{c}}, \\ \end{array}\end{eqnarray} \tag{ C.1 }$$

where the subscripts s and c denote operators acting in the state space of the system and clock respectively. The first line in equation (C.1) contains terms that are at most two-local, acting on two clock qubits. The second line contains a sum of terms that are at most (k + 2)-local, acting on k system qubits and 2 clock qubits. As written, the locality of the operator in the third line depends on the initial state $|\psi (0)\rangle$. However, it is always possible to encode the preparation of the state $|\psi (0)\rangle$ as evolution from the state ${{[\Pi _{i=0}^{N}|0\rangle }_{i}}{{]}_{s}}\otimes |0{{\rangle }_{c}}$. This can be enforced with a term in the Hamiltonian, ${{\hat{H}}_{{\rm init}}}=[\sum _{i=1}^{N}|1\rangle {{\langle 1{{|}_{i}}]}_{s}}\otimes |0\rangle \langle 0{{|}_{c}}$, which is two-local acting on a single system-qubit together with the initial clock-qubit. Therefore, examining equation (C.1) we see that the non-Hermitian clock Hamiltonian acts in a Hilbert space of dimension ${{2}^{N+\frac{T}{\delta t}}}\times {{2}^{N+\frac{T}{\delta t}t}}$, but it can be written as a sum of terms, which in the worst case act in a Hilbert space of dimension ${{2}^{k+2}}\times {{2}^{k+2}}$, where in general ${{2}^{k+2}}\ll {{2}^{N+T/\delta t}}$ (see examples below). Furthermore, for an arbitrary clock-Hamiltonian $\mathcal{H}$ containing at worst (k + 2)-local terms, it is always possible to construct a new Hamiltonian with the same ground state that contains only $2+2=4$-local terms, using the method of perturbative gadgets [22, 27].

As an example, consider the Haken–Strobl model of pure dephasing, which describes the transfer of a localized excitation (exciton) through a system of coupled molecules such as chromophores in the presence of an environment [28] . The Haken–Strobl–Lindblad equation is given by equation (2), with ${{\hat{H}}_{s}}=\sum _{m=1}^{N}{{\epsilon }_{m}}|m\rangle \langle m|+\sum _{m=1}^{N}\sum _{m^{\prime} =1}^{N}{{V}_{m,m^{\prime} }}(|m\rangle \langle m^{\prime} |+|m^{\prime} \rangle \langle m|)$ and ${{\hat{C}}_{m}}=\sqrt{\Gamma }|m\rangle \langle m|$, where N is the number of molecules. Here, $|m\rangle$ describes an excitation localized on the mth molecule with energy ${{\epsilon }_{m}}$, ${{V}_{m,m^{\prime} }}$ describes energy transfer between molecules and ${{\hat{C}}_{m}}$ describes fast vibrations due to the environment, which introduce decoherence without energy exchange. We can easily encode this Hamiltonian with qubits, if $|1\rangle$ describes an excited molecule and $|0\rangle$ describes a molecule in its ground state. Then ${{\hat{H}}_{s}}=\sum _{m=1}^{N}\frac{{{\epsilon }_{m}}}{2}(1+\hat{\sigma }_{m}^{z})+\sum _{m=1}^{N}\sum _{m^{\prime} =1}^{N}{{V}_{m,m^{\prime} }}(\hat{\sigma }_{m}^{+}\hat{\sigma }_{m^{\prime} }^{-}+\hat{\sigma }_{m}^{-}\hat{\sigma }_{m^{\prime} }^{+})$ and ${{\hat{C}}_{m}}=\frac{\sqrt{\Gamma }}{2}(1+\hat{\sigma }_{m}^{z})$. If we further assume that the excitation begins localized on the jth molecule, $|\psi (0)\rangle =|j\rangle$ and equation (C.2) becomes

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{H} & = & {{I}_{s}}\otimes {{\left[ \mathop{\sum }\limits_{t=0}^{T-\delta t}-|t+\delta t\rangle \langle t|-|t\rangle \langle t+\delta t|+|t\rangle \langle t|+|t+\delta t\rangle \langle t+\delta t| \right]}_{c}} \\ {} & {} & +\delta t{{\left[ {\rm i}\mathop{\sum }\limits_{m=1}^{N}\frac{{{\epsilon }_{m}}-{\rm i}\Gamma }{2}(1+\hat{\sigma }_{m}^{z})+\mathop{\sum }\limits_{m=1}^{N}\mathop{\sum }\limits_{m^{\prime} =1}^{N}{{V}_{m,m^{\prime} }}(\hat{\sigma }_{m}^{+}\hat{\sigma }_{m^{\prime} }^{-}+\hat{\sigma }_{m}^{-}\hat{\sigma }_{m^{\prime} }^{+}) \right]}_{s}} \\ {} & {} & \otimes {{\left[ \mathop{\sum }\limits_{t=0}^{T-\delta t}|t+\delta t\rangle \langle t|-|t\rangle \langle t+\delta t| \right]}_{c}} \\ {} & {} & +{{(1-\hat{\sigma }_{j}^{z})}_{s}}\otimes |0\rangle \langle 0{{|}_{c}}. \\ \end{array}\end{eqnarray} \tag{ C.2 }$$

We see that this Hamiltonian can be engineered with at most four-local couplings, which can be further reduced to two-local couplings using gadgets [22, 27].

The above Hamiltonian is non-interacting. In the case of an interacting system, there will be 4-local terms within the system, which lead to 6-local terms when coupled to the clock. When the interacting system is fermionic, fermion operators must be encoded into qubits using the Jordan–Wigner transformation, which introduces additional non-locality [29]. Nevertheless, with gadgets even fermionic systems can be efficiently reduced to two-local interactions [30].