Confinement and entanglement dynamics on a digital quantum computer

Confinement describes the phenomenon when the attraction between two particles grows with their distance, most prominently found in quantum chromodynamics (QCD) between quarks. In condensed matter physics, confinement can appear in quantum spin chains, for example, in the one dimensional transverse field Ising model (TFIM) with an additional longitudinal field, famously observed in the quantum material cobalt niobate or in optical lattices. Here, we establish that state-of-the-art quantum computers have reached capabilities to simulate confinement physics in spin chains. We report quantitative confinement signatures of the TFIM on an IBM quantum computer observed via two distinct velocities for information propagation from domain walls and their mesonic bound states. We also find the confinement induced slow down of entanglement spreading by implementing randomized measurement protocols for the second order Rényi entanglement entropy. Our results are a crucial step for probing non-perturbative interacting quantum phenomena on digital quantum computers beyond the capabilities of classical hardware.

H = k,n V (n) |k, n k, n| + 2h x cos k 2 |k, n k, n − 1| + |k, n k, n + 1| (1) which, for an infinite chain, can be diagonalised using the transformation |k, α = n C α J n−ν k,α (x k ) |k, n . (2) in which ν k,α = , J is the Bessel function of the first kind and the coefficient C α is used for normalisation (1). Bessels functions can also diagonalise the finite system, (2), but as the boundary effects for our protocols are very small, the analytics for an infinite chain are sufficient. The energy levels E k,α can be computed via the boundary condition that J −ν k,α (x k ) = 0. Continuing with the analogy from QCD, the masses of the mesons formed by the domain wall pairs can be calculated via the difference between energy levels and the ground state. As well as these, the allowed velocities of the mesons can be calculated as the maximal gradient of each energy level, see Fig.1(b) in the main text. Note that, as we only consider the motion of very small mesons (n = 1), we are not troubled by finite size affects when measuring their velocities and can directly compare them to the theoretical predictions which we have also checked via ED (2).
The observation of these quantities then provides direct evidence of the underlying confinement dynamics. For example, in Fig 1 (a) in the main text, with h z = 0 one velocity is seen and corresponds to the free kink motion described by the TFIM. However, if h z = 0.2 there are two velocities. The initial velocity from t = 0 until t ∼ 4J again shows free kink motion, while for times t > 15J the slower velocity of the meson governs the dynamics.
Clearly, the Suzki-Trotter decomposition will only give reliable results for small times.
Thus, in order to simulate long time periods on must use multiple trotter steps, U (t) ∼ U (∆t) U (∆t) U (∆t)..., here U (∆t) is given by the approximation in Eq.(8) The number of trotter steps required depends on the length of time that is being simulated. There are extensions to this approximation to further reduce the resulting error. Building on our previous work, we use the symmetric decomposition (3) given by: Note that e −ibt is not symmetrised in this expression as [b, c] = 0. A schematic of the gate sequence required to implement this is given in Fig.1. Here Velocity extraction.-In order to extract the initial free kink and later meson velocities from the mitigated data obtained by the IBM device the gradient of the light cone formed in ∆ zz i data is computed. Due to the inherent error in the NISQ device, there is a level of uncertainty in the saturation levels that should be used. Thus, a range of velocities are computed and the averages and standard deviations are presented in Fig.1 in the main text.
To attain the initial velocities, data for quench dynamics up to t = 4J, a time before boundary effects present on the kinks, was simulated using four trotter steps. In these results there are initialisation errors that would cover up the free kink velocity unless they are removed. Fig.2 shows ∆ zz i data collected for sites i ∈ {0, 1, 2} after forcing the minimum of each to be zero, this removes the effect of the initialisation error. To calculate a velocity one can choose a saturation level, S, and find the times at which ∆ zz i surpasses this number for each site i ∈ {0, 1, 2}.
This gives three (x, t) coordinates for which the gradient of a line of best fit gives the velocity. When computing the light cone gradients for initial velocities the same range of saturation levels were used for each value of h z for consistency.
It turns out that the meson velocities are much easier to obtain because the meson velocities themselves are observed at a larger scale. Data was collected up to t = 8J, a time such that boundary effects are not present, using seven trotter steps. An initial time of t = 4J, roughly the time at which mesons form, was used as the starting point of the light cone to measure the meson velocities. The light cones produced for a given saturation level, S = 0.21, and the corresponding velocities are shown in Fig.3. Again, the same range of saturation levels was used for each value of h z to obtain an average velocity.  defined as Here, |ψ(t) is the the full many-body wave function of the spin chain at the time t after the global quantum quench and |i are the basis states of the two kink subspace. From this definition of ψ, φ(t) can be seen as a measure of the percentage of the full wave function, |ψ(t) , that remains in the two kink subspace. In Fig.4 it is clear that the two kink subspace is responsible for > 80% of the dynamics in the protocols used in this paper. Hence, by post selecting states measured in the two kink subspace the correct states will be selected up to reasonable error. Beyond this, in Fig.5 a comparison of the results taken from the IBM device before and after error mitigation are presented. This highlights the effect of the subspace projection as well as enforcing the inversion symmetry of the initial state in order to project away background noise and extract the desired physical results. The power of this projection can be understood via the following argument. Let δ be the probability of measurement error for a qubit. For a two kink state there are just four possible erroneous spin flips that do not result in the measurement to be outside the subspace. Therefore, to first order approximation the error that will not be mitigated is just 4δ which does not scale with system size. In turn, the probability of error in simulations that can be mitigated via a projection into the two kink subspace is (N − 4)δ. At second order, the probability of two consecutive errors occurring that take the result out of and then back into the subspace is of order δ 2 . With a small δ this second order process is much less likely, and thus, this projection allows large error mitigation. Here h x = 0.5 and L = 9. As only initial velocities were calculated for h z < 0.3, φ(t) is presented in the range tJ < 4 for these longitudinal field strengths, this is highlighted by the horizontal dashed line. Even in the free kink case, h z = 0, the majority of the dynamics, > 80%, is within the two kink subspace for the times used to compute the velocities in Fig.1 (b) in the main text. Figure 5: A comparison of result collected from the IBM device with and without error mitigation. The IBM results for∆ zz i presented in Fig.2 in the main text with h x = 0.5, h z = 0.5 and L = 9. The bare result are shown first followed by the same data after post selection of the two kink subspace states, and finally the symmerised results. The background noise is considerably reduced by the two kink subspace projection and the light cones are restored by the process of forcing the inversion symmetry of the initial state. This allows the confinement physics to be observed and the gradients of light cones to be measured.
Rényi entropy measurements.-Although the qualitative structure of Rényi entanglement entropy obtained on the IBM machine is correct, in order to obtain quantitative agreement with ED results a constant shift is needed. As the circuit depths used are too low for coherence times to be the source of this problem, the two possible causes are initialisation and measurement errors or gate errors. To understand this we performed two separate tests.
Firstly, we ran a simple circuit in which we directly measure the entropy of the system initialised in the |↑↑↓↓↑↑ state with the expected result of zero. On the IBM machine this is equivalent to starting a run, initialising the state and adding the random single qubit gates on the half chain in consideration and measuring. As errors in single qubit gates are very low, any error in the entropy calculation can be assigned to initialisation and measurement error with a high degree of certainty. Using the built in measurement mitigation tool in Qiskit Ignis the entropy recorded is ∼ 0 as seen in Fig.6. Thus, we conclude that this shift is not an effect of initialisation and measurement error.
Secondly, we ran circuits with a varying number of trotter steps that only evolve the state to very short times, tJ = 0.01. This will also have an expected Rényi entropy of zero to the degree of accuracy which is obtainable by the IBM device. This protocol will include gate errors as well as initialisation and measurement errors. From Fig.6 it can be seen that the error in performing gates on the IBM device has the effect of a constant entropy shift observed for all times measured on the device that grows with the number of trotter steps used. As this shift is not time dependent it is reasonable to assume it is due to errors in implementation only and can be safely removed to observe the true physics. In fact, a basic error model can be derived to account for this shift in random measurement protocols which will be presented elsewhere (4).
As well as removing this shift, in order to obtain clear results for entropy, a system size of six sites was used. It is thus important to rule out finite size effect as the cause of the suppression in entropy. Fig.7 shows the finite size scaling of Rényi entropy. For a system size of six times up to tJ ∼ 2.5 are free of finite size effect. Although for longer times slight effects of a short chain are visible when considering longitudinal field strengths of h z = 0 and 0.5, they are small enough for the halting of entropy growth due to confinement to be visible. Furthermore, there are no finite size effects for a longitudinal field of strength h z = 0.75. Hence, the suppression seen in the IBM device can be confidently assigned to confinement dynamics.  Here, h x = 0.5. Blue results correspond to h z = 0, red results correspond to h z = 0.5 and gold results correspond to h z = 0.75. Although finite size effects are present for a chain of length L = 6 they only have a slight effect for times tJ > 2.5 and are not large enough to hide the halting effect of entropy growth due to confinement.