Floquet quantum simulation with superconducting qubits

TL;DR

This paper introduces a Floquet-based quantum simulation method using superconducting qubits to emulate spin models, offering a resource-efficient approach for programmable quantum simulators.

Contribution

It develops a Floquet theory-based algorithm for simulating XYZ spin models with superconducting qubits, including a practical recipe for implementing the transverse Ising model.

Findings

Numerical demonstration of simulating the transverse Ising model.

Benchmarking shows advantages over digital simulation with limited resources.

Protocol is feasible with current superconducting chip technology.

Abstract

We propose a quantum algorithm for simulating spin models based on periodic modulation of transmon qubits. Using Floquet theory we derive an effective time-averaged Hamiltonian, which is of the general XYZ class, different from the isotropic XY Hamiltonian typically realised by the physical setup. As an example, we provide a simple recipe to construct a transverse Ising Hamiltonian in the Floquet basis. For a 1D system we demonstrate numerically the dynamical simulation of the transverse Ising Hamiltonian and quantum annealing to its ground state. We benchmark the Floquet approach with a digital simulation procedure, and demonstrate that it is advantageous for limited resources and finite anharmonicity of the transmons. The described protocol can serve as a simple yet reliable path towards configurable quantum simulators with currently existing superconducting chips.

Equations

$${\hat{\mathcal{H}}}_{\mathrm{F}}=\frac{1}{T}{\int }_0^T{\hat{\mathcal{H}}}_0^{\prime }(t)dt,$$

$${\hat{\mathcal{H}}}_0=\sum _{j=1}^{N-1}J({\sigma }_j^x{\sigma }_{j+1}^x+{\sigma }_j^y{\sigma }_{j+1}^y).$$

$${\hat{\mathcal{H}}}_1(t)=\sum _{j=1}^{\lceil N/2\rceil }{\mathbf{h}}_{\mathrm{odd}}(t)\cdot {\mathbf{\sigma }}_{2j-1}+\sum _{j=1}^{\lfloor N/2\rfloor }{\mathbf{h}}_{\mathrm{even}}(t)\cdot {\mathbf{\sigma }}_{2j},$$

$${\hat{\mathcal{H}}}_{\mathrm{F}}=\sum _{\alpha ,{\alpha }^{\prime }=x,y,z}\sum _{j=1}^{\lfloor N/2\rfloor }\overline{{\xi }_{\alpha {\alpha }^{\prime }}}{\sigma }_{2j}^{\alpha }({\sigma }_{2j-1}^{{\alpha }^{\prime }}+{\sigma }_{2j+1}^{{\alpha }^{\prime }}),$$

$${\hat{\mathcal{H}}}_F=J_{\mathrm{sim}}\sum _{j=1}^{N-1}{\sigma }_j^x{\sigma }_{j+1}^x+h^z\sum _{j=1}^N{\sigma }_j^z\equiv {\hat{\mathcal{H}}}_{\mathrm{tIsing}},$$

$$\hat{\mathcal{H}}(t)={\hat{\mathcal{H}}}_0+{\hat{\mathcal{H}}}_1(t)=J\sum _{j=1}^{N-1}({\sigma }_j^x{\sigma }_{j+1}^x+{\sigma }_j^y{\sigma }_{j+1}^y)+\sum _{j=1}^{\lceil N/2\rceil }{\mathbf{h}}_{\mathrm{odd}}(t)\cdot {\mathbf{\sigma }}_{2j-1}+\sum _{j=1}^{\lfloor N/2\rfloor }{\mathbf{h}}_{\mathrm{even}}(t)\cdot {\mathbf{\sigma }}_{2j},$$

$$\hat{\mathcal{U}}(t)=\mathcal{T}\exp \left[-i\underset{t_0}{\overset{t}{\int }}d{t}^{\prime }{f}_{\mathrm{o}}(t^{\prime })\sum _{j=1}^{\lceil N/2\rceil }(h_{\mathrm{o}}^x{\sigma }_{2j-1}^x+h_{\mathrm{o}}^y{\sigma }_{2j-1}^y+h_{\mathrm{o}}^z{\sigma }_{2j-1}^z)+f_{\mathrm{e}}(t^{\prime })\sum _{j=1}^{\lfloor N/2\rfloor }(h_{\mathrm{e}}^x{\sigma }_{2j}^x+h_{\mathrm{e}}^y{\sigma }_{2j}^y+h_{\mathrm{e}}^z{\sigma }_{2j}^z)\right],$$

$${\hat{\mathcal{U}}}_{\mathrm{p}}(t)=\exp \left[-i\frac{g_{\mathrm{p}}(t)}{2}\sum _{j=1}^{N_{\mathrm{p}}}(n_{\mathrm{p}}^x{\sigma }_P^x+n_{\mathrm{p}}^y{\sigma }_P^y+n_{\mathrm{p}}^z{\sigma }_P^z)\right],$$

$${\hat{\mathcal{U}}}^{†}(t)\left[J\sum _{j=1}^{N-1}({\sigma }_j^x{\sigma }_{j+1}^x+{\sigma }_j^y{\sigma }_{j+1}^y)\right]\hat{\mathcal{U}}(t)=e^{\left[i\frac{g_{\mathrm{e}}(t)}{2}\sum _{j=1}^{N_{\mathrm{e}}}(n_{\mathrm{e}}^x{\sigma }_{2j}^x+n_{\mathrm{e}}^y{\sigma }_{2j}^y+n_{\mathrm{e}}^z{\sigma }_{2j}^z)\right]}{e}^{\left[i\frac{g_{\mathrm{o}}(t)}{2}\sum _{j=1}^{N_{\mathrm{o}}}(n_{\mathrm{o}}^x{\sigma }_{2j-1}^x+n_{\mathrm{o}}^y{\sigma }_{2j-1}^y+n_{\mathrm{o}}^z{\sigma }_{2j-1}^z)\right]}\times$$

$$\times \left[J\sum _{j=1}^{N-1}({\sigma }_j^x{\sigma }_{j+1}^x+{\sigma }_j^y{\sigma }_{j+1}^y)\right]e^{\left[-i\frac{g_{\mathrm{o}}(t)}{2}\sum _{j=1}^{N_{\mathrm{o}}}(n_{\mathrm{o}}^x{\sigma }_{2j-1}^x+n_{\mathrm{o}}^y{\sigma }_{2j-1}^y+n_{\mathrm{o}}^z{\sigma }_{2j-1}^z)\right]}{e}^{\left[-i\frac{g_{\mathrm{e}}(t)}{2}\sum _{j=1}^{N_{\mathrm{e}}}(n_{\mathrm{e}}^x{\sigma }_{2j}^x+n_{\mathrm{e}}^y{\sigma }_{2j}^y+n_{\mathrm{e}}^z{\sigma }_{2j}^z)\right]}.$$

$${\hat{\mathcal{R}}}_j^{\mathbf{n}}(g)=e^{-i\frac{g}{2}\mathbf{n}\cdot {\sigma }_j}={\hat{\mathcal{R}}}_j^z(\varphi ){\hat{\mathcal{R}}}_j^y(\theta ){\hat{\mathcal{R}}}_j^z(g){\hat{\mathcal{R}}}_j^y(\theta {)}^{†}{\hat{\mathcal{R}}}_j^z(\varphi {)}^{†},$$

$$\hat{\mathcal{U}}^{\dagger}(t)\hat{\mathcal{H}}_{0}\hat{\mathcal{U}}(t)=\hat{\mathcal{U}}_{\mathrm{e}}^{\dagger}\big{[}\hat{\mathcal{R}}^{z}_{\mathrm{o}}(\phi_{\mathrm{o}})\hat{\mathcal{R}}^{y}_{\mathrm{o}}(\theta_{\mathrm{o}})\hat{\mathcal{R}}^{z}_{\mathrm{o}}(-g_{\mathrm{o}})\hat{\mathcal{R}}^{y}_{\mathrm{o}}(-\theta_{\mathrm{o}})\hat{\mathcal{R}}^{z}_{\mathrm{o}}(-\phi_{\mathrm{o}})\hat{\mathcal{H}}_{0}\hat{\mathcal{R}}^{z}_{\mathrm{o}}(\phi_{\mathrm{o}})\hat{\mathcal{R}}^{y}_{\mathrm{o}}(\theta_{\mathrm{o}})\hat{\mathcal{R}}^{z}_{\mathrm{o}}(g_{\mathrm{o}})\hat{\mathcal{R}}^{y}_{\mathrm{o}}(-\theta_{\mathrm{o}})\hat{\mathcal{R}}^{z}_{\mathrm{o}}(-\phi_{\mathrm{o}})\big{]}\hat{\mathcal{U}}_{\mathrm{e}},$$

$$e^{\hat{\mathcal{M}}}{\hat{\mathcal{H}}}_0{e}^{-\hat{\mathcal{M}}}={\hat{\mathcal{H}}}_0+[\hat{\mathcal{M}},{\hat{\mathcal{H}}}_0]+\frac{1}{2!}[\hat{\mathcal{M}},[\hat{\mathcal{M}},{\hat{\mathcal{H}}}_0]]+...=\sum _{k=0}^{\infty }\frac{1}{k!}[\hat{\mathcal{M}},{\hat{\mathcal{H}}}_0{]}_k,$$

$${\hat{\mathcal{H}}}^{\prime }(t)=\sum _{\alpha ,{\alpha }^{\prime }=x,y,z}\sum _{j=1}^{N/2}{\xi }_{\alpha {\alpha }^{\prime }}{\sigma }_{2j}^{\alpha }({\sigma }_{2j-1}^{{\alpha }^{\prime }}+{\sigma }_{2j+1}^{{\alpha }^{\prime }}),$$

$$\displaystyle\xi_{xx}(t)=\sin^{2}(\theta_{e})\cos(\phi_{e})\Big{[}\cos(g_{o}[t])\cos^{2}(\theta_{o})\cos(\phi_{o})\cos(\phi_{e}-\phi_{o})-\sin(\phi_{e})\sin(g_{o}[t])\cos(\theta_{o})-\cos(g_{o}[t])\sin(\phi_{o})\sin(\phi_{e}-\phi_{o})$$

$$\displaystyle+\sin^{2}(\theta_{o})\cos(\phi_{o})\cos(\phi_{e}-\phi_{o})\Big{]}+\cos(g_{e}[t])\Big{\{}\cos(g_{o}[t])\Big{[}\cos^{2}(\theta_{o})\cos(\phi_{o})\big{(}\cos^{2}(\theta_{e})\cos(\phi_{e})\cos(\phi_{e}-\phi_{o})+\sin(\phi_{e})\times$$

$$\times$$

$$\displaystyle+\cos^{2}(\theta_{e})\cos(\phi_{e})\sin^{2}(\theta_{o})\cos(\phi_{o})\cos(\phi_{e}-\phi_{o})+\sin(\phi_{e})\sin^{2}(\theta_{o})\cos(\phi_{o})\sin(\phi_{e}-\phi_{o})\Big{\}}+\sin(g_{e}[t])\cos(\theta_{e})\times$$

$$\displaystyle\times\Big{[}\sin(g_{o}[t])\cos(\theta_{o})-\sin^{2}\left(\frac{g_{o}[t]}{2}\right)\sin^{2}(\theta_{o})\sin(2\phi_{o})\Big{]},$$

$$\displaystyle\xi_{yy}(t)=\sin^{2}(\theta_{e})\sin(\phi_{e})\Big{[}\sin(\phi_{o})\cos(\phi_{e}-\phi_{o})\big{(}\cos(g_{o}[t])\cos^{2}(\theta_{o})+\sin^{2}(\theta_{o})\big{)}+\cos(\phi_{e})\sin(g_{o}[t])\cos(\theta_{o})$$

$$\displaystyle+\cos(g_{o}[t])\cos(\phi_{o})\sin(\phi_{e}-\phi_{o})\Big{]}+\cos(g_{e}[t])\Big{\{}\cos(g_{o}[t])\Big{[}\cos^{2}(\phi_{e})\big{(}\cos^{2}(\theta_{o})\sin^{2}(\phi_{o})+\cos^{2}(\phi_{o})\big{)}+\cos^{2}(\theta_{e})\sin(\phi_{e})\times$$

$$\displaystyle\times\big{(}\cos^{2}(\theta_{o})\sin(\phi_{o})\cos(\phi_{e}-\phi_{o})+\cos(\phi_{o})\sin(\phi_{e}-\phi_{o})\big{)}+\sin(\phi_{e})\cos(\phi_{e})\sin^{2}(\theta_{o})\sin(\phi_{o})\cos(\phi_{o})\Big{]}-\sin^{2}(\theta_{e})\sin(\phi_{e})\times$$

$$\displaystyle\times\cos(\phi_{e})\sin(g_{o}[t])\cos(\theta_{o})+\sin^{2}(\theta_{o})\sin(\phi_{o})\big{(}\cos^{2}(\theta_{e})\sin(\phi_{e})\cos(\phi_{e}-\phi_{o})-\cos(\phi_{e})\sin(\phi_{e}-\phi_{o})\big{)}\Big{\}}+\sin(g_{e}[t])\cos(\theta_{e})\times$$

$$\displaystyle\times\Big{[}\sin(g_{o}[t])\cos(\theta_{o})+\sin^{2}\left(\frac{g_{o}[t]}{2}\right)\sin^{2}(\theta_{o})\sin(2\phi_{o})\Big{]},$$

$$\displaystyle\xi_{zz}(t)=\sin(\theta_{e})\sin(\theta_{o})\bigg{\{}\cos(\phi_{e}-\phi_{o})\left[\sin(g_{e}[t])\sin(g_{o}[t])+2\big{(}1-\cos(g_{e}[t])\big{)}\cos(\theta_{e})\sin^{2}\left(\frac{g_{o}[t]}{2}\right)\cos(\theta_{o})\right]$$

$$\displaystyle+\sin(\phi_{e}-\phi_{o})\left[\big{(}1-\cos(g_{e}[t])\big{)}\cos(\theta_{e})\sin(g_{o}[t])-2\sin(g_{e}[t])\sin^{2}\left(\frac{g_{o}[t]}{2}\right)\cos(\theta_{o})\right]\bigg{\}},$$

$$\displaystyle\xi_{xy}(t)=\sin^{2}(\theta_{o})\sin(\phi_{o})\Big{[}\cos(g_{e}[t])\Big{(}\cos^{2}(\theta_{e})\cos(\phi_{e})\cos(\phi_{e}-\phi_{o})+\sin(\phi_{e})\sin(\phi_{e}-\phi_{o})\Big{)}+\sin^{2}(\theta_{e})\cos^{2}(\phi_{e})\cos(\phi_{o})\Big{]}$$

$$\displaystyle+\sin^{2}(\phi_{o})\Big{[}\sin^{2}(\theta_{e})\cos^{2}(\phi_{e})\sin(g_{o}[t])\cos(\theta_{o})+\frac{1}{2}\sin^{2}(\theta_{o})\Big{(}\sin^{2}(\theta_{e})\sin(2\phi_{e})-2\sin(g_{e}[t])\cos(\theta_{e})\Big{)}\Big{]}+\sin(g_{o}[t])\cos(\theta_{o})\times$$

$$\displaystyle\times\Big{[}\cos^{2}(\phi_{e})\Big{(}\cos(g_{e}[t])\cos^{2}(\theta_{e})+\sin^{2}(\theta_{e})\cos^{2}(\phi_{o})\Big{)}+\cos(g_{e}[t])\sin^{2}(\phi_{e})\Big{]}+\frac{1}{4}\cos(g_{o}[t])\Bigg{[}4\cos(g_{e}[t])\cos^{2}(\theta_{e})\cos(\phi_{e})\cos(\phi_{o})\times$$

$$\displaystyle\times\Big{(}\sin(\phi_{e})\cos(\phi_{o})-\cos(\phi_{e})\sin^{2}(\theta_{o})\sin(\phi_{o})\Big{)}+4\sin^{2}\left(\frac{g_{e}[t]}{2}\right)\sin^{2}(\theta_{e})\sin(2\phi_{e})\cos^{2}(\theta_{o})\sin^{2}(\phi_{o})+4\sin(\phi_{e})\cos(\phi_{e})\times$$

$$\displaystyle\times\cos^{2}(\phi_{o})\Big{(}\sin^{2}(\theta_{e})-\cos(g_{e}[t])\Big{)}-4\sin^{2}(\theta_{o})\sin(\phi_{o})\cos(\phi_{o})\Big{(}\cos(g_{e}[t])\sin^{2}(\phi_{e})+\sin^{2}(\theta_{e})\cos^{2}(\phi_{e})\Big{)}-\sin(g_{e}[t])\cos(\theta_{e})\times$$

$$\displaystyle\times\Big{(}2\sin^{2}(\theta_{o})\cos(2\phi_{o})+\cos(2\theta_{o})+3\Big{)}\Bigg{]},$$

$$\displaystyle\xi_{yx}(t)=\frac{1}{8}\Bigg{\{}8\sin(g_{e}[t])\cos(\theta_{e})\sin^{2}(\theta_{o})\cos^{2}(\phi_{o})+8\cos(g_{e}[t])\sin^{2}(\theta_{o})\sin(\phi_{o})\cos(\phi_{o})\Big{(}\cos^{2}(\theta_{e})\sin^{2}(\phi_{e})+\cos^{2}(\phi_{e})\Big{)}$$

$$\displaystyle+2\cos(g_{o}[t])\bigg{[}4\cos^{2}(\theta_{o})\cos(\phi_{o})\Big{(}\sin^{2}(\theta_{e})\sin(\phi_{e})\cos(\phi_{e}-\phi_{o})-\cos(g_{e}[t])\cos(\phi_{e})\sin(\phi_{e}-\phi_{o})\Big{)}+4\cos(g_{e}[t])\cos^{2}(\theta_{e})\times$$

$$\displaystyle\times\sin(\phi_{e})\Big{(}\cos^{2}(\theta_{o})\cos(\phi_{o})\cos(\phi_{e}-\phi_{o})-\sin(\phi_{o})\sin(\phi_{e}-\phi_{o})\Big{)}+\sin(g_{e}[t])\cos(\theta_{e})\Big{(}-2\sin^{2}(\theta_{o})\cos(2\phi_{o})+\cos(2\theta_{o})+3\Big{)}$$

$$\displaystyle-4\sin(\phi_{o})\Big{(}\cos(g_{e}[t])\cos(\phi_{e})\cos(\phi_{e}-\phi_{o})+\sin^{2}(\theta_{e})\sin(\phi_{e})\sin(\phi_{e}-\phi_{o})\Big{)}\bigg{]}+8\sin^{2}(\theta_{e})\sin(\phi_{e})\sin^{2}(\theta_{o})\cos(\phi_{o})\times$$

$$\displaystyle\times\Big{(}\cos(\phi_{e}-\phi_{o})-\cos(g_{e}[t])\cos(\phi_{e})\cos(\phi_{o})\Big{)}-8\sin^{2}(\theta_{e})\sin^{2}(\phi_{e})\sin(g_{o}[t])\cos(\theta_{o})-2\cos(g_{e}[t])\sin(g_{o}[t])\cos(\theta_{o})\times$$

$$\displaystyle\times\Big{(}2\sin^{2}(\theta_{e})\cos(2\phi_{e})+\cos(2\theta_{e})+3\Big{)}\Bigg{\}},$$

$$\displaystyle\xi_{yz}(t)=\sin(\theta_{o})\Bigg{\{}\cos(g_{e}[t])\cos^{2}(\theta_{e})\sin(\phi_{e})\Big{(}2\sin^{2}\left(\frac{g_{o}[t]}{2}\right)\cos(\theta_{o})\cos(\phi_{e}-\phi_{o})+\sin(g_{o}[t])\sin(\phi_{e}-\phi_{o})\Big{)}$$

$$\displaystyle+\sin(g_{e}[t])\cos(\theta_{e})\Big{(}2\sin^{2}\left(\frac{g_{o}[t]}{2}\right)\cos(\theta_{o})\cos(\phi_{o})-\sin(g_{o}[t])\sin(\phi_{o})\Big{)}+\cos(\phi_{e}-\phi_{o})\Big{(}2\sin^{2}(\theta_{e})\sin(\phi_{e})\sin^{2}\left(\frac{g_{o}[t]}{2}\right)\cos(\theta_{o})$$

$$\displaystyle+\cos(g_{e}[t])\cos(\phi_{e})\sin(g_{o}[t])\Big{)}+\sin(\phi_{e}-\phi_{o})\Big{(}\sin^{2}(\theta_{e})\sin(\phi_{e})\sin(g_{o}[t])-2\cos(g_{e}[t])\cos(\phi_{e})\sin^{2}\left(\frac{g_{o}[t]}{2}\right)\cos(\theta_{o})\Big{)}\Bigg{\}},$$

Full text

Quantum simulation generally relies on exploiting a controllable quantum system to imitate a complex quantum system of interest Lloyd1996 . It provides an efficient way to solve classically inaccessible problems of material science Wecker2015 and quantum chemistry Wecker2014 ; Reiher2016 . Ultimately quantum simulation may give access to properties of complex quantum systems with exponential speed-up as compared to classical algorithms NielsenChuang ; Georgescu2014 . Superconducting circuits (SC) Wallraff2004 ; Devoret2013 have recently emerged as one of the prime candidates for realizing full scale quantum computers, with operations close to the fault tolerant threshold Barends2014 ; Kelly2015 ; Riste2015 ; Corcoles2015 ; Ofek2016 ; Reed2012 and the possibility of scaling these advances to larger systems Bejanin2016 . From the simulation point of view, various protocols were proposed and realized experimentally Houck2012 , including simulation of spin systems LasHeras2014 ; Salathe2015 ; Barends2016 , fermionic models Barends2015 , quantum chemistry OMalley2015 , thermalization Neill2016 , synthetic magnetic fields Roushan2016 , ultrastrong coupling Ballester2012 ; Langford2016 , and gauge field theories Mezzacapo2015 . To demonstrate the full potential of quantum simulation with SC, however, there is a need for protocols which can outperform classical protocols for realistic near term medium size systems.

Typically protocols for quantum simulation can be divided into digital and analog (or emulator) types. While these techniques are similar, they exploit different methods to achieve a quantum speed-up. The digital approach relies on discretizing the Hamiltonian evolution using a set of quantum gates. A protocol for simulating an arbitrary unitary $\hat{U}(t)=\exp(-i\hat{\mathcal{H}}t)$ governed by a Hamiltonian $\hat{\mathcal{H}}$ not available in a physical setup, exploits the sequential implementation of the available unitaries $\hat{U}_{\ell}$ represented by gates. The corresponding unitary of a single digital step $j$ of duration $\delta t$ can be constructed as $\hat{U}_{j}(\delta t)=\prod_{\ell}\hat{U}_{\ell}$. This string of unitaries can be recast in terms of Hamiltonians $\hat{U}_{j}(\delta t)=\prod_{m}e^{-i\hat{\mathcal{H}}_{m}\delta t}$. An implementation of $N_{\mathrm{Tr}}\rightarrow\infty$ of these Trotter steps combines into a unitary $\hat{U}(t)=\lim_{N_{\mathrm{Tr}}\rightarrow\infty}\hat{U}_{j}(\delta t)^{N_{\mathrm{Tr}}}\approx e^{-i\hat{\mathcal{H}}t}$, where $\hat{H}=\sum_{m}\hat{\mathcal{H}}_{m}$ and $t=\delta tN_{\mathrm{Tr}}$ NielsenChuang . Given a universal set of gates any required Hamiltonian can in principle be simulated. This poses the challenge of implementing several quantum gates for successful simulation, leading to errors if they have limited fidelity. This issue can possibly be overcome by quantum error correction but this requires substantial overhead in resources. The digital approach is widely used for quantum simulation with SC LasHeras2014 ; Salathe2015 ; Barends2016 ; Barends2015 ; OMalley2015 ; Ballester2012 ; Langford2016 , as it is tunable and does not require changing sample layout to simulate different models.

Analog quantum simulation relies on the actual physical implementation of the required Hamiltonian, corresponding to the emulation of a targeted real system. This was realised on various platforms, including cold atoms in optical lattices Bloch2012 ; Simon2011 and trapped ions Kim2010 ; Cenko2015 . By exploiting continuous time dynamics, Trotterization errors are minimized and analog protocols can therefore have superior simulation fidelities compared to digital approaches. However, they are restricted to the types of Hamiltonian physically realizable in the system. In particular, this is the case for SC quantum systems. Current highly coherent chains of transmons, sketched in Fig. 1, are limited to nearest neighbour flip-flop type of interqubit interaction provided by the capacitive coupling between them. Thus, they are confined to simulate isotropic XY type spin-1/2 model Wendin2016 . The accessible models can be enriched by implementing different connections between the qubits to engineer various nonlinear couplings Neumeier2013 ; Dumur2015 ; Sameti2016 , or allowing for modulation of interqubit interactions to break the rotating wave approximation Bertet2006 ; Mezzacapo2014 ; Kapit2015 . However, this adds extra complications to the experiments and potentially introduces additional errors.

Here, we propose an alternative analog-like simulation strategy, which can be performed without modifications of the system. It is based on using a Floquet basis to perform quantum simulation with superconducting circuits, and can be extended to ground state preparation via quantum annealing. The idea relies on the time-dependent modulation of the Hamiltonian $\hat{\mathcal{H}}(t)=\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{1}(t)$, where $\hat{\mathcal{H}}_{0}$ is a time-independent part, and $\hat{\mathcal{H}}_{1}(t)=\hat{\mathcal{H}}_{1}(t+T)$ is periodic with a period $T=2\pi/\omega$. The corresponding unitary operator for a single period can be rewritten as an effective evolution with a time-independent Hamiltonian given by the Magnus expansion Bukov2015 ; Iadecola2015 . When the frequency of modulation $\omega$ is much bigger than the coupling in the static Hamiltonian, $\omega/\|\hat{\mathcal{H}}_{0}\|\gg 1$, the dynamics of the system can be conveniently represented in terms of a period averaged Floquet Hamiltonian,

[TABLE]

where $\hat{\mathcal{H}}^{\prime}_{0}(t)$ denotes the static Hamiltonian $\hat{\mathcal{H}}_{0}$ rewritten in the interaction picture with respect to the oscillating part. The resulting Hamiltonian $\hat{\mathcal{H}}_{F}$ may possess qualitatively different behavior compared to $\hat{\mathcal{H}}_{0}$. Such Floquet type simulation recently gained attention in the cold atom Goldman2014 and condensed matter physics communities Rudner2013 ; Jiang2011 , where Floquet topological insulators and gauge fields were introduced. Finally, in circuit QED the Floquet quasienergies of a single qubit were studied recently Deng2015 .

In this paper we consider superconducting qubits with isotropic XY coupling, and show that by exploiting fast driving of each site we can efficiently tailor the effective Hamiltonian of the system. For concreteness we here focus on one dimensional chains, but the method can easily be extended to more dimensions. First, we describe how the approach can be used to simulate the dynamics of the transverse Ising model showing that as compared to digital protocols higher fidelity can be attained. Next, we simulate quantum annealing to the ground state of the transverse Ising model and find that the Floquet approach outperforms the digital for restricted resources, when limited by the finite anharmonicity of the transmons. Finally, we describe an algorithm for simulating the spin-1/2 XYZ model with the Floquet approach.

*System and Hamiltonian.—*As a physical realization we consider a capacitively coupled chain of qubits (Fig. 1), where periodically oscillating effective magnetic fields are applied at chosen lattice sites. The qubits can be of transmon Riste2015 , xmon Barends2013 , gmon Chen2014 , or gatemon-type Casparis2016 , with the main requirements being tunability and high fidelity operation. The time independent Hamiltonian in the rotating frame $\hat{\mathcal{H}}_{0}$ contains a nearest-neighbour flip-flop interaction with bare coupling $J$, corresponding to the isotropic XY spin-1/2 model

[TABLE]

Here $\sigma^{\alpha}_{j}$ ($\alpha=x,y,z$) are spin-1/2 Pauli operators at lattice site $j$, $N$ denotes total number of qubits in the chain, and we consider open boundary conditions. The time dependent Hamiltonian $\hat{\mathcal{H}}_{1}(t)$ contains a periodic magnetic field $\mathbf{h}_{j}(t)$ which rapidly oscillates along arbitrary axes, and we assume that it differs between even and odd sites,

[TABLE]

where $\lfloor x\rfloor$ and $\lceil x\rceil$ denote floor and ceiling functions, respectively. We assume that the magnetic field is sharply turned on at time $t=0$, and thus explicitly account for the kick operator contribution, which is relevant in the current setting targeting realistic quantum simulation Goldman2014 . Going to the rotating frame with respect to $\hat{\mathcal{H}}_{1}$ and integrating over a period as in Eq. (1) we get the reduced Floquet Hamiltonian in the general form [SM, A]:

[TABLE]

where the time-averaged coefficients $\overline{\xi_{\alpha\alpha^{\prime}}}$ are defined in the Supplemental Material SM , and are controlled by the amplitude and alignment of the effective magnetic fields.

*Transverse Ising model: quantum dynamics.—*We first consider the Floquet simulation of the tranverse Ising model, which represents a particular case of the more general Floquet Hamiltonian (4). It can be realized with a drive of the form $\mathbf{h}_{\mathrm{even}}(t)=\lambda\omega\cos(\omega t)\mathbf{e}^{x}+2h^{z}(1+\mathcal{J}_{0}[4\lambda])^{-1}\cos(2\lambda\sin[\omega t])\mathbf{e}^{z}$ and $\mathbf{h}_{\mathrm{odd}}(t)=h^{z}\mathbf{e}^{z}$, where $\lambda$ is a drive parameter, and $\mathcal{J}_{0}[x]$ denotes the zeroth order Bessel function of the first kind. The $z$-directed terms additionally introduce an effective transverse magnetic field $h^{z}$. In the infinite frequency limit $|J|/\omega\rightarrow 0$ and for $\lambda=1.20241$ ($\mathcal{J}_{0}[2\lambda]=0$) this leads to an effective Hamiltonian of the form [SM, B]

[TABLE]

where for later reference $J_{\mathrm{sim}}$ describes the effective simulated coupling of the model. In the transverse Ising case it corresponds to the original bare coupling, $J_{\mathrm{sim}}=J$, and we note that the ratio between the effective magnetic field in the Floquet basis and the Ising interaction, $h^{z}/J_{\mathrm{sim}}$, can be controlled by the drive parameters. To verify the approach, we numerically calculate the full dynamics of the system with time periodic driving of frequency $\omega/|J|=50$ and access the Floquet dynamics by looking at stroboscopic times $nT$, where $n$ is an integer [Fig. 2 (a)]. The results are compared to an ideal simulation of the transverse Ising Hamiltonian (5), with the initial state $|\psi_{\mathrm{in}}\rangle=\bigotimes_{j=1}^{N/2}(|\uparrow\rangle_{2j-1}-i|\downarrow\rangle_{2j-1})/\sqrt{2}\otimes|\uparrow\rangle_{2j}$. Additionally, we benchmark the Floquet simulation with a digital protocol [SM, C]. It uses the isotropic XY interaction and its rotated version, and approaches the Ising model in the limit of large $N_{\mathrm{Tr}}$ LasHeras2014 . The results in Fig. 2(b) show that the Floquet simulation closely follows the exact dynamics at short time, but deviates at later stages due to the finite Floquet frequency. On the contrary, the digital approach shows substantial deviations for this limited number of Trotter periods even at short times, but will have rapid convergence with more Trotter steps (see below). We note that a comparison of the number of Floques periods and Trotter steps may not be a fair comparison, since the latter involves multiple gates. A more detailed comparison is done below.

*Transverse Ising model: quantum annealing.—*Next, we study the ground state preparation of the simulated model. To access a ground state, we perform the quantum annealing procedure Kadowaki1998 , which also serves as a basis for adiabatic quantum computing Fahri2000 , and may solve NP-complete problems Fahri2001 . Using the Floquet basis, we design the Hamiltonian $\hat{\mathcal{H}}_{F}(t)=J_{\mathrm{sim}}\sum_{j=1}^{N-1}\sigma_{j}^{x}\sigma_{j+1}^{x}+(1-t/t_{f})h^{z}\sum_{j=1}^{N}\sigma_{j}^{z}$, where the effective magnetic field is linearly turned off during a total annealing time $t_{f}$. Here we consider $J_{\mathrm{sim}}<0$, which allows to achieve the ground state of the ferromagnetic $x$-Ising Hamiltonian. The ideal target state is an entangled GHZ-like state $|\psi_{\mathrm{T}}\rangle=(|+\rangle^{\otimes N}+|-\rangle^{\otimes N})/\sqrt{2}$, where $|\pm\rangle={(|\downarrow\rangle\pm|\uparrow\rangle)}/\sqrt{2}$, and we start from the trivial initial state $|\psi_{\mathrm{in}}\rangle=|\downarrow\rangle^{\otimes N}$.

The results of the annealing procedure are shown in Fig. 3. The dynamics of the system, quantified as the fidelity of the instantaneous wavefunction of the system with the ideal target state, $F=|\langle\psi_{\mathrm{T}}|\psi(t)\rangle|^{2}$, is shown in Fig. 3(a) for a 4 qubits chain. Blue dots correspond to Floquet simulation at stroboscopic times, which closely follows the red solid curve of the ideal continuous annealer. The blue oscillatory curve corresponds to the full dynamics, with fast oscillation arising from the drive term. To study the scaling with the system size, we perform fixed time ($t_{f}=15|J_{\mathrm{sim}}|^{-1}$) annealing for chains of various length [Fig. 3(b)]. Fixing the drive frequency to moderate values, we observe that the final infidelities of the Floquet simulator, $1-F_{\mathrm{Floquet}}(t_{f})$, and digital simulator, $1-F_{\mathrm{digital}}(t_{f})$, have similar scalings with the system size, both adding extra infidelity on top of the continuous evolution, and largely dependent on $\omega$ and $N_{\mathrm{Tr}}$ as described below. In Fig. 3(c) we show the dependence on the frequency of the periodic drive of the Floquet infidelity, measured with respect to the finite time annealing state. Here the frequency is rescaled by the total annealing time, such that $\omega t_{f}/2\pi$ shows the number of stroboscopic periods. The infidelity can be reduced for large $\omega$, with results converging to continuous annealing infidelity for $\omega/|J|\rightarrow\infty$. The analogous behavior for the digital approach corresponds to the variation of the number of Trotter steps, and is shown in Fig. 3(d). While a direct comparison between two approaches is complicated, the general tendency can be deduced: the Floquet approach has smaller infidelity for small number of steps and limited resources, while the digital approach has better scaling if a large number of Trotter steps $N_{\mathrm{Tr}}$ can be implemented.

*Imperfections.—*To describe a realistic quantum simulator, we study the influence of a finite anharmonicity $A$ of the transmon qubits, which will be a major limitation to our approach. Driving a transmon with a finite anharmonicity leads to leakage of information from the logical subspace. To account for this we consider a full Hamiltonian of a SC chain [SM, D], and perform numerical simulations including doubly occupied states of the transmons. As an example we use annealing of the transverse Ising model with $N=4$. The resulting infidelity of the simulation is shown in Fig. 4. First, we fix the value of the anharmonicity, and calculate the infidelity as a function of the Floquet frequency [Fig. 4(a)]. We observe that contrary to the ideal circuit, the infidelity is minimized for a finite (optimal) drive frequency $\omega^{\mathrm{opt}}$. We note that the window of frequencies in which the infidelity stays close to minimal is typically broad. In Fig. 4(b) we show the optimized infidelity of the simulation as a function of $A$ (blue curve).

To benchmark the results of the Floquet simulation we compare it to the digital simulation. Assuming the same value $\epsilon$ for the single- and two-qubit gate errors, the error for a single Trotter step is equal to $\epsilon_{\mathrm{Tr}}=(5N-4)\epsilon$. We optimize the total infidelity $1-F_{\mathrm{tot}}=1-(1-\epsilon_{\mathrm{Tr}})^{N_{\mathrm{Tr}}}(1-\epsilon_{\mathrm{dig}})$ with respect to the number of Trotter steps. Here $\epsilon_{\mathrm{dig}}$ is the error coming from digitization of the dynamics [as shown in Fig. 3(d)]. To highlight the relevant parameter range, we assign a cut-off based on the single gate time $t_{\mathrm{gate}}\geq 35A^{-1}$. The results are plotted as horizontal lines in Fig. 4(b), and show that the Floquet approach can outperform the digital approach unless very high fidelity gates with $\epsilon<10^{-4}$ are used. Furthermore, the Floquet approach is highly advantageous for small values of $A/J_{\mathrm{sim}}$, which for a given transmon anharmonicity $A$ is the regime where the simulation is finished the fastest and thus has the least influence of decoherence. Thus for shorter time of simulation/higher error rates the Floquet is advantageous compared to the digital approach. We highlight that while digital approaches typically exploit DRAG techniques to remove leakage Motzoi2009 ; Martinis2014 , the presented Floquet approach is not specifically designed to work for small $A$, and it may be possible to improve on this issue using few-tone drives.

*Estimates.—*To quantify the performance we consider numbers which can be achieved with currently available transmon setups Barends2013 ; Chen2014 . Taking the anharmonicity to be $A=2\pi\times 300$ MHz, drive frequency $\omega=2\pi\times 9.8$ MHz, nearest-neighbour coupling $J=2\pi\times 1$ MHz (reduced compared to most setups), and $t_{f}=2.4~{}\mu$s, the four qubit chain can be annealed to the ground state of the Ising model with $1-F_{\mathrm{Flq}}=0.037$ (ideal continuous annealing gives $1-F_{\mathrm{cont}}=0.00616$). An additional error will arise from dephasing, but its effect will be small for highly coherent qubits ($T_{2}>10~{}\mu$s).

Reaching a similar performance with the digital strategy is highly challenging and would require single and two qubit gate operation times of $18$ ns $=33A^{-1}$ and $\epsilon=10^{-4}$ accuracy. The single Trotter step duration for $N=4$ is $t_{\mathrm{Tr}}=0.162~{}\mu$s, and with $t_{f}=2.4~{}\mu$s this can allow for 14 Trotter steps. The corresponding digitization error is $1-F_{\mathrm{dig}}=0.041$ (dephasing should be added separately).

*Generic XYZ Hamiltonians.—*The Floquet approach may be extended to simulate generic spin-1/2 models represented by XYZ type spin Hamiltonian, $\hat{\mathcal{H}}_{\mathrm{XYZ}}=\sum_{j=1}^{N-1}(J^{x}\sigma_{j}^{x}\sigma_{j+1}^{x}+J^{y}\sigma_{j}^{y}\sigma_{j+1}^{y}+J^{z}\sigma_{j}^{z}\sigma_{j+1}^{z})$. These are of the so called non-stoquatic type, where recent results have suggested that they can give enhanced computational powers Nishimori2016 . We consider couplings $J_{y}=2J^{x}/3$, $J^{z}=J^{x}/3$, and $J^{x}=J<0$. This configuration can be realized by a uniform periodic magnetic field in the $x$ direction at all lattice sites [SM, E]. As compared to the transverse Ising case this Hamiltonian possesses small energy gaps, and in the absence of the additional transverse field is difficult to anneal even with the ideal continuous Hamiltonian. The results are qualitatively similar to the Ising case, showing that the Floquet approach outperforms digital simulation unless ultrahigh fidelity gates are implemented, although high-quality simulation is very challenging in both cases.

*Conclusion.—*We have presented a scheme for a reconfigurable and tunable superconducting quantum simulator based on the transmon qubits. Utilizing the Floquet approach, we show that the originally limited isotropic XY interaction can be transformed into transverse Ising or XYZ type spin-1/2 Hamiltonian. The approach was used for simulation of multiqubit system dynamics and preparation of non-trivial groundstates. The Floquet simulation was shown to perform better than a digital scheme for restricted resources, and represents a realistic path for modern SC quantum simulators.

*Acknowledgements.—*The research was funded by the European Union Seventh Framework Programme through ERC Grant QIOS (Grant No. 306576). O.K. acknowledges useful discussions with Pedram Roushan, Charles Neill, Guanyu Zhu, Mohammad Hafezi, Göran Johansson, and Stefan Filipp.

.1 Derivation of the generic Floquet Hamiltonian

Here we present the derivation of Eq. (4) in the main text, and simultaneously describe the general Floquet Hamiltonian originating from arbitrary axis magnetic field oscillation for the two sublattices.

We start with the transmon Hamiltonian written in the form

[TABLE]

and consider the effective magnetic fields to be different for even and odd sublattices. Here $\lfloor x\rfloor$ and $\lceil x\rceil$ denote floor and ceiling functions, correspondingly. To have an explicit form of $\mathbf{h}_{\mathrm{even/odd}}(t)$ we decompose it into $\mathbf{h}_{\mathrm{even/odd}}(t)=f_{\mathrm{e/o}}(t)(\mathbf{e}^{x}h^{x}_{\mathrm{e/o}}+\mathbf{e}^{y}h^{y}_{\mathrm{e/o}}+\mathbf{e}^{z}h^{z}_{\mathrm{e/o}})$, where $f_{\mathrm{e/o}}(t)=f_{\mathrm{e/o}}(t+T)$ is some periodic function, $\mathbf{e}^{x,y,z}$ form a Cartesian basis, and $h^{\alpha}_{\mathrm{e/o}}$ are components of an effective magnetic field for even (indexed by e) and odd (indexed by o) sites. This corresponds to the magnetic field oscillating along a certain axis, which is different for even and odd sublattices. Here we consider the time-dependence to be the same for all spin components, such that the axis of magnetic field does not precess. This restriction decreases the number of independent tuning parameters. Note that the effective magnetic field used for the simulation of the transverse Ising model, deviates slightly from this form. This deviation is, however, a perturbation and will be dealt with below. If the fields do not have this form, the time-dependent Hamiltonian $\hat{\mathcal{H}}_{1}(t)$ does not commute with itself at different times, thus largely complicating the solution. This more general setting will be considered in future works.

We perform a unitary transformation with respect to the rapidly oscillating time-dependent part $\hat{\mathcal{H}}_{1}(t)$. This is done using the unitary operator:

[TABLE]

where $\hat{\mathcal{T}}\{...\}$ is the time-ordering operator and $t_{0}$ is an initial switch-on time. The unitary (S2) can be factorized into odd and even sublattices parts, and time-ordering disappears as the magnetic field oscillates along a fixed axis. This yields $\hat{\mathcal{U}}(t)=\hat{\mathcal{U}}_{\mathrm{o}}(t)\hat{\mathcal{U}}_{\mathrm{e}}(t)$ with

[TABLE]

where indices $\mathrm{p}=\mathrm{o,e}$ and $\mathrm{P}=\mathrm{2j-1,2j}$ help denote sublattice parity, and we define the time-integrated functions $g_{\mathrm{p}}(t)\equiv 2\int\limits_{t_{0}}^{t}dt^{\prime}f_{\mathrm{p}}(t^{\prime})h_{\mathrm{p}}$. Here $h_{\mathrm{p}}=\sqrt{(h^{x}_{\mathrm{p}})^{2}+(h^{y}_{\mathrm{p}})^{2}+(h^{z}_{\mathrm{p}})^{2}}$ defines the absolute value of the effective magnetic field vector, and $n^{x,y,z}_{\mathrm{p}}=h^{x,y,z}_{\mathrm{p}}/h_{\mathrm{p}}$ correspond to the normalized Cartesian components. The rotated Hamiltonian then reads $\hat{\mathcal{H}}^{\prime}(t)=\hat{\mathcal{U}}^{\dagger}(t)\hat{\mathcal{H}}(t)\hat{\mathcal{U}}(t)-i\hat{\mathcal{U}}^{\dagger}(t)\partial_{t}\hat{\mathcal{U}}(t)$. The second term is divided into two parts and each of these can be rewritten using the relation for the derivative of an arbitrary time-dependent matrix $A(t)$, being $e^{-A(t)}\partial_{t}e^{A(t)}=\dot{A}(t)-[A(t),\dot{A}(t)]/2!+[A(t),[A(t),\dot{A}(t)]]/3!-...$. For the case of $A(t)=-i\frac{g_{\mathrm{p}}(t)}{2}\sum\limits_{j=1}^{N_{\mathrm{p}}}(n^{x}_{\mathrm{p}}\sigma^{x}_{P}+n^{y}_{\mathrm{p}}\sigma^{y}_{P}+n^{z}_{\mathrm{p}}\sigma^{z}_{P})$ considered here the commutators vanish, and in total the derivative term gives $-\hat{\mathcal{H}}_{1}(t)$. However, we emphasize that this conclusion would not be true for general time-dependence of the Cartesian components of an effective magnetic field, where additional derivative-dependent terms appear.

Next, we need to calculate the matrix product terms of the form

[TABLE]

Let us perform the unitary rotation for each sublattice consecutively. For this, it is convenient to go from Cartesian to the spherical coordinate frame $\mathbf{n}_{\mathrm{p}}=\{n^{x}_{\mathrm{p}},n^{y}_{\mathrm{p}},n^{z}_{\mathrm{p}}\}\leftrightarrow\{\cos\phi_{\mathrm{p}}\sin\theta_{\mathrm{p}},\sin\phi_{\mathrm{p}}\sin\theta_{\mathrm{p}},\cos\theta_{\mathrm{p}}\}$. Next, noting that the unitary operation is defined as a spin rotation with respect to a fixed axis $\mathbf{n}_{\mathrm{p}}$, it can be decomposed into Cartesian rotations as:

[TABLE]

where the rotation operator is defined as $\hat{\mathcal{R}}^{\alpha}(\varphi)=e^{-i\frac{\varphi}{2}\sigma^{\alpha}}$ ($\alpha=x,y,z$). Then, the Hamiltonian after the odd sublattice transformation can be obtained through rotations

[TABLE]

and the subsequent transformation for the even sublattice can be performed in a similar fashion.

Finally, to get a closed expression for the transformed Hamiltonian we use the Baker-Campbell-Hausdorff formula

[TABLE]

where $[\hat{\mathcal{M}},\hat{\mathcal{H}}_{0}]_{k}$ denotes $k$-th order nested commutator. We proceed with calculating the commutators and resumming the series. After straightforward but tedious algebra we can get the Hamiltonian in a rotating frame

[TABLE]

with the coefficients

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we state explicitly the time-dependence of even and odd integral terms $g_{e}[t]$ and $g_{o}[t]$.

Once we have rotated the Hamiltonian into a suitable frame, the corresponding unitary operator for the evolution during a period $T$ can be rewritten using the Magnus expansion Bukov2015SM :

[TABLE]

with the Floquet Hamiltonian $\hat{\mathcal{H}}_{\mathrm{F}}=\hat{\mathcal{H}}^{(0)}_{\mathrm{F}}+\hat{\mathcal{H}}^{(1)}_{\mathrm{F}}+\hat{\mathcal{H}}^{(2)}_{\mathrm{F}}+...$, which consists of corrections in $\mathcal{O}[1/\omega]$. They can be written as

[TABLE]

and higher-order terms can be written in a similar way using nested commutators. We notice that $||\hat{\mathcal{H}}_{\mathrm{F}}^{(k)}||\sim(1/\omega)^{k}$, and thus for very small time intervals, given by the period of the oscillating term $T=2\pi/\omega\rightarrow 0$ (infinite frequency limit), the effective Floquet Hamiltonian is represented by the period-average of the time-dependent Hamiltonian written in Eq. (S10). In this work we consider the Floquet Hailtonian only to lowest order, $\hat{\mathcal{H}}_{\mathrm{F}}=\hat{\mathcal{H}}_{\mathrm{F}}^{(0)}$, while higher order corrections $\sim(1/\omega)^{k}$ ($k>0$) are accounted for in the numerical integration of the full time-dependent Hamiltonian.

Finally, let us choose the form of the oscillatoric magnetic field and find the period-averaged Floquet Hamiltonian for the transmon chain. This can be chosen in the form

[TABLE]

where $\lambda_{\mathrm{e/o}}$ are constants of order unity, and $\varphi_{\mathrm{e/o}}$ are initial phases for the modulation. The latter is of high importance in the Floquet formalism, as it leads to kick-operator terms which change the basis of the system, but do not enter the effective time-independent Hamiltonian Goldman2014_SM . However, in the current study we are interested in actual protocols with Floquet simulation, where the drive term is abruptly turned on at time point $t_{0}=0$, and the initial phase of the drive may be important. Here we consider zero initial phases $\varphi_{\mathrm{e/o}}=0$, such that the kick operator is unity.

The integral functions $g_{\mathrm{e/o}}[t]$ are given by

[TABLE]

Then, the period averaged coefficients $\overline{\xi_{\alpha\alpha^{\prime}}}=(1/2\pi)\int_{0}^{2\pi}d\tau\xi_{\alpha\alpha^{\prime}}(\tau)$ can be written in the form:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The above equations define the exact form of the generic Hamiltonian (4) from the main text, and thus describe the possible Hamiltonians accessible for the Floquet quantum simulation with this method, assuming different even/odd periodic cosine modulation.

.2 Transverse Ising model derivation and intuitive explanation

In this section we provide a detailed procedure to engineer the transverse Ising Hamiltonian as an effective Floquet Hamiltonian of the isotropic XY model with transverse and longitudinal fields. For this we consider a system with two (odd and even) sublattices, where only one of the sublattices experience fast oscillations of the magnetic field (see sketch in Fig. S1).

The intuitive way to describe the realization of the Ising model ($\propto\sigma_{1}^{x}\sigma_{2}^{x}$) from the isotropic XY case (with interaction type $\propto\sigma_{1}^{x}\sigma_{2}^{x}+\sigma_{1}^{y}\sigma_{2}^{y}$) is to take two spins and imagine one to be rotating in a magnetic field. Starting from the flip-flop interaction, if we choose the axis of the magnetic field to be in the $\mathbf{e}^{x}$ direction, nothing will happen to the first interaction term. At the same time, the rotation of the second spin leads to oscillations of the second term between $\pm\sigma_{1}^{y}\sigma_{2}^{y}$, and also induce a $\propto\sigma_{1}^{y}\sigma_{2}^{z}$ interaction component. For large frequencies of oscillation and carefully chosen drive amplitude, the plus and minus components will cancel each other, as well as the cross-interaction components, ultimately leaving the Ising term $\propto\sigma_{1}^{x}\sigma_{2}^{x}$ only.

Having in mind the aforementioned intuitive explanation, and taking the full solution from section A, the starting Hamiltonian for the simulation of the Ising model in the transmon chain reads:

[TABLE]

where $\lambda(\equiv\lambda_{\mathrm{even}})$ is a drive parameter, and $\mathcal{J}_{0}[x]$ denotes the zeroth order Bessel function of the first kind. In the infinite frequency limit $|J|/\omega\rightarrow 0$ and for $\lambda=1.20241$ (such that $\mathcal{J}_{0}[2\lambda]=0$) the interaction term can be reduced to the Ising type. Additionally, the last term in Eq. (S15) is designed to introduce a transverse effective magnetic field and can be written as:

[TABLE]

The first term in Eq. (S16) is a static magnetic field on the odd sublattice, and commutes trivially with the fast oscillation part. However, the magnetic field on the even sublattice can be modified by the drive. The second term representing the magnetic field, deviates from the general form considered in Sec. .1, since it does not have the same time dependence. As opposed to the other field, however, the magnitude of this field does not increase with increasing $\omega$ and can thus be treated as a perturbation. Since it doesn’t commute with the main driving field, it will be strongly modified by it. Going to the rotating frame with the unitary operator $\hat{\mathcal{U}}_{R}(t)=\exp\left\{i\lambda\sin(\omega t)\right\}$, the magnetic term becomes

[TABLE]

Finally, performing the period averaging, the $\sigma^{y}_{2j}$ term vanishes trivially since the sine function oscillates between positive and negative values. At the same time, given that $\lambda$ is fixed by the condition $\mathcal{J}_{0}[2\lambda]=0$, the integral $\int_{0}^{2\pi}dx\cos^{2}(2\lambda\sin[x])=\pi(1+\mathcal{J}_{0}[4\lambda])$ gives a finite result

[TABLE]

thus allowing to introduce an effective transverse field $h^{z}$.

.3 Digital simulation of transverse Ising model

In this section we describe the digital simulation protocol, which we use to benchmark the performance of the Floquet quantum simulator. The transverse Ising model simulation, which we will consider, was theoretically described in Ref. LasHeras2014_SM and experimentally realized in Ref. Salathe2015_SM . The circuit scheme is shown in Fig. S2.

The algorithm relies on the preparation of a unitary transformation with the effective Hamiltonian of interest using repetitions of a small step, corresponding to the Trotterization procedure. The protocol for the simulation of an arbitrary k-local Hamiltonian $\hat{\mathcal{H}}$ (generally not available in the physical setup) relies on the sequential implementation of the available parts of a Hamiltonian $\hat{\mathcal{H}}_{k}$ (constructed from gates) such that $\sum_{k}\hat{\mathcal{H}}_{k}=\hat{\mathcal{H}}$. The corresponding unitary of a single digital step $j$ of duration $\delta t$ reads

[TABLE]

and the implementation of $N_{\mathrm{Tr}}\rightarrow\infty$ Trotter steps combine into the unitary $\hat{U}(t)=\lim_{N_{\mathrm{Tr}}\rightarrow\infty}\hat{U}_{j}(\delta t)^{N_{\mathrm{Tr}}}\approx e^{-i\hat{\mathcal{H}}t}$.

In this spirit, the implementation of the transverse Ising model model was proposed to rely on multiple applications of the Trotter step graphically shown in Fig. S2. It starts with the implementation of $U_{\mathrm{XY}}=\exp\{-i\delta t\sum\limits_{j=1}^{N/2}\hat{\mathcal{H}}^{(+)}_{2j-1,2j}\}$, where by $\hat{\mathcal{H}}^{(+)}_{j,j^{\prime}}=\frac{J}{2}(\sigma^{x}_{j}\sigma^{x}_{j^{\prime}}+\sigma^{y}_{j}\sigma^{y}_{j^{\prime}})$ we define the simple application of XY interaction for each pair of qubits. Next, this unitary can be rotated by applying $\pi$ X rotations at every second qubit, $\hat{R}_{x}=\exp\{-i\frac{\pi}{2}\sum_{j=1}^{N/2}\sigma^{x}_{2j}\}$, which leads to $\hat{R}_{x}^{\dagger}U_{\mathrm{XY}}\hat{R}_{x}=\exp\{-i\delta t\sum\limits_{j=1}^{N/2}\hat{\mathcal{H}}^{(-)}_{2j-1,2j}\}$, where we define $\hat{\mathcal{H}}^{(-)}_{j,j^{\prime}}=\frac{J}{2}(\sigma^{x}_{j}\sigma^{x}_{j^{\prime}}-\sigma^{y}_{j}\sigma^{y}_{j^{\prime}})$ as the XY Hamiltonian with the YY term flipped by rotation. Finally, the last layer in the Trotter step implements the transverse fields with $U_{\mathrm{Z}}=\exp\{-i\delta th^{z}\sum\limits_{j=1}^{N}\sigma_{j}^{z}\}$. Once the Trotter step is repeated many times, the non-commuting Hamiltonian parts can be added, thus implementing the transverse Ising model digitally.

The same considerations can be repeated for the digitized annealing procedure Barends2016_SM . Here the important part is to keep the phase applied by $U_{\mathrm{Z}}$ gates consistent with the adiabatic evolution.

.4 Accounting for the finite anharmonicity

In this section we consider the transmon chain Hamiltonian and account for finite anharmonicity. This can be described by the Hamiltonian

[TABLE]

where $\hat{a}_{j}^{\dagger}$ ($\hat{a}_{j}$) corresponds to the creation (annihilation) operator for excitations of the $j$-th transmon circuit. The first term in Eq. (S20) corresponds to nearest-neighbour capacitive coupling for transmons. The second term in curly brackets denotes an effective magnetic field in the $z$ direction given by the flux-bias dependent detuning $\Delta_{j}$ and the microwave drive terms $\Omega_{j}$ corresponding to an effective magnetic field in $xy$ plane. The last term of Eq. (S20) corresponds to the anharmonicity of the circuit $A$ provided by Josephson junctions. In the case of infinitely large anharmonicity the Hamiltonian (S20) can be projected onto the lowest occupation subspace for each qubit $\{|0\rangle,|1\rangle\}_{j}$, accounting only for singly excited circuits. This allows for a spin-1/2 description of the chain, and subsequent simulation of quantum magnetism. However, in realistic transmon samples the anharmonicity is typically small, and higher states of the circuit must be accounted for [see sketch in Fig. S3(a)]. In particular, this is important for the case of a strong microwave drive $\Omega$, as it leads to non-zero occupation of higher lying states, corresponding to the leakage of information out of the logical subspace. This can largely decrease the fidelity, and typically is the bottleneck for fast digital computation.

In this study we consider the effects of finite $A$ by expanding the Hilbert space for each site to have doubly occupied states, $\{|0\rangle,|1\rangle,|2\rangle\}_{j}$. As a test case we take the ground state preparation of the transverse Ising model, studied for the case of infinite anharmonicity in the preceding sections. The annealing schedule is chosen in the same form, making use of the $h^{x}\leftrightarrow\Omega(\Omega^{*})$ and $h^{z}\leftrightarrow 2\Delta$ correspondence.

The main results of the finite $A$ scaling for the Floquet quantum simulation is presented in Fig. 4 and the corresponding section of the main text. Here we provide more details of the Floquet to digital benchmarking procedure. The estimate of the digital protocol infidelity accounts for several contributions. The first comes from the Trotterization procedure $\epsilon_{\mathrm{dig}}(N_{\mathrm{Tr}})$ and depends strongly on the number of Trotter steps, favouring long sequences. The second contribution is a total infidelity from gate operations $\epsilon_{\mathrm{gates}}(N_{\mathrm{Tr}})=1-\left[1-(5N-4)\epsilon\right]^{N_{\mathrm{Tr}}}$, which increases with the number of Trotter steps. The optimization procedure is performed for different values of the gate errors $\epsilon$. The results are plotted in Fig. S3(b), and the optimal Trotter step numbers $N_{\mathrm{Tr}}^{\mathrm{opt}}$ is shown to decrease in the case of large gate errors. The corresponding optimal infidelity [Fig. S3(c)] shows a significant increase for $\epsilon>10^{-3}$.

To compare the digital and Floquet approaches, we should compare how each of the two approaches could be implemented on comparable physical systems. To this end we consider transmons with the same anharmonicity $A$, and assume that they also have comparable decoherence rates (but note that we assume that the physical coupling $J$ can be different in the two scenarios). Therefore, to have a similar influence of decoherence in the two approaches we assume that the simulations need to be completed in the same time. To simulate the same evolution this requires that the two approaches have the same $J_{\mathrm{sim}}$ and thus the same values of $A/J_{\mathrm{sim}}$. In the Floquet case this is defined by the $A/\omega$ and $\omega/J$ ratios. For the digital simulation, the relation is more subtle and relies on the scaling of the gate time with $A$ for a fixed error value. The full discussion of this complex subject lies beyond the scope of present study, and for simplicity we just assume that each gate can be implemented in a time $t_{\mathrm{gate}}=c/A$, where $c$ is a constant which controls the quality of the gate. Taking the existing studies Rol2016_SM ; Motzoi2009_SM ; Martinis2014_SM , and considering a best case scenario, we set $c=35$. For realistic devices with $A=2\pi\times 300$ MHz this will correspond to 18 ns gates, and we will consider a low error of $\epsilon=10^{-4}$. This will be extremely challenging to achieve with current technologies, but for the example considered in the main text this is what is required to achieve a performance comparable to what is achievable with the Floquet approach. With more realistic numbers the performance of the digital approach will be less ideal and thus not comparable to the result which can be achieved with the Floquet simulation.

.5 Floquet simulation of spin-1/2 XYZ model

To simulate the XYZ model in the Floquet basis, we start with the time-dependent Hamiltonian in the form

[TABLE]

where we consider the oscillating effective magnetic field to be homogeneous for all sites. $\hat{\mathcal{H}}_{z}$ describes the part of the Hamiltonian responsible for implementing the static $z$-oriented magnetic field for annealing. Considering $\mathbf{h}(t)=\lambda\omega\cos(\omega t)\mathbf{e}^{x}$ with $\omega$ being the largest energy scale, and setting $\lambda=1.81144$, we eliminate the cross-terms and are left with a Floquet Hamiltonian of the form

[TABLE]

where the couplings are $J^{y}=2J^{x}/3$, $J^{z}=J^{x}/3$, and $J^{x}=J<0$. Here the simulated coupling changes for the YY and ZZ interaction components, and we consider $J^{x}=J_{\mathrm{sim}}$ as a reference. We note that as compared to the transverse Ising case this Hamiltonian possesses a small energy gap, and in the absence of an additional transverse field it is difficult to anneal even with the ideal continuous Hamiltonian.

To characterize the Floquet and digital simulation procedures we plot the instantaneous infidelity with respect to the continuous annealing case, observing how closely one can follow the ground state [Fig. S4(a)]. The blue curve for the Floquet simulation at stroboscopic times shows that the deviation begins to grow once we approach the critical point. To compare with the digital procedure, we plot the infidelity for the Trotterization approach. We assume the number of Trotter steps is equal to half of number of Floquet periods, $N_{\mathrm{Tr}}=1/2(t_{f}/T)=477$ (see below for the details of the circuit). This will be an upper bound for the number of Trotter steps for the digital simulator for the same resources. This can be deduced from the digital simulation protocol assuming that the time required to implement the two-qubit gate is inversely proportional to the coupling strength, $t_{\mathrm{gate}}\sim|J|^{-1}$, and in addition the gates on the two sublattices have to be applied separately. We note that considering different ordering of the gates results in largely different results for the final state infidelity. While we have not performed a full optimization of this ordering, the results presented in the figure is the result of the optimization over 24 different possibilities for a Trotter step composition. This suggests that the digital procedure is strongly model dependent, and extra resources are required for sequence optimization Wecker2015_SM ; Reiher2016_SM . In Fig. S4(b) we plot the dependence of the Floquet XYZ annealing on the drive frequency, showing the resources necessary for high fidelity annealing as a function of the number of stroboscopic periods. A similar analysis is performed for the digital procedure, where the dependence of Trotter step number is considered [Fig. S4(c)]. Akin to the transverse Ising annealing case described in the main text, the Floquet approach shows smaller infidelity for limited number of steps, in particular for $N=4$. Unlike the transverse Ising model, however even for $N_{\mathrm{Tr}}$ as large as thousand steps, the digital approach fails provide smaller infidelity as compared to Floquet approach, but this may change if even higher number of Trotter steps are considered.

*Digital XYZ model.—*The considered digital simulation protocol, originally described in Ref. LasHeras2014_SM , is sketched in Fig. S5.

It relies on the sequential rotation of the basis for nearest-neighbour interaction, such that in the limit of a large number of Trotter steps it sums up to $\sum_{j=1}^{N-1}(J^{x}\sigma_{j}^{x}\sigma_{j+1}^{x}+J^{y}\sigma_{j}^{y}\sigma_{j+1}^{y}+J^{z}\sigma_{j}^{z}\sigma_{j+1}^{z})$. First, the XY unitary is performed, implementing $U_{\mathrm{XY}}=\exp\{-i\delta t\sum\limits_{j=1}^{N-1}\alpha_{xy}J(\sigma^{x}_{j}\sigma^{x}_{j+1}+\sigma^{y}_{j}\sigma^{y}_{j+1})\}$, where $\alpha_{xy}$ is some constant. Next, applying $\pi/2$ rotation around the $x$ axis for each qubit, $\hat{R}_{x}=\exp\{-i\frac{\pi}{4}\sum_{j=1}^{N}\sigma^{x}_{j}\}$, the two-qubit unitary can be transformed to $U_{\mathrm{XZ}}=\exp\{-i\delta t\sum\limits_{j=1}^{N-1}\alpha_{xz}J(\sigma^{x}_{j}\sigma^{x}_{j+1}+\sigma^{z}_{j}\sigma^{z}_{j+1})\}$. Subsequent $\pi/2$ rotation around the $y$ axis will implement the $U_{\mathrm{YZ}}=\exp\{-i\delta t\sum\limits_{j=1}^{N/2}\alpha_{yz}J(\sigma^{y}_{j}\sigma^{y}_{j+1}+\sigma^{z}_{j}\sigma^{z}_{j+1})\}$ interaction. For instance, the final chosen configuration of $J_{y}=2J^{x}/3$, $J^{z}=J^{x}/3$ can be achieved by choosing $J_{x}=J$ (can be different from $J_{\mathrm{sim}}$), $\alpha_{xy}=2/3$, $\alpha_{xz}=1/3$, $\alpha_{yz}=0$, simplifying the gate sequence. Finally, the $U_{\mathrm{Z}}$ operation introduces an effective magnetic field in the $z$ direction which allows for annealing to the ground state of the XYZ model. The linear schedule can then be designed similarly to the non-stoquastic case considered in Ref. Barends2016_SM . While the sequence represented in Fig. S5 will work perfectly in the $N_{\mathrm{Tr}}\rightarrow\infty$ limit, we note that the order of the unitaries $\{\mathcal{S}\}=\{U_{\mathrm{XY}},U_{\mathrm{XZ}},U_{\mathrm{YZ}},U_{\mathrm{Z}}\}$, which form a Trotter step, will alter the final infidelity for the annealed state. Thus, for the digital simulation procedure we consider 24 permutations of unitaries for the set $\mathcal{S}$, and choose the sequence of the step which yields minimal infidelity.

We note that alternatively the XYZ model can be simulated with controlled-phase (ZZ) gates as described by Barends et al. Barends2016_SM .