Superconducting Circuit Architecture for Digital-Analog Quantum Computing

We propose a superconducting circuit architecture suitable for digital-analog quantum computing (DAQC) based on an enhanced NISQ family of nearest-neighbor interactions. DAQC makes a smart use of digital steps (single qubit rotations) and analog blocks (parametrized multiqubit operations) to outperform digital quantum computing algorithms. Our design comprises a chain of superconducting charge qubits coupled by superconducting quantum interference devices (SQUIDs). Using magnetic flux control, we can activate/deactivate exchange interactions, double excitation/de-excitations, and others. As a paradigmatic example, we present an efficient simulation of an $\ell\times h$ fermion lattice (with $2<\ell \leq h$), using only $2(2\ell+1)^2+24$ analog blocks. The proposed architecture design is feasible in current experimental setups for quantum computing with superconducting circuits, opening the door to useful quantum advantage with fewer resources.


Introduction
It is known that the calculation of the exact dynamics of a quantum many-body system is in general a challenging task. When the system is complex enough, analytical and numerical solutions are not possible. However, quantum simulation (QS) allows us to overcome this difficulty by using controllable and manipulable quantum systems, known as quantum simulators, to study another non-controllable one [1,2]. QS can be classified into three different groups: analog quantum simulations (AQS), digital quantum simulations (DQS), and digital-analog quantum simulations (DAQS). Formally, DQS and Digital-Quantum Computing (DQC) are equivalent, and we will claim the same for DAQS and Digital-Analog Quantum Computing (DAQC).
In AQS, we may reproduce a given target Hamiltonian for some parameter regimes of the simulated model and of the quantum simulator, which is not universal [3,4]. In DQS, we can perform a sequence of quantum gates in the quantum simulator, which is similar to what happens in DQC. In this case, we can approximate any unitary evolution, specifically, the unitary evolution of the target Hamiltonian model or the given quantum algorithm [5]. Even if DQC is universal, it is less accurate than AQS and Analog Quantum Computing, requiring quantum error correction to scale up and being impractical for current platforms. On the other hand, the recently proposed DAQS aims at getting the best of analog and digital paradigms, with more versatility in the target models or algorithms, higher accuracy, efficient administration of coherence time, and more suitable for current noise intermediate scale quantum (NISQ) architectures [6]. In DAQS, we use a continuous set of complex many-body interactions as a resource (analog blocks), offered naturally by the architecture possibilities of the quantum platform used as quantum computer. We complement it with arXiv:2103.15696v2 [quant-ph] 14 May 2022 accessible continuous sets of single-qubit operations (digital steps), providing versatility of target Hamiltonian models or algorithms [7][8][9][10].
Several physical platforms have been used as quantum simulators, such as optical lattices [11], trapped ions [12] and superconducting circuits [13][14][15]. The latter has particular features allowing for current scalability and design flexibility, and produced the recent claim of quantum supremacy in quantum computing [16,17]. At the same time, diverse target models have been proposed and implemented, principally for DQS and AQS, such as quantum chemistry systems [18][19][20][21], high-energy physics [22][23][24][25], and condensed matter physics [26][27][28]. In the long run to reach fault-tolerant quantum computers in the future, proposals for DAQS appear promising for this and next generation of co-design quantum computers [29][30][31][32]. To achieve that, it would be desirable to widely enhance the DAQC algorithmic mappings as well as the variety of accessible quantum computer geometries and topologies.
In this work, we propose a superconducting circuit design for DAQC. It is composed of a chain of charge qubits coupled through grounded superconducting quantum interference devices (SQUIDs), where the nearest-neighbor qubits are off resonance (different energy gap), similar ideas have been proposed using variable electrostatic fields [33,34]. In this proposal, the SQUIDs modify the resonant condition among nearest-neighbor qubits, which allows us to produce different and independent interactions like an exchange or double excitation/de-excitation term. Also, as the SQUIDs are physically distant, we could manipulate them individually, activating/deactivating several interactions for obtaining parametrized multiqubit gates, also proving a large family of analog multibody Hamiltonians, being suitable for efficient implementations of DAQC protocols. Finally, we test our architecture performing the simulation of the Fermi-Hubbard model, where we need only 2(2 + 1) 2 + 24 analog blocks for the simulation of a × h ( ≤ h) Fermi-Hubbard lattice.

The model
2.1 Two-qubit model First, let us consider the two-qubit system showed in Fig. 1. It consists of two charge qubits coupled through a grounded symmetric SQUID. The Lagrangian of the circuit is given by where f = (d/dt) f (t), represent the time derivate of a function f , E eff J s = 2E J s cos (ϕ ext ) is the effective Josephson energy of the SQUID, and ϕ j = 2πΦ j /Φ 0 is the the superconducting phase, with the superconducting flux quantum Φ 0 = h/2e, and 2e is the electrical charge of a Cooper pair. Moreover Φ 1 , Φ 2 , and Φ s are node fluxes defined in Fig. 1. The effective inductance of the SQUID L J s (ϕ ext ) = (Φ 0 /2π) 2 /E eff J s can be tuned by the external magnetic flux ϕ ext (t), providing a tunable boundary condition [35,36].
By applying the Legendre transformation, we obtain the Hamiltonian (see appendix A) Qubit 1 Qubit 2 SQUID Figure 1 Circuit diagram of two charged qubits (green) coupled through a grounded SQUID (blue). E J 1 (2) , C J 1 (2) , C g 1 (2) , and V g 1 (2) are the Josephson energy, Josephson capacitance, gate capacitance, and gate voltage of the qubit 1 (2) respectively. E Js and C s are the Josephson energy and effective capacitance of the SQUID. Moreover, C c is the coupling capacitance, and Φ 1 , Φ 2 , and Φ s are node fluxes that define the degrees of freedom of the circuit. with H coupling = g 12 Q 1 Q 2 + g 1s Q 1 Q s + g 2s Q 2 Q s , where Q j = ∂L/∂Φ j is the charge (conjugate momenta) of the jth node given by Q 1(2) = (C 1(2) + C c )Φ 1(2) − C g 1(2) V g 1(2) − C c Φ s , Q s = (C s + 2C c )Φ s − C c (Φ 1 + Φ 2 ), (4) and the effective Josephson capacitances are defined as with ). Moreover, the gate-charge numbers read and couplings strengths are given by

Figure 3
Ratio between the SQUID impedance Z s and the qubit 1(2) impedance Z 1(2) as a function of the external magnetic flux ϕ ext .
Here, we consider the regime of high plasma frequency for the SQUID, where the charge energy is small compared to the Josephson energy, and the plasma frequency of the SQUID is far exceeding the frequency of the qubits (see Fig. 2), then we can consider Φ s Φ 1 (Φ 2 ) and Φ s Φ 1 (Φ 2 ) [35]. In addition, we also consider the low impedance for the SQUID (see Fig. 3), which allow us consider Φ s Φ 1 (Φ 2 ). Based on the above conditions, we obtain the next relation for Q s (see appendix A) Now, using the Euler-Lagrange equations we obtain (see appendix A) using the same above conditions we get the relation for ϕ s as where we approximate sin (ϕ s ) ≈ ϕ s .
We note that, we can write the charge in the node j as Q j = 2en j . Promoting the classical variables {n j , ϕ j } to quantum operators {n j ,φ j } with the commutation relation [e iφ j ,n j ] = e iφ j [37], and applying Eqs. (8) and (10) to Eq. (2), we obtain the quantum mechanical Hamiltonian describing our circuit aŝ where the effective coupling strength andĤ j sub is the Hamiltonian of the jth subsystem, given bŷ with E C j = e 2 /2(C j + C c ),n g j = −C g j V g j /2e, and In the following discussion, we considern g 1 =n g 2 = 0.5, and = 1. It is convenient to write the circuit Hamiltonian in the charge basis, it meansn i = n j n j |n j n j | and cos (φ j ) = 1/2( n j |n j n j + 1| + n j |n j + 1 n j |) [37]. Due to the anharmonicity of H j sub (see appendix A), we can perform the two-level approximation in order to obtain the effective Hamiltonian where ω 1 = E J 1 , ω 2 = E J 2 and σ α j is Pauli matrix element of the jth charge qubit and I is the identity operator. Now, we will consider the external flux ϕ ext to be composed by a DC signal and a small AC signal as with |A 1 |, |A 2 | |ϕ DC |, with which we can approximate whereĒ J s = 2E J s cos (ϕ DC ). Then, we can rewrite the Hamiltonian in Eq. (15) aŝ where Now, we write the Hamiltonian of Eq. (18) in the interaction picture with respect tô H 0 = 2 i= j ω j σ z j /2 and perform the rotating wave approximation (RWA), obtaininĝ Here, we make use of ∆ 12 = ω 1 −ω 2 , µ 12 = ω 1 +ω 2 and we neglect the fast oscillating terms proportional to exp(±i(∆ 12 + ν 1(2) )t), exp(±i(µ 12 + ν 1(2) )t), exp(±i∆ 12 t), and exp(±iµ 12 t). As the qubits are far from resonance and considering {g 0 , A 1 g 1 /2, A 2 g 1 /2} {∆ 12 , µ 12 , ν 1 , ν 2 }, the RWA is justified (for more details see appendix A). Considering ν 1 = ∆ 12 and ν 2 = µ 12 , the Hamiltonian in Eq. (20) turnŝ where we neglect the fast oscillating terms proportional to exp(±i(∆ 12 − ν 2 )t) and exp(±i(µ 12 − ν 1 )t). We recall that, for a proper choice of the phasesφ 1 andφ 2 in Eq. (21), we can engineer different interactions as those in Tab. 1.

Three-qubit model
In the three-qubit model, we consider the circuit given by Fig. 4. It is composed of a chain of three charge qubits coupled through grounded SQUIDs. As in the previous case, we consider far off-resonance nearest-neighbor qubits. Following the same procedure of the two-qubits model, we get the next effective Hamiltonian (see appendix B) where ω 3 = ω 1 = E J 1 , ω 2 = E J 2 , and the time-dependent signal reads Moreover the coupling strength g ( j) 0 and g ( j) 1 ( j = {1, 2}) are given by To visualize the dynamics of the system, we write the Hamiltonian in the interaction picture. After we consider the resonant conditions ν (1) 1 = ν (2) 1 = ∆ 12 and ν (1) 2 = ν (2) 2 = µ 12 and neglect the fast oscillating terms, the Hamiltonian in the interaction picture readŝ is the interaction Hamiltonian between jth and ( j + 1)th qubit. By choosing proper phase parameters, we can engineer different interaction operators between adjacent qubits, like in the previous case. It is possible to generalize this expression for a chain of qubits coupled through grounded SQUIDs (see Fig. 5), where we define the qubits in odd positions as qubit 1 with frequency ω 1 and the qubits in even positions as qubit 2 with frequency ω 2 . In the following discussion, we consider the amplitude of the two harmonic signals to be the same, that is A ( j) 1 = A ( j) 2 = A and the coupling strength g ( j) 0 = g 0 and g ( j) 1 = g 1 . By considering the resonant conditions ν ( j) 1 = ∆ 12 , ν ( j) 2 = µ 12 and choosing proper phase parametersφ ( j) 1 andφ ( j) 2 , we can engineer again a family of interactions between nearestneighbor qubits as is shown in Tab. 2. Note that the phaseφ j 1 required to achieve ±σ x j σ y j+1 and ±σ y j σ x j+1 are different for odd and even j.
The controllability and flexibility of the interactions that our proposal offers, give us the possibility to the implement of a large variety of Hamiltonians in an analog way, such as Dzyaloshinskii-Moriya, XY, homogeneous and inhomogeneous spin chains. Such analog Hamiltonians could be very useful for DAQS and DAQC, where we can use such analog Hamiltonians like a resource (complex multibody gate) for simulating more complex systems, like quantum chemistry physics, condensed matter phenomena in spin lattices [38,39], and shortcuts to adiabaticity in digitezed adiabatic quantum computing [40].
The approach we presented in this work is intimately linked to the nature of the Jordan-Wigner mapping that requires the quantum simulation algorithm follows a linear sorting of the lattice sites to reproduce the fermionic anti-commutation relation, being a qubit-chain a natural simulator of the fermion models. Naturally, the search for an experimentally feasible fermion to qubit mapping approaching a two-dimensional lattice beyond the Jordan-Wigner transformation is an open question that deserves further investigation. There is a recent work that proposed a novel mapping in this direction, but the experimental realization is still an open question [41]. Also, the use of qubit lattices for the Jordan-Wigner transformation has been proposed; nevertheless, the gates number scaling is the same as that in the qubit-chain case [8].
The current proposal could be extended to a two-dimensional (2D) array of charged qubits coupled through grounded SQUIDs generating a more complex family of Hamiltonians, opening the door to more efficient simulations. Nevertheless, for 2D structures, we can have the non-trivial problem of cross-talk between the different SQUID and loops in the circuit, requiring a deep feasibility study which is not the scope of this article. In the next section, we show a particular example about the efficient DAQC of a complex system, the Fermi-Hubbard model. This example will illustrate the versatility of our design.

Digital-Analog Quantum Computation
Hubbard model represents the interaction between the neighboring sites, which is defined by hopping element and Coulombic interaction on the same site, called on-site interaction [42]. In this section, we are interested in the simulation of the hopping terms of a × h fermion-lattice (with ≤ h). The Hamiltonian of a × h fermion-lattice (see Fig. 6 (a)) reads where c † j,α (c j,α ) are the creation (annihilation) operators of the jth site, with the number operator n j,↑(↓) = c † j,↑(↓) c j,↑(↓) , and spin-α, with α = {↑, ↓}. To suppress the index α, we map the × h lattice to a equivalent 2 × h spin-less lattice by where b † k are the creation operation over the site k for the lattice given by Fig. 6 (b). Using these operators the Hamiltonian Eq. (27) can be rewrite as Finally, we map the 2 × h fermion lattice to a spin−1/2 chain using the Wigner-Jordan transformation (see Fig. 6 (c)), where we represent b j and b † j as a combination of Pauli matrices Before to write the equivalent spin chain Hamiltonian, we define the operator where j k, and U α j = e −i π 4 σ α j σ α j+1 . After some algebraic manipulation, we obtain (for details see appendix C) where H hori , H verti and H coul are given by which correspond to the Hamiltonian of the horizontal hopping, vertical hopping and coulomb interaction respectively, with and We notice that each term of H hori in Eq. (33) correspond to the set of horizontal hopping interactions which do not share sites in the lattice, these eight interactions are represented by arrows with different colors in Fig. 7 (blue, red, green and brown), and different textures (solid and dashed). Also, we highlight that each interactions given by Eq. (34) can be performed in an analog way since the sub-gates involved are applied in different qubits and can be done together without interfering with each other (see appendix C), giving us the analog resource for the digital-analog simulation. The number m 1 is the number of the hopping terms corresponding to the blue(red) solid/dash arrows, and m 2 is the number of the hopping terms corresponding to the green(brown) solid/dash arrows (see Fig. 7 Now, we approximate the time evolution of our system using the Trotter expansion [43] as follows where U word (t) = e −iH word t . Using Eq. (33) we have First, for U hori we simulate eight types of interactions (see Fig. 7 and Eq. (33)), each of these interactions need three gates, as we mention above each of this gate can be simulated in an analog way. Therefore, to simulate U hori (the horizontal hopping) we need 24 gates, as we can see in Eq. (34).
Second, for U verti , i.e. all the vertical hopping terms, we need 2(2 + 1) types of interactions, they are shown in Fig. 8 with different colors and different textures. We note that all the interactions of the same type can be performed at the same time, due to they do not share sites during its implementation, i.e. all the interactions with the same color and same texture can be simulated in parallel. Now, to simulate each type of interaction we need 2 + 1 analog gates (see Eqs. (33) and (34)). Then, to simulate all the vertical hopping terms, we need 2(2 + 1) 2 analog gates. Therefore to simulate all the hopping terms in a × h fermion lattice with ≤ h we need only 2(2 + 1) 2 + 24 = 8 2 + 8 + 26 gates. It is necessary to highlight that to our knowledge, the more efficient proposal for the quantum simulation of the Hubbard model is a trapped ion one [44], which uses a multibody entangling gate, needing 8(2h − − h) + 20 gates for a × h fermion lattice. It means that for a square lattice ( = h), the trapped ion proposal needs 2(8 2 − 8 + 10) gates, which for large need almost two times more gates than our proposal. It means that even if superconducting circuits cannot use multi-body entangling gates like trapped ions, it still is useful by a suitable design as in this work.
On the other hand, U coul (t/n), correspond to the free energy of the original Hubbard model (Eq. (27)) and we will not consider for the Hopping dynamics simulation. Nevertheless, it can be simulated using three gates, i.e. the analog interaction j σ x 2 j−1 σ x 2 j plus two rotations in the y-axis (the local terms σ j z can be mapped correspond to the free energy of our simulator).
The simulation time can be easily derivate as follow. Each type of interaction involve unitary gate of the form U a = exp(− iA 2 t nÔ ) and of the form U b = exp(−i π 4Ô ). From Eq. (37) we obtain that the time for each kind of gates (a and b) is respectively. From the simulation of U hori we have 8 gates of the class a, and from U verti we have 2(2 + 1) gates of the class a. Then the time necessary to perform all these gates is (4 + 10)τ a . As we have a total of 2(2 + 1) 2 + 24 gates, the number of type b gates is 8 2 + 4 + 16. Therefore, the total time for the simulation is We note that for the case = 2, we need fewer gates, in particular, to simulate U hori we need half of the gates, it means 12, for this case, the simulation time also decreases and is given by Finally, the character digital-analog of our simulation is given by the use of analog gates in each digital step, it means gates that act over several qubits simultaneously. In the next section, we present the numerical results of a quantum simulation of the hopping interaction of a 2 × 3 fermion-lattice.

Numerical results: 2 × 3 fermion lattice
As we mention above, for the case of = 2 we only need 12 gates to simulate U hori , then for 2 × h lattice, we need 62 gates per Trotter step  Figure 9 shows the types of hopping interactions to simulate for a 2 × 3 fermion-lattice. In Fig. 9 (a), we can see the four interactions to describe the horizontal hopping, where the solid arrows correspond to the forward hopping, and the dashed arrows correspond to the backward hopping. For vertical hopping, it requires ten types of interactions, as is shown in Fig. 9 (b). The sequence of the gates for different hopping interactions is shown in appendix D and appendix E.
For the simulation, we map the 2×3 fermion lattice into a 12 qubit chain described by the Hamiltonian Eq. (26). The parameters that we consider for the simulation are summarized in the table 3, where t represent the simulated time (evolution time of the system to be simulated). Figure 10 show the fidelity | ψ sim (t)|ψ(t) | 2 of our simulation for different initial states and 10, 20 and 30 Trotter steps, where |ψ (t) is the state at time t of the real model, and |ψ sim (t) is the state given by the simulation, which simulate the evolution at a time t. If we think the fermion lattice as a 3×2 matrix, where each element can be ↑, ↓ or vacuum, the initials states are: Fig. 10 (a), ↑ for the sites (1, 1), (2,2) and (3,1), and the rest in vacuum; Fig. 10 (b), ↑ for the sites (1, 1) and the rest in vacuum; Fig. 10 (c), quantum superposition between the state ↑↓ for all sites and vacuum for all sites. Finally Fig. 10 (d), shows the mean fidelity over 1000 random states. cQED model 0.05(GHz) Ag 1 /2π 0.08(GHz)

Conclusion
We have designed a superconducting circuit architecture suitable for DAQC, in the sense of providing a wide family of analog Hamiltonians as source of analog blocks and more flexibility for this pragmatic quantum computing paradigm. We test our design by the numerical calculation of the quantum simulation of a 3 × 2 fermion lattice described by the Hubbard model. We find that for a × h lattice ( ≤ h), we need 2(2 + 1) 2 + 24 analog blocks, which depend only on one of the dimensions of the lattice and improves. To our knowledge, this result would improve previous achievements for simulating the Hubbard model by a factor 2 for large square lattice ( = h) [44]. Moreover, the total simulation time for 30 Trotter steps is less than 0.2 µs with an ideal fidelity around 0.97 (only digital error), which makes our proposal experimentally feasible. Finally, we consider this work provides an important boost to the DAQC paradigm, paving the way for computing and simulating complex systems in quantum platforms, while approaching us to useful quantum advantage with fewer algorithmic and hardware resources.

Appendix Appendix A: Simplified two-qubit Hamiltonian
For completeness, we describe the derivation of the simplified Hamiltonian of the two-qubit system, which consists of two charge qubits coupled through a grounded SQUID as shown in Fig. 1 of the main text. The Lagrangian of the circuit reads where ϕ j = 2πΦ j /Φ 0 , with Φ 0 = h/2e is the superconducting flux quantum and 2e is the electrical charge of a Cooper pair. Moreover E eff J s = 2E J s cos (ϕ ext ) is the effective Josephson energy of the SQUID. Now, we calculate the conjugate momenta (node charge) By applying the Legendre transformation H(Φ k , Q k ) = k Q k Φ k − L, we obtain the Hamiltonian where the effective Josephson capacitances, the gate-charge numbers and the coupling strength are defined as follows with C j = C g j +C J j ( j = {1, 2}), and C 3 = C c (C 1 +C 2 )(C s +C c )+C 2 c C s +C 1 C 2 (2C c +C s ). In the following discussion, we calculate the simplified Hamiltonian between the two CPBs by applying the two approximations Φ s Φ 1(2) (Φ s Φ 1(2) ), and Φ s Φ 1 (2) i.e. we consider the SQUID in phase regime with high plasma frequency, and meanwhile the low impedance. With the first approximation Φ s Φ 1(2) , we can neglect the terms proportional to Φ s in Eq. (42) obtaining the relation between nodes charge as follows Next, we derive the relation between the nodes flux, by calculating the Euler-Lagrange (E-L) equations ∂L/∂Φ j − d(∂L/∂Φ j )/dt = 0, which govern the dynamics of our system and by applying the two approximations, we can neglect the terms proportional to Φ s in Eq. (46), and meanwhile approximate sin ϕ s = ϕ s obtaining Moreover, with the approximation Φ s Φ 1(2) , we approximate cos ϕ s ≈ (1 − ϕ 2 s /2), with which we can keep the potential energy of the SQUID up to the second-order, and by replacing Eq. (45), and Eq. (47) in the Hamiltonian in Eq. (43), we obtain where the effective Josephson capacitance and gate-charge number which shows that the replacement of the Q s in terms of Q 1 (2) in the Hamiltonian in Eq. (43) corrects the effective Josephson capacitances and gate-charge numbers of the simplified model. Moreover, we define which depend on the external flux ϕ ext . Now by promoting the classical variables to quantum operators i.e. Q j →Q j = 2en j and ϕ j →φ j with the commutation relation [e iφ j ,n j ] = e iφ j , we obtain the following quantum Hamiltonian is the Hamiltonian of the jth subsystem with the charge energy E C j = e 2 /2C J i . In the following discussion, we considern g 1 =n g 2 = 0.5 and = 1. Note that the free Hamiltonian of the subsystem H j sub in Eq. (52) includes both the bare CPB Hamiltonian and the nonlinear term proportional to sin (φ j ) 2 . To study how this extra term affects the anharmonicity of the subsystem, we plot the energy spectrum of the free HamiltonianĤ j sub as a function of offset chargen g j in Fig. 11 for different E J j /E C j and γ j (ϕ ext )/E C j . It shows that as long as we keep E J j /E C j in charge regime, the increase of γ j (ϕ ext ) will not destroy the anharmonicity of the subsystem, and level of anharmonicity still depends on ratio E J j /E C j . Thus, in the following discussion, we assume both CPBs working in charge regime, and write the Hamiltonian in Eq. (48) in the number basis withn j = m m|m j m j |, cos (φ j ) = 1 2 ( m |m j m j +1|+H.C) and sin (φ j ) = − i 2 ( m |m j m j + 1| − H.C), where |m j is the mth exited state of the jth subsystem. After we perform the two-level approximation, the operator sin (φ j ) and the nonlinear term γ j (ϕ ext ) sin (φ j ) 2 in the subsystem basis read where σ α j is the Pauli matrix, I j is the identity operator. Thus, we can neglect the term proportional to sin (φ j ) 2 , as it only provides a shift to the qubit frequency. Finally, we obtain the simplified Hamiltonian as followŝ Next, we consider the external flux ϕ ext to be composed of a DC signal and a small AC signal as ϕ ext = ϕ ext (t) = ϕ DC + ϕ AC (t), where ϕ AC (t) = A 1 cos (ν 1 t +φ 1 ) + A 2 cos (ν 2 t +φ 2 ) and |A 1 |, |A 2 | |ϕ DC |, with which we can approximate Here,Ē J s = 2E J s cos ϕ DC . By replacing Eq. (55) in the Hamiltonian Eq. (54), we obtain where the coupling strength To visualize the dynamics of our system, we go to the interaction picture characterized by the free HamiltonianĤ 0 = ω 1 2 σ z 1 + ω 2 2 σ z 2 and perform the rotating wave approximation (RWA) obtaininĝ where ∆ 12 = ω 1 −ω 2 , µ 12 = ω 1 +ω 2 , and we neglected the fast oscillating terms proportional to exp(±i(∆ 12 + ν 1(2) )t), exp(±i(µ 12 + ν 1(2) )t), exp(±i∆ 12 t), and exp(±iµ 12 t), since we consider the qubits are far from resonance and the coupling strength {g 0 , g 1 A 1 /2, g 1 A 2 /2} {∆ 12 , µ 12 , ν 1 , ν 2 }. Next, we assume that ν 1 = ∆ 12 and ν 2 = µ 12 , with which we can perform the second RWA neglecting the fast oscillating terms proportional to exp(±i(∆ 12 − ν 2 )t) and exp(±i(µ 12 − ν 1 )t) in Eq. (58) obtaininĝ To prove the justification of the RWA we applied, we plot the population inversion between the states |01 and |10 in Fig. 12. It shows that the numerical results calculated from the Hamiltonian in Eq. (59) (red line and orange line) coincide well with the results calculated Eq. (56) (green line and blue line), which proves the validity of the approximation we applied. Finally, by replacing σ + j = (σ x j + iσ y j )/2 and σ − j = (σ x j − iσ y j )/2 in Eq. (59), we obtain the interaction Hamiltonian in terms of Pauli matrices as followŝ and the interactions we can engineer with different phaseφ 1 andφ 2 are shown in the Tab. 1 of the main text.

Appendix B: Simplified three-qubit Hamiltonian
In this section, we derive the simplified Hamiltonian of the three-qubit system, where the three charge qubits coupled with each other through the grounded SQUIDs as shown in Fig. 4 of the main text. The Lagrangian of the circuit reads where E eff J s j = 2E J s cos ϕ ( j) ext is the effective Josephson energy of the jth SQUID and we assume that the third CPB to be the same as the first one i.e. C g 3 = C g 1 , C J 3 = C J 1 , and V g 3 = V g 1 . Now we calculate the conjugate momenta (node charge) Q j = ∂L/∂Φ j By applying the Legendre transformation H(Φ k , Q k ) = k Q k Φ k − L, we obtain the Hamiltonian where the effective Josephson capacitances , the coupling strength and the gate-charge numbers with ], and C j = C J j + C g j ( j = {1, 2}). Next, we calculate the effective Hamiltonian of this three-qubit model, by applying the approximations Φ s 1(2) and Φ s 1(2) Φ 1 (2,3) , where we assume the both SQUIDs in phase regime with high plasma frequency and low impedance. With the first approximation, we can neglect the terms proportional to Φ s 1 (2) in Eq. (62) obtaining the relation between nodes charge as follows Moreover, to obtain the relation between nodes flux, we calculate the E-L equations, and by applying the two approximations, we can neglect the terms proportional to Φ s 1 (2) and approximate sin (ϕ s 1(2) ) = ϕ s 1(2) obtaining Meanwhile, with the condition Φ s 1(2) Φ 1(2,3) , we can approximate cos ϕ s 1(2) ≈ (1−ϕ 2 s 1(2) /2), where we can keep the potential energy of the SQUID up to the second-order and by replacing Eq. (67) and Eq. (68) in the Hamiltonian in Eq. (63), we obtain By promoting the classical variables to quantum operators, i.e. Q j →Q j = 2en j and ϕ j →φ j with the commutation relation [e iφ j ,n j ] = e iφ j in Eq. (69), we obtain the quantum Hamiltonian where the Hamiltonian of the subsystem readŝ and the charge energy E C j = e 2 /(2C J j ). In the following discussion, we considern g 1 = n g 2 = 0.5 and = 1. As mentioned in appendix A, the term proportional to (sinφ j ) 2 in the subsystem Hamiltonian in Eq. (72) does not destroy the anharmonicity of the system. Thus, in charge regime, we can safely perform the two-level approximation and write the Hamiltonian in Eq. (71) in the subsystem basis, where the operator sin (φ j ) = σ y j /2, and the nonlinear term proportional to sin (φ j ) 2 can be regarded as a shift to the qubit frequency obtaininĝ where ω 1 = E J 1 , ω 2 = E J 2 , and the coupling strength γ 12 ϕ (1) ext /4, γ 23 ϕ (2) ext /4 depend on the external flux through the first, and the second SQUID, respectively. Here we consider the external flux ϕ ( j) ext ( j = {1, 2}) to be composed of a DC signal and a small AC signal as ϕ ( j) whereĒ ( j) J s = 2E J s cos ϕ ( j) DC . By replacing Eq. (74) in the Hamiltonian in Eq. (73), we obtainĤ where the coupling strength To visualize the dynamics of the system, we go to the interaction picture characterized by the free HamiltonianĤ 0 = ω 1 σ z 1 /2 + ω 2 σ z 2 /2 + ω 1 σ z 3 /2. Moreover, we consider the resonant conditions ν ( j) 1 = ∆ 12 = ω 1 − ω 2 and ν ( j) 2 = µ 12 = ω 1 + ω 2 and perform the RWA obtaininĝ where we neglected the fast oscillating terms proportional to exp ±i(∆ 12 + ν ( j) 1(2) )t , exp ±i(µ 12 + ν ( j) 1(2) )t , exp(±i∆ 12 t), exp(±iµ 12 t), exp ±i(∆ 12 − ν ( j) 2 )t , and exp ±i(µ 12 − ν ( j) 1 )t , as we consider the qubits are far from resonance, and the coupling strength Notice that the phaseφ ( j) 1 required to achieve ±σ x j σ y j+1 and ±σ y j σ x j+1 is different for odd and even j, and the interactions we can engineer are shown in Tab. 2.

Appendix C: Mapping Fermion Hubbard model to spin model
For completeness, we derive the Hamiltonian of an h × fermion-lattice (see Fig. 6 (a)) in terms of spin operators by applying the Wigner-Jordan transformation. The Hamiltonian can be expressed as where A is the kinetic energy, B is the on-site repulsion, c † j,α (c j,α ) are the creation (annihilation) operators that act over the jth site, n j,↑(↓) = c † j,↑(↓) c j,↑(↓) is the number operator, and α =↑, ↓ is the spin component. To suppress the index α, we map this lattice to an equivalent 2 × h lattice as shown in Fig. 6 are the creation (annihilation) operation over the site k for the lattice. Now the Hubbard Hamiltonian can be written in terms of where the three terms correspond to the horizontal hopping Hamiltonian, the vertical hopping Hamiltonian, and the Coulomb interaction, respectively. By applying the Wigner-Jordan transformation, we map the operator b † j (b j ) to the combination of Pauli matrices as follows where σ k j is the k-Pauli-matrix associated with the spin−1/2 of the jth position of the chain as shown in Fig. 6 (c). Using the Eq. (82) and with k > j, we obtain where Z k j = ⊗ k = j σ z . And we can see from Fig. 14 that k − j is an even number for all j and k involved in the interaction, therefore Now, we simulate this multi-body interaction using only two-body gates. Before everything, we note that , which allow us construct all multi-body operators of the form given by Eq. (84). In the following discussion, we derive the Hamiltonian corresponding to the horizontal hopping H hori , vertical hopping H verti , and Coulomb interaction H coul in terms of Pauli matrices, respectively.

C.1 Horizontal hopping
We first calculate the horizontal hopping Hamiltonian, which involves eight types of interactions. And all these interactions can be divided into two directions (see Fig. 7), i.e. the forward hopping (solid arrows) b † j b j+2 and the backward hopping (dashed arrows) b † j+2 b j , where we can write both of them in terms of Pauli matrices with Eq. (84) and Eq. (85a), Now we calculate the Hamiltonian corresponding to the blue solid arrows (see Fig. 7), which contains the hopping terms where p = 2 − 4 − (−1) +1 . By replacing Eq. (86) in Eq. (87), we obtain where with m 1 = 2 − 1 − (−1) +1 /4, and m 2 = 2 − 3 + (−1) +1 /4, which corresponds to the number of the hopping terms of the blue(red) solid/dashed arrows, and green(brown) solid/dashed arrows respectively (see Fig. 7), where m 1 + m 2 = − 1. As for the blue dashed arrows, we have Thus for the horizontal hopping corresponding to the blue arrows, we have Following the previous procedure, we obtain the Hamiltonian corresponding to the red arrows, green arrows, and brown arrows (see Fig. 7) as follows 2,4) , where the first term and the second term correspond to the forward hopping and the backward hopping, respectively. Finally, the Hamiltonian of the horizontal hopping in terms of Pauli matrices read where the first term corresponding to downward hopping (dashed arrows), and the second term corresponding to upward hopping (solid arrows) as shown in Fig. 8. To construct the term σ x j Z 2 + j−1 j+1 σ x 2 + j in Eq. (94), we insert Eq. (85b) on both sides of σ x j+ −1 σ y j+ as follows And for σ y j Z 2 + j−1 j+1 σ y 2 + j , it can be written in terms of By replacing Eq. (95) and Eq. (96) in Eq. (94), we obtain the vertical hopping Hamiltonian between jth and ( j + 2 )th fermion as follows with which we can write the hopping terms in each column together, obtaining the vertical hopping hamiltonian as follows where we define (99)

C.3 Coulomb interaction
For the Coulomb interaction, each term in the last sum in the Hamiltonian in Eq. (81) reads where I j is the identity operator and the Hamiltonian of the coulomb interaction reads Finally, we write the Hamiltonian of the Hubbard model in terms of spin-1/2 operators

Appendix D: Digital decomposition of the hopping Hamiltonian for a 2 × 3 Fermion Hubbard model
In this section, we decompose the exact evolution of the horizontal hopping and vertical hopping of a 2 × 3 Fermion Hubbard model respectively into a sequence of discrete gates by applying Trotter expansion.

D.1 Horizontal hopping
For a 2 × 3 fermion lattice, we first consider the horizontal hopping Hamiltonian H * hori as shown in Fig. 9(a) 10 , corresponding to the horizontal hopping for the spin-up fermion (blue solid/dashed arrows) and spin-down fermion (red solid/dashed arrows), respectively (see Fig. 9(a)). By applying the JW transformation, we map the fermonic creation and annihilation operators onto spin operators, and finally obtain the horizontal Hamiltonian as follows where the four terms correspond to the blue solid arrows, blue dashed arrows, red solid arrows, and red dashed arrows, respectively in Fig. 9(a). Now we approximate the time evolution of horizontal hopping Hamiltonian by applying the first-order Trotter expansion where U * word (t) = e −iH * word t .

D.2 Vertical hopping
The Hamiloinian of the vertical hopping can be written in terms of where represents the hopping between the th qubit and ( + 4)th qubit.
To avoid the sub-gates required in the same interaction sharing the qubits, we define the eight terms in Eq. (106) into five groups {h 1 , h 6 }, {h 2 , h 7 }, {h 3 , h 8 }, {h 4 }, and {h 5 }, where each group includes both upward (solid arrows) and downward (dashed arrows) hopping as shown in Fig. 15. Furthermore, all the interactions with the same color and same texture (solid/dashed) can be simulated at the same time i.e. to simulate the vertical hopping of a 2 × 3 fermion lattice, it requires ten types of interactions, and each interaction needs with the interaction Hamiltonian between the jth and ( j + 1)th qubit Here, we neglected the fast oscillating terms proportional to exp(±i(∆ 12 + ν ( j) 1(2) )t), exp(±i(µ 12 + ν ( j) 1(2) )t), exp(±i∆ 12 t), exp(±iµ 12 t), as we assume that the adjacent qubits are far from resonance and the coupling strength Fig. 16, we show that we can activate coupling terms    Table 6 Phase parameters required to simulate the time evolution of the vertical hopping h 5 , see Eq. (108e).    Table 8 Phase parameters required to simulate the time evolution of the vertical hopping h 3 and h 8 , see

Vertical Hopping Operatorφ
Eq. (108c).  Table 9 Phase parameters required to simulate the time evolution of the vertical hopping h 2 and h 7 , see

Vertical Hopping Operatorφ
Eq. (108b).  Table 10 Phase parameters required to simulate the time evolution of the vertical hopping h 1 and h 6 , see