Dynamical Quantum Phase Transitions of the Schwinger Model: Real-Time Dynamics on IBM Quantum

Abstract

Simulating the real-time dynamics of gauge theories represents a paradigmatic use case to test the hardware capabilities of a quantum computer, since it can involve non-trivial input states’ preparation, discretized time evolution, long-distance entanglement, and measurement in a noisy environment. We implemented an algorithm to simulate the real-time dynamics of a few-qubit system that approximates the Schwinger model in the framework of lattice gauge theories, with specific attention to the occurrence of a dynamical quantum phase transition. Limitations in the simulation capabilities on IBM Quantum were imposed by noise affecting the application of single-qubit and two-qubit gates, which combine in the decomposition of Trotter evolution. The experimental results collected in quantum algorithm runs on IBM Quantum were compared with noise models to characterize the performance in the absence of error mitigation.

1. Introduction

The availability of noisy intermediate-scale quantum (NISQ) devices in cloud access platforms is a fundamental step towards the quantum computing era. Nonetheless, the limited number of available qubits and the absence of controllable error probabilities prevents these systems from actually outperforming current classical computing capabilities. Multiple hardware setups have been engineered for quantum computing purposes, with different advantages regarding gates’ fidelity and experimental realization. Examples of NISQ devices are represented by circuits with superconducting transmon qubits [1,2,3], ion traps or optical lattices hosting Rydberg atoms [4,5,6,7,8], and qubits encoded in photonic modes of optical setups [9,10], this list being far from exhaustive.

In this framework, high-energy physics represents an interesting testbed for quantum devices. On the one hand, quantum computation can be applied “downstream”, to optimize data analysis and event reconstruction from experiments [11,12,13,14,15]. On the other hand, the “upstream” investigation of gauge theories, especially in their lattice formulation [16,17,18,19,20], can benefit from the possibility of performing quantum simulations of regimes not achievable with perturbative techniques. Long-standing questions related to low-energy processes in quantum chromodynamics (QCD) are still far from the current capabilities of Monte Carlo techniques, due to the sign problem in fermionic amplitudes [21,22,23]. To overcome some of these limitations, research at the interface of quantum information, condensed-matter, and high-energy physics is targeting the adoption of new theoretical and computational tools. The state-of-the-art in the field is represented by tensor network methods, able to reduce the exponential complexity to a polynomial one for states characterized by short-range entanglement [24,25,26,27,28,29,30]. These methods are suitable candidates to obtain a breakthrough in non-perturbative regimes: some preliminary studies about quantum electrodynamics (QED) in one spatial dimension proved the ability of tensor networks in describing a wide phenomenology, such as vacuum phase transition, the string-breaking mechanism, and scattering processes [31,32,33]. On the other hand, highly entangled quantum systems must be studied by means of specifically designed setups, since their complexity cannot be managed with current classical computing capabilities.

A major role in making quantum computation and simulation effective to solve practical problems in NISQ devices is played by error correction and mitigation [34,35,36]. In the context of digital real-time evolution [37,38,39,40,41,42], such procedures should keep the quantum state of the system in the physical subspaces allowed by the gauge constraint. Unfortunately, these techniques are expected to limit the quantum advantage due to the required classical information processing. Nevertheless, current simulations on NISQ devices would benefit from the aforementioned techniques in the description of targeted physical phenomena. The hardest operation for any device is represented by a controlled gate, whose implementation varies in each platform, as with ion traps [4,41,42] or superconducting circuits [37,38,39,40]. These setups elaborate single-qubit states by means of laser pulses, tuned with the targeted transitions. Instead, entangling gates are specifically engineered by exploiting the properties of the hardware, thus leading to crucially different fidelities and times required for the physical implementation. Alternatively, analog simulations, e.g., in optical lattices, can adopt different strategies, such as a periodic drive to obtain energy terms endowed with the same symmetry characterizing lattice QED [43,44]. In this case, errors can induce gauge-invariance-breaking terms, which can lead to an emergent prethermal behavior [45,46,47,48].

This paper aimed at testing the superconducting qubit systems available in the IBM Quantum platform [1] in a simple lattice gauge theory application. We implemented digital evolutions generated by the QED Hamiltonian in $1+1$ dimensions, consisting of non-commuting local contributions [4,31,32,33,37,41,49,50,51,52,53,54,55,56] and displaying the occurrence of dynamical quantum phase transitions (DQPTs) [57,58,59] in specific cases of quantum quenches. A discretization of the electromagnetic field is required to implement gauge degrees of freedom in the quantum simulator: accomplishing such a task by replacing the continuous $\mathrm{U}\left(1\right)$ gauge group with one of the cyclic groups ${\mathbb{Z}}_n$ guarantees a unitary implementation of the gauge connections [49]; another possible choice is represented by the family of quantum link models [60]. In order to assign each link of the lattice to a single qubit, we chose a ${\mathbb{Z}}_2$ discretization, which ensures a minimal circuit depth [52], thus reducing error sources. The study of real-time dynamics involves three crucial stages: the preparation of initial states, the evolution, and the final measurements. Gate errors affecting a digital simulation accumulate in a more or less coherent manner, which is affected by more variables and more error sources than the stand-alone characterization of gates. Time evolution is partitioned into steps to monitor the measurements statistics’ variation without the inclusion of error correction and mitigation. The ground state preparation required by the chosen quench protocol was specifically designed in order to minimize errors and characterize the first stage’s output statistics. Then, we characterized the effectiveness of time evolution, choosing to analyze the system in proximity of a DQPT, where the dynamics is particularly sensitive to noise [51], thus framing the simulation in most unsafe conditions. A concluding estimation of the amount of error probability reduction required for a partial observation of the targeted DQPT was made by analyzing the statistics of the collected results.

The paper is structured as follows. In Section 2, we introduce the lattice Schwinger model and the Jordan–Wigner transformation that maps it into a qubit system and describe the quench protocol and the DQPTs expected in the model. In Section 3, we describe the experimental scheme composed of ground state preparation and the subsequent Trotter evolution. The collected results are compared with the simulated evolution affected by the error probabilities of noisy gates. In Section 4, we discuss our outcomes in light of the existing literature and present a possible outlook for our research.

2. The Lattice Schwinger Model

QED in $1+1$ dimensions, also known as the Schwinger model, is a $\mathrm{U}\left(1\right)$ gauge theory describing the interaction of the electromagnetic field, consisting of only an electric component, and a fermionic particle with mass m and charge g. The model can be discretized on a one-dimensional lattice with spacing a, by associating with each lattice site x an anticommuting field ${\psi }_x$, which represents a spinless fermion, while links between each pair of neighboring sites host the gauge degrees of freedom, described by the electric field $E_{x,x+1}$ and the vector potential $A_{x,x+1}$. The latter determines the gauge connection $U_{x,x+1}={\mathrm{e}}^{\mathrm{i}a{A}_{x,x+1}}$, characterized by the generalized canonical commutation relation $[E_{x^{\prime },x^{\prime }+1},U_{x,x+1}]={\delta }_{x,x^{\prime }}{U}_{x,x+1}$. The lattice model Hamiltonian for a finite lattice with N sites reads [4,31,32,33,49]

$$H=-\frac{\mathrm{i}}{2a}\sum _{x=0}^{N-1}\left({\psi }_x^{†}{U}_{x,x+1}{\psi }_{x+1}-\mathrm{H}.\mathrm{c}.\right)+m\sum _{x=0}^{N-1}{(-1)}^x{\psi }_x^{†}{\psi }_x+\frac{g^{2}a}{2}\sum _{x=0}^{N-1}{E}_{x,x+1}^2,$$

(1)

where periodic boundary conditions [50,51] require the identification $N\equiv 0$. The model involves staggered (Kogut–Susskind) fermions [61], described by single-component spinors ${\psi }_x$, with negative-mass components encoded in odd-x sites. The physical subspace ${\mathcal{H}}_G$ is spanned by states $|\varphi \rangle$ satisfying the Gauss law constraint $G_x|\varphi \rangle =0$ at all sites x, where, for a ${\mathbb{Z}}_n$ gauge group,

$$G_x=\sqrt{\frac{n}{2\pi }}\left(E_{x,x+1}-E_{x-1,x}\right)-{\psi }_x^{†}{\psi }_x-\frac{{(-1)}^x-1}{2}.$$

(2)

The electric field was simulated in the following through a ${\mathbb{Z}}_n$ discretization of U$\left(1\right)$ [31,32,33,49,52,53] with $n=2$. Unlike in the quantum link models [50,51], where the electric field is replaced by a spin operator, the ${\mathbb{Z}}_n$ model is based on replacing gauge connections with permutation matrices [49]. In the case of ${\mathbb{Z}}_2$, the electric field in each link $(x,x+1)$ can have two eigenstates, which will be labeled as ${\left\{\right|\ell =\frac{1}{2}\rangle }_x,|\ell =-\frac{1}{2}{\rangle }_x\}$, with

$$E_{x,x+1}{\left|\pm \frac{1}{2}\right.〉}_x=\pm \frac{\sqrt{\pi }}{2}{\left|\pm \frac{1}{2}\right.〉}_x,$$

(3)

while the gauge connections act as [31,32,33]

$$U_{x,x+1}{\left|\pm \frac{1}{2}\right.〉}_x={\left|\mp \frac{1}{2}\right.〉}_x.$$

(4)

An immediate implication of the ${\mathbb{Z}}_2$ model is the irrelevance in the Hamiltonian (1) of the electric field energy, which becomes a constant.

The simplest nontrivial periodic lattice is composed by $N=2$ sites: the states spanning the physical subspace ${\mathcal{H}}_G$ of the ${\mathbb{Z}}_2$ model in this simple case are shown in Figure 1. The states in Panels (a) and (d) represent two “Dirac vacua”, with a filled negative-mass and an empty positive-mass site. In these states, the total electric field is constant and nonvanishing. These observations motivated the notation ${|vac\rangle }_{\pm }$ for these two states, where the index is related to the sign of the background electric field. Particle hopping leads to the remaining “mesonic” basis states $|e^{+}{e}^{-}{\rangle }_L$ and $|e^{+}{e}^{-}{\rangle }_R$, represented in Panels (b) and (c), respectively, where the index refers to the counterclockwise (L) or clockwise (R) hopping of the fermion from the negative- to the positive-mass site.

The Jordan–Wigner transformation maps the spinor field into a spin system [4,52], which corresponds to our qubit register, as

$${\psi }_x=\prod _{\ell <x}\left(\mathrm{i}{Z}_{\ell }\right)\frac{X_x+\mathrm{i}{Y}_x}{2},{\psi }_x^{†}=\prod _{\ell <x}(-\mathrm{i}{Z}_{\ell })\frac{X_x-\mathrm{i}{Y}_x}{2},$$

(5)

where X, Y, and Z are Pauli matrices, ${\sigma }^{\pm }=\frac{X\pm \mathrm{i}Y}{2}$, and occupied sites correspond to qubit states $|\downarrow \rangle$. The Hamiltonian of the resulting spin system is

$$H=H_J+H_m=\frac{J}{2}\sum _{x=0}^{N-1}\left({\sigma }_x^{-}{U}_{x,x+1}{\sigma }_{x+1}^{+}+\mathrm{H}.\mathrm{c}.\right)-\frac{m}{2}\sum _{x=0}^{N-1}{(-1)}^x{Z}_x=\sum _{x=0}^{N-1}{h}_x,$$

(6)

where the free parameter corresponds to a coupling constant $J=\frac{1}{a}$, once energy is scaled in units of mass m.

2.1. Dynamical Quantum Phase Transitions

We aimed at studying the non-equilibrium dynamics of the described lattice Schwinger model following a quantum quench [50,51,57]. Generally, in this protocol, one considers a family of Hamiltonians $H\left(\gamma \right)$ that depends on a tunable parameter and prepares an initial state coinciding with the ground state $|{\psi }_g\rangle$ of $H_0=H\left({\gamma }_0\right)$. At $t=0$, the Hamiltonian suddenly switches to $H=H\left({\gamma }_f\right)$, determining the evolution $|\psi \left(t\right)\rangle ={\mathrm{e}}^{-\mathrm{i}tH}|{\psi }_g\rangle$, characterized by the survival (or Loschmidt) amplitude:

$$\mathcal{G}\left(t\right)=\langle {\psi }_g|\psi \left(t\right)\rangle .$$

(7)

To identify a possible DQFT, we search for the zeros of the Loschmidt echo:

$$\mathcal{L}\left(t\right)={\left|\mathcal{G}\left(t\right)\right|}^2={\mathrm{e}}^{-N\lambda \left(t\right)},$$

(8)

which depends on the number N of degrees of freedom and on the rate function $\lambda \left(t\right)$, which becomes divergent in correspondence with the aforementioned zeros.

The targeted evolution generated by the Hamiltonian (6) is determined by a single free parameter, as the Loschmidt amplitude phase $\phi (J,t)=arg\mathcal{G}\left(t\right)$ is undefined in correspondence with the critical points. Nonetheless, in their neighborhood, $\phi$ is expected to be smooth up to a discontinuity line of $2\pi$ starting from the criticality. This characterization corresponds to a vortex, with a winding number:

$$\nu =\frac{1}{2\pi }{\oint }_{\mathcal{C}}\mathrm{d}\mathit{s}·\nabla \phi ,$$

(9)

where $\mathcal{C}$ is a loop in the $(J,t)$ plane [50].

The adopted protocol quenches the Kogut–Susskind staggered fermions at $t=0$ by inverting the mass sign: $H(m,J)⟶H(-m,J)$. The quenched Hamiltonian can be decomposed into parity sectors, as described in Appendix A:

$$H(-m,J)=H^{(-)}\oplus H^{(+)}.$$

(10)

In the even sector, the evolution is generated by

$$H^{(+)}=\left(\begin{array}{cc}-\frac{J^2-m^2}{\sqrt{m^2+J^2}}& \frac{2Jm}{\sqrt{m^2+J^2}}\\ \frac{2Jm}{\sqrt{m^2+J^2}}& \frac{J^2-m^2}{\sqrt{m^2+J^2}}\end{array}\right),$$

(11)

in the subspace spanned by the basis $\left\{\right|{\psi }_g\rangle ,|{\psi }_{\overline{g}}\rangle \}$, made of the even eigenstates of the initial Hamiltonian $H_0=H(m,J)$, which are associated with the lowest and highest eigenvalue $E_g=-\sqrt{m^2+J^2}$ and $E_{\overline{g}}=+\sqrt{m^2+J^2}$, respectively. The odd parity sector involves the eigenstates $|{\psi }_e\rangle$ and $|{\psi }_{\overline{e}}\rangle$ of $H(m,J)$, which are independent of J and characterized by the eigenvalues $\pm m$; these two states are still eigenstates of the quenched Hamiltonian $H(-m,J)$, which only inverts their eigenvalues.

The Loschmidt amplitude for the initial ground state:

$$|{\psi }_g\rangle =a_g{\left(\right|vac\rangle }_{+}+{|vac\rangle }_{-})+b_g{\left(\right|e^{+}{e}^{-}\rangle }_L+|e^{+}{e}^{-}{\rangle }_R),$$

(12)

derived in Appendix A, reads

$$\mathcal{G}\left(t\right)={(2a_g^2-2b_g^2)}^2{\mathrm{e}}^{-\mathrm{i}{E}_{g}t}\left(1+\frac{J^2}{m^2}{\mathrm{e}}^{-\mathrm{i}(E_{\overline{g}}-E_g)t}\right),$$

(13)

with

$$a_g=\frac{1}{\sqrt{2\left(1+p_g^2\right)}},b_g=\frac{p_g}{\sqrt{2\left(1+p_g^2\right)}},\mathrm{with}{p}_g=\frac{m}{J}-\sqrt{\frac{m^2}{J^2}+1}.$$

(14)

DQPTs are observed for $J=m$ at times

$$t_j=\frac{(2j+1)\pi }{2E_{\overline{g}}}=\frac{(2j+1)\pi }{2\sqrt{2}m},$$

(15)

yielding the Rabi oscillations between $|{\psi }_g\rangle$ and $|{\psi }_{\overline{g}}\rangle$ expected from Equation (11), as shown in Figure 2. The behavior of the phase, reported in Figure 2b, features vortices corresponding to Loschmidt echo nodes, while the remaining discontinuities in survival maximum values compensate each other.

2.2. Ground State Preparation

The protocol presented in Section 2.1 requires the preparation of the input state $|{\psi }_g\rangle$, namely the ground state of $H(m,J)$. Based on reasons clarified in Section 2.4, the degrees of freedom of the lattice are assigned to the four qubits of the ibmq_manila circuit $|q_0{q}_1{q}_2{q}_3\rangle$:

• $q_0$ and $q_3$ host the “electric field” states of the ${\mathbb{Z}}_2$ links;
• The staggered spinless fermions are encoded in $q_1$ and $q_2$.

The four physical states refer to the following computational basis states: ${|vac\rangle }_{-}=|1011\rangle$, $|e^{+}{e}^{-}{\rangle }_L=|0101\rangle$, $|e^{+}{e}^{-}{\rangle }_R=|1100\rangle$, and ${|vac\rangle }_{+}=|0010\rangle$, expressed according to the IBM qiskit notation $|\uparrow \rangle =|0\rangle$ and $|\downarrow \rangle =|1\rangle$ for matter sites and $|\frac{1}{2}\rangle =|0\rangle$ and $|-\frac{1}{2}\rangle =|1\rangle$ for links. Since each state can be unambiguously identified by the first two qubits $|q_0{q}_1\rangle$, one can associate with the ground state $|{\psi }_g\rangle$ an auxiliary product state of two qubits:

$$|{\psi }_g^{\prime }\rangle =a_g\left(\right|10\rangle +|00\rangle )+b_g\left(\right|01\rangle +|11\rangle )=\frac{1}{\sqrt{2}}\left(\right|0\rangle +|1\rangle )\otimes \sqrt{2}\left(a_g\right|0\rangle +b_g|1\rangle ),$$

(16)

with the amplitudes corresponding, in the DQPT condition $J/m=1$, to $a_g=0.653$ and $b_g=-0.271$.

The ground state for the complete four-qubit system is obtained by acting with $CNOT$ two-qubit gates:

$$CNO{T}_{ij}|q_i{q}_j\rangle =|q_i,q_i\oplus q_j\rangle ,$$

(17)

which increase the amount of entanglement in the system. For this reason, containing the error probability entailed by these gates is essential to guarantee an effective quantum computation, which cannot be efficiently simulated by classical computers. The circuit chosen for the ground state preparation reads

$$|{\psi }_g\rangle =CNO{T}_{32}CNO{T}_{03}CNO{T}_{13}CNO{T}_{02}{X}_2|{\psi }_g^{\prime }\rangle \otimes |00\rangle ,$$

(18)

and is pictorially represented in Figure 3.

2.3. Noise Models

The simulation of the circuits included an error probability entailed by each gate application, generally described by the bit flip and phase flip error channel $\rho \mapsto \mathcal{D}\left[\rho \right]={\sum }_{i=0}^3{K}_i\rho K_i^{†}$, with

$$K_0=\sqrt{1-p_x-p_y-p_z}𝟙,K_1=\sqrt{p_x}X,K_2=\sqrt{p_y}Y,K_3=\sqrt{p_z}Z.$$

(19)

Such quantum channels, explicitly analyzed in Appendix B, are associated with every single-qubit gate employed both in the state preparation and in its time evolution, while for the two-qubit gates, the independent error probabilities $(p_x,p_y,p_z)$ can be varied in Equation (19) to define

$$\rho \mapsto \tilde{\mathcal{D}}\left[\rho \right]=\sum _{i,j=0}^3{\tilde{K}}_{ij}\rho {\tilde{K}}_{ij}^{†}with{\tilde{K}}_{ij}=K_i\otimes K_j$$

(20)

for a two-qubit density matrix $\rho$. Each circuit includes also reset and measurement gates, which are affected in the simulations only by bit flips [62,63], implemented by a single noise contribution $K_1$, thus corresponding to $p_z=p_y=0$.

The comparison of the simulations with the outputs of IBM Quantum was evaluated in terms of the trace distance:

$$T({\rho }_{\mathrm{ibmq}},{\rho }_{\mathrm{sim}})=\frac{1}{2}\left|\right|{\rho }_{\mathrm{ibmq}}-{\rho }_{\mathrm{sim}}{\left|\right|}_1=\frac{1}{2}\mathrm{Tr}\left[\sqrt{{({\rho }_{\mathrm{ibmq}}-{\rho }_{\mathrm{sim}})}^{†}({\rho }_{\mathrm{ibmq}}-{\rho }_{\mathrm{sim}})}\right],$$

(21)

which quantifies the similarity between the output state of simulated state ${\rho }_{\mathrm{sim}}$ and the actual output of the IBM hardware ${\rho }_{\mathrm{ibmq}}$.

2.4. Trotter Evolution

The evolution determined by the Hamiltonian (6), composed of non-commuting local terms $h_x$, can be approximated by a Trotter decomposition basis on local unitary operators:

$${\mathrm{e}}^{-\mathrm{i}Ht}={\mathrm{e}}^{-\mathrm{i}{\sum }_x{h}_{x}t}={\left({\mathrm{e}}^{-\mathrm{i}{h}_{N-1}\Delta t}{\mathrm{e}}^{-\mathrm{i}{h}_{N-2}\Delta t}\dots {\mathrm{e}}^{-\mathrm{i}{h}_0\Delta t}\right)}^{\frac{t}{\Delta t}}+\mathcal{O}(\Delta t).$$

(22)

The improved approximation that would in principle be provided by the Suzuki–Trotter formula [52] is not well suited in this framework, because it would require a larger number of gates for circuit implementation.

A decomposition of each term in Equation (22) according to the available set of gates was formulated in [52]. Here, we present its specific application to the ${\mathbb{Z}}_2$ gauge group [31,49], where $U_{x,x+1}=U_{x,x+1}^{†}=X_{x,x+1}$. The fermionic hopping contribution can be equivalently expressed as

$$H_J=\frac{J}{4}\sum _{x=0}^{N-1}\left[X_{x,x+1}\left(X_x{X}_{x+1}+Y_x{Y}_{x+1}\right)\right]=\sum _{x=0}^{N-1}{h}_{J,x}.$$

(23)

The evolution related to the Trotter time steps generated by the three qubits’ interaction in Equation (23) was implemented according to the Cartan decomposition [64,65,66]. Concerning the periodic lattice with $N=2$ sites, we considered for clarity the hopping term $h_{J,0}$, acting on the subsystem $|q_0{q}_1{q}_2\rangle$:

$$\begin{array}{cc}\hfill {\mathrm{e}}^{-\mathrm{i}{h}_{J,0}\Delta t}& =K^{†}AK,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill K& =CNO{T}_{12}CNO{T}_{01}{H}_1{H}_{0}CNO{T}_{12},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill A& =𝟙\otimes R_z(J\Delta t/2)\otimes R_z(-J\Delta t/2),\hfill \end{array}$$

(26)

where $H_i$ is the Hadamard gate acting on $q_i$ and $R_z\left(\alpha \right)={\mathrm{e}}^{-\mathrm{i}Z\alpha /2}$. The remaining term $h_{J,1}$ acts in an analogous way on the subsystem $|q_1{q}_2{q}_3\rangle$, as represented in Figure 4.

The decomposition first rotates the product basis states:

with ${|\pm \rangle =|(-1)}^{q_0}\rangle =\frac{1}{\sqrt{2}}\left(\right|0\rangle \pm |1\rangle )$ the X eigenstates. The action of Hadamard gate $H_0$ is required to entangle the state of the associated link with matter sites in the following steps: where the bar stands for the logical NOT. A further elaboration of the above $|Z{Z}_{q_1{q}_2}\rangle$ states expression simplifies the application of Cartan decomposition in the evolution with diagonal operators:

$$\begin{array}{cc}\hfill |Z{Z}_{q_1{q}_2}(\Delta t)\rangle =& A|Z{Z}_{q_1{q}_2}\rangle =𝟙 \otimes R_z(J\Delta t/2)\otimes R_z(-J\Delta t/2)\left[{|(-1)}^{q_1}\pm \rangle \otimes \frac{1}{\sqrt{2}}\left(|0,q_1\oplus q_2{\rangle +(-1)}^{q_1}|1,\overline{q_1\oplus q_2}\rangle \right)\right]\hfill \\ \hfill =& {|(-1)}^{q_1}\pm \rangle \otimes \frac{1}{\sqrt{2}}\left({\mathrm{e}}^{-\mathrm{i}\frac{J\Delta t}{4}(1-{(-1)}^{q_1\oplus q_2})}|0,q_1\oplus q_2{\rangle +(-1)}^{q_1}{\mathrm{e}}^{\mathrm{i}\frac{J\Delta t}{4}(1+{(-1)}^{\overline{q_1\oplus q_2}})}|1,\overline{q_1\oplus q_2}\rangle \right),\hfill \end{array}$$

(27)

such that states satisfying $q_1\oplus q_2=1$ acquire a time-dependent phase, while there is no evolution outside the physical subspace ${\mathcal{H}}_G$ with $q_1\oplus q_2=0$, as predicted by Equation (6).

Focusing on states with $q_1\oplus q_2=1$, the next circuit steps related to $K^{†}$ yield

$$\begin{array}{cc}\hfill |Q_{q_1{q}_2}(\Delta t)\rangle =& CNO{T}_{01}CNO{T}_{12}|Z{Z}_{q_1{q}_2}(\Delta t)\rangle \hfill \\ \hfill =& \frac{1}{2}\left({\mathrm{e}}^{-\mathrm{i}\frac{J\Delta t}{2}}|00\rangle \pm {(-1)}^{q_1}{\mathrm{e}}^{-\mathrm{i}\frac{J\Delta t}{2}}|11\rangle +{(-1)}^{q_1}{\mathrm{e}}^{\mathrm{i}\frac{J\Delta t}{2}}|01\rangle \pm {\mathrm{e}}^{\mathrm{i}\frac{J\Delta t}{2}}|10\rangle \right)\otimes |1\rangle ,\hfill \end{array}$$

(28)

followed by the last part of the decomposition:

$$CNO{T}_{12}{H}_1{H}_0|Q_{q_1{q}_2}(\Delta t)\rangle =\left\{\begin{array}{c}cos\frac{J\Delta t}{2}|001\rangle -\mathrm{i}sin\frac{J\Delta t}{2}|110\rangle ,\mathrm{if}{q}_0=0,q_1=0,\hfill \\ cos\frac{J\Delta t}{2}|101\rangle -\mathrm{i}sin\frac{J\Delta t}{2}|010\rangle ,\mathrm{if}{q}_0=1,q_1=0,\hfill \\ cos\frac{J\Delta t}{2}|010\rangle -\mathrm{i}sin\frac{J\Delta t}{2}|101\rangle ,\mathrm{if}{q}_0=0,q_1=1,\hfill \\ cos\frac{J\Delta t}{2}|110\rangle -\mathrm{i}sin\frac{J\Delta t}{2}|001\rangle ,\mathrm{if}{q}_0=1,q_1=1,\hfill \end{array}\right.$$

(29)

as expected by the action of ${\mathrm{e}}^{-\mathrm{i}{H}_{J,0}\Delta t}$. Depending on states of the remaining link, encoded in $|q_3\rangle$, there are states not belonging to ${\mathcal{H}}_G$ that show a time evolution: they correspond to the ones shown in Figure 1 with a reversed matter site occupation. The remaining contributions in the Trotter expansion of Equation (22) are the diagonal mass terms $H_m={\sum }_{x=0}^{N-1}{h}_{m,x}$ of Equation (6), expressed by

$${\mathrm{e}}^{-\mathrm{i}{h}_{m,x}\Delta t}=R_z(-{(-1)}^{x}m\Delta t)={\mathrm{e}}^{\mathrm{i}{(-1)}^{x}m{Z}_x\Delta t/2},$$

(30)

as reported in Figure 4.

The topology of ibmq_manila in Figure 5a is well suited for the implementation of the Trotter evolution, since every CNOT involves nearest-neighbor qubits. In a noiseless scenario, the presented Trotter evolution would yield the Rabi oscillations in Figure 2c, corresponding to the analytical solution of the evolution by the quenched Hamiltonian.

3. Simulations of Real-Time Dynamics

The experimental results presented in the following were collected in the IBM Quantum platform [1]. The circuit test was based on the simplest periodic lattice required for the implementation of the Schwinger model described in Section 2, composed of $N=2$ sites for the matter field and an equal number of links endowed with the ${\mathbb{Z}}_2$ gauge group. This choice allows for the optimization of the number of gates involved in each Trotter time step, compared to higher-dimensional discretizations of the $\mathrm{U}\left(1\right)$ gauge group [52]. The simulations include a noise model referring to an error probability affecting each gate, but do not take into account effects related to coherence times, as discussed in Section 4.

The ground state preparation procedure presented in Section 2.2 performs better in terms of fidelity than the Python package qiskit built-in command QuantumCircuit. initialize, as shown in Figure 5. Errors due to the use of $CNOT$ gates between non-neighboring qubits were investigated by implementing the same preparation in two different topologies, shown in Figure 5a,b: in ibm_nairobi, the qubits are encoded as follows (see the qubit labels in Panel (b)): $q_0\to “2”$, $q_1\to “0”$, $q_2\to “3”$, $q_3\to “1”$. Actually, Figure 3 shows the presence of three $CNOT$ gates involving $q_3$, together with each one of the remaining qubits. For this reason, this is convenient to encoding it into the highest-degree node of Figure 5b, while the ibm_nairobi circuit is limited to the first four qubits. Despite their structural differences, the median values of the fidelities obtained with the ibmq_manila and ibm_nairobi topologies are essentially the same and are about $0.7$, showing in both cases a much higher efficiency with respect to the implementation of QuantumCircuit.initialize. However, the collected statistics in Figure 5c, referring to 80 runs for each ground state preparation modality, shows that the interquartile range obtained with ibm_nairobi is smaller than the one provided by ibmq_manila. Moreover, fluctuations towards low values in the former case are much less relevant.

The readout of the output states is based on state_tomography_circuits, which exploits for our four-qubit circuit the Pauli basis, resulting in $3^4$ circuits required by the related orthogonal measurements [67]. Simulated noise models include the effects of bit flips in the last measurement part of the circuit [62,63].

The ground state preparation was simulated by defining a noise model in AerSimulator. Three different models were compared with the ibmq_manila output through the trace distance defined in Equation (21) [68]. Each gate appearing in Figure 3 is affected by error probabilities expressed by the error channels (19)–(20). The simulations in Figure 6 use the following models of probability assignment:

(a)

Single- and two-qubit gates share the same probability parameters $(p_x,p_y,p_z)$, generally different along the three axes;

(b)

Single-qubit gates have the same error probability along each noise direction $p_1=p_x=p_y=p_z$; two-qubit gates have an analogous property, but are characterized by an independent probability $p_2$;

(c)

Two parameters $p_1$ and $p_2$ quantify the error probability along both X and Z for single- and two-qubit gates, respectively, while errors along Y are neglected.

Figure 6. Trace distance $T({\rho }_{\mathrm{ibmq}},{\rho }_{\mathrm{sim}})$ of the simulated noise models affecting circuits’ output with respect to the averaged ground state of ibmq_manila by varying the error probabilities. In Panel (a), contour plots in the plane $(p_x,p_z)$ refer to different values of $p_y$. Panels (b,c) evaluate the variation in terms of single- and double-qubit gate error $(p_1,p_2)$ with the inclusion or exclusion of Y errors, respectively.

Figure 6. Trace distance $T({\rho }_{\mathrm{ibmq}},{\rho }_{\mathrm{sim}})$ of the simulated noise models affecting circuits’ output with respect to the averaged ground state of ibmq_manila by varying the error probabilities. In Panel (a), contour plots in the plane $(p_x,p_z)$ refer to different values of $p_y$. Panels (b,c) evaluate the variation in terms of single- and double-qubit gate error $(p_1,p_2)$ with the inclusion or exclusion of Y errors, respectively.

The simulations reported in Figure 6a were averaged over 20 realizations of the noise models, while those in Panels (b)–(c) were obtained from 50 realizations. Panel (a) shows a contour plot corresponding to each value of $p_y$, obtained by interpolating the trace distance evaluated over a grid consisting of $21\times 21$ points with spacing $1\times {10}^{-3}$ in the $(p_x,p_z)$ plane. The considered output of ibmq_manila is a ground state averaged over 80 realizations. The mean value of the trace distance minimum in the variation of $p_y$ is $10.81\times {10}^{-2}\pm 3.3\times {10}^{-3}$, attained at the averaged coordinates in the plane $(p_x,p_z)=((10.7\pm 1.1)\times {10}^{-3},(6.3\pm 0.9)\times {10}^{-3})$. These results are stable with respect to variations in $p_y$. Indeed, a Y error is provided by a sequence of simultaneous bit flips and phase flips, whose probability is much smaller than that of each single error. Panels (b)–(c) show the trace distance of the two noise models with independent parameters for single- and two-qubit gates with respect to the ground state ${\rho }_{\mathrm{ibmq}}$ experimentally prepared 80 times, with no error along Y in the case of Panel (c). The minimum values of the trace distance in the two cases can be considered equal within the statistical fluctuations. Moreover, the minima are found in correspondence with $(p_{1,xyz},p_{2,xyz})\simeq (5\times {10}^{-3},1\times {10}^{-2})$ in Panel (b) and $(p_{1,xz},p_{2,xz})\simeq (7.5\times {10}^{-3},1.5\times {10}^{-2})$ in Panel (c), highlighting a scale factor of $3/2$ due to the absence of the error along Y in the latter case. The minimum values of the trace distance, which equal $9.77\times {10}^{-2}$ in Panel (b) and $11.16\times {10}^{-2}$ in Panel (c), express the ability to distinguish the output ${\rho }_{\mathrm{sim}}$ from the experimental one ${\rho }_{\mathrm{ibmq}}$ approximately once out of ten times.

The ground state in the input evolves according to the quench protocol presented in Section 2.1. Non-commuting local terms in the Hamiltonian, determined by the Jordan–Wigner transformation (6), are circumvented by means of the Trotter evolution described in Section 2.4, which was implemented using the ibmq_manila topology. The Loschmidt echo and the overlaps with the remaining physical states are shown in Figure 7, where the different curves were obtained by varying the time step length $\Delta t$ and by averaging 10 experimental realizations of evolution for each case.

The high error probability translates into a fast convergence towards the maximally mixed state ${\rho }_{\infty }={\left(\dim \mathcal{H}\right)}^{-1}𝟙$. This trend shows a striking deviation for the probability ${\left|\langle {\psi }_e|\psi \left(t\right)\rangle \right|}^2$ during the first two time steps, probably driven by a coherent error accumulation. To focus on this behavior, as well as to limit the computational resources required in simulations, the trace distance is averaged over the first three time steps $t_i=(i-1)\Delta t$, using

$$\overline{T({\rho }_{\mathrm{ibmq}},{\rho }_{\mathrm{sim}})}=\frac{1}{3\Delta t}\sum _{i=1}^{3}T({\rho }_{\mathrm{ibmq}}\left(t_i\right),{\rho }_{\mathrm{sim}}\left(t_i\right))\Delta t,$$

(31)

expressing the mean probability of distinguishing the evolution outputs of ibmq_manila from the simulated ones.

The optimal error probabilities for the ground state preparation are identified in Figure 6c. Simulations of the evolution including noise models fix these parameters assigned to the gates in Figure 3. The Trotter evolution given in Figure 4 includes gates affected by error probabilities, which we implemented according to the following two models:

• Single- and two-qubit gates are characterized by the same arbitrary parameters $(p_x,p_z)$;
• Equal error probabilities in the X and Z directions, but they take generally different values for single-qubit gates ($p_1=p_x=p_z$) and two-qubit gates ($p_2=p_x=p_z$).

The time step that determines the best resolution trade-off in view of investigating the DQPT point is $\Delta t=0.1$. In the averaged trace distance (31), the argument ${\rho }_{\mathrm{ibmq}}$ is averaged over ten experimental realizations of the density matrix evolution, as represented in Figure 7, while ${\rho }_{\mathrm{sim}}$ is averaged over 20 noise realizations for each time step of the evolution. The contour plots in Figure 8 interpolate the evaluation over a grid composed by $31\times 31$ points with spacing ${10}^{-3}$. They describe how distinguishable both aforementioned noise models are with respect to the evolution in ibmq_manila. The minimum value for both models in Panels (a)–(b) is obtained by a further evaluation along the elongated direction, and it is approximately equal to $23.5\times {10}^{-2}$, corresponding to $(p_x,p_z)\simeq (1.1\times {10}^{-2},1.5\times {10}^{-2})$ in the case of Panel (a) and $(p_1,p_2)\simeq (1\times {10}^{-2},1.6\times {10}^{-2})$ in the case of Panel (b). These values are used in Figure 9 to compare the evolution of physical states. In Panels (a)–(d), the overall behavior characterized by a convergence towards the maximally mixed state is captured by the noise models, but non-negligible deviations corresponding to the probability ${\left|\langle {\psi }_e|\psi \left(t\right)\rangle \right|}^2$ in Panel (c) are still not captured by the model, probably because they are driven by coherently accumulated errors, not included for stand-alone gates in current noise models’ implementation. A comparable high number of single- and two-qubit gates makes the two noise models overlap, as represented in Figure 8, so we can focus on the simpler one described by the probabilities $(p_x,p_z)$. Such a probability array is denoted as $\mathbf{p}$ in Figure 9e–h, where the comparison of the reduced noise regimes aimed at estimating a threshold such that a revival is observed after the first DQPT. This translates into the research of a non-monotonic behavior of the Loschmidt echo, which requires an overall error probability per gate ten-times lower than in current implementations.

4. Discussion

Error correction and mitigation techniques are crucial for NISQ devices in order to efficiently simulate the targeted dynamics [34,35,36]. An example was given in [37], concerning the real-time dynamics on IBM Quantum of a periodic lattice model for a $1+1$ QED model with $N=4$ matter sites. There, the discretization of the gauge degrees of freedom was based on a different (non-unitary) truncation. Moreover, the targeted quantities were the vacuum energy and pair production, without a focus on the observation of DQPTs. The exploitation of parity sectors and allowed momenta entails a large reduction of the required qubit number, yielding a scheme able to constrain the evolution in the physical subspace. The reduced Hamiltonian for the matter degrees of freedom in the targeted sector generates a Trotter dynamics implemented in circuits through the Cartan decomposition. Aiming at the zero noise extrapolation, a procedure adopting repeated application of noisy $CNOT$ gates, was implemented. The circuit depth for each time step increases with respect to the decomposition described in this work, thus allowing us to estimate our targeted first DQPT (with time step $\Delta t=0.1$) corresponding to the $T_2$ coherence time, because of the 10-time-step limit mentioned in [37]. Indeed, the total number of circuit moments in our Trotter evolution with 10 steps is equal to 200, as shown in Figure 4. We have to include the ground state preparation depth in Figure 3, because of the $SWAP$ gates related to three $CNOT$ between non-nearest neighbor qubits. The maximal gate temporal extent allowed by the coherence time is slightly lower than the effective one, thus signaling an overestimation of our error probabilities.

The comparison of our results with an ion trap simulation of the lattice Schwinger model with $N=4$ matter sites in [4] has to take into account the much lower number of Mølmer–Sørensen gates, determined by the higher value of the time step. This is related to the different purpose of the aforementioned work, which aimed at characterizing the pair production mechanism.

The evolution in the proximity of a DQPT is considerably more affected by noise [50,51], as simulated for a transverse field Ising model [40] with error rates comparable to those obtained in our study. Nonetheless, the Hamiltonian terms of an Ising model concern at most spin pairs, thus reducing the circuit complexity for Trotter product formulas. In the case of commuting Hamiltonian terms, at the basis of plaquette dynamics without matter degrees of freedom [38], the Trotter product is not required, thus yielding a further reduction of the circuit depth. Concerning the estimated error probabilities in our analysis, their magnitudes were confirmed in the study of scalar Yukawa coupling [69].

The current experimental realization of ion traps and Rydberg atoms in optical lattices shows a higher value for the average gate fidelity of entangling gates [5,70,71]. The introduction of thermal effects, as well as an increased number of parameters for gates’ errors must be considered in order to improve the proposed noise models. Such a detailed description is motivated by the determined “optimal” error probabilities for ground state preparation, which do not coincide with those characterizing time evolution. Moreover, the coherent error accumulation would require more sophisticated noise models, such as correlated dissipation in subsequent quantum channels. The inclusion of error correction and mitigation [37,38,39,40] will be investigated in future research to keep the dynamics in the physical subspace and to balance the noise affecting DQPTs’ observation.

5. Conclusions

We studied the possibility of simulating the real-time dynamics of a model of QED in $1+1$ dimensions, on an elementary lattice composed of two fermionic sites, implemented on IBM Quantum [1]. More specifically, we analyzed the dynamics after a mass quench, close to a dynamical quantum phase transition. The considered quench protocol requires ground state preparation, based on an optimized circuit able to outperform built-in functions, as measured by fidelity with the ideal state. Limitations in observing DQPTs were described in terms of the error probabilities associated with each gate. Different noise models were simulated and compared to capture the main features of the measured evolution, thus determining a marginal contribution of noise along the Y direction. These minimal models revealed the partial observation of the targeted DQPTs’ phenomena in circuit implementations with a reduced error probability. The estimated error rate also indicated a promising implementation on ion traps, such as those available in the IonQ platform [72].