Digital-analog quantum computing of fermion-boson models in superconducting circuits

In solid-state physics, the Hubbard–Holstein (HH) model is commonly used to describe fermionic lattices interacting with phonons¹. This model captures a range of phenomena, such as the formation of polarons—quasiparticles consisting of an electron or hole dressed by a phonon cloud—and their bound states^2,3,4. It plays a vital role in the determination of the transport properties of biomolecules^5,6, as well as the correlation effects and localization phenomena in materials such as cuprates, fullerides, and manganites^7,8,9,10. Consequently, the HH model and similar systems, describing fermionic lattices coupled to bosons, are fundamental in various research fields and technologies in chemistry, materials, and high-energy physics.

Strongly correlated fermionic systems, including those involving fermion-phonon interactions, pose significant challenges due to their complexity. Achieving exact solutions for general fermionic systems becomes impractical beyond a few particles, leading to the need for approximations and symmetry conditions. Various theoretical and computational techniques have been developed to address these challenges, including mean-field theories^11,12,13, quantum Monte Carlo methods^14,15, density functional theory¹⁶, tensor network methods^17,18, and the density matrix renormalization group algorithm^19,20. These approaches provide ways to approximate the behavior of many-body systems and gain insights into their properties, although mean-field theories overlook important correlation effects. These techniques play a crucial role in advancing technology in the era of the second quantum revolution.

Quantum simulation and quantum computing offer advanced approaches for studying the dynamics of strongly correlated systems²¹. By encoding the targeted model into controllable quantum processors, the quantum evolution can be mimicked in a digital, analog, or digital-analog manner. In digital simulations, a specific sequence of single-qubit and two-qubit gates, forming a universal set of quantum logic gates, is used to emulate any quantum evolution. However, the fidelity of single and two-qubit gates, as well as qubit coherences, limit the usefulness of this approach in our noisy intermediate-scale quantum (NISQ) era. In contrast, analog simulations naturally perform the desired quantum evolution within specific parameter ranges with a high accuracy, but it only works for the proposed specific model. Examples include quantum simulations of light–matter interactions in trapped ions^22,23 and superconducting circuits²⁴. A merged approach called digital-analog quantum computing (DAQC)²⁵ bridges the gap between these two paradigms. DAQC leverages analog dynamics from the simulator Hamiltonian, in the spirit of analog blocks, which are combined with digital steps, using the Trotter-Suzuki decomposition in an algorithmic manner. In this co-design approach, the quantum algorithms and hardware architecture are developed together, tailoring each to optimize the overall system performance²⁶.

Quantum simulations have been proposed for studying strongly correlated fermionic systems, including e–p interactions, on various platforms such as trapped ions²⁷, cold atoms^28,29,30, and superconducting circuits^31,32. For qubit-based platforms, the key challenge lies in mapping fermionic degrees of freedom onto spin degrees of freedom, achieved, for example, through the Jordan–Wigner transformation. Trapped ions benefit from the Mølmer–Sørensen gate, whether as a well-behaved two-qubit gate or a multiqubit quantum operation. Alternatively, superconducting architectures have also been advanced, demonstrating comparable scalability features through parallel interactions and geometric considerations³³. Moreover, superconducting circuits offer an ideal platform for compact digital-analog quantum simulations of strongly correlated fermionic systems with multiple bosonic modes, like the Hubbard–Holstein model. Such versatility is due to the possible addition of multiple waveguide cavities or LC oscillators.

We present a superconducting circuit design for compact simulation of the HH model within the DAQC paradigm, applicable to one- and two-dimensional arrays. Our design comprises qubits, resonators, and superconducting quantum interference devices (SQUIDs) arranged along a chain. This architecture enables natural qubit-resonator and tunable qubit-qubit interactions, facilitating the implementation of analog blocks to encode the evolution under the HH model. Such hybrid continuous-variable and discrete-variable quantum computing hardware combines the strengths of both systems, enabling advanced quantum computations and error correction, and is supported by novel compilation techniques and formal architectures to bridge applications and hardware³⁴. We employ the DAQC protocol on a toy model of a half-filling two-site HH Hamiltonian, numerically calculating the fidelity with respect to the exact evolution. Comparing our results to a purely digital approach, we obtain enhanced accuracy. Furthermore, by evaluating resource scaling in terms of gate count, we surpass the performance of existing trapped-ion platforms. Our proposal is thus relevant for exploring phase diagrams, polaron formation, and applications in materials, with implications in chemistry and high-energy physics.

The HH model is a minimal description of how electrons on a lattice interact due to their electric charge and couple to phononic vibrations. The model includes e–p interactions, e–e repulsion, and electron hopping between lattice sites. For an N-site lattice, it is described by the Hamiltonian

where ({c}_{jsigma }^{dagger }({c}_{jsigma })) is the fermionic creation (annihilation) operator for site j and spin σ, ({a}_{j}^{dagger }({a}_{j})) is the creation (annihilation) operator for the jth bosonic mode, and ({hat{n}}_{jsigma }={c}_{jsigma }^{dagger }{c}_{jsigma }). The parameters ω₀ and U are the bosonic frequency and on-site interaction, respectively. From the coupling terms, k is the hopping energy in the lattice, 〈j, l〉 refers to nearest-neighbor sites, and g represents the e–p coupling strength. We note that if ω₀ = g = 0 (without bosonic modes), the resulting Hamiltonian corresponds to the Hubbard model. On the other hand, if U = 0 (fermions without charge), we obtain the Holstein model.

To simulate the HH model, we need to encode the bosonic and fermionic degrees of freedom. For the first one, in superconducting circuits, we can use transmission lines or LC oscillators, obtaining an analog to the phonons in the HH model. On the other hand, the fermion operators (with spin) can be mapped to spinless fermion operators in a lattice of double the size, that is c_j,↑ → b_2j and c_j,↓ → b_2j+1. Finally, to encode spinless fermionic operators, we perform the Jordan–Wigner (JW) transformation^35,36, ({b}_{j}to mathop{prod }nolimits_{j = 1}^{i-1}(-{sigma }_{j}^{z}){sigma }_{i}^{-}). In this manner, we map the HH lattice to a spin chain interacting with bosonic modes (see Fig. 1). The latter can be described by the Hamiltonian

where the term (g({a}_{j}+{a}_{j}^{dagger })) was absorbed by the change of variables a_j → a_j + g/ω₀ and therefore (bar{U}=U-4{g}^{2}/{omega }_{0}). We can see that the first and second terms in Eq. (2) correspond to the free energy of bosonic and spin degrees of freedom, easily implementable by a set of LC oscillators (or transmission lines) and qubits. The next two terms represent the qubit-qubit and qubit-boson interactions, which can be done in superconducting circuits using capacitive coupling and local rotations in the qubits. The last terms correspond to the fermionic tunneling terms, which can be implemented by using consecutive two-body rotations, as shown in ref. ³³.

a The HH lattice to be simulated. The bosonic frequency is ω₀, the hopping energy between lattice sites is − k, and g is the e–p coupling strength. b Mapping of HH lattice to an equivalent spinless fermion lattice using the JW transformation. Each site (shaded) is represented by two fermionic pairs, and the on-site repulsion is given by U. c Representation of the mapped 1D HH model (each site with a different color). d The corresponding superconducting circuit architecture representing the spin-boson chain using a 1D linear chain of transmon qubits (({{rm{Q}}}_{j}^{L(R)})) coupled via resonators (R_j).

Full size image

Digital-analog quantum simulator

We describe here an architecture that allows to implement the Hamiltonian of Eq. (2) efficiently: a chain of qubits coupled to a set of resonators. To simulate this Hamiltonian, we propose the superconducting circuit architecture shown in Fig. 1d. Such a superconducting design consists of a linear array of transmon qubits (({{rm{Q}}}_{j}^{L(R)}), where j is the site index and L(R) denotes the left or right qubit at the same site) coupled through grounded superconducting quantum interference devices (SQUIDs). In addition, the chain of qubits has LC resonators (R_j) coupled at each site by grounded SQUIDs. For the transmon array, we intersperse devices with different spectra. It means that nearest-neighbor qubits are off-resonance, i.e., qubits in odd positions have a frequency different from qubits in even positions, and the resonators also have a different frequency than the qubits. Therefore, our architecture can be seen as a chain of building blocks, enclosed by a dashed square in Fig. 1, composed of three different devices: a left qubit with frequency ω_L, a right qubit with frequency ω_R, and a resonator with frequency ω. A similar architecture without resonators was proposed in ref. ³³, and another one including experimentally performed resonators in ref. ³⁷, both in the context of quantum simulations.

The interactions between qubits and between each qubit and the resonator at each site are controlled through a single SQUID. This is achieved by applying an external flux to the SQUID, consisting of a DC component and a small AC signal. The flux is expressed as (see Supplementary Note A)

where

Here, α = {+, −} and ({varphi }_{ab}^{alpha }={A}_{ab}cos ({nu }_{ab}^{alpha }t+{tilde{varphi }}_{ab}^{alpha })) with the amplitude A_ab, frequency ν_ab and phases ({tilde{varphi }}_{ab}^{alpha }). With this form of the signal, we specify which part of the subsystem is controlled by ({tilde{varphi }}_{ab}^{alpha }). We can then turn on (off) the specific phases ({tilde{varphi }}_{ab}^{alpha }) to activate (deactivate) the couplings. The full evolution of the qubit-qubit and qubit-resonator interactions is achieved in a sequence of Trotter steps.

After the high-plasma-frequency and low-impedance-SQUID approximations, our architecture is described by the interaction-picture Hamiltonian

where ({H}_{s}^{(j)}) and ({H}_{s-s}^{(j)}) are the jth building-block Hamiltonian and the interaction between building blocks j and j + 1, respectively. Such Hamiltonians have the form (for a detailed calculation, see Supplementary Note A)

where ({C}_{ab}^{pm }=cos {tilde{varphi }}_{ab}^{+}pm cos {tilde{varphi }}_{ab}^{-}), ({tilde{C}}_{ab}^{pm }=cos {tilde{phi }}_{ab}^{+}pm cos {tilde{phi }}_{ab}^{-}), and similar definitions for ({S}_{ab}^{pm }) and ({S}_{ab}^{pm }), replacing (cos) by (sin). The phases ({tilde{varphi }}_{ab}^{pm }) are associated with the magnetic flux through the SQUIDs inside each building block, and ({tilde{phi }}_{ab}^{pm }) are the phases of the magnetic flux through the SQUIDs that couple different building blocks. The magnetic flux parameters required to activate or deactivate specific interactions are shown in Supplementary Table 1. We note that by adjusting the different phases, we can engineer a wide range of qubit-qubit and qubit-boson interactions.

Simulations

For simulation, we consider the dimensionless Hamiltonian H/k and choose the timescale in units of 1/k. For different materials, the range of parameters of the HH model of Eq. (1) is such that k and U lie between 0.2 and 10 eV, g is between 0 and 100 eV, and ω lies in the range 0.01 – 1 eV³⁸. The corresponding physical time scales of the system are between fs and ps. Considering a microwave driving with qubit and resonator frequencies in the range 5 – 10 GHz, the Trotter time step for the digital blocks comprising the single and two qubit gates are in the range 10 –500 ns³⁹, while for the analog blocks comprising the bosonic and qubit-bosonic gates, they are around 72 ns³⁷. The initial states of the system are defined using the ordering of the Hilbert space as (leftvert {Q}_{1},{Q}_{2},{Q}_{3},{Q}_{4},re{s}_{1},re{s}_{2}rightrangle), where Q_i and res_i denote the qubit and resonator subspaces, respectively.

In Fig. 2, we show the time-dependent state fidelity ∣〈ψ_exact∣ψ_sim〉∣² of the noise-based simulation as a function of time for 50 Trotter steps. The time step δt is kept constant, so that the total evolution consists of 50 discrete steps. The noisy DAQC simulation incorporates amplitude-damping errors for both the resonators and the qubits. The damping factors are determined based on the duration of the analog and digital blocks T_block divided by their relaxation times T₁ of the resonators and qubits, such that γ = T_block/T₁. The duration of the analog block is 50 ns, the duration of the digital block is 200 ns, and the relaxation time for both qubits and resonators was set to T₁ = 80 μs. This leads to calculated damping factors of ({gamma }_{{rm{res}}}=0.000625) for the resonators and γ_qubit = 0.0025 for the qubits.

The initial state is half-filling, ({psi }_{i}=leftvert 0,1,1,0,0,0rightrangle). Each resonator was truncated to eight levels. a weak coupling (g = 0.1k), b strong coupling (g = k) and c ultra-strong coupling (g = 5k).

Full size image

Unlike digital blocks, which accumulate errors from both amplitude damping and gate imperfections, analog blocks are only affected by amplitude damping without additional noise sources introduced by discrete gate operations. Noise is applied using Kraus operators to model amplitude damping, where each Kraus operator A_k(t) is constructed for a given time t, decay rate γ, and system dimension. These operators modify the quantum state iteratively by applying them in time steps, ensuring that the normalization condition ({sum }_{k}{A}_{k}^{dagger }(t){A}_{k}(t)=I) is satisfied to conserve total probability. For the reservoir components, noise is applied using a set of tensor product Kraus operators specific to the subsystem dimensions, while for the qubit components, similar noise channels target individual qubits. In the noisy digital simulation, depolarizing errors are combined with thermal relaxation errors for each gate, with gate fidelities derived from IBM devices. The digital noise model uses the Qiskit noise model to construct these errors, handling both single-qubit and two-qubit operations by combining depolarizing and thermal relaxation components to reflect realistic noise conditions in digital quantum simulations.

The oscillations in the fidelity curves of the ideal cases, shown by the orange triangles (Digital) in Fig. 2, and the red dots (DAQC) in the same figure, appear from the non-commuting terms in the Hamiltonian, for which the Trotter decomposition of the unitary time evolution operator introduces a phase factor. Such oscillations have been documented in previous studies^27,40. In the noisy simulations, it is observed that noise affects these oscillations, and they are canceled out. The reduced oscillation amplitude in the DAQC approach is expected as fewer Trotter decompositions are performed. The DAQC approach, therefore, achieves more stable dynamics with fewer artifacts, making it a powerful method for simulating complex quantum systems.

The digital simulations were carried out using Bosonic Qiskit⁴¹, where each resonator was truncated to n levels. This requires an order of ({log }_{2}n) number of digital gates each, apart from the gates required for the qubits. The terms a^†a, a^† + a, and σ_z(a^† + a) can be implemented using cv_r(), cv_d(), and cv_c_d() gates, respectively. This is provided by Bosonic Qiskit and corresponds to the phase space rotation, displacement, and controlled displacement operators, respectively. When decomposed into digital gates, these bosonic gates can be realized using combinations of single-qubit rotations U3(θ, ϕ, λ) and CNOTs.

Since entangling gates are common in NISQ processors, we estimate the number of CNOTs required for each of these bosonic gates for a single Trotter step. We plot the CNOTs required per bosonic gate with increasing resonator levels (n) in Fig. 3. The CNOTs required for the gates, cv_r() and cv_d(), grow equally. In case of cv_c_d(), it grows with a very high margin, requiring a total of 67 CNOTs for n = 8.

Number of CNOTs (N_CNOT) required to perform each bosonic gate as a function of the number of levels (n) for each bosonic mode for a pure digital approach.

Full size image

From a dynamic point of view, our DAQC protocol can produce very high fidelities. For the two-site HH model, we compute the overlap between the exact solution and the numerical simulation using the DAQC protocol for different regimes, as shown in Fig. 4. Based on the observation that a high-fidelity simulation of the HH model is feasible, as seen in Fig. 2, we investigate which coupling parameters (k, g, and U) can be realized by the superconducting architecture to accurately simulate the model in these regimes. The coupling regimes suitable for a two-site simulation are presented in Fig. 4 (ideal simulation). In Fig. 4a, the red marker indicates the final fidelity of the ideal DAQC simulation shown in Fig. 2. The range of fidelities in Fig. 2 is above 0.988, which is consistent with the ideal DAQC curve in Fig. 2. The architecture simulates the two-site HH model efficiently within the given parameter regime and can therefore be extended to a linear chain or a 2D architecture. For the latter, there is a possibility to demonstrate a quantum advantage for interacting fermion-boson systems. We observe that the fidelity in Fig. 4 is higher than 0.98, making our proposal suitable for studying the dynamical properties of condensed matter systems. Moreover, since the scaling of our proposal depends only on the smaller dimension of the lattice, we anticipate that these results can be maintained for larger systems.

DAQC fidelity ∣〈ψ_exact∣ψ_sim〉∣² as a function of a the on-site interaction U and the e–p coupling g in units of the hopping energy k. The red dot represents the last point of the ideal DAQC simulation curve in Fig. 3 (red). b U and k in units of g, and c k and g in units of U for 50 Trotter steps. The initial state is half-filling, ({psi }_{i}=leftvert 0,1,1,0,0,0rightrangle), where each resonator was truncated to eight levels.

Full size image

Quantum simulation of condensed-matter physics

The formation and decay of bound states due to electron–electron and electron–phonon interactions have profound effects on the electrical and thermal transport properties of materials⁴². Despite its simplicity, the HH model captures some of the important non-trivial phenomena arising from these interactions. The DAQC protocol proposed here can be used to simulate the real-time dynamics of this model, not easily accessible in condensed-matter experiments. For instance, spatially-resolved and temporally-resolved density-density correlation functions, accessible here, were used in previous numerical studies to characterize polaron formation^43,44. Exotic bound states such as repulsively bound pairs could also be simulated—these states can occur in lattice systems when the interaction energy is larger than the kinetic energy bandwidth. They are expected to decay quickly due to dissipation but were recently observed in antiferromagnetic compounds⁴⁵ and in cold atoms in optical lattices⁴⁶.

As an example of the physical quantities accessible with the proposed hardware and algorithm, we compute the electron double occupation (langle {hat{n}}_{juparrow }{hat{n}}_{jdownarrow }rangle) and the total number of phonons (langle {hat{n}}_{ph}rangle =langle {hat{a}}_{2}^{dagger }{hat{a}}_{2}rangle +langle {hat{a}}_{2}^{dagger }{hat{a}}_{2}rangle) as functions of time. These quantities characterize the dynamics of a repulsively bound pair (with U ≫ k) initially localized on site j = 1. We consider a weak and a strong electron–phonon coupling in the regime g ≤ 5k.

For weak coupling g = 0.1k, shown in Fig. 5a, the total electron double occupation 〈n_1↑n_1↓〉 + 〈n_2↑n_2↓〉 stays close to the initial value 1, indicating that the electrons remain in a bound state. The high-frequency oscillation (with frequency close to U) results from a virtual pair-breaking and reassociation process where one electron tunnels to the second site and back. The site-resolved double occupation also shows a slower oscillation associated with the tunneling of the bound pair between the two sites. The total number of phonons in Fig. 5d is close to zero in this case. For g = 5k [Fig. 5b], the oscillations become irregular as the strong coupling to phonons introduces more degrees of freedom and more eigenfrequencies contribute to the dynamics. The total double occupation is slightly reduced, and the number of phonons is nonzero, which suggests that the pair can be broken when the interaction energy is dissipated by the creation of phonons. These observations illustrate the possibilities of the DAQC algorithm for the quantum simulation of condensed-matter physics relevant to material properties. They also highlight its advantage over a fully digital algorithm, since a high fidelity is crucial for the measurement of high-frequency oscillations.

a The electron double occupation (langle {hat{n}}_{juparrow }{hat{n}}_{jdownarrow }rangle) on sites j = 1 and j = 2 and the total double occupation (langle {hat{n}}_{1uparrow }{hat{n}}_{1downarrow }rangle +langle {hat{n}}_{2uparrow }{hat{n}}_{2downarrow }rangle) computed with the DAQC algorithm. These quantities show coherent oscillations as functions of time for weak coupling g = 0.1k. b For strong coupling g = 5k, the oscillations become irregular and the average total double occupation is reduced. c The total number of phonons (langle {hat{n}}_{ph}rangle =langle {hat{a}}_{1}^{dagger }{hat{a}}_{1}rangle +langle {hat{a}}_{2}^{dagger }{hat{a}}_{2}rangle) is close to zero for weak coupling. For strong coupling, it is nonzero with irregular oscillations due to the large number of degrees of freedom involved in the dynamics. Here, the on-site interaction and phonon frequency are U = ω₀ = 8k.

Full size image

We developed a co-design DAQC approach for the quantum simulation of strongly correlated fermion-boson models in superconducting circuits. The proposed digital-analog encoding allows us to simulate the problem Hamiltonian with high fidelity and fewer resources as compared to the purely digital case. We obtain a circuit-depth scaling proportional to ℓ² for a ℓ × h lattice where ℓ ≤ h. Such a compact scaling cannot be obtained with purely digital quantum computing, where the number of controlled gates and qubits grows with both dimensions in the lattice. Conversely, the scaling provided by other platforms like trapped ions grows with both dimensions, ℓ ⋅ h, and does not have the versatility to encode several bosonic modes⁴⁷. This implies that for square lattices where ℓ = h, both platforms exhibit similar scaling. However, in elongated lattices, the SC circuit architecture maintains a lower circuit depth, which can lead to improved performance and reduced error accumulation.

To demonstrate the performance of our model, we consider a two-site HH model at half-filling. We perform noise-based simulations considering the current state-of-the-art, where the DAQC approach exhibits superior performance in comparison to the digital. For various coupling regimes that can be explored with the device, the model can achieve high fidelity values, larger than 0.98 in all cases.

Additionally, to show the relevance of our model for condensed matter systems and its applicability in simulation of materials, we explored the dynamics of a repulsively bound electron pair in a two-site system, where the coupling to phonons introduces a dissipation channel and leads to pair breaking.

Our encouraging results suggest that these types of co-design encoding are suitable to outperform classical computation capabilities in the NISQ era, making our current quantum technology useful for industrial problems related to condensed-matter physics.

Activating digital and analog blocks

We now focus on mimicking each interaction term of the target Hamiltonian of Eq. (2). The on-site repulsion ({sigma }_{2j}^{z}{sigma }_{2j-1}^{z}={sigma }_{L,j}^{z}{sigma }_{R,j}^{z}) and the e–p terms can be realized simultaneously by turning off all the magnetic fluxes (({A}_{QQ}={A}_{QR}^{(j)}={bar{A}}_{QQ}=0)), applying a Hadamard gate over each qubit, and then turning on the signal QQ and QR in each building block (A_QQ ≠ 0 and A_QR ≠ 0), while adjusting ({tilde{varphi }}_{QQ}^{+}={tilde{varphi }}_{QQ}^{-}=0). This allows the system to evolve according to

with A_zz ~ A_QQ and A_e−ph ~ A_QR. The operators ({sigma }_{R(L),j}^{x}) act on the right (left) qubit at the same site j. Finally, we again apply Hadamard gates over all qubits, turning off all the signals. Accordingly, the implementation of the on-site interaction and electron–phonon coupling terms can be implemented as

where H^⊗2n are Hadamard gates applied over each qubit.

For the fermionic part, we consider an HH model with a ℓ × h lattice with ℓ ⋅ h = N and h > ℓ/2, implying a spinless fermion lattice of 2ℓ × h sites. We note that in such a geometry, there are horizontal and vertical nearest-neighbor sites. For the horizontal ones, l = j + 1. In Eq. (2), we then have the fermionic term

Here, the first part can be written as

where the gate is ({U}_{xx}={e}^{-ifrac{pi }{4}{sum }_{j}{sigma }_{j,R}^{x}{sigma }_{j+1,L}^{x}}). Such nearest-neighbor Ising interactions are natively available in superconducting hardwares⁴⁸. Moreover, these interactions can be implemented via the evolution generated by the device Hamiltonian in Eq. (4), setting ({bar{A}}_{QQ},ne, 0), ({tilde{phi }}_{QQ}^{+}=0), and ({tilde{phi }}_{QQ}^{-}=pi) for a duration of (t=pi /({m}_{1}{tilde{A}}_{QQ})). Consequently, the terms from the first line of Eq. (7) can be simulated by the gate sequence ({U}_{xx},{U}_{xy}(t),{U}_{xx}^{dagger }), where the operation U_xy(t) can be generated from the device Hamiltonian by setting ({A}_{QQ}={bar{A}}_{QQ},ne, 0) and ({tilde{varphi }}_{QQ}^{+}={tilde{varphi }}_{QQ}^{-}=-pi /2). The last gate of the sequence, ({U}_{xx}^{dagger }), is implemented by setting ({tilde{phi }}_{QQ}^{+}=-{tilde{phi }}_{QQ}^{-}=pi /2).

In the same way, the second line in Eq. (7) can be implemented as

where ({U}_{yy}={e}^{-ifrac{pi }{4}{sum }_{j}{sigma }_{j,R}^{y}{sigma }_{j+1,L}^{y}}), which can be implemented by the evolution generated by the device Hamiltonian of Eq. (4) with ({bar{A}}_{QQ}ne 0) and ({tilde{phi }}_{QQ}^{+}={tilde{phi }}_{QQ}^{-}=0) over a duration (t=pi /({m}_{1}{tilde{A}}_{QQ})). Thus, the terms present in the second line of Eq. (7) can be simulated by the gate sequence ({U}_{yy},{U}_{yx}(t),{U}_{yy}^{dagger }), where the operation U_yx(t) can be implemented by setting ({A}_{QQ}={bar{A}}_{QQ}ne 0), where ({tilde{varphi }}_{QQ}^{+}=-{tilde{varphi }}_{QQ}^{-}=pi /2) and ({tilde{phi }}_{QQ}^{+}={tilde{phi }}_{QQ}^{-}=-pi /2).

For the vertical fermionic terms l = ℓ + j, we have

This term is more complex, and to implement it, we can follow a similar approach to ref. ³³. Without loss of generality, we consider odd ℓ and write the first term of Eq. (10) as

Here, s_j = (ℓ − 1)/2 + j,

and

Circuit depth

To implement terms of the form given by Eq. (11), we need to perform 2ℓ rotations plus the interaction σ^xσ^y, obtaining a circuit depth of 2ℓ + 1. As we have shown, all these two-body rotations can be done by adjusting the phases of the different signals (see Supplementary Table 1). In a similar way, we can write the second term of Eq. (10) as

which adds 2ℓ + 1 layers to the circuit depth. We note that with 4ℓ + 2 layers, we can do one vertical interaction. Nevertheless, it can be shown that all the vertical interactions can be done at the same time for the same column. This is because all the interactions commute, meaning that they can be done together as in the proposal of ref. ³³. Therefore, we need 2ℓ(2ℓ + 1) layers to simulate all the vertical interactions. It is important to highlight that this is the only term that provides scaling to our algorithm.

The circuit depth of our algorithm is now three for all the e–p and on-site interaction terms, six for all the horizontal fermionic terms, and 2ℓ(2ℓ + 1) for all the vertical ones. This means we need a total circuit depth of 2ℓ(2ℓ + 1) + 9 to simulate a h × ℓ HH lattice, where ℓ ≤ h and ℓ > 1. For the case of a chain (ℓ = 1), the circuit depth is 9 and does not scale with the chain size. As an example, we show the circuit for our DAQC encoding for a two-site HH model in Fig. 6. Therefore, if we want to simulate a quantum evolution of the HH model using our DAQC approach, we have, in the interaction picture,

Here, N is the number of Trotter steps, and each exponential can be implemented as mentioned before. For such evolutions, our circuit depth is given by 2Nℓ(2ℓ + 1) + 9. It is particularly interesting if we compare our findings with the digital approach, where the circuit depth always scales with both dimensions. We also note that in a digital quantum computer, the implementation of bosonic modes requires ({log }_{2}(N)) qubits, where N is the maximum energy level considered in the bosonic dynamics. This makes our DAQC approach also efficient in terms of hardware resources. Another platform that offers similar possibilities is trapped ions²⁷. Nevertheless, it fails in the ability to implement multiple bosonic modes and, as it uses Mølmer and Sørensen gates, the circuit depth also scales with the two dimensions of a lattice.

The gates in green are the analog blocks, while the blue ones are the digital steps. The circuit depth is 9 which does not scale as we add more sites in the 1D array.

Full size image