Implementation: trotterized evolution

For our global quench protocol, we first need to prepare the initial state of our system. Since the IBM quantum computers are initialized in the state \(\left|\,\uparrow \uparrow \,\cdots \ \right\rangle\) by default, both of our choices of initial states are tensor product states in the \(z\)-basis and thus can be created by applications of the Pauli \(X\) gate. Next, we need to implement the time evolution, which proceeds by three main steps that we will briefly outline here.

First, we discretize time, that is, we split the time evolution operator, \(\hat{U}(t)={e}^{-i\hat{H}t}\), into a sequence of discrete operators, i.e., \(\hat{U}(t)=\hat{U}(\Delta t)\cdots \hat{U}(\Delta t)\) with fixed \(\Delta t\). Each application of the discrete operator \(\hat{U}(\Delta t)\) is called a Trotter step. This is illustrated in Fig. 1a where we apply an increasing number of Trotter steps to reach later times.

In our simulations we fix \(\Delta t\) to fix the accuracy of our approximation. However, since we can only implement up to \(5\) Trotter steps, we add additional data points by varying \(\Delta t\) in the final Trotter step. More explictly, consider Trotter steps \(M\) and \(M+1\), we can add additional data points by using the evolution operator

$${e}^{-i\hat{H}t}\approx \hat{U}{(\Delta t)}^{M}\hat{U}(\delta t),$$ (10)

where \(\delta t=\Delta t/r\), where \(r\) is the number of data points we want between \(t=M\Delta t\) and \(t=(M+1)\Delta t\). Since the accuracy of the trotter decomposition \(U(\delta t)\) is better than that of \(U(\Delta t)\) (see below), these extra data points will have errors intermediate between that of the \({M}^{\text{th}}\) and \({(M+1)}^{\text{th}}\) steps.

Second, we perform a Trotter decomposition of (trotterize) the unitary evolutions, that is, we approximately decompose the operator \(\hat{U}(\Delta t)\) into a sequence of unitaries that act on at most two neighbouring qubits. In the following we use either a basic or symmetric Trotter decomposition, shown in Fig. 1b, c, respectively. For \(m\) Trotter steps of length \(\Delta t\), the error of the symmetric decomposition is of order \({\mathcal{O}}(m{(\Delta t)}^{3})\), compared with \({\mathcal{O}}(m{(\Delta t)}^{2})\) for the basic decomposition. Note, however, that due to the symmetric structure, when we apply several Trotter steps we can combine several layers of gates. This means that we only need an extra two layers of gates compared with the basic decomposition regardless of the number of Trotter steps.

Third, we must efficiently decompose these two qubit operators into the gates that can be directly implemented on the quantum device, which are the CNOT gate and arbitrary single qubit unitaries. An efficient decomposition is found in ref. 92, which we summarise below. The result is that if \(U\,{

e}\, 0\) then \(\hat{B}\) and \(\hat{C}\) (defined in the figure caption) can be implemented using three CNOTs and with \(U=0\) this can be reduced to only two CNOTs.

Trotter Decomposition

In this paper, we use a Trotter decomposition (commonly known as a trotterization) of the unitary time evolution operator \(\hat{U}(t)={e}^{-i\hat{H}t}\). That is, we want to approximate these operators by a sequence of more easily implemented operators, namely those that act on at most two qubits.

As a starting point, consider a Hamiltonian of the form \(\hat{H}=\hat{A}+\hat{B}\), where \([\hat{A},\hat{B}]\,{

e}\, 0\). Then, since these operators do not commute in general, we have that

$${e}^{-i\hat{H}t}={e}^{-i\hat{A}t}{e}^{-i\hat{B}t}+{\mathcal{O}}({t}^{2}),$$ (11)

which can be naturally extended to Hamiltonians that are sums of more than two terms. To use this fact for trotterized evolution we can use the two steps outlined in the main text. We will now go through the details of these steps in more detail, for the Hamiltonian Eq. (3):

$$\hat{H}=-J\sum _{j=1}^{N-1}\left({\hat{\sigma }}_{j}^{x}{\hat{\sigma }}_{j+1}^{x}+{\hat{\sigma }}_{j}^{y}{\hat{\sigma }}_{j+1}^{y}\right)+U\sum _{j=1}^{N-1}{\hat{\sigma }}_{j}^{z}{\hat{\sigma }}_{j+1}^{z}+\sum _{j=1}^{N}{h}_{j}{\hat{\sigma }}_{j}^{z},$$ (12)

which we rewrite here for convenience.

First, we discretize time and split the evolution operator \(\hat{U}(t)\) into a product of discrete evolution operators \(\hat{U}(\Delta t)\), i.e.,

$${e}^{-i\hat{H}t}={\left({e}^{-i\hat{H}\Delta t}\right)}^{M},$$ (13)

where \(\Delta t=t/M\) and \(M\) is the number of "Trotter steps". Since the Hamiltonian commutes with itself, Eq. (13) is exact. To perform time evolution, we will typically fix \(\Delta t\) and increase the number of Trotter steps \(M\) to reach later times.

The second step is to approximate each of these discrete time operators in a similar manner to Eq. (11). For notational simplicity, let us first define the operators

$${\hat{A}}_{j}={e}^{-i{h}_{j}{\hat{\sigma }}_{j}^{z}\Delta t},\qquad {\hat{B}}_{j}={e}^{-i(U{\hat{\sigma }}_{j}^{z}{\hat{\sigma }}_{j+1}^{z}-J({\hat{\sigma }}_{j}^{x}{\hat{\sigma }}_{j+1}^{x}+{\hat{\sigma }}_{j}^{y}{\hat{\sigma }}_{j+1}^{y}))\Delta t}.$$ (14)

Using these operators, we can make the approximation

$${e}^{-i\hat{H}\Delta t}=\left(\prod _{j}{\hat{A}}_{j}\right)\left(\prod _{j\,\text{even}\,}{\hat{B}}_{j}\right)\left(\prod _{j\,\text{odd}\,}{\hat{B}}_{j}\right)+{\mathcal{O}}({(\Delta t)}^{2}),$$ (15)

which corresponds to the schematic quantum circuit show in Fig. 1b in the main text, and which we will refer to as the basic trotterization. If we wish to evolve to time \(t=M\Delta t\), then we find that the error is \({\mathcal{O}}(\Delta t)\), which is controlled by the size of the Trotter step \(\Delta t\). We can, therefore, improve the accuracy of the approximation by decreasing \(\Delta t\), however, this must be balanced against the cost of needing more Trotter steps, as explained in the main text.

To improve the accuracy of our simulations, we can use better approximations to the discrete evolution operators by way of higher-order Trotter decompositions.93 The leading error term in Eq. (11) is due to the non-zero commutator \([\hat{A},\hat{B}]\). By compensating for this error, we can increase the order of the leading error term.

The only higher-order decomposition that we will consider is the symmetrized Trotter step. Let us again start with the simple case of \(\hat{H}=\hat{A}+\hat{B}\). The symmetric decomposition would then be

$${e}^{-i\hat{H}t}={e}^{-i\frac{t}{2}\hat{A}}{e}^{-it\hat{B}}{e}^{-i\frac{t}{2}\hat{A}}+{\mathcal{O}}({t}^{3}).$$ (16)

The error in this decomposition is of higher-order due to the symmetry which ensures that \(U(-t)=U{(t)}^{\dagger }={U}^{-1}(t)\). This means that the even-order error terms vanish and the leading error is \(\sim \!{(\Delta t)}^{3}\). See ref. 93 for more details and for an iterative method for constructing higher-order decompositions.

For the Hamiltonian Eq. (3) in question, let us again make some definitions to simplify notation

$$\begin{array}{c}{\hat{A}}_{j}={e}^{-i{h}_{j}{\hat{\sigma }}_{j}^{z}\frac{\Delta t}{2}},\quad {\hat{B}}_{j}={e}^{-i(U{\hat{\sigma }}_{j}^{z}{\hat{\sigma }}_{j+1}^{z}-J({\hat{\sigma }}_{j}^{x}{\hat{\sigma }}_{j+1}^{x}+{\hat{\sigma }}_{j}^{y}{\hat{\sigma }}_{j+1}^{y}))\frac{\Delta t}{2}},\\ {\hat{C}}_{j}={e}^{-i(U{\hat{\sigma }}_{j}^{z}{\hat{\sigma }}_{j+1}^{z}-J({\hat{\sigma }}_{j}^{x}{\hat{\sigma }}_{j+1}^{x}+{\hat{\sigma }}_{j}^{y}{\hat{\sigma }}_{j+1}^{y}))\Delta t},\end{array}$$ (17)

then the symmetric Trotter decomposition is

$${e}^{-i\hat{H}\Delta t}=\left(\prod _{j}{\hat{A}}_{j}\right)\left(\prod _{j\,\text{even}\,}{\hat{B}}_{j}\right)\left(\prod _{j\,\text{odd}\,}{\hat{C}}_{j}\right)\left(\prod _{j\,\text{even}\,}{\hat{B}}_{j}\right)\left(\prod _{j}{\hat{A}}_{j}\right)+{\mathcal{O}}({(\Delta t)}^{3}),$$ (18)

which is shown schematically in Fig. 1c of the main text.

Measurement

When making a measurement we will find the system in one of the many-body states \(\left|\lambda \right\rangle\) in this basis, e.g., \(\left|\,\uparrow \downarrow \downarrow \cdots \ \right\rangle\) or \(\left|\,\downarrow \uparrow \downarrow \cdots \ \right\rangle\). By performing multiple runs and measurements, we can approximate the probability of measuring each of the many-body states, that is, we can extract \({\left|{\alpha }_{\lambda }\right|}^{2}\), where \({\alpha }_{\lambda }\) is the probability amplitude for the state \(\left|\lambda \right\rangle\). These probabilities can then be used to construct the observables. To be more concrete, consider the expectation value of the operator \({\hat{\sigma }}_{j}^{z}\) on site \(j\). This is computed as follows,

$$\left\langle \psi (t)\right|{\hat{\sigma }}_{j}^{z}\left|\psi (t)\right\rangle =\sum _{\lambda :j=\uparrow }| {\alpha }_{\lambda }{| }^{2}-\sum _{\lambda ^{\prime} :j=\downarrow }| {\alpha }_{\lambda ^{\prime} }{| }^{2},$$ (19)

where the sums are over all states with the \({j}^{\text{th}}\) spin up or down respectively, and \(| {\alpha }_{\lambda }{| }^{2}\) is the proportion of measurements for which we found the state \(\lambda\). In all the following experiments we will use 8192 measurements per data point, which means that the statistical error for these local correlators is \(\sim 0.01\), which is too small to be included in our figures.

Error mitigation: physical Hilbert space

As we noted in the previous section, due to the presence of conserved quantities in the Hamiltonian evolution, the Hilbert space of states splits into those that are physically allowed by the evolution and those that are not. In particular, the models we consider have conserved net magnetization \({S}_{z}={\sum }_{j}{\hat{\sigma }}_{j}^{z}\). This fact turns out to be advantageous, and allows us to perform rudimentary error mitigation – a term that we use to distinguish it from scalable error correction methods, since in this case we deal only with the lowest order errors. The idea is to simply disregard any measurements for states outside of the physical Hilbert space. Let use present a simplified argument for why this might be a good thing to do, where we will first assume that only bit-flips can occur.

A single bit-flip in the course of the evolution will take us outside of the Hilbert space, and let us denote the probability of this happening as \(\Delta\). However, the lowest order of errors within the physical Hilbert space is \({\Delta }^{2}\), i.e., we need at least two bit-flip errors to get back to the same total magnetization. Hence, by discarding counts outside of the physical Hilbert space we reduce the leading order error to \({\Delta }^{2}\). If the probability, \(\Delta\), is sufficiently small, then we can rely on these perturbative arguments. However, if the error rate is large enough, then multiple bit-flip errors can become significant and the error mitigation will be ineffective.

In Fig. 6b, we show the percentage of measurements that are within the physical Hilbert space, that is, the percentage of measurements that are kept. For the first data points, we are simply measuring the initial state, which has an average measurement fidelity per qubit of \(\sim\!\! 95 \%\) leading to a value of approximately \({(95 \% )}^{6}\approx 70 \%\). At long times the number of retained states stabilises and approaches a value that corresponds to the percentage of states in the total Hilbert space that are physical – in other words, the percentage probability that a randomly selected state is in the physical subspace. Once this number of discarded measurements is reached, the errors are very large and the error mitigation is no longer effective.

Bit-flips are not the only type of errors that could occur. As an example, there are also phase errors, which do not necessarily change the net magnetization. However, we note that the constrained data does typically have improved accuracy, and so we use the restriction to the physical Hilbert space for all our subsequent data.

Quantum Circuits

Here, we will go over some of the details necessary to implement the trotterized evolution operators on the IBM devices. We will only cover those elements of direct relevance to this paper and refer the reader to ref. 1 for an introduction to quantum circuits and quantum computation. We will first introduce the quantum gates that can be implemented on the IBM quantum computers, and then decompose the two qubit unitary operations that appear in our trotterized evolution operator in terms of the elementary one and two qubit gates.

A quantum circuit consists of an array of quantum channels – which represent the physical qubits—and a series of quantum gates that are applied to them. These quantum gates are unitary operators that can be applied to one or more of the qubits. The IBM quantum devices can implement the CNOT gate along with an arbitrary single qubit gate, parametrized by three phases. The cobination of these gates forms a universal set that is, any \(N\) qubit gate can be implemented using a combination of these gates, and in principle an arbitrary quantum computation can be performed.

There is a collection of single particle gates that are useful for writing quantum circuits. We write down a list of the most frequently used gates and how they are implemented on the IBM machines. Consider the computational basis to be the tensor product of single qubit states in the \(z\)-basis, i.e., \(\{\left|\,\uparrow \right\rangle ,\left|\,\downarrow \right\rangle \}\). All matrix forms of the gates are given in this basis and all measurements are made in this \(z\)-basis. It is important to note that gate multiplication reads left to right, whereas matrix multiplication is right to left, i.e. .

Firstly, we have the Pauli matrices, which in the standard quantum information notation are

(20)

Secondly, we have the Hadamard and the \(S\) and \(T\) phase gates,

(21)

And finally, we have the \(X,Y\) and \(Z\) rotation gates,

(22)

which correspond to rotations of the qubit around the \(x,y\) and \(z\) axes respectively. All single qubit gates can be written as a product of these rotation gates, up to a phase. This phase is global and is not measurable and can therefore be omitted. In the IBM machine, all single qubit gates can be directly implemented using

(23)

For slightly less general gates the IBM computer implements either \({U}_{2}(\phi ,\lambda )={U}_{3}(0,\phi ,\lambda )\) or \({U}_{1}(\lambda )={U}_{3}(0,0,\lambda )\), which use fewer physical operations and shorter real time. Before running the circuits, we can combine all strings of single qubit gates into a single one of these three single qubit gates, using the functions available in qiskit.42

The most important two qubit gate for our purposes is the CNOT gate

(24)

This gate flips the second qubit depending on the state of the first. This gate allows the two qubits to become entangled, and combined with general single qubit gates forms a set capable of universal quantum computation, see ref. 1 for a proof. The CNOT is the only multi-qubit gate currently that can be directly implemented on the IBM quantum machines.

Also of interest to us is the reversed CNOT gate

(25)

where we differentiate the reversed CNOT from the CNOT because of the directionality of the IBM machines, i.e., only CNOTs in a given direction can be implemented along the qubit connections. If a CNOT is applied (programmatically) in the wrong directly, the above transformation using Hadamard gates will be applied by qiskit implicitly. Since the single qubit gate fidelities are an order of magnitude better than that of the CNOTs this transformation is not costly, and these additional gates will often be incorporated into other strings of single qubit gates.

Change of basis

When we make a measurement on the quantum machine it is with respect to a given basis, which we take to be the \(z\)-basis. However, we may choose to change the basis for several reasons such as: to measure different operators, to prepare an initial state, or to apply a gate which is more efficiently implemented in a different basis.

We will consider only local changes of basis, i.e. a change of basis for the individual qubits. While a general local change of basis can be implemented using the general single qubit gates above, the most frequently used will be those that change from the Z basis to the X or Y basis.

To change to the X basis, we use the Hadamard gate, \(H\). This implements the transformation

$$Z\to X,\qquad X\to Z,\qquad Y\to -Y.$$ (26)

Note that this mapping is its own inverse, which is a result of the Hadamard gate being both unitary and Hermitian.

To change to the Y basis, we use a combination of the Hadamard and \(S\) gates. We can implement the basis change with the combination \(HSH\), which maps

$$Z\to Y,\qquad Y\to Z,\qquad X\to -X.$$ (27)

Note that this combination of gates is not its own inverse but instead is \(H{S}^{\dagger }H\).

Two qubit gates

The Trotter decomposition allows us to write the general unitary time evolution operator approximately as a product of single and two qubit unitary operators. We, therefore, want to find a way to write a general two qubit unitary in terms of the CNOT gate and single qubit gates that can be applied directly on the IBM devices. The optimal decomposition is found in ref. 92. We briefly review the main results of this paper that are of direct relevance to us.

The optimal decomposition uses the fact that a general matrix in \(U(4)\) can be decomposed as \(U=({A}_{1}\otimes {A}_{2})\cdot N(\alpha ,\beta ,\gamma )\cdot ({A}_{3}\otimes {A}_{4})\),94 where

$$N(\alpha ,\beta ,\gamma )=\exp \left[i\left(\alpha {\sigma }^{x}\otimes {\sigma }^{x}+\beta {\sigma }^{y}\otimes {\sigma }^{y}+\gamma {\sigma }^{z}\otimes {\sigma }^{z}\right)\right].$$ (28)

As a quantum circuit, this can be written as

(29)

This operator \(N(\alpha ,\beta ,\gamma )\) is of direct interest to us for quantum dynamics since it is already of the form required for our Trotter decomposition. This gate can be constructed using a minimum of three CNOTs. The optimal decomposition for \(N(\alpha ,\beta ,\gamma )\) is given by the quantum circuit

(30)

where \(\theta =\frac{\pi }{2}-2\gamma ,\phi =2\alpha -\frac{\pi }{2}\), and \(\lambda =\frac{\pi }{2}-2\beta\). Note that in ref. 92 they use a different sign convention for the rotation gates. Despite the apparent asymmetry of the decomposition, this sequence of gates is symmetric with respect to swapping the two qubits.

For certain cases of our Hamiltonian, namely when \(U=0\), the \(N(\alpha ,\beta ,\gamma )\) gate is more general than we need. By restricting ourselves to \(N(\alpha ,0,\gamma )\) (plus single qubit basis changes) we can reduce the number of CNOTs required to two. This gives us access to matrices of the form

$$N(\alpha ,0,\gamma )=\exp [i(\alpha {\sigma }^{x}\otimes {\sigma }^{x}+\gamma {\sigma }^{z}\otimes {\sigma }^{z})].$$ (31)

We can proceed with the help of the so-called Magic Matrix92,94

(32)

Using this matrix we find \({{\mathcal{M}}}^{\dagger }N(\alpha ,0,\gamma ){\mathcal{M}}={e}^{i\gamma {\sigma }^{z}}\otimes {e}^{i\alpha {\sigma }^{z}}\), which in turn implies that \(N(\alpha ,0,\gamma )={\mathcal{M}}({e}^{i\gamma {\sigma }^{z}}\otimes {e}^{i\alpha {\sigma }^{z}}){{\mathcal{M}}}^{\dagger }\), since \({\mathcal{M}}\) is unitary. As a quantum circuit this can be written as

(33)

This gate can be further simplified by noting that a product of single qubit gates is another single qubit gate. Furthermore, \([S,{R}_{z}(\theta )]=0\), and \(H{R}_{z}(\theta )H={R}_{x}(\theta )\) which gives

(34)

where we have arbitrarily flipped the circuit with respect to the two qubits.

Choosing the best qubits

In all our numerics we used between 6 and 10 of the qubits of the IBM machines, which is only a subset of the available 20 qubits. Hence, we wish to find the best such subset so that we get the most accurate results from the machine. We do this by using a simple iterative algorithm, which we will outline here. We note that "best" is a matter of definition involving the balance of many different parameters. We define best to mean the set of qubits that has the lowest average CNOT errors. We do, however, impose restrictions on the allowed measurement error and T2 times.

The algorithm consists of the following steps:

1. Restrict the set of allowed qubits to \(\tilde{N}\) by discarding any for which the measurement error exceeds a measurement threshold or has a T2 time lower than a given threshold. 2. List all CNOTs that connect the allowed qubits. 3. Set value \(M=N-1\), where \(N\) is the number of qubits. 4. Construct a restricted list of \(M\) CNOTs that have the lowest errors. 5. From this restricted list of CNOTs, iteratively, construct all chains that connect \(N\) qubits. 6. If no such chain is found and \(M\,{<}\,\tilde{N}\), then set \(M=M+1\) and repeat from step 4. If no chain is found and \(M=\tilde{N}\), then add an extra qubit to the allowed list, increasing \(\tilde{N}\to \tilde{N}+1\) by increasing the measurement threshold and repeat from step 2. 7. Once a set of possible chains is found, pick the one with the lowest average CNOT error.

We note that this algorithm is not necessarily efficient or scalable, but it works for the current size of the quantum computers. See refs. 95,96 for a alternative approaches. For NISQ devices, we expect that we will similarly want to pick the best subset of qubits, rather than all of them, to achieve greater accuracy simulations and computations. Therefore, an efficient algorithm for large systems is likely to be important.

Data across different days

One practical issue with using the current IBM machines, is the fact that they are recalibrated on approximately a daily basis. Due to the large errors inherent in the devices, this can have a large effect on the results of a simulation and even mean that a different set of qubits is the "best" on different days. This has the unfortunate consequence that the data obtained from the machine will depend on the day on which it was computed. In this section we show a comparison of the results across three days, demonstrating this variability. Thankfully, we also show that the conclusions that we have made in the previous section remain valid independent of the day on which we performed the computations.

First, in Fig. 7a we show the constrained IBM data for the local density of the end spin—corresponding to that shown in Fig. 2b—obtained across the three days. This demonstrates the large fluctuations in our simulations that are due to the different calibrations of the device. For comparison, we show the results from two runs from the same day separated by approximately 10 h in Fig. 7a. Unlike in Fig. 7a the IBM device was not recalibrated between runs. The difference between the runs (blue squares) is of the order of 5%, which is comparable to the measurement and statistical errors.

Fig. 7 Comparing data from different days. a Simulation of the magnetization on the end spin compared across three consecutive days. b Comparison of two separate runs for the local magnetization on the end site after a quench from a domain wall with \(N=6\) and Hamiltonian (case (i)). The data sets were obtained on 6 May 2019 with approximately 10 h between runs. The blue squares show the difference between runs. In both figures we impose the conservation laws, see main text. Full size image

Fortunately, the day-to-day fluctuations do not affect the qualitative behaviour that is simulated by the machine. For example, in Fig. 8 we compare the data for the spin correlator, corresponding to Fig. 4c of the main text. While some of the behaviour is reproduced in all three subfigures—such as the linear light-cone spreading at short-times—the quality of the simulations in Fig. 8a is quite clearly better. For instance, we see a much clearer linear spread of the correlations, and it is the only subfigure 8a where the correlations convincingly reach the boundary of the system.

Fig. 8 Data computed across three consecutive days for the spin correlator. Corresponding to Fig. 4 of the main text. The subfigures are labelled by the dates on which the simulation was run. Full size image

We also consistently see the correct qualitative behaviour for the spread of particles, \({N}_{\text{half}}\), in Fig. 9. Here, we show the data corresponding to Fig. 3a. While the behaviour differs quantitatively across the days, the qualitative physical behaviour is clear in all plots. This verifies the robustness of the results that we presented above.

Fig. 9 Data computed across three consecutive days. Corresponding to Fig. 3a. The date for each subfigure is shown at the top. Full size image

Finally, in Table 1, we show the different qubits that were chosen on three consecutive days. We also show the minimum, average and maximum values for the readout errors, CNOT errors and T2 decoherence times for the 6-qubit chain. These values were computed using the calibration data provided by IBM. This table shows that across the three day period, the optimal set of qubits can change and the corresponding values can fluctuate significantly.

Table 1 Values for IBM machine for the three consecutive days of simulations. Full size table

While there is clear variation across the three subsequent days that we studied, we point out that the parameters provided by IBM, shown in Table 1, are not necessarily good indicators of the quality of the simulation. For instance, the numbers in the table might suggest that on 12 March 2019 the accuracy of the quantum computer was the worst of the three, yet we find the opposite when considering the data for our simulations, see e.g., Fig. 8. While it is true that the results obtained from the other IBM devices are generally worse, which also reflects a significant difference in gate errors etc., it seems that the small set of numbers in Table 1 is not able to fully characterize the machine.