Parallel in time dynamics with quantum annealers

Recent years have witnessed an unprecedented increase in experiments and hybrid simulations involving quantum computers. In particular, quantum annealers. There exist a plethora of algorithms promising to outperform classical computers in the near-term future. Here, we propose a parallel in time approach to simulate dynamical systems designed to be executed already on present-day quantum annealers. In essence, purely classical methods for solving dynamics systems are serial. Therefore, their parallelization is substantially limited. In the presented approach, however, the time evolution is rephrased as a ground-state search of a classical Ising model. Such a problem is solved intrinsically in parallel by quantum computers. The main idea is exemplified by simulating the Rabi oscillations generated by a two-level quantum system (i.e. qubit) experimentally.

It is needless to say that simulating dynamical systems with near-term quantum technology poses one of the most difficult and technologically challenging endeavor 1 . Various computations of certain aspects of manybody quantum physics can already be assisted by the existing hardware [2][3][4] . For instance, recent experiments have demonstrated that quantum annealers 5 can be turned into neural networks that can learn the ground state energy of a physical system 6 . A similar task can also be accomplished with fewer qubits using quantum gates 7,8 .
The aforementioned examples characterize static processes where there is no real-time dynamics being simulated directly. Noticeably, near-term quantum annealers do simulate quantum annealing, which is a timedependent phenomenon. However, the optimization problem itself, i.e., the one to be solved by the annealer, exhibits no time dependence 9 . Thus, following the time evolution, even of a single qubit on a quantum annealer is a challenging task for the current technology. This should, nonetheless, be possible at least in principle. Indeed, a time-dependent quantum problem can be (re)formulated as a static one, defined on an appropriately enlarged Hilbert space 10 . This is realized using the Feynman's clock operator 11, 12 . This observation naturally encapsulates a family of powerful algorithms referred to as parallel in time or parareal methods, often invoked to simulate the system's dynamics on heterogeneous classical hardware 13,14 . The latter techniques effectively take advantage of the fact that a part of the evolution can be distributed and carried out in parallel. Nevertheless, with such an approach, one can never reach full parallelism on any classical hardware (of the Turing type) due to the communication bottlenecks 15 . Nonetheless, these limitations do not apply to the quantum hardware. Quite the contrary, quantum computers operate in parallel and any algorithm (cf. Refs. [16][17][18] ) they execute needs to be carefully designed from scratch to utilize their intrinsic parallelism fully.
As a proof of concept, in this article, we demonstrate how present-day quantum annealers may be programmed to simulate dynamical systems in parallel (due to their noisiness only in the specific regime of the problem's parameters). In particular, we determine the time evolution of a single qubit (Rabi oscillations) solely from experiments conducted on the newest D-Wave 2000Q quantum chip [19][20][21][22] . At the same time, due to the underlying connectivity (all-to-all) and the extensive amount of qubits it requires, the proposed algorithm constitutes a natural test which can determine the usefulness of various annealing technology realized by e.g. the Floquet annealer 23 , the large-scale (photonic 24 ) Ising machines [25][26][27][28] , and the Fujitsu digital annealer 29

parallel in time dynamics
Consider a dynamical system (e.g. a quantum system isolated from its environment 30 ) whose behavior can be described by a L dimensional and possibly time-dependent, Kamiltonian K(t) ("Kamiltonian" refers to a Hamiltonian-like operator, K 31 .). The system dynamics is encoded, at all times, in a (quantum) state, |ψ(t)� , whose evolution is governed by a Schrödinger like equation 32 , This first order differential equation admits a unique solution |ψ(t)� := U(t, t 0 )|ψ(t 0 )� , where propagates an arbitrary initial state, |ψ(t 0 )� , from t 0 to t ≥ t 0 whereas T denotes the time-ordering operator 33 .
Such an ordering can be omitted whenever [K(t), K(t ′ )] = 0 . In particular, for time independent systems, is a Hamiltonian, the evolution operator (2) is unitary and the dynamics (1) is norm preserving and reversible, i.e.
To solve Eq. (1), one usually discretizies the time interval [t 0 , t] selecting N distinct moments, i.e. t := t N−1 > · · · > t n+1 > t n > · · · > t 0 . The dynamics can then be formulated as a sequence of unitary gates, acting on an initial state. Note, each U n := U(t n+1 , t n ) can also be formally expressed using Eq. (2). Practically, however, for small time steps, all gates U n are approximated using variety of methods 32 . Those include exact diagonalization for small system 34 , Suzuki-Trotter decomposition 35 , commutator-free expansion 36 or sophisticated tensor networks techniques 37 .
The latter equation provides a starting point for various sequential numerical schemes for solving differential equations on classical computers 38 . In principle, however, those gates could also be realized on a quantum computer, which could then resolve the unitary dynamics efficiently 32 . Unfortunately, current quantum hardware does not allow for such gates to be constructed yet. Nevertheless, the underlying idea behind decomposition (3) can be harnessed to formulate an optimization problem that can be solved by present-day quantum annealers 4 . This is the main idea we put forward in this work.
Indeed, consider a superposition of quantum states in different moments of time t n , where the clock states are orthonormal, �t n |t m � = δ nm . With the corresponding clock operator, one obtains C | � = 0| � . Thus, | � is the ground state of C . Obviously, this state is not unique since we have specified neither initial nor boundary condition. However, introducing a penalty, say C 0 , allows one to provide additional constrains. In particular, specifying that C 0 = |t 0 ��t 0 | ⊗ (I − |ψ 0 ��ψ 0 |) , the following linear system encodes Eq. (1) subjected to |ψ(t 0 )� = |ψ 0 � . For hermitian systems, the above complex linear system of N × L equations expresses the reversible dynamics of the system in terms of a sequence of unitary gates (3). The hermitian clock operator can also be derived from e.g. time-embedded discrete variational principle (that is, the principle of least action) 10 . The idea can be further extended to open quantum systems 39 .
To solve the dynamics expressed in Eq. (6) on a quantum annealer one needs to formulate it as an optimization problem 12,40,41 . Moreover, such an optimization needs to be encoded via the Ising spin-glass Hamiltonian 5 (or QUBO 42 ) defined on a particular sparse graph called chimera 43 (or pegasus 44 ). We stress that these particular graphs are specific for the D-Wave hardware, and other architectures allow for a different connectivity between qubits. For instance, all-to-all (complete graph) in case of the (classical) Fujitsu Digital annealer 29 . Furthermore, at least complex fixed-point arithmetic is also required to express quantum states in consecutive moments of time 45 . Here, we incorporate a strategy introduced only recently in Ref. 46 , cf. also Ref. 45 for real matrices. To this end, we employ a natural correspondence between complex numbers and real 2 × 2 matrices, namely a + bi � → aÎ + ibσ y , to represent A using only real entries.
We further rely on a straightforward observation that the solution to Eq. (6), expanded in the standard basis as |x� = x i |i� , also minimizes the following functional h(x) = �A |x� − |��� 2 and vice versa. That is, a global minimum of h, i.e. x 0 is a solution (6) as h(x 0 ) = 0 . Moreover, when the simulated system is hermitian then A is positive definite. Therefore, x 0 is also a minimum of A |�� = |t 0 � ⊗ |ψ 0 �, A = C + |t 0 ��t 0 | ⊗ I, Since variables x i are real, the objective functions f (x) can not be programmed directly to be optimized on a quantum annealer. Nevertheless, one can obtain the so called fixed-point representation for each x i as a linear combination of binary variables q α i 45 The above correspondence is constructed with the assumption that R bits of binary representation are used for every real number in the solution vector. In our approach, the order of magnitude of the solution's coefficients is also assumed, i.e. x i ∈ [−2 D , 2 D ] for a fixed D ∈ N.
Therefore, the minimization problem to be solved on a quantum annealer can finally be formulated as where The constant energy contribution, f 0 , can be omitted as both f (q) and f (q) − f 0 have the same optimal solution q 0 . Since f (q 0 ) = 0 , one can easily asses the quality of the solution found by any heuristic approach. For small N, QUBO (9) is defined on a complete graph [cf. Fig. 1b] with |V | = R × N × (2L) vertices. In contrast, when N ≫ 2 the number of edges is equal to the number of nonzero elements of A which is sparse. Currently, the biggest complete graph that can be embedded on the 2000Q chip has |V | = 65 vertices ( |V | = 180 for the Pegasus topology 47 ), cf. Fig. 1. It is worth mentioning that classical solvers (hardware-based or otherwise) usually offer better connectivity and thus can realize much denser graphs without the need for embedding. For example, the so-called coherent Ising machines (among others) can incorporate complete graphs consisting of the order of 10 3 vertices 26 . Therefore, QUBO generated from the dynamics provide a natural "stress" test for those machines which can asses their usefulness in simulating physics.

Quantum annealing
Adiabatic quantum computing can be seen as an alternative paradigm of computation 5 . Essentially, it is equivalent to the gate model of quantum computation that uses logical gates operating on quantum states to implement quantum algorithms 32 . The main idea is based on the quantum adiabatic theorem 48 . When a system starting www.nature.com/scientificreports/ from its ground state is driven slowly enough, it has time to adjust to any change, and thus it can remain in the ground state during the entire evolution. Assume a quantum system is prepared in the ground state of an initial ("simple") Hamiltonian H 0 . Then, it will slowly evolve to the ground state of the final ("complex") Hamiltonian H p that one can harness to encode the solution to an optimization problem. In particular, the dynamics of the current D-Wave 2000Q quantum annealer is supposed to be governed by the following time-dependent Hamiltonian (cf. Ref. 9 ) where the problem Hamiltonian H p realizes the spin-glass Ising model defined on the chimera graph, (E , V ) , specified by its edges and vertices, The annealing time τ varies from microseconds to milliseconds depending on the programmable schedule 9 . Typically, during the evolution g(s) varies from g(0) ≫ 0 [i.e. all spins point in the . Note, the Hamiltonian H p is classical in a sense that all its terms commute. Thus, its eigenstates translate directly to classical optimization variables, q α i , which we introduced to encode the time evolution (6) as QUBO (9). The Pauli operators σ z i , σ x i describe the spin degrees of freedom in the z-and x-direction respectively.
Dimensionless real couplers, J ij ∈ [−1, 1] , and magnetic fields, h i ∈ [−2, 2] , are programmable. In practice, the actual values of those parameters that are sent to the quantum processing unit differ from the ones specified by the user by a small amount δJ ij , δh i 49 . This is due to various reasons including noise effects which we will neglect in this work (cf. Ref. 4,50 ).
Most practical optimization problems are defined on dense graphs which can be embedded onto the chimera graph 51 . There is, however, a substantial overhead that effectively limits the size of problems that can be solved with current quantum annealers. This is, nonetheless, an engineering issue that will most likely be overcome in the near future 23,47 .

Results
To exemplify the main idea we consider a two-level quantum system (qubit) whose Hamiltonian reads where σ y is the Pauli spin matrix in the y-direction. For the sake of simplicity, we further set ω = π/2 . Moreover, due to the limited number of qubits and sparse connectivity of D-Wave quantum annealers, we mainly consider the system's evolution at six distinct integer time points, starting from |ψ 0 � = |0� . This ensures that the dynamics can be captured precisely with two bits of precision per component of the state vector, thus allowing one to run experiments on the D-Wave 2000Q annealer. For the illustrative purposes we reconstruct �σ z �(t).
As depicted in Fig. 1, the low noise D-Wave 2000Q annealer was able to capture the dynamics faithfully [cf. Fig. 1d,f], for τ = 200 µs, 2000 µs. Therein, probability distributions, ρ , were constructed from 10 4 anneals. There were no post-processing involved and the Boltzmann temperature, β , was set to its default value. Furthermore, to construct the dynamics (Fig. 1d,f) only the lowest energy state, reported by the device, was utilized. For this particular problem excited states (also returned by the annealer) are not, a priori useful. This experiments demonstrate an improvement in comparison to the (not that) older generation, results for which are shown in Fig. 1c and e.
In contrast, results obtained from an emulation of the D-Wave output with tensor networks (cf. Ref. 52 ) are presented in Fig. 2. As a reference point, we have also included solutions found by the CPLEX optimizer 53 . Both these solvers, being purely classical, exhibit superior performance in comparison to the D-Wave quantum annealers (Assuming sufficiently large precision of all A ij .). This is noticeable especially for problems that require bigger graphs resulting from higher precision-(R ≥ 3 , N = 6 ), cf. Fig. 2a-or extra time points ( N > 6 , R = 2 ), cf. Fig. 2b. Similar degradation of the solution quality with the increasing problem size has been observed, e.g., in Ref. 54,55 in the context of problems requiring complete graphs, cf. Fig. 1. The behavior, as mentioned above, is expected from an early stage device which is prone to errors. Their origins, however, are anything but straightforward to pinpoint precisely. In stark contrast, there is yet another source of errors that is related to the precision of J ij , and h i 56 . Those errors are believed to be predominant for the type of simulations introduced in this work. Indeed, Fig. 2c-h shows the destructive (above all not monotonic) effect of the limited precision-r, of the problem coefficients A ij -on the solution. Beyond a certain threshold, neither the D-Wave annealer nor the aforementioned classical heuristics can reproduce the dynamics (i.e., oscillations) accurately.
As a final note, we stress that the precision R (the only one that influences the graph size) is the number of bits required for representing the discretized continuous variables, x [cf. Eq. (8)], which is different from r. The latter is an inherent characteristic of the hardware which stems from the DAC quantization step (for more details cf. D-Wave's technical notes 56 ).

conclusions
In this article, we have proposed a parallel in time approach to simulate dynamical systems with the quantum annealing technology. While one should not expect current quantum annealers to be faster/better than classical computers, our results constitute, first and foremost, a proof of concept demonstrating how the technology can www.nature.com/scientificreports/ be employed to simulate the time evolution of simple (e.g. two-level) quantum systems (in certain regime of the corresponding parameters). This task is a priori difficult for the current prototypical quantum computers which are prone to errors and has been designed mostly to simulate static phenomena. Furthermore, not only the Ising instances we have generated can be executed on the commercially available D-Wave annealers, but they can also be tested on: coherent Ising machines [24][25][26][27][28] , the Floquet annealer 23 , and the Fujitsu digital annealer 29 that celebrate all-to-all connectivity. This provides a practical "stress" test for those machines which can determine their usefulness in simulating various time-dependent properties of physical systems.