Lattice gauge theory and dynamical quantum phase transitions using noisy intermediate scale quantum devices
LLattice gauge theory and dynamical quantum phase transitions using noisyintermediate scale quantum devices
Simon Panyella Pedersen ∗ Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark
N. T. Zinner † Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark andAarhus Institute of Advanced Studies, Aarhus University, DK-8000 Aarhus C, Denmark (Dated: August 21, 2020)With noisy intermediate scale quantum devices in mind, we propose a class of superconductingcircuits for the general implementation of U(1) lattice gauge models via the formalism of quantumlink models. The circuit can be modularly scaled to any lattice configuration. Simulating the circuitdynamics with realistic circuit parameters we find that it implements the target dynamics with asteady average fidelity of 99 .
5% or higher. The principles of these circuits can be generalized toimplement other, more complicated gauge symmetries. Finally, we consider readout of the circuitusing a method that yields information about all the degrees of freedom with resonators coupleddispersively to only a subset of them. As an example of interesting physics to study with thiswe show how such a quantum link model on a periodic chain exhibits dynamical quantum phasetransitions by studying the Loschmidt amplitude and a novel gauge invariant string order parameter.The zeros of the Loschmidt amplitude as well as the zeros of our order parameter are revealed byvortices in their phases, which can be counted by a topologically invariant winding number.
I. INTRODUCTION
There is an increasing study of non-equilibrium quan-tum dynamics as improving experimental quantum con-trol makes it accessible [1]. Quantum simulators havebeen realized with cold atoms in optical lattices, ions, andsuperconducting quantum circuits (SQCs) among others,and have already been used to study exciting dynamicalphenomena like time crystals [2, 3], many-body localiza-tion [4, 5], prethermalization and thermalization [6–8],and particle-antiparticle creation and annihilation [9]. Anewly emerging addition to these is dynamical quantumphase transitions (DQPTs) [10–22]. These phase transi-tions have been studied experimentally [23–26], and offera broad spectrum of fascinating physics, like a connec-tion to topology [27–30], allowing for the definition ofdynamical topological order parameters [26, 31], vortexdynamics [32], scaling and universality [11], and a show-ing both a connection to underlying equilibrium phasetransitions [10, 28, 33], as well as being completely inde-pendent of them [34–36], the latter showing their trulynon-equilibrium nature. A particularly interesting typeof system for the study of dynamics is gauge theories,specifically lattice gauge theories (LGTs) [37–39]. Gaugetheories are at the basis of our understanding of particlephysics, and are notoriously difficult to handle both an-alytically and numerically. They are thus ideally suitedfor analog simulation [40–47]. DQPTs in gauge theorieshave been studied numerically [17, 30], as well as ana-lytically in the non-interacting limit in Ref. [30], but ∗ [email protected] † [email protected] have yet to be observed experimentally to the best of ourknowledge. In this work we show how to obtain LGTs, inthe form of quantum link models (QLM) [17, 40, 48–51],in SQCs in a fully consistent way, showing that we getthe desired Hamiltonians with very high fidelity. We usethe example of U(1) to demonstrate this, and we showthat exploring the new field of DQPTs is possible withNISQ-era devices [52]. Through use of the eigenmodesof the capacitive network (also known as normal or elec-trical modes) [53–57] the circuit implements three spin-1 / .
5% or higher,with most of the loss caused by leakage to higher levels,which could be further suppressed at the cost of slowerdynamics. The circuit can be scaled in a modular way toconstruct any desired spin lattice configuration. We thenstudy a particular U(1) symmetric system, the massiveSchwinger model [9, 30, 58, 59], exhibiting DQPTs af-ter a quench for all system sizes considered, and we showhow the smallest, non-trivial version can be realized withour SQC. Furthermore, we provide a readout scheme forhow to observe this in concrete setups, inspired by thatof Refs. [55, 56]. This makes it possible now to use NISQdevices to do precision studies of LGTs and DQPTs. a r X i v : . [ qu a n t - ph ] A ug II. RESULTSA. Target QLM Hamiltonian
The system we wish to realize with our circuit is thefollowing example of an interacting LGT: the (1+1)DU(1) invariant theory of fermions on a periodic latticeinteracting via an electric-like field. We represent thefermionic field with spinless, staggered mass fermionson the sites of the lattice, and transform these viathe Jordan-Wigner transformation [51, 60–62] into spin-1 / / / H = N − X n =0 (cid:20) − ( − n m σ zn + J (cid:0) σ + n S + n,n +1 σ − n +1 + H.c. (cid:1)(cid:21) (1)where N is the number of matter sites, which must beeven to conserve the symmetry between particles and an-tiparticles, m is the staggered mass of the fermions, and J is the matter-gauge coupling strength. σ αn and S αn,n +1 with α = z, + , − are Pauli-Z, step-up and -down matricespertaining to the matter field spin at site n , and the gaugefield spin on the link connecting site n and site n + 1, re-spectively. After going through our proposal for a circuitrealization of this system, we will consider quenches ofthe sign of the mass m → − m and the resulting DQPTs. B. Circuit realization
We present a circuit that implements two matter sitespins and a gauge link spin interacting via a direct three-body XXX-coupling, which through appropriate tuningyields the desired U(1) interaction seen in Eq. (1) in therotating wave approximation (RWA). The circuit scalesnaturally in a modular fashion, and could be used tocreate large 1D chains, 2D lattices, or any other config-uration of matter sites interacting through gauge links.Hence, the circuit could be used to experimentally imple-ment dynamics in 1D or 2D models, to study for examplevacuum quenches similar to what we will look at below,or strong CP-breaking in gauge models. Fig. 1c shows adiagrammatic implementation of a plaquette of four sites,hinting how a 2D configuration would have to be made.Early work was done to indicate that circuits could beused for simulating LGTs [62–71], but this work did notconsider concrete cases in detail, nor any checks whetherthe circuits actually realize the right dynamics with highfidelity. Here we do those things for the first time. We di-rectly compare the time evolution operator implementedby our circuit with the desired time evolution operatorof the U(1) spin-1 / K . Ideally this capacitance iszero, but has been included to study the effects of cou-pling to ground as well as capacitive coupling to controlor readout devices. There are four external fluxes, Φ i for i = 1 , , ,
4, threaded through loops that each con-sist of two Josephson junctions in parallel, in the purplebranches. The Josephson junctions themselves are imag-ined as implemented with SQUIDs such that the Joseph-son energy of each has an increased interval of possiblevalues and high tunability. In Fig. 1b the layout of thecircuit is chosen such that the Φ i only pass through thecircuit loops with the Josephson junctions pertaining to E s , using airbridges [75–80]. The modular scalability ofthe circuit has also been made explicit in this diagram.The idea is to make copies of the circuit in sequence,while using the same branches for the matter site spins.The circuit can thus be quite intuitively scaled to a chainof matter sites interacting through gauge links. A mat-ter site branch could potentially also be shared by morethan two copies of the circuit, making it possible to real-ize more complicated configurations. This would howeverresult in many wires connecting to the same branch andrequire many bridges. The ability to cross conductors viaairbridges makes SQCs a suitable platform to implementperiodic boundary conditions. Bridges make it possibleto access all points in a complicated circuit, while keep-ing it planar. In order to simulate the periodic systemdescribed by Eq. (1) for the simple case of N = 4 (whichwe show below hosts DQPTs), we would have four copiesof the circuit put together in this way, forming a square.A simple diagram of this can be seen in Fig. 1c. In thenext section we will consider readout of the circuit bydispersively coupling resonators to just the matter sitebranches, i.e. readout of the square would be done bycoupling resonators to its corner branches (the red andblue).In the following we will consider only a single copy ofthe circuit, using circuit parameters optimized for thatcase. Putting several together to form a square or somelarger system, the circuit parameters would have to be re-tuned. However, only the modes on the branches whichare shared with the new copies will be affected by them,i.e. the matter field modes. These do, however, havetheir own circuit parameters, which affect only the mat- φ φ φ φ C E CE c E s (cid:12) Φ C E c E s (cid:12) Φ CE c E s (cid:12) Φ C E C E c E s (cid:12) Φ K KKK φ φ φ φ (cid:12) Φ (cid:12) Φ (cid:12) Φ (cid:12) Φ a ) b ) c ) Figure 1. (a) The diagram of the circuit used to implement three spins interacting via a direct three-body XXX-coupling.Notice the identical circuit parameters. The identical grounding of each node has been included only for completeness ofthe study. We implement the spins via the eigenmodes of the capacitive network. The branches pertaining to each of thesemodes have been coloured separately. The blue and red branches will implement matter sites spins, while the purple branchesimplement the gauge link spin. (b) The same circuit but now folded to make the modular scalability of the circuit completelyclear. Multiple copies of the circuit, sharing the matter site branches pairwise, will implement a chain of matter sites coupledvia gauge links. (c) A simple diagram of how four copies of the circuit could be put together to implement the periodic N = 4version of Eq. (1). ter modes. These are the parameters with a subscript inFig. 1a on the red and blue branches. Hence, the effectof the new copies could be compensated for by changingthese parameters correspondingly. Hence, we are assuredthat if there is a good set of parameters for the case ofjust one copy, we could easily find parameters for severalcopies put together. In our study we have found thatit is not difficult to optimize the parameters for the cir-cuit, and in the Supplemental Information we show anexample of these that we use in the following.Truncating the anharmonic modes to their two lowestlevels, we find that this circuit implements the followingHamiltonian (see the Methods section) H s = −
12 Ω σ z −
12 Ω g σ zg −
12 Ω σ z + J z g σ z σ zg + J z σ z σ z + J zg σ zg σ z + J z g σ z σ zg σ z + J x g σ x σ xg σ x (2)where the Ω’s and J ’s are spin model parameters. Tocalculate the parameters we use a method introducedby the authors in Ref. [81], which avoids approximat-ing the trigonometric functions via a Taylor expansion,but instead takes their full effect into account. Thisgives more accurate parameters, when truncating the fluxHamiltonian, and can be used for any sine or cosine ofa linear combination of the flux coordinates. The ex- act dependence of the spin model parameters on the cir-cuit parameters, and details on their derivation can beseen in the Supplemental Information. The circuit hasthus resulted in three spins interacting through severalZ-type couplings, and a direct XXX-coupling. Throughappropriate tuning the XXX-coupling can be reduced to σ +0 σ + g σ − +H.c. in a RWA, which is exactly the U(1) gaugecoupling in Eq. (1), σ +0 S +0 , σ − + H.c.. The main effect ofthe principally undesirable Z-type couplings is to shift en-ergy levels, which we can compensate for. More impor-tantly there are interactions between the spin-1 / / / H s = H + H int ,where H contains all the Z-type terms and H int is justthe XXX-coupling. Consider H s in a frame rotating withrespect to H + m ( σ − σ ) H R = e iH t ( H − H ) e − iH t = − mσ z + 12 mσ z + J x g X p,r,s ∈{ + , −} e − iω prs t σ p σ rg σ s (3)where the sum is over all eight combinations of the three σ ± i , and the frequency of their phase is given by ω prs = p (Ω − m ) + r Ω g + s (Ω + m ) + 2 prsJ z g (4)If the system is now to tuned such that for example ω ++ − = − ω −− + = 0, then the operator σ +0 σ + g σ − + H.c.will be resonant, as desired. All other combinations willbe off-resonant, and would disappear in a RWA, as longas the spin transition frequencies and their differencesare much larger than the J z , which is already somethingwe must fulfil to justify the truncation to the spin- sub-space, see the Methods section. Furthermore, we have re-covered the staggered mass via the terms − mσ z + mσ z .Thus in an appropriately rotating frame, the Hamiltonianin Eq. (2) implemented by the circuit indeed recreates theone-dimensional U(1) quantum link model of Eq. (1) fortwo matter sites and the link between them.In finding appropriate circuit parameters, we optimizeaccording to the effective spin model parameters in orderto compensate for the effects of the higher levels. Wethen do a check of the overall behaviour of the circuit,ensuring that it works as intended, including no signifi-cant surviving interactions with higher levels. To do thiswe use average fidelity [72]. In the Methods section weexplain how our use of the average fidelity gives a directcomparison of the time evolution operators of circuit andthe target system, described in the previous, while takinginto account the higher levels of the circuit.While we have explained how the XXX-coupling in aRWA yields the desired U(1) interaction term, we do notactually use this approximation, but instead retain allterms to show that they do indeed not disturb the de-sired dynamics significantly. As mentioned, we presentthe specific set of optimized circuit parameters which isused in these simulations in the Supplemental Informa-tion. All simulations are performed without includingnoise. However, with present superconducting qubit lifetimes [74, 89, 90] we do not believe noise would signifi-cantly disturb the results presented here. In Fig. 2 thecalculated average fidelity of the circuit’s implementa-tion of the target dynamics can be seen in black. The Time [ns]0.990.9951 A v e r a g e fi d e li t y Figure 2. In black: The average fidelity of the circuit’simplementation of the target dynamics in a rotating frame.The fidelities are close to or above 99 .
5% and keep steady overlong times, with an oscillation on a short time scale. In green:The fidelity without taking leakage to the higher levels intoaccount. From this we see that the largest part, 0 . − . . fidelity is about or above 99 .
5% at all times, and whileit oscillates on a short timescale, it seems to keep steadyover the plotted interval. Hence, the implementation ofthe desired dynamics is good, and stable in the sensethat we are not accumulating error or continuously los-ing population to the higher levels. We seem to lose asmall fraction of the population immediately, which thenpartly oscillates back and forth. In green is plotted thesame average fidelity plus the leakage to higher levels, i.e.this plot shows the fidelity if we do not take leakage intoaccount. Hence, we can see that about 0 . − .
45% fidelityis lost because of population leaking to the higher levelsof the circuit, while about 0 .
05% is lost due to effectiveinteractions induced by virtual processes involving thehigher levels. These high and steady fidelities show di-rectly how our superconducting circuit truly implementsthe desired dynamics, with circuit parameters availableto experiments. Hence, the circuit is a strong candidatefor studying the U(1) QLM with present, NISQ-era de-vices.The circuit design principles we have used here, i.e.looking at the eigenmodes of the capacitive network ina symmetric circuit to achieve multi-body couplings andsuppressing as many undesired interactions as possible,could be used to achieve other interesting gauge invari-ant systems. It would be an obvious next step to worktowards higher gauge symmetries, like SU(2), or to at-tempt to implement gauge link operators with three lev-els. The latter would allow for the study of confinement,and might be implemented by using two spin-1 / C. Readout for state tomography
To perform readout of the circuit we would use amethod inspired by Refs. [55, 56]. Here they performquantum state tomography of two qubits by measuringthe dispersive shift of a resonator coupled to just one ofthem. The idea is that strong ZZ-couplings shift the en-ergy of one qubit conditioned on the state of the othersufficiently such that it can be seen in the dispersive shiftof the resonator. Hence, where normally one observes twoshifts of the resonator corresponding to the two eigenval-ues of σ z , one would see four shifts corresponding to thefour combinations of eigenvalues from the two qubits.Similarly, we imagine doing readout of just the mattersite spins, but still gain information about the gauge linkspins. A resonator coupled through identical capacitorsto the two nodes pertaining to a matter site mode (thered and blue branches in Fig. 1) will couple to only thatmode. Hence, usual dispersive readout of the spin can beperformed. The ZZ-couplings between this mode and itsneighbouring gauge and matter modes will make it pos-sible to derive some information about them as well. Itis easier to couple to the matter modes, as they live on asingle branch between two nodes, while the gauge modeslive on four branches between four nodes. We thereforepropose measuring on all matter modes and comparingthe data to extract information about the whole system.In the next section we consider quench dynamics ofEq. (1), finding DQPTs using the Loschmidt amplitude G ( t ) = h ψ (0) | ψ ( t ) i (i.e. the overlap between the initialand time-evolved state). An important thing to note isthat this quantity in the circuit’s own frame will notbe the same as in the rotating frame. This is essen-tially because it involves an odd number of time evo-lution operators. Transforming to a rotation frame via | ψ ( t ) i → U ( t ) | ψ ( t ) i , for some unitary, time dependentoperator U ( t ), the Loschmidt amplitude transforms as h ψ (0) | ψ ( t ) i → h ψ (0) | U ( t ) | ψ ( t ) i . These two are not iden-tical. However, the Loschmidt amplitude is often mea-sured by performing state tomography [91–95] and thencalculating G ( t ) from the results [21, 23–26, 32]. With fullinformation about the state of the circuit at any time,the Loschmidt amplitude can easily be calculated. Inour case we have U ( t ) = e iH , eff t , where H , eff is an ef-fective version of H , taking renormalization from higherlevel interactions into account. H , eff , which describesthe effective energies, could be determined in a separateexperiment, by setting E s = 0 via flux tuning, thus turn-ing off the XXX-coupling, and then initializing in eachof the spin-1 / ω r = ω r − g r ∆ + α − (cid:18) g r ∆ − g r ∆ + α (cid:19) σ z (5)Here ω r is the bare resonance frequency, σ z pertains tothe qubit, ∆ = ω q − ω r is the detuning between thetransition frequency ω q of the qubit and ω r , α is the an-harmonicity of the qubit, and finally g r is the strengthof the dispersive coupling. The transition frequency ω q is in our case a combination of the bare spin transitionfrequency and ZZ-coupling strengths. If we consider cou-pling a resonator to spin 0 in H s of Eq. (2), then we cansee that ω q = Ω − (cid:0) J z g σ zg + J z σ z + J z g σ zg σ z (cid:1) where the operators are to be understood as the specificeigenvalue they take on in the state that the circuit hascollapsed to as we measured it. We can now consider thedispersive shift χ ijk = ω r − ω r of the resonator frequencyas a function of the bare detuning ∆ = Ω − ω r , where i, j, k = 0 , , g, | i or | i . In Fig. 3 the eight shifts, corresponding to theeight spin-1 / and ∆ (with ∆ the equivalent of ∆ for the secondspin and resonator) such that the eight shifts are as dis-tinct as possible, and where comparing shifts from bothof the spins helps to determine the state of the gaugelink spin in addition to the two matter site spins. Inorder to remain in the dispersive region we must sat-isfy g r / | ∆ | , g r / | ∆ + α | (cid:28) / ω q . Looking at Eq. (5), we can see that χ ijk ∼ g r / | ∆ | , g r / | ∆ + α | , and thus the conditions forthe dispersive regime can be written as χ ijk (cid:28) g r . Thisessentially means we must stay well away from the re-gions where χ ijk diverges. These are marked with verticaldashed lines in Fig. 3. If for the sake of example we con-sider a resonator coupling strength of g r = 20 × π MHz,we get the energy scale shown on the right y -axis of Fig. 3.In the insets of Fig. 3 can be seen zoomed in regions whichare between χ ijk = ± × π MHz. If we choose a bare de-tuning within these regions, we could use the shift of thefirst resonator (the left inset) to distinguish between thedashed and solid lines, i.e. the state of spin 0, and use thesecond resonator (the right inset) to distinguish betweenthe blue/navy and the red/brown lines, i.e. the stateof spin 1. This is similar to dispersive measurement ofqubits, but we can then use this information to excludesome of the theoretical predictions of the measurements,making it easier to distinguish between the shifts caused − −
100 0 100 200 300∆ , [2 π MHz] − . − . . . . . . . χ i j k , [ g / π M H z ] − −
200 0 200 400 600 800∆ , [2 π MHz] ijk 000001010011100101110111 − − − − − − [ π M H z ] Figure 3. Left and right, the shift χ ijk of the resonance frequency of a resonator dispersively coupled to spin 0 or spin 1,respectively, for all eight spin-1 / n = Ω n − ω r for n = 0 , ijk = 0 , | i or | i , of each of the three spins. Assuming a resonatorcoupling strength of g r = 20 × π MHz, we get the energy scale shown on the right y -axis. The insets show intervals of the baredetunings where the dispersive shifts of each state are distinct enough to distinguish between the eight spin state with just thesetwo measurements, particularly when the information from each measurement is compared. To be in the dispersive regime wemust have χ ijk (cid:28) g r , i.e. we must stay well away from the points where χ ijk diverges, which are marked with vertical dashedlines. by the state of the gauge link spin. With a resolutionof 1 × π MHz in a measurement of the dispersive shift,which is experimentally feasible [24, 55, 56, 101], it wouldbe possible to distinguish the states of the circuit withthis or even a smaller choice of g r .The above analysis is approximate, as it uses only thebare spin model parameters, and the formula Eq. (5) isderived for a single qubit with some specific bare tran-sition frequency coupled to a resonator. In our case itis clear that the higher levels of the circuit would af-fect these calculations, and it is in fact a resonator, orseveral, coupled to a larger circuit, a system of multipleinteracting qubits or spins. A more accurate analysis us-ing numerical methods could be carried out to find theactual shifts of the resonance frequency of the resonatordependent on the circuit state. However, the quantita-tive results would be the same, namely that the differentstates would result in different shifts. It would then be amatter of determining whether those shifts would be suf-ficient to distinguish the states in a measurement, usingthe comparative method outlined above. Whether or notthis is the case is in the end a consequence of the cho-sen circuit parameters, so one could optimize the circuitparameters with respect to these considerations. Thesedispersive readouts could then be used to perform a fullquantum state tomography, yielding all information nec-essary to study the dynamics of the system. D. DQPTs in QLM chain
A recent study by Zache et al. has considered DQPTsin continuum and LGT models. In Ref. [30], DQPTswere found through the study of vorticity in an appropri-ate order parameter, implying that the transitions havea topological nature. Here, we consider the QLM equiv-alent of their system, described by
Eq. (1), over a largerrange of parameters and write down a gauge invariantstring order parameter as an example of what the aboveproposed realization of a spin-1 / G ( t ), fol-lowing a quench of the sign of the mass term, m → − m ,can be found by looking for vortices in their respectivephases in the Methods section. Hence, there is a topolog-ical aspect to the DQPTs we observe, both through theorder parameter, but also directly in the Loschmidt am-plitude. The vortices of the order parameter are dynam-ical in the sense that as a function of the matter-gaugeinteraction strength, they are created, move around, andcan annihilate with each other [32]. We find that theLoschmidt echo L ( t ) = |G ( t ) | and its zeros, correspond-ing to the zeros of the Loschmidt amplitude G ( t ) itself,converge to a certain structure for larger system sizes,particularly for small times. The zeros appear to be con-verging to lines in the thermodynamic limit. Similar re-sults have been found in other systems with DQPTs, andis similar to how the zeros of a partition function areknown to converge to lines in the thermodynamic limit[10, 15, 28, 102].This quench corresponds to a maximal quench of thevacuum angle. The vacuum angle is a parameter thatmay be included in Quantum Chromodynamics as wellas the Schwinger model, relating to the non-trivial struc-ture of their vacua [103, 104], and quantifying a CP-violating term. For more information see [59, 103–110].We will, however, not be considering the vacuum angleitself, but will rather focus on the quench, and the sub-sequent DQPTs. Explicitly we will be initializing thesystem in the ground state of the pre-quench Hamilto-nian H i = H ( m, J ) at time t = 0, and then perform uni-tary time evolution according to the post-quench Hamil-tonian H f = H ( − m, J ). We then study the Loschmidtamplitude, and Loschmidt echo L ( t ) = |G ( t ) | (the ze-ros of these quantities define the occurrence DQPTs), aswell as a both spatially and temporally non-local orderparameter g ( k, t ) = h ψ (0) | g ( k ) | ψ ( t ) i that we introduce,where g ( k ) is the Fourier transform of the gauge invariantstring order parameters of the system, which we discussfurther in the Methods section. We study the dynamicsof the Loschmidt amplitude and the order parameter fordifferent system sizes, N = 2 , , ..., L ( t ), for all system sizes. Asmentioned these zeros are accompanied by vortices in thephase of the order parameter, which can be counted by awinding number (see the Methods section). Hence, evenin the smallest system this order parameter reveals thestructure of L ( t ), which repeats itself as the system size isincreased, and exhibits non-trivial vortex dynamics. Theorder parameter that we employ here is related to thegauge-invariant time-ordered Green’s function computedin the context of lattice gauge theory in [30] and our workdemonstrates how to transfer this to QLM models.In Fig. 4 can be seen L ( t ) for the considered range of J/m and t , for system sizes N = 4 , ,
16. The zeros of L ( t ) are marked with circles colour-coded according tothe orientation of the vortices, i.e. whether they windupwards counter-clockwise or clockwise, correspondingto being right- or left-handed. It is unclear whetherthe orientation of the vortices in the phase of G ( t ) hasany significance. Most of them are left-winding with afew right-winding at late times in the middle and lowerframes of Fig. 4. The zeros of L ( t ) are found by apply-ing the method described above to the phase, φ G ( t ) , of G ( t ). This phase has two components φ G ( t ) = φ dyn + φ P ,where φ dyn is the dynamical phase defined by φ dyn = − R t dt h ψ ( t ) | H f | ψ ( t ) i = − t h ψ (0) | H f | ψ (0) i , while the φ P is a purely geometric phase called the Pancharatnamgeometric phase [111], which is an extension of the con-cept of Berry’s phase [21, 111, 112]. It is the Pancharat-nam phase which contains the vortices, while the dynam-ical varies smoothly as a function of J/m and is linear in t . The zeros of the order parameter are indicated withtriangles and a different colour code for distinguishabil-ity. We can see how the zeros of the order parameter liealong the troughs of L ( t ). Intuitively, one might think of Figure 4. Contour plots of L ( t ) for the considered range of J/m and t , for system sizes N = 4 , ,
16 from top to bottom.The zeros, as found by considering vortices in the phase of G ( t ), are marked with circles coloured blue for right-windingvortices, and yellow for left-winding. A representative setof the vortices of the order parameter are plotted with arrowheads. White arrow heads pointing to the right indicate right-winding vortices, and red arrow heads pointing left indicateleft-winding vortices. These trace out the continuous lines inthe ( J, t )-plane where the order parameter vanishes. Partic-ularly for N = 4 in the top panel it can be clearly seen howvortices of opposite orientation move around the ( J, t )-planeand annihilate. The structure of L ( t ) and its zeros has a clearpattern that is present in all three plots. For N = 16 L ( t )is very close to zero in large areas, and its zeros, particularlyat early times, trace out a curve following the center of theselobes. the vortices as charges, and when two charges of oppo-site sign are near each other, they screen each other off,allowing them to be in areas where L ( t ) takes on largervalues. For N = 4 in the upper panel of Fig. 4 it is partic-ularly clear in for example the region 1 < J/m <
2, howvortices in the phase of the order parameter of oppositeorientation are created as a pair, move in the
J, t -plane,and annihilate with each other. Furthermore, we see howthe order parameter does not have zeros for very smalland very large
J/m . Intuitively it is reasonable that forboth
J/m (cid:28)
J/m (cid:29) L ( t ) being small in largelobes that start at around J/m = 0 . J/m while L ( t ) increases. The zeros of L ( t ) occurnear the center of these lobes. This patterns appearsalready at N = 4 and repeats itself for larger systemsizes. We thus see how even the small system with N = 4reproduces features of the much larger N = 16 system,where the behaviour of L ( t ) has started to converge. Thismeans that even a small experimental realization of thissystem could yield interesting results. While the evensmaller system with N = 2 does show zeros both in g and L ( t ) the behaviour of these is considerably differentand simpler than for larger N . We explain this differentbehaviour by its small size. Its dynamics are necessarilyquite different simply because it is so limited. Hence,while that system may be interesting for experimentalstudies in itself, it does not reflect the behaviour of large N systems as well as N = 4 does. In a future work wewill elaborate on our results regarding this quench of aU(1) QLM.From the above we see how this systems hosts non-trivial post-quench dynamics, and how even a small ver-sion of the system reflects the behaviour of a much largerversions. This makes it relevant for experimental studyusing the SQC we proposed above. Using quantum statetomography with the proposed method of measurement,all information necessary for calculating the Loschmidtamplitude and our order parameter could be extractedto study the small version of the system. III. DISCUSSION
We have shown how to realize lattice gauge theoriesin superconducting quantum circuits. Specifically, wehave provided a method for general circuits to implementquantum link models with a high average fidelity. Thisopens up the possibility for experimental study of quan-tum link models in NISQ-era devices. We have proposeda superconducting circuit, which realizes three spin-1 / .
5% or above, using realistic circuitparameters. We have proposed an approach for read-out of the circuit, using a method of resonators coupleddispersively to a subset of the circuit spins, but whichnonetheless gives information about all the spins, by ex- ploiting their pairwise ZZ-couplings. Using the same setof circuit parameters again we showed how this methodcan principally reduce the complexity of measurement ofthe circuit’s state.As a demonstration of the principles in our work wehave studied a periodic (1+1)D spin-1 / Z n would beinteresting to develop. Further work could also be doneon the specific system studied here. It would be inter-esting to link the dynamical quantum phase transitionsfound here to a potential underlying equilibrium phasetransition or entropy production [19, 21, 23]. IV. METHODSA. Circuit Hamiltonian
We will here give a few details on the derivation of thespin Hamiltonian of our circuit, as shown in Fig. 1. Wedefine the node fluxes of the circuit φ = ( φ , φ , φ , φ ) T [57], but will be working in the eigenmodes of the ca-pacitive network [53–56, 81], ψ = ( ψ CM , ψ , ψ g , ψ ) T ,defined through φ = −
01 0 −
11 0 − − ψ This results in no interactions through the capacitors,greatly reducing the complexity of the interactions in thesystem. We will furthermore introduce a new set of ex-ternal fluxes, Ψ , Ψ g and Ψ , of which the Φ i are certainsimple, linear combinations. We set these new externalfluxes to be constant Ψ j = − π/ j = 0 , g,
1. In thesecoordinates and with these choices of external fluxes thecircuit Hamiltonian becomes H c = K − q + K − gg q g + K − q − E cos ψ − E cos ψ − E c cos ψ cos ψ g cos ψ − E s sin ψ sin ψ g sin ψ (6)where the q j are momentum variables conjugate to the ψ j , and K − jj are the diagonal entries of the inverse capac-itance matrix in the ψ j coordinates. This Hamiltoniandescribes three transmon-like [73] anharmonic oscillatormodes, interacting only through the interesting triple co-sine and sine interactions. These are direct, completelyeven and completely odd, three-body interactions. Thesine functions come about as a consequence of settingΨ j = − π/
2. Recasting each of the ψ j and q j variablesin terms of harmonic oscillator operators, i.e. bosoniccreation and annihilation operators, a † j and a j , and trun-cating the system to the two lowest levels of each, yieldsthe spin Hamiltonian in Eq. (2). This shows how the ψ - and ψ -modes will represent matter site spins, whilethe ψ g -mode has the role of gauge link spin. For a fullderivation of the Hamiltonian, and explicit expressionsfor the spin model parameters, see the Supplemental In-formation. B. Z-type couplings and working regime
Though the Z-type couplings in Eq. (2) are principallyundesirable, their effect is mainly to make the tuningof the circuit slightly more complicated. They do notdisturb the main and interesting feature of the system.Looking at the Hamiltonian, the Z-type couplings essen-tially just shift the energy levels of the system, but itturns out that this does not affect the detuning of thestates we are interested in, except for the ZZZ-couplingwhich, however, is a weak coupling. More importantly,however, is the fact that the coupling strengths of theZ-type couplings are representative for even interactionsbetween the spin-1 / a † j a † j a j a j + H.c. or a † j a † j a † j a † j + H.c., where there is an even number of cre-ation and annihilation operators, a † j and a j , for each an-harmonic mode. When working in the spin-1 / one excitation, these inter-action will either simply have no effect because they cannot remove two excitations from the modes, or they willbe suppressed if the coupling strengths are much smallerthan the spin transition energies. However, higher ordercontributions from such interactions will turn out to af-fect both the detuning of | g i and | g i , and thestrength of the XXX-coupling. That is, the system will undergo virtual excitations and de-excitations, which ef-fectively change the spin model parameters. Again, thisdoes not disturb the main feature of the circuit, but sim-ply renormalizes the parameters, making the tuning ofthe circuit more complicated, as we must now tune effec-tive parameters and not the bare ones which appear inEq. (2). Specifically, we are interested in tuning | g i and | g i into resonance, with all other states de-tuned. Because of the staggered mass in Eq. (1), the firstof these two states represents a particle-antiparticle pairconnected by a gauge field flux tube, and the second rep-resents two empty sites with the gauge field pointing theopposite direction. If these are in resonance while all oth-ers are not, the XXX-coupling will implement the propergauge invariant interaction, reducing to σ +0 σ + g σ − + H.c..This corresponds to ω ++ − = − ω −− + = 0, as discussedbelow Eq. (4).Though the Z-type couplings would all disappear bysetting E c = 0, i.e. removing the junctions pertain-ing to E c , this would result in the anharmonicity of thegauge link mode vanishing. The anharmonicities α j ofthe three modes, which justify the truncation to the twolowest level of each anharmonic oscillator, can be seen inthe Supplemental Information. We note that α g is in-deed proportional to E c , because the gauge link mode isonly affected by the Josephson junctions on the purplebranches. The anharmonicities of the matter site modes, α and α , each have contributions from both E c , and E or E respectively. The truncation is justified if interac-tions between the spin- subspace and higher levels of theHilbert space of Eq. (2) are suppressed. This will, brieflyput, be the case if even interaction strengths are muchsmaller than spin transition frequencies and their differ-ences, and if odd interaction strengths are much smallerthan the anharmonicities, i.e. for our system if J z (cid:28) Ω j and J x g (cid:28) α j . C. Comparison of circuit and target dynamics
To compare the dynamics of our circuit, as describedby Eq. (6), with the target dynamics, described byEq. (1), we use average fidelity. Average fidelity is a mea-sure of how well a certain process implements a desiredoperation. In our case the process is the time evolutionof the circuit, determined by the full circuit Hamiltonianin Eq. (6), in a frame rotating such that the resultingHamiltonian is H R in Eq. (3), where m is set to half theeffective detuning between | g i and | g i (here 0and 1 refer to the ground and excited state of the spin,and the subscripts denote which spin it is), and the barecoupling strength J x g is replaced with the effective cou-pling between these two states. To take contributionsbeyond renormalized parameters from higher level inter-actions into account, we truncate the anharmonic modesof the circuit to the lowest four levels, when simulatingdynamics. We only rotate the spin-1 / H + m ( σ − σ ), where all en-0tries pertaining to levels higher than the spin-1 / H from Eq. (1) with two matter sites and a gauge linkbetween them. We compare only the dynamics of thespin-1 / / / H . For details on how the averagefidelity is calculated, using the clever formula of Ref. [72],see the Supplemental Information. The average fidelitythus compares the very time evolution operators them-selves for the circuit and the target system. Followingour analysis, and as the mass and coupling strength ofthe target Hamiltonian are chosen to be half the effectivedetuning and the effective coupling strength respectively,the fidelity is a priori expected to be quite high. Sincethis is all done with four levels included in each anhar-monic mode, the fidelity will, however, be a measure ofhow much the higher levels affect the dynamics of thecircuit beyond just the renormalization of the mass andcoupling strength. In particular, some population willbe lost to the higher levels, and just as virtual processescontribute to the strength of the XXX-coupling, they willalso to some extent induce other effective interactions.These will disturb the desired dynamics and might begauge variant, resulting in population moving outside ofthe chosen gauge sector of the spin-1 / D. Order parameter
Our order parameter is g ( k, t ) = h ψ (0) | g ( k ) | ψ ( t ) i ,where | ψ ( t ) i is the state at time t and g ( k ) = X m =0 , N − X n =0 e − ikd m ( n ) σ − m n − Y i = m S α m ( n ) i,i +1 σ + n This is a sum over two representative sites m = 0 , n = 0 , ..., N −
1. The summand is a Fouriercoefficient times a string operator, consisting of two mat-ter site operators, one at the representative site m andthe other at n , connected by the gauge link operatorsbetween the two sites, making the total operator gaugeinvariant. Here the products of link operators are overthe shortest path between site m and site n , i.e. eithercounter-clockwise or clockwise along the circular lattice,see Fig. 5. Likewise d m ( n ) is the distance from site m to n along the shortest path, with d m ( n ) being positive forclockwise paths and negative for counter-clockwise paths.Similarly, α m ( n ) = − for clockwise paths and α m ( n ) = +for counter-clockwise paths, ensuring the gauge invari-ance of the summands. For sites on the exact oppositeside of the circular lattice, i.e. m − n = N/
2, the two paths around the lattice are equidistant and so both areincluded, see Fig. 5b. Thus the operator g ( k ) is essen-tially the Fourier transform of the gauge invariant stringoperators connecting the sites 0 and 1 with all sites of thelattice. The order parameter g ( k, t ) is then the Fouriertransform of the amplitudes of a matter excitation mov-ing from either site 0 or 1 to site n , via the shortest path,in the time between initialization and t , in a gauge in-variant manner. The reason we have both a term for site0 and one for site 1 is to make the operator symmetricwith respect to particles and antiparticles, as site 0 is aparticle site while site 1 is an antiparticle site. E. Counting vortices
Considering G ( t ) as a complex function in the ( J, t )-plane and g ( k, t ) a complex function in the ( k, t )-plane(with J fixed), we can find their zeros in a similar fashion,namely by looking at their phase. Zeros of a smooth com-plex function of two variables are accompanied by vor-tices in the phase of that function. This leads to a methodfor numerically finding the vortices as detailed below,which is an adaptation of the work in Ref. [113], devel-oped for computing Chern numbers in momentum space.When a complex function, say f = re iϕ , becomes zero atsome critical point, its phase ϕ is undefined at that point.For a smooth function of two variables, f = f ( x, y ), thisresults in a vortex in ϕ surrounding the critical point.The intuition of this is that while ϕ is undefined at thecritical point it is otherwise smooth, up to discontinuitiesof 2 π . If there were no discontinuity around the criticalpoint there would obviously be a meaningful, smooth ex-tension of ϕ at the undefined point, which contradictsthe fact that f goes to zero. The phase must thus have aline of 2 π discontinuity extending from the critical point,which we will refer to as tears (as in torn fabric). Startinga this discontinuity and going around the critical point, ϕ must then attain all possible values between − π and π in a smooth way, as there would otherwise again bea meaningful, smooth extension of ϕ , or discontinuitieswith a different value than 2 π . Hence, going around thecritical point in a closed curve, ϕ will go through all val-ues from − π to π , and have a discontinuity of 2 π betweenthe extremal values. Such vortices can be counted by awinding number ν = 12 π I C d l · ∇ ϕ where C is a closed curve. This number may then be con-sidered a dynamical topological order parameter [26, 31],as it is a parameter changing its value with time, takingon discrete values which only depend on the topology of ϕ , i.e. its vortices, and whether C encloses these vortices.Such a winding number essentially detects and counts thenumber of times ϕ has discontinuously jumped from π to − π along C , minus the number of times ϕ has jumpedfrom − π to π , i.e. how many times in total ϕ increases1 a ) d ( ) = − , α ( ) = + d ( ) = , α ( ) = − b ) d ( ) = − , α ( ) =+ d ( ) = , α ( ) = − Figure 5. Some examples of the paths taken by the string operators which are summed over in the order parameter for asystem size of N = 8 matter sites. a) Two examples, one originating from particle site 0 and the other from the antiparticlesite 1, showing the sign conventions for clockwise and counter-clockwise paths. b) An example of the case where m − n = N/ by 2 π along the closed curve. Hence, closing such a curvetightly around a vortex, the winding number essentiallyjust detects that a tear enters the area enclosed by thecurve without exiting again. Together with the image ofthe vortices always being tailed by these tears, we see thatto find the vortices we must simply be able to identifythe tears and their ends. The tears can only end eitherat the edge of the considered parameter space or at avortex. Hence, one only needs sufficient resolution (in asimulation or data) to distinguish discontinuous jumps of2 π from the jumps between the data points on a coarsegrid in order to find the zeros of the function f . Thismakes it possible to find and study zeros and vorticeswith a minimal numerical computational effort, and wehave used an algorithm based on this idea to do so in oursystem. Furthermore, one can be certain that these willbe true zeros of the function, and not points where thecomplex function merely has a very small modulus. Thisis otherwise principally quite hard to do, as unavoidablenumerical imprecision would usually make it necessaryto set an arbitrary limit on when the modulus is small enough to indicate that the function has actually becomezero. Essentially, the vortices are easy to find numericallyeven with low resolution because they are extended struc-tures (as opposed to the single point where the functionvanishes). This extended nature is also the reason theycan be counted by a winding number. V. DATA AVAILABILITY
The data that support the findings in this study areavailable from the corresponding author upon reasonablerequest.
VI. CODE AVAILABILITY
All code used to generate the presented data in thispaper is available upon reasonable request. [1] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum sim-ulation, Reviews of Modern Physics , 153 (2014).[2] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee,J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter,A. Vishwanath, N. Y. Yao, and C. Monroe, Observationof a discrete time crystal, Nature , 217 (2017).[3] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou,J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani,C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Observation of discrete time-crystalline order ina disordered dipolar many-body system, Nature ,221 (2017).[4] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Luschen,M. H. Fischer, R. Vosk, E. Altman, U. Schneider, andI. Bloch, Observation of many-body localization of inter-acting fermions in a quasirandom optical lattice, Science , 842 (2015). [5] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W.Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe,Many-body localization in a quantum simulator withprogrammable random disorder, Nature Physics , 907(2016).[6] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa,B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Dem-ler, and J. Schmiedmayer, Relaxation and Prethermal-ization in an Isolated Quantum System, Science ,1318 (2012).[7] B. Neyenhuis, J. Zhang, P. W. Hess, J. Smith, A. C. Lee,P. Richerme, Z.-X. Gong, A. V. Gorshkov, and C. Mon-roe, Observation of prethermalization in long-range in-teracting spin chains, Science Advances , e1700672(2017).[8] C. Neill, P. Roushan, M. Fang, Y. Chen, M. Kolodru-betz, Z. Chen, A. Megrant, R. Barends, B. Campbell,B. Chiaro, A. Dunsworth, E. Jeffrey, J. Kelly, J. Mutus,P. J. J. O’Malley, C. Quintana, D. Sank, A. Vainsencher,J. Wenner, T. C. White, A. Polkovnikov, and J. M. Mar-tinis, Ergodic dynamics and thermalization in an iso-lated quantum system, Nature Physics , 1037 (2016).[9] E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg,A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz,P. Zoller, and R. Blatt, Real-time dynamics of latticegauge theories with a few-qubit quantum computer, Na-ture , 516 (2016), arXiv:arXiv:1605.04570v1.[10] M. Heyl, A. Polkovnikov, and S. Kehrein, Dynami-cal quantum phase transitions in the transverse-fieldising model, Physical Review Letters , 1 (2013),arXiv:arXiv:1206.2505v2.[11] M. Heyl, Scaling and Universality at Dynamical Quan-tum Phase Transitions, Physical Review Letters , 1(2015), arXiv:arXiv:1505.02352v2.[12] E. Canovi, P. Werner, and M. Eckstein, First-OrderDynamical Phase Transitions, Physical Review Letters , 265702 (2014), arXiv:1408.1795.[13] D. Pekker, G. Refael, E. Altman, E. Demler, andV. Oganesyan, Hilbert-Glass Transition: New Uni-versality of Temperature-Tuned Many-Body Dynami-cal Quantum Criticality, Physical Review X , 011052(2014).[14] R. Vosk and E. Altman, Dynamical Quantum PhaseTransitions in Random Spin Chains, Physical ReviewLetters , 217204 (2014).[15] M. Schmitt and S. Kehrein, Dynamical quantum phasetransitions in the Kitaev honeycomb model, PhysicalReview B , 075114 (2015).[16] A. A. Zvyagin, Nonequilibrium dynamics of a systemwith two kinds of fermions after a pulse, Physical Re-view B , 075122 (2017).[17] Y.-P. Huang, D. Banerjee, and M. Heyl, DynamicalQuantum Phase Transitions in U(1) Quantum LinkModels, Physical Review Letters , 250401 (2019),arXiv:1808.07874.[18] M. Lacki and M. Heyl, Dynamical quantum phasetransitions in collapse and revival oscillations of aquenched superfluid, Physical Review B , 1 (2019),arXiv:1812.02209.[19] B. O. Goes, G. T. Landi, E. Solano, M. Sanz, and L. C.C´eleri, Wehrl entropy production rate across a dynami-cal quantum phase transition (2020), arXiv:2004.01126.[20] A. A. Zvyagin, Dynamical quantum phase transitions(Review Article), Low Temperature Physics , 971 (2016), arXiv:1701.08851.[21] M. Heyl, Dynamical quantum phase transitions: A re-view, Reports on Progress in Physics , 10.1088/1361-6633/aaaf9a (2018), arXiv:1709.07461.[22] M. Heyl, Dynamical quantum phase transitions: A briefsurvey, Epl , 25 (2019), arXiv:arXiv:1811.02575v1.[23] P. Jurcevic, H. Shen, P. Hauke, C. Maier, T. Bry-dges, C. Hempel, B. P. Lanyon, M. Heyl, R. Blatt, andC. F. Roos, Direct Observation of Dynamical Quan-tum Phase Transitions in an Interacting Many-BodySystem, Physical Review Letters , 080501 (2017),arXiv:1612.06902.[24] X. Y. Guo, C. Yang, Y. Zeng, Y. Peng, H. K. Li,H. Deng, Y. R. Jin, S. Chen, D. Zheng, and H. Fan,Observation of a Dynamical Quantum Phase Transitionby a Superconducting Qubit Simulation, Physical Re-view Applied , 1 (2019), arXiv:1806.09269.[25] K. Xu, Z.-H. Sun, W. Liu, Y.-R. Zhang, H. Li,H. Dong, W. Ren, P. Zhang, F. Nori, D. Zheng, H. Fan,and H. Wang, Probing the dynamical phase transi-tion with a superconducting quantum simulator (2019),arXiv:1912.05150.[26] X.-Y. Xu, Q.-Q. Wang, M. Heyl, J. C. Budich, W.-W.Pan, Z. Chen, M. Jan, K. Sun, J.-S. Xu, Y.-J. Han, C.-F.Li, and G.-C. Guo, Measuring a dynamical topologicalorder parameter in quantum walks, Light: Science &Applications , 7 (2020), arXiv:1808.03930.[27] H. Hu and E. Zhao, Dynamical topology of quantumquenches in two dimensions (2019), arXiv:1911.02211.[28] S. Vajna and B. D´ora, Topological classification of dy-namical phase transitions, Physical Review B - Con-densed Matter and Materials Physics , 10.1103/Phys-RevB.91.155127 (2015), arXiv:arXiv:1409.7019v1.[29] I. Hagym´asi, C. Hubig, Legeza, and U. Schollw¨ock,Dynamical Topological Quantum Phase Transitions inNonintegrable Models, Physical Review Letters , 1(2019), arXiv:1904.00867.[30] T. V. Zache, N. Mueller, J. T. Schneider, F. Jendrze-jewski, J. Berges, and P. Hauke, Dynamical Topolog-ical Transitions in the Massive Schwinger Model witha θ Term, Physical Review Letters , 1 (2019),arXiv:arXiv:1808.07885v1.[31] J. C. Budich and M. Heyl, Dynamical topologi-cal order parameters far from equilibrium, PhysicalReview B , 10.1103/PhysRevB.93.085416 (2016),arXiv:arXiv:1504.05599v3.[32] N. Fl¨aschner, D. Vogel, M. Tarnowski, B. S. Rem, D.-S.L¨uhmann, M. Heyl, J. C. Budich, L. Mathey, K. Sen-gstock, and C. Weitenberg, Observation of dynamicalvortices after quenches in a system with topology, Na-ture Physics , 265 (2018).[33] Z. Huang and A. V. Balatsky, Dynamical Quan-tum Phase Transitions: Role of Topological Nodes inWave Function Overlaps, Physical Review Letters ,086802 (2016).[34] M. Fagotti, Dynamical Phase Transitions as Proper-ties of the Stationary State: Analytic Results afterQuantum Quenches in the Spin-1/2 XXZ Chain (2013),arXiv:1308.0277.[35] S. Vajna and B. D´ora, Disentangling dynamical phasetransitions from equilibrium phase transitions, PhysicalReview B , 161105 (2014).[36] F. Andraschko and J. Sirker, Dynamical quantum phasetransitions and the Loschmidt echo: A transfer matrix approach, Physical Review B , 125120 (2014).[37] K. G. Wilson, Confinement of quarks, Phys. Rev. D ,2445 (1974).[38] J. Kogut and L. Susskind, Hamiltonian formulation ofWilson’s lattice gauge theories, Physical Review D ,395 (1975).[39] J. Smit, Introduction to Quantum Fields on a Lattice ,Cambridge Lecture Notes in Physics (Cambridge Uni-versity Press, 2002) pp. 1–284.[40] D. Banerjee, M. Dalmonte, M. M¨uller, E. Rico, P. Ste-bler, U.-J. Wiese, and P. Zoller, Atomic Quantum Simu-lation of Dynamical Gauge Fields Coupled to FermionicMatter: From String Breaking to Evolution after aQuench, Physical Review Letters , 175302 (2012),arXiv:arXiv:1205.6366v2.[41] S. K¨uhn, J. I. Cirac, and M. C. Ba˜nuls, Quantum sim-ulation of the Schwinger model: A study of feasibility,Physical Review A - Atomic, Molecular, and OpticalPhysics , 1 (2014).[42] V. Kasper, F. Hebenstreit, M. Oberthaler, andJ. Berges, Schwinger pair production with ultracoldatoms, Physics Letters B , 742 (2016).[43] A. S. Dehkharghani, E. Rico, N. T. Zinner, and A. Ne-gretti, Quantum simulation of Abelian lattice gauge the-ories via state-dependent hopping, Physical Review A , 043611 (2017).[44] A. Mil, T. V. Zache, A. Hegde, A. Xia, R. P. Bhatt,M. K. Oberthaler, P. Hauke, J. Berges, and F. Jen-drzejewski, Realizing a scalable building block of aU(1) gauge theory with cold atomic mixtures (2019),arXiv:1909.07641.[45] L. Sanchez-Palencia, Constructing Field Theories UsingQuantum Simulators, Physics , 10 (2020).[46] T. V. Zache, T. Schweigler, S. Erne, J. Schmiedmayer,and J. Berges, Extracting the Field Theory Descrip-tion of a Quantum Many-Body System from Exper-imental Data, Physical Review X , 011020 (2020),arXiv:1909.12815.[47] B. Yang, H. Sun, R. Ott, H.-Y. Wang, T. V. Zache,J. C. Halimeh, Z.-S. Yuan, P. Hauke, and J.-W. Pan,Observation of gauge invariance in a 71-site quantumsimulator (2020), arXiv:2003.08945.[48] D. Horn, Finite matrix models with continuous localgauge invariance, Physics Letters B , 149 (1981).[49] P. Orland and D. Rohrlich, Lattice gauge magnets:Local isospin from spin, Nuclear Physics B , 647(1990).[50] S. Chandrasekharan and U.-J. Wiese, Quantum linkmodels: A discrete approach to gauge theories, NuclearPhysics B , 455 (1997), arXiv:9609042v2 [arXiv:hep-lat].[51] P. Hauke, D. Marcos, M. Dalmonte, and P. Zoller,Quantum Simulation of a Lattice Schwinger Model ina Chain of Trapped Ions, Physical Review X , 041018(2013), arXiv:arXiv:1306.2162v3.[52] J. Preskill, Quantum Computing in the NISQ era andbeyond, Quantum , 79 (2018).[53] N. Bergeal, R. Vijay, V. E. Manucharyan, I. Siddiqi,R. J. Schoelkopf, S. M. Girvin, and M. H. Devoret,Analog information processing at the quantum limitwith a Josephson ring modulator, Nature Physics , 296(2010).[54] M. Kounalakis, C. Dickel, A. Bruno, N. K. Lang-ford, and G. A. Steele, Tuneable hopping and nonlinear cross-Kerr interactions in a high-coherence supercon-ducting circuit, npj Quantum Information , 38 (2018),arXiv:1802.10037.[55] T. Roy, S. Kundu, M. Chand, S. Hazra, N. Nehra,R. Cosmic, A. Ranadive, M. P. Patankar, K. Damle, andR. Vijay, Implementation of Pairwise Longitudinal Cou-pling in a Three-Qubit Superconducting Circuit, Physi-cal Review Applied , 054025 (2017), arXiv:1610.07915.[56] T. Roy, S. Hazra, S. Kundu, M. Chand, M. P.Patankar, and R. Vijay, A programmable three-qubitsuperconducting processor with all-to-all connectivity,arXiv:1809.00668 (2018).[57] U. Vool and M. H. Devoret, Introduction to quan-tum electromagnetic circuits, International Journalof Circuit Theory and Applications , 897 (2017),arXiv:arXiv:1610.03438v2.[58] J. Schwinger, Gauge Invariance and Mass. II, Phys. Rev. , 2425 (1962).[59] S. Coleman, R. Jackiw, and L. Susskind, Charge Shield-ing and Quark Confinement in the Massive SchwingerModel*, Annals of Physics , 267 (1975).[60] P. Jordan and E. Wigner, ¨Uber das Paulische ¨Aquiv-alenzverbot, Zeitschrift f¨ur Physik , 631 (1928).[61] L. Susskind, Lattice fermions, Physical Review D ,3031 (1977).[62] D. Marcos, P. Rabl, E. Rico, and P. Zoller,Superconducting Circuits for Quantum Simulationof Dynamical Gauge Fields, Physical Review Let-ters , 10.1103/PhysRevLett.111.110504 (2013),arXiv:arXiv:1306.1674v2.[63] D. Marcos, P. Widmer, E. Rico, M. Hafezi, P. Rabl, U.-J. Wiese, and P. Zoller, Two-dimensional lattice gaugetheories with superconducting quantum circuits, Annalsof Physics , 634 (2014), arXiv:arXiv:1407.6066v2.[64] A. Mezzacapo, E. Rico, C. Sab´ın, I. L. Egusquiza,L. Lamata, and E. Solano, Non-Abelian SU(2)Lattice Gauge Theories in Superconducting Cir-cuits, Phys. Rev. Lett. , 240502 (2015),arXiv:arXiv:1505.04720v2.[65] K. Le Hur, L. Henriet, A. Petrescu, K. Plekhanov,G. Roux, and M. Schir´o, Many-body quantum electro-dynamics networks: Non-equilibrium condensed matterphysics with light, Comptes Rendus Physique , 808(2016).[66] M. J. Hartmann, Quantum simulation with interactingphotons, Journal of Optics , 104005 (2016).[67] G. K. Brennen, G. Pupillo, E. Rico, T. M. Stace, andD. Vodola, Loops and Strings in a Superconducting Lat-tice Gauge Simulator, Physical Review Letters , 1(2016), arXiv:arXiv:1512.06565v2.[68] H. Alaeian, C. W. S. Chang, M. V. Moghaddam, C. M.Wilson, E. Solano, and E. Rico, Lattice gauge fields viamodulation in circuit QED: The bosonic Creutz ladder(2018), arXiv:1805.12410.[69] N. Klco, E. F. Dumitrescu, A. J. McCaskey, T. D. Mor-ris, R. C. Pooser, M. Sanz, E. Solano, P. Lougovski,and M. J. Savage, Quantum-classical computation ofSchwinger model dynamics using quantum computers,Physical Review A , 032331 (2018).[70] P. Roushan, C. Neill, A. Megrant, Y. Chen, R. Bab-bush, R. Barends, B. Campbell, Z. Chen, B. Chiaro,A. Dunsworth, A. Fowler, E. Jeffrey, J. Kelly, E. Lucero,J. Mutus, P. J. O’Malley, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White,E. Kapit, H. Neven, and J. M. Martinis, Chiralground-state currents of interacting photons in a syn-thetic magnetic field, Nature Physics , 146 (2017),arXiv:arXiv:1606.00077v3.[71] Z. H. Yang, Y. P. Wang, Z. Y. Xue, W. L. Yang, Y. Hu,J. H. Gao, and Y. Wu, Circuit quantum electrodynamicssimulator of flat band physics in a Lieb lattice, PhysicalReview A , 1 (2016), arXiv:arXiv:1603.04686v2.[72] M. A. Nielsen, A simple formula for the average gatefidelity of a quantum dynamical operation, Physics Let-ters A , 249 (2002).[73] J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I.Schuster, J. Majer, A. Blais, M. H. Devoret, S. M.Girvin, and R. J. Schoelkopf, Charge-insensitive qubitdesign derived from the Cooper pair box, Phys. Rev. A , 42319 (2007), arXiv:0703002v2 [arXiv:cond-mat].[74] A. Gyenis, P. S. Mundada, A. Di Paolo, T. M. Haz-ard, X. You, D. I. Schuster, J. Koch, A. Blais, andA. A. Houck, Experimental realization of an intrin-sically error-protected superconducting qubit (2019),arXiv:1910.07542.[75] M. Abuwasib, P. Krantz, and P. Delsing, Fabricationof large dimension aluminum air-bridges for supercon-ducting quantum circuits, Journal of Vacuum Science& Technology B, Nanotechnology and Microelectronics:Materials, Processing, Measurement, and Phenomena , 031601 (2013).[76] Z. Chen, A. Megrant, J. Kelly, R. Barends,J. Bochmann, Y. Chen, B. Chiaro, A. Dunsworth,E. Jeffrey, J. Y. Mutus, P. J. J. O’Malley, C. Neill,P. Roushan, D. Sank, A. Vainsencher, J. Wenner, T. C.White, A. N. Cleland, and J. M. Martinis, Fabricationand characterization of aluminum airbridges for super-conducting microwave circuits, Applied Physics Letters , 052602 (2014).[77] Z.-C. Jin, H.-T. Wu, H.-F. Yu, and Y. Yu, Fabricationof Al air-bridge on coplanar waveguide, Chinese PhysicsB , 100310 (2018).[78] A. Dunsworth, R. Barends, Y. Chen, Z. Chen,B. Chiaro, A. Fowler, B. Foxen, E. Jeffrey, J. Kelly, P. V.Klimov, E. Lucero, J. Y. Mutus, M. Neeley, C. Neill,C. Quintana, P. Roushan, D. Sank, A. Vainsencher,J. Wenner, T. C. White, H. Neven, J. M. Martinis, andA. Megrant, A method for building low loss multi-layerwiring for superconducting microwave devices, AppliedPhysics Letters , 063502 (2018).[79] H. Mukai, K. Sakata, S. J. Devitt, R. Wang, Y. Zhou,Y. Nakajima, and J. S. Tsai, Pseudo-2D supercon-ducting quantum computing circuit for the surfacecode: the proposal and preliminary tests (2019),arXiv:1902.07911.[80] V. Schmitt, Design, fabrication and test of a four su-perconducting quantum-bit processor , Theses, Universit´ePierre et Marie Curie - Paris VI (2015).[81] S. P. Pedersen, K. S. Christensen, and N. T. Zinner, Na-tive three-body interaction in superconducting circuits,Physical Review Research , 033123 (2019).[82] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando,S. Gustavsson, and W. D. Oliver, A Quantum En-gineer’s Guide to Superconducting Qubits (2019),arXiv:1904.06560.[83] M. A. Rol, F. Battistel, F. K. Malinowski, C. C. Bultink,B. M. Tarasinski, R. Vollmer, N. Haider, N. Muthusub- ramanian, A. Bruno, B. M. Terhal, and L. DiCarlo, Fast,High-Fidelity Conditional-Phase Gate Exploiting Leak-age Interference in Weakly Anharmonic Superconduct-ing Qubits, Physical Review Letters , 120502 (2019).[84] P. Zhao, P. Xu, D. Lan, X. Tan, H. Yu, and Y. Yu, High-contrast ZZ interaction using multi-type superconduct-ing qubits (2020), arXiv:2002.07560.[85] N. J. S. Loft, M. Kjaergaard, L. B. Kristensen, C. K.Andersen, T. W. Larsen, S. Gustavsson, W. D. Oliver,and N. T. Zinner, Quantum interference device for con-trolled two-qubit operations (2018), arXiv:1809.09049.[86] T. Bækkegaard, L. B. Kristensen, N. J. S. Loft, C. K.Andersen, D. Petrosyan, and N. T. Zinner, Realizationof efficient quantum gates with a superconducting qubit-qutrit circuit, Scientific Reports , 13389 (2019).[87] S. E. Rasmussen, K. S. Christensen, and N. T. Zinner,Controllable two-qubit swapping gate using supercon-ducting circuits, Physical Review B , 134508 (2019).[88] R. E. Barfknecht, S. E. Rasmussen, A. Foerster, andN. T. Zinner, Realizing time crystals in discrete quan-tum few-body systems, Physical Review B , 144304(2019).[89] Z. Wang, S. Shankar, Z. K. Minev, P. Campagne-Ibarcq,A. Narla, and M. H. Devoret, Cavity Attenuators forSuperconducting Qubits, Phys. Rev. Applied , 14031(2019).[90] S. Touzard, A. Kou, N. E. Frattini, V. V. Sivak, S. Puri,A. Grimm, L. Frunzio, S. Shankar, and M. H. Devoret,Gated Conditional Displacement Readout of Supercon-ducting Qubits, Phys. Rev. Lett. , 80502 (2019).[91] L. Qin, L. Xu, W. Feng, and X.-Q. Li, Qubit state to-mography in a superconducting circuit via weak mea-surements, New Journal of Physics , 033036 (2017).[92] M. Baur, A. Fedorov, L. Steffen, S. Filipp, M. P.da Silva, and A. Wallraff, Benchmarking a QuantumTeleportation Protocol in Superconducting Circuits Us-ing Tomography and an Entanglement Witness, Physi-cal Review Letters , 040502 (2012).[93] S. Filipp, P. Maurer, P. J. Leek, M. Baur, R. Bianchetti,J. M. Fink, M. G¨oppl, L. Steffen, J. M. Gambetta,A. Blais, and A. Wallraff, Two-Qubit State Tomogra-phy Using a Joint Dispersive Readout, Phys. Rev. Lett. , 200402 (2009).[94] L. Dicarlo, J. M. Chow, J. M. Gambetta, L. S.Bishop, B. R. Johnson, D. I. Schuster, J. Majer,A. Blais, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf,Demonstration of two-qubit algorithms with a super-conducting quantum processor, Nature , 240 (2009),arXiv:0903.2030.[95] M. Steffen, M. Ansmann, R. C. Bialczak, N. Katz,E. Lucero, R. McDermott, M. Neeley, E. M. Weig, a. N.Cleland, and J. M. Martinis, Measurement of the En-tanglement of Two Superconducting Qubits via StateTomography, Science , 1423 (2006).[96] A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller,Measuring entanglement growth in quench dynamics ofbosons in an optical lattice, Physical Review Letters , 1 (2012), arXiv:1205.1521.[97] H. Pichler, L. Bonnes, A. J. Daley, A. M. L¨auchli,and P. Zoller, Thermal versus entanglement en-tropy: A measurement protocol for fermionic atomswith a quantum gas microscope, New Journal ofPhysics , 10.1088/1367-2630/15/6/063003 (2013),arXiv:1302.1187. [98] R. Islam, R. Ma, P. M. Preiss, M. E. Tai, A. Lukin,M. Rispoli, and M. Greiner, Measuring entanglemententropy through the interference of quantum many-body twins (2015), arXiv:1509.01160.[99] N. M. Linke, S. Johri, C. Figgatt, K. A. Lands-man, A. Y. Matsuura, and C. Monroe, Measuring theR´enyi entropy of a two-site Fermi-Hubbard model on atrapped ion quantum computer, Physical Review A ,052334 (2018).[100] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch,C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F.Roos, Probing R´enyi entanglement entropy via random-ized measurements, Science , 260 (2019).[101] J. Heinsoo, C. K. Andersen, A. Remm, S. Krinner,T. Walter, Y. Salath´e, S. Gasparinetti, J.-C. Besse,A. Poto \ ifmmode \ checkc \ else ˇc \ finik, A. Wallraff, andC. Eichler, Rapid High-fidelity Multiplexed Readout ofSuperconducting Qubits, Phys. Rev. Applied , 34040(2018).[102] G. E. Fisher, Lectures in Theoretical Physics , Lecturesin Theoretical Physics: Lectures Delivered at the Sum-mer Institute for Theoretical Physics, University of Col-orado, Boulder, Vol. 7 (University of Colorado Press,Boulder, 1965).[103] C. G. Callan, R. F. Dashen, and D. J. Gross, The struc-ture of the gauge theory vacuum, Physics Letters B ,334 (1976).[104] R. Jackiw and C. Rebbi, Vacuum periodicity in a Yang-Mills quantum theory, Physical Review Letters , 172(1976).[105] S. Coleman, More about the massive Schwinger model,Annals of Physics , 239 (1976).[106] T. Sch¨afer and E. V. Shuryak, Instantons in QCD, Re-views of Modern Physics , 323 (1998).[107] G. Gabadadze and M. Shifman, QCD vacuum andAxions: What’s Happening?, International Journal ofModern Physics A , 3689 (2002), arXiv:0206123v1[arXiv:hep-ph].[108] R. D. Peccei, The strong CP problem and axions, Lec-ture Notes in Physics , 3 (2008), arXiv:0607268v1[arXiv:hep-ph].[109] R. D. Peccei, D. B. Tanner, and K. A. van Bibber, WhyPQ?, AIP Conference Proceedings , 7 (2010).[110] J. E. Kim and G. Carosi, Axions and the strong CPproblem, Reviews of Modern Physics , 557 (2010),arXiv:arXiv:0807.3125v2.[111] S. Pancharatnam, Generalized theory of interferenceand its applications, Proceedings of the Indian Academy of Sciences - Section A , 398 (1956).[112] J. Samuel and R. Bhandari, General Setting for Berry’sPhase, Physical Review Letters , 2339 (1988).[113] T. Fukui, Y. Hatsugai, and H. Suzuki, Chern numbersin discretized Brillouin zone: Efficient method of com-puting (spin) Hall conductances, Journal of the PhysicalSociety of Japan , 1674 (2005). VII. ACKNOWLEDGEMENTS
The authors would like to thank Torsten Zache for in-sightful discussions. The authors would also like to thankK. S. Christensen, N. J. S. Loft, S. E. Rasmussen, L. B.Kristensen and T. Bækkegaard for general discussionspertaining to this work. The authors acknowledge sup-port from the Independent Research Fund Denmark, theCarlsberg Foundation, and AUFF through the Jens Chr.Skou fellowship program.
VIII. AUTHOR CONTRIBUTIONS
S.P.P. and N.T.Z. proposed the system, designed thestudy and methods. S.P.P. conducted the simulationsand analysed the results. S.P.P. wrote the initial draftof the manuscript. S.P.P. and N.T.Z. contributed to theediting towards the final version for publication.
IX. COMPETING INTERESTS
The authors declare no competing interests.
X. ADDITIONAL INFORMATION
Supplementary information is available for this pa-per at
Correspondence and requests for materials shouldbe addressed to S.P.P. upplemental Information: Lattice gauge theory and dynamical quantum phasetransitions using noisy intermediate scale quantum devices
Simon Panyella Pedersen ∗ Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark
N. T. Zinner † Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark andAarhus Institute of Advanced Studies, Aarhus University, DK-8000 Aarhus C, Denmark (Dated: August 21, 2020)
SUPPLEMENTARY METHODS
In the five sections below we give additional details on the following: First a full derivation of the Hamiltonian,anharmonicities and further details of our proposed circuit, using a method for truncation devised by the authors thatdoes not implement a Taylor expansion, yielding more exact spin model parameters. Secondly, we briefly detail howwe numerically found the necessary effective spin model parameters, which take renormalization from interactionswith levels outside the spin-1 / A. Deriving the circuit Hamiltonian
In a superconducting circuit of transmon qubits there is natively two-body interactions, as an interaction is betweennodes and is mediated by a branch, which could of course never connect more than two nodes at a time. To get multi-body interactions we therefore need to change coordinates. Below we change coordinates to the eigenmodes of thecapacitance network, which simultaneously removes all interactions from capacitors (making the system simpler), andallows for multi-body interactions, in particular the three-body XXX-coupling, we are interested in. The intuitionis that these modes are more spread out on the circuit and so more than two at a time can ”be at the same place”and interact. In particular multi-body interactions come from Josephson junctions terms, as these are non-linear.Capacitance and linear inductive terms are square in the coordinates, and so could never couple more than two at atime.The only way (seemingly) to get a direct, triple, odd (i.e. originates from a product of creation/annihilationoperators which all have an odd power) interaction like the XXX-coupling using capacitors, linear inductors, andJosephson junctions in a superconducting circuit, is to shift the position of the minimum of the potential away fromthe origin. Otherwise, all circuit elements contribute with terms that are even in the number of flux coordinates. Thisshifting is done by introducing (constant) drives of the potential terms, such that the different functions and theirminima are shifted relative to each other. Shifting the minimum and expanding around it makes no difference forcapacitive and linear inductive terms. Capacitive terms involve the derivative and so do not care about a constantshift, while linear inductive terms can only give first and second order terms (in the fluxes), and all first order termsdisappear when we expand around an extremal point, like the potential minimum. Hence, we only get new andinteresting terms from the Josephson junction. This in turn give us a myriad of messy terms. Consider the followingcosine term, in the context of a larger system, where the coordinates are chosen such that the minimum of the fullpotential is at the origin, φ = 0, but such that the argument of the cosine is shiftedcos( φ + α ) = cos( α ) − sin( α ) φ −
12 cos( α ) φ + 16 sin( α ) φ + 124 cos( α ) φ + O ( φ ) ∗ [email protected] † [email protected] a r X i v : . [ qu a n t - ph ] A ug φ φ φ φ C E L CLE c E s (cid:12) Φ C L E c E s (cid:12) Φ CLE c E s (cid:12) Φ C E L C L E c E s (cid:12) Φ K KKK
Figure 1. The circuit for implementing the XXX-coupling in a more general form than the used in the main text. We havehere included linear inductors on all branches for the sake of generality.
Here φ represents any combination of fluxes and α is a shift combining different (constant) drives in the cosine, andshifts of the coordinates/potential. We are interested in the φ term for our XXX interaction. We seek a circuit whosesymmetry ensures that all the undesired terms cancel. As we seek a three-body interaction we will work with a circuitof four nodes, all-to-all connected, and in the basis of the eigenvectors of the capacitance matrix, i.e. the eigenmodesof the capacitive network. Working in this basis removes interactions from capacitors, reducing the complexity by alot, and also allows for more than two coordinates to be involved in the same term, yielding multi-body interactions.In Fig. 1 can be seen our circuit in a more general form than was used in the main text. We have here includedlinear inductors on all branches for the sake of studying the most general case of the circuit. Notice the symmetry ofthe circuit parameters. Defining the original node flux coordinates ~φ = ( φ , φ , φ , φ ), the circuit has the followingcapacitance matrix C = C + C + K − C − C − C − C C + C + K − C − C − C − C C + C + K − C − C − C − C C + C + K We define new coordinates ~ψ = ( ψ CM , ψ , ψ g , ψ ), and make a change of coordinates according to ~φ = −
01 0 −
11 0 − − ~ψ Notice that the coordinates ψ and ψ are simply φ − φ and φ − φ , making the terms corresponding to thesebranches affect only these modes individually, while φ g will be affected by the four branches with identical circuitparameters. This transformation diagonalizes the capacitance matrix to K = K C + 4 C + 2 K C + K
00 0 0 4 C + 4 C + 2 K It is important here that the capacitance to ground is identical for all nodes, as it will otherwise introduce interactions.Now let us consider the Josephson junctions on the branches pertaining to the ψ g mode. Taking either of the fourjunctions with subscript c or s , we find that they can be neatly summed as followscos( φ − φ ) + cos( φ − φ ) + cos( φ − φ ) + cos( φ − φ )= cos( − ψ + ψ g + ψ ) + cos( ψ + ψ g + ψ ) + cos( ψ − ψ g + ψ ) + cos( − ψ − ψ g + ψ )= 2 cos( ψ ) cos( ψ g + ψ ) + 2 cos( ψ ) cos( ψ g − ψ )= 4 cos( ψ ) cos( ψ g ) cos( ψ )Hence, these four junctions result in some very neat interactions between the modes and contributions to their energiesand anharmonicities. We choose the spanning tree and the exact layout of the circuit such that the external fluxes,Φ i , i = 1 , , ,
4, enter in the junctions with subscript s . We then choose the Φ i ’s as sums of three other fluxes Ψ j , j = 1 , g,
2, where the sign of each Ψ j in each cosine term must be the same as the corresponding ψ j . For example, ifwe include Φ in the first cosine in the above calculation, cos( φ − φ + Φ ) = cos( − ψ + ψ g + ψ + Φ ), then we cansee that choosing Φ = − Ψ + Ψ g + Ψ results in a Ψ j -term for each ψ j with the same sign. Doing the same with allfour cosines, we essentially shift ψ j → ψ j + Ψ j . The result of the above calculation would then be4 cos( ψ + Ψ ) cos( ψ g + Ψ g ) cos( ψ g + Ψ )Hence, we can use these to tune the interactions pertaining to E s . In particular choosing Ψ j = − π for all j , we geta triple sine interaction. The triple sine will have its minimum at ψ j = 0, as all other terms already did, and so thisshift will only affect this particular term of the Hamiltonian. Expanding to fourth order, it yields exactly a ψ ψ g ψ interaction corresponding to the desired XXX coupling. Likewise for the four linear inductors we find( φ − φ ) + ( φ − φ ) + ( φ − φ ) + ( φ − φ ) = ( − ψ + ψ g + ψ ) + ( ψ + ψ g + ψ ) + ( ψ − ψ g + ψ ) + ( − ψ − ψ g + ψ ) = 4( ψ + ψ g + ψ )Thus they do no not introduce any interaction, as mentioned above, but merely add to the harmonicity and energyof the modes. The cosines pertaining to E c will yield ZZ interaction, which we would rather be without, but theyalone contribute to ψ g ’s anharmonicity. If we instead had only the identical linear inductors, we would get similarcontributions to the mode energies, but no interactions and also no contribution to anharmonicities. In this case ψ g would be a pure harmonic oscillator, but we would have fewer interactions. The other two modes get theiranharmonicities from the E and E junctions. The inductors turn out to not be necessary but could be used to tunethe energies of ψ and ψ if needed. In order for the XXX interaction to work the modes need to be detuned in aproper way, as described below.With these calculations and in these coordinates, the Hamiltonian of the circuit becomes H = 12 14 C + 4 C + 2 K p + 12 14 C + K p + 12 14 C + 4 C + 2 K p + (4 E L + 4 E L ) ψ − E cos(2 ψ ) + 4 E L ψ g + (4 E L + 4 E L ) ψ − E cos(2 ψ ) − E c cos ψ cos ψ g cos ψ − E s sin ψ sin ψ g sin ψ where E L = L and likewise for the others. We will not be performing the usual expansion to fourth order beforerecasting and truncation, but will rather use a method developed by the authoer for in a previous paper [1], withwhich the trigonometric functions can be truncated without performing any approximation. To do so we performthe recasting in terms of harmonic oscillator modes first. We introduce bosonic creation and annihilation operators, a † j , a j , via ψ j = r j √ a † j + a j ) p j = i √ r j ( a † j − a j )Here r j is parameter which quantifies the ”size” of the ψ j coordinate. When performing the usual expansion to fourthorder, one is in fact expanding to fourth order in r j / √
2, and so the usual transmon regime (of circuit parameters) iswhen r j / √ r j are defined as r j = (cid:18) K p j K ψ j (cid:19) / where K p j and K ψ j are the coefficients of the p j and the ψ j terms, respectively, in the Hamiltonian, when it isexpanded in powers of ψ j . If we expand H to just second order to find these terms, we get H = 12 14 C + 4 C + 2 K p + 12 14 C + K p + 12 14 C + 4 C + 2 K p + (4 E L + 4 E L ) ψ − E (cid:0) ψ (cid:1) + 4 E L ψ g + (4 E L + 4 E L ) ψ − E (cid:0) ψ (cid:1) − E c (cid:18) ψ (cid:19) (cid:18) ψ g (cid:19) (cid:18) ψ (cid:19) − E s ψ ψ g ψ + O ( ψ j )= 12 14 C + 4 C + 2 K p + 12 14 C + K p + 12 14 C + 4 C + 2 K p + (4 E L + 4 E L − E − E c ) ψ + (4 E L − E c ) ψ g + (4 E L + 4 E L − E − E c ) ψ + O ( ψ j )where we have ignored constant terms and all terms involving more than two ψ j are hidden in O ( ψ j ). From this wecan see that for example K p = 12 14 C + 4 C + 2 KK φ = 4 E L + 4 E L − E − E c and similarly for the other coordinates. With these we can define r = [8(2 C + 2 C + K ) (2 E L + 2 E L + E + E c )] − / r g = [4(4 C + K )(2 E L + E c )] − / r = [8(2 C + 2 C + K ) (2 E L + 2 E L + E + E c )] − / We will be writing the trigonometric functions in terms of complex exponentials, as this helps us truncate them lateron. Thus, we can write the Hamiltonian in terms of the creation and annihilation operators as H = −
12 14 C + 4 C + 2 K r ( a † − a ) −
12 14 C + K r g ( a † g − a g ) −
12 14 C + 4 C + 2 K r ( a † − a ) + 4 ( E L + E L ) r a † + a ) − E (cid:16) e i √ r ( a † + a ) + e − i √ r ( a † + a ) (cid:17) + 4 E L r g a † g + a g ) + 4 ( E L + E L ) r a † + a ) − E (cid:16) e i √ r ( a † + a ) + e − i √ r ( a † + a ) (cid:17) − E c (cid:16) e ir ( a † + a ) / √ + e − ir ( a † + a ) / √ (cid:17) (cid:16) e ir g ( a † g + a g ) / √ + e − ir g ( a † g + a g ) / √ (cid:17) (cid:16) e ir ( a † + a ) / √ + e − ir ( a † + a ) / √ (cid:17) − i E s (cid:16) e ir ( a † + a ) / √ − e − ir ( a † + a ) / √ (cid:17) (cid:16) e ir g ( a † g + a g ) / √ − e − ir g ( a † g + a g ) / √ (cid:17) (cid:16) e ir ( a † + a ) / √ − e − ir ( a † + a ) / √ (cid:17) In order to truncate this Hamiltonian to any number of desired levels in each mode, we must calculate the correspondingmatrix elements of each operator in the Hamiltonian. The matrix elements are calculated in the Fock basis, and sooperators like ( a † j ± a j ) are easily dealt with. Indeed, it is the trigonometric functions, now in terms of exponentials,that require some work. Here we recognize that the exponentials are on the form of a displacement operator, knownfrom harmonic oscillators, D ( ξ ) = e ξa † − ξ ∗ a Here ξ is some complex number, and the fundamental property of the displacement operator is that D ( ξ ) | i = | ξ i ,where | i is the ground state of the harmonic oscillator, a | i = 0, and | ξ i is a coherent state. The coherent states canin turn be written out as a series, specifically D ( ξ ) | i = e −| ξ | / ∞ X j =0 ξ j √ j ! | j i Furthermore, we know the commutations relations of the displacement operator D ( ξ ) a † = ( a † − ξ ∗ ) D ( ξ )With these relations and the definition | n i = √ n ! ( a † ) n | i , we can calculate the needed matrix elements. Looking atthe Hamiltonian we see that we will only need to consider ξ = ik for some real number k . Let us therefore considerthe general case, and find the matrix elements of D ( ik ) = e ik ( a † + a ) . Let us also note that since the operator a † + a is symmetric, the operator D ( ik ) = e ik ( a † + a ) is also symmetric, i.e. h n | e ik ( a † + a ) | m i = h m | e ik ( a † + a ) | n i . We find withthe above relations h | e ik ( a † + a ) | i = e − k / h | e ik ( a † + a ) | i = ike − k / h | e ik ( a † + a ) | i = (1 − k ) e − k / With these we can write the operator e ik ( a † + a ) truncated to two levels in terms of Pauli matrices and the identity M h e ik ( a † + a ) i = (cid:20)(cid:18) − k (cid:19) + k σ z + ikσ x (cid:21) e − k / where M n [ · ] is the n × n matrix representation of an operator. It is now just a matter of replacing exponentials withthis expression and reducing the result, i.e. we don’t need to perform the truncation calculation more than once, fromnow on it is just a mechanical method of inserting the right expressions into our equations. Likewise we can find thefollowing standard two-level expressions for the remaining operators M [ a † a ] = 12 (1 − σ z ) M [ a † − a ] = − iσ y M [( a † − a ) ] = − σ z M [ a † + a ] = σ x M [( a † + a ) ] = 2 − σ z This procedure is easily expanded to include more levels in the truncation. Performing this exact truncation of thetrigonometric functions can reduce the numerical difficulty of simulating the dynamics of a circuit while includinghigher levels. Such a study is often relevant to check the actual effect of the higher levels of a circuit rather than onlystudying the approximate dynamics of system truncated to two levels in each mode.We can thus calculate the two level Hamiltonian exactly. We do so H = 12 14 C + 4 C + 2 K p + 12 14 C + K p + 12 14 C + 4 C + 2 K p + (cid:18) L + 2 L (cid:19) ψ − E cos(2 ψ ) + 2 L ψ g + (cid:18) L + 2 L (cid:19) ψ − E cos(2 ψ ) − E c cos ψ cos ψ g cos ψ − E s sin ψ sin ψ g sin ψ = −
12 14 C + 4 C + 2 K r ( a † − a ) −
12 14 C + K r g ( a † g − a g ) −
12 14 C + 4 C + 2 K r ( a † − a ) + 4 ( E L + E L ) r a † + a ) − E (cid:16) e i √ r ( a † + a ) + e − i √ r ( a † + a ) (cid:17) + 4 E L r g a † g + a g ) + 4 ( E L + E L ) r a † + a ) − E (cid:16) e i √ r ( a † + a ) + e − i √ r ( a † + a ) (cid:17) − E c (cid:16) e ir ( a † + a ) / √ + e − ir ( a † + a ) / √ (cid:17) (cid:16) e ir g ( a † g + a g ) / √ + e − ir g ( a † g + a g ) / √ (cid:17) (cid:16) e ir ( a † + a ) / √ + e − ir ( a † + a ) / √ (cid:17) − i E s (cid:16) e ir ( a † + a ) / √ − e − ir ( a † + a ) / √ (cid:17) (cid:16) e ir g ( a † g + a g ) / √ − e − ir g ( a † g + a g ) / √ (cid:17) (cid:16) e ir ( a † + a ) / √ − e − ir ( a † + a ) / √ (cid:17) = 12 14 C + 4 C + 2 K r (2 − σ z ) + 12 14 C + K r g (2 − σ zg ) + 12 14 C + 4 C + 2 K r (2 − σ z )+ 2 ( E L + E L ) r (2 − σ z ) − E (cid:2)(cid:0) − r (cid:1) + r σ z (cid:3) e − r + 4 E L r g − σ zg )+ 2 ( E L + E L ) r (2 − σ z ) − E (cid:2)(cid:0) − r (cid:1) + r σ z (cid:3) e − r − E c (cid:20)(cid:18) − r (cid:19) + r σ z (cid:21) " − r g ! + r g σ zg − r (cid:19) + r σ z (cid:21) e − ( r + r g + r ) / − i E s [ √ ir σ x ][ √ ir g σ xg ][ √ ir σ x ] e − ( r + r g + r ) / = − " E + 4 E L + 4 E L + E c ) r + 2 E r e − r + E c r − r g ! (cid:18) − r (cid:19) e − r r g + r σ z − " E L + E c ) r g + E c r g (cid:18) − r (cid:19) (cid:18) − r (cid:19) e − r r g + r σ zg − " E + 4 E L + 4 E L + E c ) r + 2 E r e − r + E c r (cid:18) − r (cid:19) − r g ! e − r r g + r σ z − E c r r g (cid:18) − r (cid:19) e − r r g + r σ z σ zg − E c r r − r g ! e − r r g + r σ z σ z − E c r g r (cid:18) − r (cid:19) e − r r g + r σ zg σ z − E c r r g r e − r r g + r σ z σ zg σ z − √ E s r r g r e − r r g + r σ x σ xg σ x Defining the spin-model parametersΩ = 2( E + 4 E L + 4 E L + E c ) r + 2 E r e − r + E c r − r g ! (cid:18) − r (cid:19) e − r r g + r Ω g = 2(4 E L + E c ) r g + E c r g (cid:18) − r (cid:19) (cid:18) − r (cid:19) e − r r g + r Ω = 2( E + 4 E L + 4 E L + E c ) r + 2 E r e − r + E c r (cid:18) − r (cid:19) − r g ! e − r r g + r J z g = − E c r r g (cid:18) − r (cid:19) e − r r g + r J z = − E c r r − r g ! e − r r g + r J zg = − E c r g r (cid:18) − r (cid:19) e − r r g + r J z g = − E c r r g r e − r r g + r J x g = −√ E s r r g r e − r r g + r we get the final Hamiltonian H = −
12 Ω σ z −
12 Ω g σ zg −
12 Ω σ z + J z g σ z σ zg + J z σ z σ z + J zg σ zg σ z + J z g σ z σ zg σ z + J x g σ x σ xg σ x Furthermore, the anharmonicities can be calculated with the same exact method, by including a third level anddefining the anharmonicities as α = E − E , where E ij is the energy gap between level i and j , and the energiesare defined as the diagonal entries of the Hamiltonian. With this we get the following anharmonicities α = − E r e − r − E c r − r g ! (cid:18) − r (cid:19) e − r r g + r ,α g = − E c r g (cid:18) − r (cid:19) (cid:18) − r (cid:19) e − r r g + r ,α = − E r e − r − E c r (cid:18) − r (cid:19) − r g ! e − r r g + r In the main text we did not include the linear inductors in the circuit. Removing them corresponds simply to setting E L , E L , E L = 0 in the above expressions. B. Derive effective parameters
In the above section we found exact spin model parameters, taking two levels in each mode into account. For thedetuning and XXX-coupling strength, however, we need to numerically find effective parameters, as the above barespin model parameters are renormalized by virtual processes involving the higher levels of the system. We consideronly the effective detuning and XXX-coupling strength this way, because they need to be accurate in order for thecircuit to result in the desired dynamics, while anharmonicities and other parameters do not have such exact demandsimposed on them. The effective detuning we are looking for is the detuning between the states | g i and | g i Time [ns] . . . . . . F i d e li t y Simulated fidelitiesFit yielding J eff = − . × π MHz,and ∆ eff = 12 . × π MHz
Figure 2. with respect to the σ +0 σ + g σ − + H.c. interaction. The effective coupling strength we are looking for is the strength ofthis interaction. To find these we simply simulate the dynamics of the circuit Hamiltonian truncated to four levels ineach mode and initialized in | g i (which is also a valid state, when each mode has four states). We then calculatethe fidelity with the state | g i at subsequent times. Ideally, these two levels behave as if they are only coupledto each other, and as such follow the evolution of a single two level system, with a certain detuning and couplingbetween the states, described by H = − ωσ z + Jσ x = (cid:18) − ω JJ ω (cid:19) (1)Initializing in one of the two states described by this Hamiltonian, the evolution of the fidelity, F ( t ), of the other isthen described by F ( t ) = J ω / J sin r ω J t ! . (2)Hence, fitting the simulated data with a sine curve gives us information about ∆ eff between | g i and | g i , andalso the effective strength, J eff , of the coupling between these states. Fig. 2 shows an example of such a fit yieldingthe effective parameters we discussed in the main text. Clearly the fidelity does follow a sine, as if these two statesare all alone in the world. This is effectively the case because all other interactions are suppressed by detuning, butas explained it is virtual processes resulting from these suppressed interactions that renormalizes the detuning andinteraction strength such that we can not directly use the bare expressions found above. The bare expressions areoff by a significant amount because of this renormalization as can be seen in the example below, Fig. 3. Using thedescribed method here to find the effective parameters, however, makes it just as easy to optimize with respect to theeffective parameters as the bare, though of course it takes more computation time. C. Optimizing circuit parameters and calculating
Higher order contributions from interactions with states outside the spin- subspace mean that we must considereffective spin model parameters when optimizing the circuit parameters. We are in particular interested in thedetuning between | g i and | g i , and the coupling strength of the XXX-coupling. Furthermore, we want thenumerical value of the anharmonicities of the modes to be about 100 × π MHz or larger in order to be able to addressthe spins for initialization and control [2]. The effective detuning, ∆ eff , and the effective coupling strength, J eff , arefound as described in the previous section.We want to show that our circuit can be used to realize the quench dynamics presented in the main text. In thiscase, the effective detuning should not be zero, but rather it defines the staggered mass m . Considering the rotatedHamiltonian H R in Eq. (3) of the main text, it can be seen that we will have m = ∆ eff /
2, as the energy differencebetween the states | g i and | g i is that of two particles, i.e. twice the mass, 2 m = ∆ eff , consistent with ∆ eff being the energy gap between the two states. Likewise, if J is the desired strength of the matter-gauge coupling inEq. (1) of the main text, then J = 2 J eff , because of the factor 1 / J/m = 4 J eff / ∆ eff to find good circuit parameters producing spin model parameters that would be interestingfor the analog simulation of the dynamics we found in the main text. The anharmonicities are obtained with the bareexpressions given above. These too will be renormalized to some extent, but since we do not need them to have somespecific value, they should simply be large, it is sufficient to calculate their bare value.Since only the Josephson energies of a SQC can be tuned in situ , it is difficult to actually perform the appropriatequench of the circuit. Instead we intend for the circuit to be constructed with the post-quench parameters. The quenchwill then be implemented by initializing the system in the ground state of the pre-quench Hamiltonian. Whether wehave the system in its pre-quench setup, go into its ground state, and then quench to the post-setup, or simplystart with the system in the post-quench setup, and then quickly initialize in the ground state of the pre-quenchHamiltonian, we will see the same resulting dynamics. This moves the difficulty from performing a fast quench toperforming a fast initialization. With this in mind we tune the circuit parameters to yield a negative J/m (as thequench is to a negative mass, m → − m ), corresponding to post-quench parameters.In Fig. 3a,b we show a set of circuit parameters, which yield J/m = − . K = 1 fF. As mentioned in the main text, K would ideally be zero, which is why we fixed it at this low value, but an experimental realization where it is completelynegligible would be even better. E s , which controls the strength of J x g , and thus J eff , is likewise fixed at the smallvalue of E s = 0 . × π GHz to ensure that odd interactions with higher levels are suppressed by the anharmonicities.That is, the higher level interactions caused by the triple sine term, which would move population outside of thespin-1 / E s also being smaller,or interactions with higher levels being less suppressed. In Fig. 3c the resulting bare spin model parameters canbe seen. Notice that some of the parameters are negative, despite being plotted with bars pointing in the positivedirection. The bars should be used as graphical comparison of absolute sizes, while the numbers indicate the exactvalues. Clearly the Z-type coupling strengths are much smaller than the spin transition frequencies, which is neededas explained in the main text. In Fig. 3d the effective detuning and XXX-coupling strength, and the anharmonicitiescan be seen. We see that J eff is noticeably smaller than the bare J x g . Furthermore, it is much smaller than theanharmonicities, as required to make the dynamics remain in the spin-1 / J/m = 4 J eff / ∆ eff = − . pre-quench J/m = 2 .
0) which accordingto Fig. 4 of the main text would result in interesting dynamics of the order parameter and Loschmidt amplitudewithin a time of tm = 2, corresponding to t = 49 . J/m considered in the main text. The tuning effectively takes place via two options. Either J eff = J/ E s , or ∆ eff = 2 m is tuned by changing all parameters but E s (which does also change J eff , but onlyslightly). While tuning J eff is easy because it is proportional to E s , ∆ eff is quite sensitive to the circuit parameters.This is essentially because in the parameter regime we are interested in, ∆ eff is about the same size as J eff , whichitself has to be much smaller than the anharmonicities, say a factor of 10, and finally the anharmonicities are about1 −
2% of the spin transition frequencies. Hence, the effective detuning ∆ eff of the spins must be about 0 . − . J/m , other nearby values could be achieved by varying just the Josephson energies, making it possible to usethe same circuit to study different values of
J/m . D. Time evolution in rotating frame
To find an expression for the time evolution operator in a rotating frame, we essentially just have to take advantageof the fact that we can transform between the non-rotating and the rotating frame at any time using the sametime-dependent operator. Imagine some state in the non-rotating frame | ψ ( t ) i = U ( t ) | ψ (0) i , where U ( t ) is the timeevolution operator in the non-rotating frame. We then define a unitary, time-dependent operator R ( t ), which we useto define the state in the rotating frame | ψ ( t ) i R = R ( t ) | ψ ( t ) i , and which satisfies that the two frames are identical at t = 0, i.e. R (0) = ∞ . Let U R ( t ) be the time evolution operator of the rotating frame that we seek. We then quite0 C C C K [fF]66.89 6.13 21.99 1.00 ( a ) E E E c E s [ π GHz]131.60 198.49 51.90 0.50 ( b ) Ω Ω g Ω J z g J z J zg J z g J x g [ π GHz]17.51 18.52 36.04 ( c ) ∆ eff J eff α α g α [ π MHz]12.85 -6.38 -165.96 -206.37 -549.79 ( d ) [ π MHz]-27.54 -10.90 -41.79 -0.26 -15.96
Figure 3. Circuit and spin model parameters. Notice that some parameters are negative, despite being plotted with positivebars. The bars are for graphical comparison of sizes, while the numbers above the bars give the exact value. (a),(b) Optimizedcircuit parameters implementing
J/m = − . K has been fixed at 1 fF as we want it to be as lowas possibly. E s has been fixed at 0 . × π GHz to ensure that the higher level odd interactions of the triple sine-term are weakenough to be suppressed by the anharmonicities. (c) The resulting bare spin model parameters. The Z-type coupling strengthsare much smaller than the spin transition frequencies as required. (d) The effective detuning and XXX-coupling strength, aswell as the anharmonicities. The anharmonicities are much larger than the effective XXX-coupling strength as required, andwe can see that
J/m = 4 J eff / ∆ eff = − . simply have U R ( t ) | ψ (0) i = U R ( t ) | ψ (0) i R = | ψ ( t ) i R = R ( t ) | ψ ( t ) i = R ( t ) U ( t ) | ψ (0) i From this we conclude that U R ( t ) = R ( t ) U ( t ). Hence, if for example we have a Hamiltonian H = H + H int in thenon-rotating frame, i.e. U ( t ) = e − iHt , and we want to look at the frame rotating with respect to H , i.e. R ( t ) = e iH t ,then the time evolution operator in the rotating frame is U R ( t ) = e iH t e − iHt . This is the result we used in the previoussection. E. Calculating the average fidelity
After finding appropriate circuit parameters we do a check of the overall behaviour of the circuit, ensuring thatit works as intended, including no disturbing interactions with higher levels. To do this we use average fidelity [3].1Average fidelity is a measure of how well a certain process implements a desired operation. In our case the processis the time evolution of the circuit, determined by H c in Eq. (6) of the main text, in a frame rotating such that theresulting Hamiltonian is H R in Eq. (3), where m = ∆ eff / J eff (rather than the bare J x g ), as explained. To take contributions from higher level interactions into account we again truncate to the lowest four levels. We only rotate the spin-1 / H + m ( σ − σ ), where all entriespertaining to levels higher than the spin-1 / H from Eq. (1) of the main text with two matter sites and thegauge link between them. We compare only the dynamics of the spin-1 / / H . The average fidelity will thus be comparing the very time evolutionoperators themselves for the circuit and the target system. The mass and coupling strength of the target Hamiltonianare chosen to be ∆ eff / J eff . The fidelity would thus a priori be expected to be quite high, but since this is alldone with four levels included in each anharmonic mode, the fidelity will be a measure of how much the higher levelsaffect the dynamics of the circuit beyond just the renormalization of the mass and coupling strength. In particular,some population will be lost to the higher levels, and just as virtual processes contribute to the strength of the XXX-coupling, they will also to some extent induce other effective interactions. These will disturb the desired dynamics andmight be gauge-variant, resulting in population moving outside of the chosen gauge sector of the spin-1 / F ( U, E ) = P j Tr (cid:16) U U † j U † E ( U J ) (cid:17) + d d ( d + 1) (3)where U is the target operation, and E is a quantum map representing the process used to implement the targetoperation. The U j are a unitary basis for operators on the relevant space of states, which has dimension d . In ourcase the target operation is time evolution U = e − iHt , with H = H (∆ eff / , J eff ) from Eq. (1). This is an operatoron the spin-1 / U j are a unitary basis for operators on that subspace. The quantum map E takes an operator ρ pertaining to the spin-1 / E ( ρ ) = P e iH , eff t e − iH c t P † ρP e iH c t e − iH , eff t P † ,where P is the projection operator from the four-level subspace to the spin-1 / ρ into the four-level space, performs time evolution according to the four-level version of thecircuit Hamiltonian H c , and then projects the result back into the spin-1 / H , eff is an effective four-level equivalent of H , taking energy shifts caused by virtual processes into account. The construction e iH , eff t e − iH c t implements time evolution according to H c in the frame rotating with respect to H , eff by first time evolving forwardswith H c and then backwards with H , eff . See the next section for details on this. We determine H , eff numerically byturning off the XXX-coupling and time evolving each of the eight spin-1 / H (expanded to the four-level Hilbert space). As the Hamiltonian only involves Z-type interactions, and other, verysuppressed interactions caused by the virtual processes, the spin-1 / e − iE eff t , where E eff is their effective energy. With these H , eff can be constructed. [1] S. P. Pedersen, K. S. Christensen, and N. T. Zinner, Physical Review Research , 033123 (2019).[2] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, “A Quantum Engineer’s Guide toSuperconducting Qubits,” (2019), arXiv:1904.06560.[3] M. A. Nielsen, Physics Letters A303