```
import qiskit
from qiskit import QuantumCircuit
from mitiq import zne
from mitiq.interface.mitiq_qiskit.qiskit_utils import initialized_depolarizing_noise
import time
import numpy as np
import matplotlib.pyplot as plt
```

# Use ZNE to simulate quantum many body scars with Qiskit on IBMQ backends#

This tutorial shows how to error mitigate a quantum simulation using ZNE, and is applied to a case where the dynamics show signs of quantum many body scars (QMBS). The example in this tutorial is based on the model studied in *Chen et al. PRR (2022)* [57]. A QMBS is a physical phenomenon where a systems dynamics gives rise to non-ergodic behaviour depending on the initial state. This is in contrast to the ergodic behaviour which is expected from the eigenstate thermalization hypothesis which predicts that a system would thermalize, i.e. a relaxation of the systems dynamics towards a thermal ensemble. For a more elaborate discussion, the interested reader is referred to *Serbyn et al. NatPhys (2021)* [58].

The Hamiltonian that is studied is the mixed-field Ising model (MFIM). Recently, in an analog quantum simulation of the MFIM the quantum many body scar phenomenon was observed experimentally using Rydberg atoms in optical tweezers (see *Bernien et al. Nat (2017)* [59]).

The MFIM Hamiltonian can be written in terms of Pauli matrices as follows

This is an Ising model with an Ising interaction strength \(V\), with longitudinal field with strength proportional to \(V\) and transverse field with strength proportional to \(\Omega\). Note that this Ising chain is defined with open boundary conditions, i.e. the strength of the field at the boundaries (\(i = 1\) and \(i = L\)) is a factor 2 smaller than at the other sites of the chain. More information on the model can be found in the article.

The dynamics of this model is governed by the Schrödinger equation

which can formally be solved as

To simulate the dynamics using a gate sequence one performs a Trotter decomposition of this unitary operator, that is

which is a product of different unitary operators. Finally, one can express each of these unitary operators as a gate sequence of single-qubit gates or two-qubit gates that are subsequently applied.

The resulting gate sequence is defined in the following function

```
def trotter_evolution_H(L: int, V: float,
Omega: float, dt: float) -> qiskit.QuantumCircuit:
'''Return the circuit that performs a time step.
Args:
L: Length of the Ising chain
V: Ising interaction strength
Omega: Transverse field strength
dt: Time step of unitary evolution
'''
cq = qiskit.QuantumCircuit(L)
# Apply Rx gates:
for ii in range(L):
cq.rx(2*Omega*dt, ii)
# Apply Rz gates:
cq.rz(-2*V*dt, 0)
for ii in range(1, L-1):
cq.rz(-4*V*dt, ii)
cq.rz(-2*V*dt, L-1)
# Mitiq ZNE raises an error for the usage of rzz.
# We will give an explicit implementation of the
# 2-CNOT implementation of the Rzz gate:
for ii in range(1, L-1, 2):
cq.cnot(ii, ii+1)
cq.rz(2*V*dt, ii + 1)
cq.cnot(ii, ii+1)
for ii in range(0, L-1, 2):
cq.cnot(ii, ii+1)
cq.rz(2*V*dt, ii + 1)
cq.cnot(ii, ii+1)
return cq
```

By subsequently applying this unitary operator (or its circuit equivalent), one (approximately) time evolves according to the system Hamiltonian.

In this example we will limit ourselves to calculating the behaviour of the staggered magnetization in the \(z\)-direction. This is defined as

The following function calculates the staggered z-magnetization for a given set of raw measurement counts

```
def staggered_mz(L: int, counts: qiskit.result.counts.Counts) -> list:
'''Calculate the staggered z-magnetization for each count
Args:
L: Length of the Ising chain
counts: raw measurement counts
'''
sz = 0
ncounts = 0
sz_list = []
# We use (L-1)-ii since the state strings start with
# the last qubit and end with the first qubit.
stag = (-1)**(1 + np.array([(L-1)-ii for ii in range(L)]))
for state, state_count in counts.items():
bit_array = np.array([-2*int(i)+1 for i in state])
ncounts += state_count
sz_val = np.sum(stag*bit_array)/L
for _ in range(state_count):
sz_list.append(sz_val)
return sz_list
```

For this tutorial we will not use real (quantum) hardware. If you do wish to do so, you can change the backends below to your desired backend.

Note: Using an IBM quantum computer requires a valid IBMQ account. See https://quantum-computing.ibm.com/ for instructions to create an account, save credentials, and see online quantum computers.

```
USE_REAL_HARDWARE = False
```

```
#from qiskit_ibm_provider import IBMProvider
#
#if USE_REAL_HARDWARE:
# provider = IBMProvider(token="MY_IBM_QUANTUM_TOKEN") # Get the API token in https://quantum-computing.ibm.com/account
# backend = provider.get_backend("ibmq_qasm_simulator") # Set quantum computer here!
#else:
# backend = qiskit.Aer.get_backend("qasm_simulator") # Default to a simulator.
```

```
backend = qiskit.Aer.get_backend("qasm_simulator") # Default to a simulator.
```

We set up an executor that simulates the desired circuit a certain amount of shots and returns the measurement statistics of our desired expectation value in the following function

```
def ibmq_executor_full(circuit: qiskit.QuantumCircuit,
shots: int = 8192) -> list:
"""Returns the expectation value of each shot
Args:
circuit: Circuit to run (can also be a list of circuits)
shots: Number of times to execute the circuit to compute
the expectation value.
"""
if NO_NOISE:
# Simulate the circuit with noise
job = qiskit.execute(
experiments=circuit,
backend=qiskit.Aer.get_backend("qasm_simulator"),
optimization_level=0, # Important to preserve folded gates.
shots=shots,
)
else:
if USE_REAL_HARDWARE:
# Run the circuit on hardware
job = qiskit.execute(
experiments=circuit,
backend=backend,
optimization_level=0, # Important to preserve folded gates.
shots=shots
)
else:
# Simulate the circuit with noise
noise_model = initialized_depolarizing_noise(noise_level=0.02)
job = qiskit.execute(
experiments=circuit,
backend=backend,
noise_model=noise_model,
basis_gates=noise_model.basis_gates,
optimization_level=0, # Important to preserve folded gates.
shots=shots,
)
# Convert from raw measurement counts to the expectation value
if type(circuit) == list:
# In case the input is a list of circuits
all_counts = [job.result().get_counts(i) for i in range(len(folded_circuits))]
sz_list = [staggered_mz(L, counts) for counts in all_counts]
else:
# In case the input is a single circuit
counts = job.result().get_counts()
sz_list = staggered_mz(L, counts)
return sz_list
def ibmq_executor(circuit: qiskit.QuantumCircuit, shots: int = 8192) -> float:
"""Returns the expectation value to be mitigated (averaged over the shots).
Args:
circuit: Circuit to run (can also be a list of circuits)
shots: Number of times to execute the circuit to compute the expectation value.
"""
sz_list = ibmq_executor_full(circuit, shots)
if type(sz_list[0]) == list:
expectation_value = [np.mean(kk) for kk in sz_list]
else:
expectation_value = np.mean(sz_list)
return expectation_value
```

We now have all necessary components to perform a time evolution of the Ising Hamiltonian. The Hamiltonian parameters used in the following are those used in the article to create Fig. 3 and 4 of Ref. [1]. However, in the following we will stick to a smaller Ising chain with 6 sites for the sake of simplicity.

As an initial test we will time evolve the system over one time step \(dt\) (you can change this by enlarging the parameter \(n_{dt}\)). The system is initialized in the Néel state, that is the \(\vert 010101\rangle\) state. We simulate results with noise, and calculate the unmitigated and mitigated result. Additionally, we also calculate the result in case no noise is present in the system, that is, an ideal Trotted decomposed time evolution.

```
# System parameters
L = 6
V = 1
Omega = 0.24
dt = 1
n_dt = 11 #number of time steps
```

```
# Initialise quantum circuit
test_circuit = qiskit.QuantumCircuit(L)
# Initialise time step quantum circuit
cc = trotter_evolution_H(L, V, Omega, dt)
# Initialise the Néel state |010101...>
for ii in range(1,L,2):
test_circuit.x(ii)
# Time evolve n_dt time steps
for _ in range(n_dt):
test_circuit = test_circuit.compose(cc)
# Measure in computational basis
test_circuit.measure_all()
```

```
NO_NOISE = False
```

```
t1 = time.time()
unmitigated = ibmq_executor(test_circuit)
t2 = time.time()
print(f"Unmitigated result {unmitigated:.3f}, after {t2-t1:.4f}s.")
```

```
Unmitigated result 0.191, after 0.1905s.
```

```
t1 = time.time()
mitigated = zne.execute_with_zne(test_circuit, ibmq_executor)
t2 = time.time()
print(f"Mitigated result {mitigated:.3f}, after {t2-t1:.4f}s.")
```

```
Mitigated result 0.392, after 1.4475s.
```

```
NO_NOISE = True
t1 = time.time()
unmitigated = ibmq_executor(test_circuit)
t2 = time.time()
print(f"No noise result {unmitigated:.3f}, after {t2-t1:.4f}s.")
```

```
No noise result 0.739, after 0.0935s.
```

The above zero noise extrapolation application uses the default scaling factors [1., 2., 3.] and applies a Richardson extrapolation.

# Low level checks of scale factors and fitting factory (i.e. extrapolation scheme/ fit)#

In this section we will study in more detail the scaling factors used in our ZNE mitigation scheme. We will have a look at the expectation values as a function of various scaling factors to get some clues to decide which noise extrapolation scheme we should use. The obtained results are noisy since there are different sources of error, eg. shot noise and gate noise. When using an extrapolation scheme for the ZNE another source of error is extrapolation uncertainty. To take this noise into account, and the corresponding uncertainty, we also calculate the standard deviation of our sample and the standard error (the uncertainty on the estimated expectation value). We continue with the earlier defined test circuit.

We start with an arbitrary choice of scaling parameters below, and we chose a large number of them to get a grasp of the scaling behavior. Furthermore, will use the random gate folding function to incorporate the scaling factors into the folded circuit. Note that the random gate folding function applies the unitary folding map \(G\rightarrow GG^\dagger G\) to a random subset of gates of the input circuit . In our specific case it is good to note that since the RZZ is explicitely given in terms of the CNOT and RZ gate, these individual gates can/ will also be folded.

The folded circuit is more sensitive to gate errors since it has a number of gates approximately equal to scale_factor * n, where n is the number of gates in the input circuit.

```
scale_factors = [1., 1.5, 2., 2.5, 3., 3.5, 4, 4.5, 5]
folded_circuits = [
zne.scaling.fold_gates_at_random(test_circuit, scale)
for scale in scale_factors
]
# Check that the circuit depth is (approximately) scaled as expected
for j, c in enumerate(folded_circuits):
print(f"Number of gates of folded circuit {j} scaled by: {len(c) / len(test_circuit):.3f}")
```

```
Number of gates of folded circuit 0 scaled by: 0.997
Number of gates of folded circuit 1 scaled by: 1.485
Number of gates of folded circuit 2 scaled by: 1.974
Number of gates of folded circuit 3 scaled by: 2.463
Number of gates of folded circuit 4 scaled by: 2.951
Number of gates of folded circuit 5 scaled by: 3.440
Number of gates of folded circuit 6 scaled by: 3.928
Number of gates of folded circuit 7 scaled by: 4.417
Number of gates of folded circuit 8 scaled by: 4.906
```

We can now calculate the expectation value, standard deviation and standard error for each of these folded circuits and show a plot that will allow us to investigate the behaviour of this circuit for the different scaling factors.

```
NO_NOISE = False
t1 = time.time()
sz_list = ibmq_executor_full(folded_circuits)
t2 = time.time()
exp_vals = [np.mean(kk) for kk in sz_list]
err = [np.std(kk)/np.sqrt(len(kk)) for kk in sz_list]
plt.figure()
plt.errorbar(scale_factors, exp_vals, yerr=err,
linestyle="", marker=".", c="k", ecolor="red")
plt.xlabel(r"scaling factor $\lambda$")
plt.ylabel(r"$<Z_\pi>/L$")
plt.xlim([0, np.max(scale_factors)+0.1])
plt.show()
```

First, we note that the overview of the behavior of the expectation value under influence of the various scaling factors need not be similar for the various circuits used in our quantum simulation. That is, the various circuits used to reach a certain time in the time evolution/ simulation of the system of interest. It can be advisable to check the behavior in more detail for the various circuits one uses to obtain the optimal extrapolation routines. For the sake of simplicity, in the following, we will limit ourselves to one kind of extrapolation scheme for a fixed set of scaling factors.

The above figure suggests that if one decides to perform the extrapolation with fewer scaling factors, e.g. only three, one has to make a proper choice of scaling factors depending on the extrapolation function one chooses. In this particular case, an exponential fit would seem reasonable over this specific range of scaling factors. If one choses for a set of scaling parameters that are all smaller than 2, a linear fit might also yield reasonable results.

Below we choose a set of scaling factors and perform several extrapolation schemes to make an initial choice of extrapolation scheme in the further simulations. First, we choose the scaling factors [1, 1.25, 1.5] and create the corresponding folded gates using random gate folding.

```
scale_factors = [1., 1.25, 1.5]
folded_circuits = [
zne.scaling.fold_gates_at_random(test_circuit, scale)
for scale in scale_factors
]
# Check that the circuit depth is (approximately) scaled as expected
for j, c in enumerate(folded_circuits):
print(f"Number of gates of folded circuit {j} scaled by: {len(c) / len(test_circuit):.3f}")
```

```
Number of gates of folded circuit 0 scaled by: 0.997
Number of gates of folded circuit 1 scaled by: 1.244
Number of gates of folded circuit 2 scaled by: 1.485
```

We can use calculate the expectation values predicted from these folded gates for each scaling factor and use various fitting factories to determine the ZNE expectation value.

```
NO_NOISE = False
t1 = time.time()
sz_vals = ibmq_executor_full(folded_circuits)
exp_vals = [np.mean(kk) for kk in sz_vals]
std_vals = [np.std(kk) for kk in sz_vals]
```

```
zero_noise_value = zne.ExpFactory.extrapolate(scale_factors, exp_vals, asymptote=0.5)
print(f"Extrapolated zero-noise value (exponential factory): {zero_noise_value:.3f}")
```

```
Extrapolated zero-noise value (exponential factory): 0.318
```

```
zero_noise_value = zne.LinearFactory.extrapolate(scale_factors, exp_vals)
print(f"Extrapolated zero-noise value (linear factory): {zero_noise_value:.3f}")
```

```
Extrapolated zero-noise value (linear factory): 0.383
```

```
zero_noise_value = zne.RichardsonFactory.extrapolate(scale_factors, exp_vals)
print(f"Extrapolated zero-noise value (Richardson factory): {zero_noise_value:.3f}")
```

```
Extrapolated zero-noise value (Richardson factory): 1.066
```

```
NO_NOISE = True
t1 = time.time()
unmitigated = ibmq_executor(test_circuit)
t2 = time.time()
print(f"No noise result {unmitigated:.3f}")
```

```
No noise result 0.739
```

These results suggest that for this particular choice of scaling factors (and this particular noisy simulation), the linear extrapolation scheme as well as the Richardson extrapolation scheme are closest to the no-noise result. In the following, we will chose the linear extrapolation scheme and peform perform the earlier quantum simulation once more with the scaling factors [1., 1.25, 1.5].

# Linear factory#

In the following we will chose the linear extrapolation as our fitting factory. The zero noise extrapolated value and its standard error can be straight forwardly calculated (see e.g. *James et al. (2021)* [60]). One can calculate the intercept \(b_0\) and its standard error of a linear fit as follows

where the slope \(b_1\) and its standard error can be calculated as

and where \(\sigma^2\) is the variance of the residual errors \(\epsilon_i = y_i - b_0 - b_1 x_i\) , one can estimate this variance by calculating the residual standard error (RSE)

where RSS is the residual sum of squares given by

Note that \((x_i, y_i)\) are all data points, that is the various predictions of each shot for our expectation values at the different scaling factors. Also note that in the above, strictly speaking one assumes that errors \(\epsilon_i\) are uncorrelated with common variance \(\sigma^2\). Even if this is not exactly valid, the estimates on the standard error are usually good approximations.

In practice, the extrapolation of the intercept is performed using the linear fitting factory. To obtain the standard error on the intercept we write the following function

```
def LinearFit(xvals: list, yvals: list) -> list:
"""Returns the intercept, slope and their standard errors
of a linear fit Y = b0 + b1*X performed on the input data.
Args:
xvals: list of x values of the data set
yvals: list of corresponding y values of the data set
"""
# Transform lists to numpy arrays
xvals = np.array(xvals)
yvals = np.array(yvals)
# Calculate the fitting parameters
b1_num = np.sum((xvals - np.mean(xvals))*(yvals - np.mean(yvals)))
b1_den = np.sum((xvals - np.mean(xvals))**2)
b1 = b1_num/b1_den
b0 = np.mean(yvals) - b1*np.mean(xvals)
# Estimate the standard deviation
e_i = yvals - b0 - b1*xvals
RSS = np.sum(e_i*e_i)
RSE = np.sqrt(RSS/(len(e_i) - 2))
ystd = np.sqrt(RSE)
se_b0 = ystd*np.sqrt(1/len(xvals) + np.mean(xvals)**2/b1_den)
se_b1 = ystd/np.sqrt(b1_den)
# Store the fitting parameters
fit_parameters = [b0, b1, se_b0, se_b1]
return fit_parameters
```

```
# Prepare the x data set and the y data set
xvals = []
yvals = []
for kk in range(len(sz_vals)):
for jj in range(len(sz_vals[kk])):
xvals.append(scale_factors[kk])
yvals.append(sz_vals[kk][jj])
# Perform the linear fit
fit_params = LinearFit(xvals, yvals)
b0, b1, se_b0, se_b1 = fit_params
```

```
print(f"The linear fit yields the intercept b0 = {b0:.3f} with standard error SE(b0) = {se_b0:.3f}.")
```

```
The linear fit yields the intercept b0 = 0.383 with standard error SE(b0) = 0.026.
```

# Create a time evolution plot up to a certain time#

Finally, we will simulate the time evolution of the system up to various times to create a figure of the behaviour of the staggered \(z\)-magnetization throughout time. The system is initialized in the Néel state and we time evolve for a maximum of \(n_{dt} = 40\) steps.

```
def simulate_circuits(circuit: qiskit.QuantumCircuit) -> list:
'''Calculate the expectation value and its standard error of
a given circuit or list of circuits.
Args:
circuit: single circuit or list of (folded) circuits
'''
sz_vals = ibmq_executor_full(circuit)
if type(sz_vals[0]) == list:
# Linear fitting of mitigated value for folded circuits:
exp_vals = [np.mean(kk) for kk in sz_vals]
std_vals = [np.std(kk) for kk in sz_vals]
expectation_value = zne.LinearFactory.extrapolate(scale_factors, exp_vals)
# Prepare the x data set and the y data set
xvals = []
yvals = []
for kk in range(len(sz_vals)):
for jj in range(len(sz_vals[kk])):
xvals.append(scale_factors[kk])
yvals.append(sz_vals[kk][jj])
# Perform the linear fit
fit_params = LinearFit(xvals, yvals)
b0, _, se_value, _ = fit_params
else:
# Calculate mitigated value and its standard
# error for single circuit
expectation_value = np.mean(sz_vals)
se_value = np.std(sz_vals)/np.sqrt(len(sz_vals))
return [expectation_value, se_value]
```

```
# System parameters
L = 6
V = 1
Omega = 0.24
dt = 1
n_dt = 40
unmitigated_list = []
mitigated_list = []
mitigated_alt_list = []
nonoise_list = []
se_unmitigated_list = []
se_mitigated_list = []
se_mitigated_alt_list = []
se_nonoise_list = []
```

```
T1 = time.time()
# Results for t = 0:
# Initialise quantum circuit
circuit = qiskit.QuantumCircuit(L)
# Initialise the Néel state |010101...>
for ii in range(1,L,2):
circuit.x(ii)
# Measure in computational basis
circuit.measure_all()
# Run circuit on the backends
NO_NOISE = False
# Unmitigated simulation
unmitigated, se_unmitigated = simulate_circuits(circuit)
# Mitigated simulation with the custom linear error bars
folded_circuits = [
zne.scaling.fold_gates_at_random(circuit, scale)
for scale in scale_factors
]
mitigated, se_mitigated = simulate_circuits(folded_circuits)
# Mitigated simulation with error bars from repeated simulations
# The advantage of this method is that it can be used for various
# fitting factories. However, one has to repeat the simulation
# multiple times, possibly with a smaller number of shots (for
# simulation time reasons we choose a small number of repetitions
# here)
LinFac = zne.inference.LinearFactory(scale_factors=[1.,1.25,1.5])
mit_list = []
n_runs = 3
for _ in range(n_runs):
mit = zne.execute_with_zne(circuit, ibmq_executor, factory=LinFac)
mit_list.append(mit)
mit_alt, se_mit_alt = [np.mean(mit_list),
np.std(mit_list, ddof=1)/np.sqrt(n_runs)]
NO_NOISE = True
# No noise simulation
nonoise, se_nonoise = simulate_circuits(circuit)
# Store results
unmitigated_list.append(unmitigated)
mitigated_list.append(mitigated)
mitigated_alt_list.append(mit_alt)
nonoise_list.append(nonoise)
se_unmitigated_list.append(se_unmitigated)
se_mitigated_list.append(se_mitigated)
se_mitigated_alt_list.append(se_mit_alt)
se_nonoise_list.append(se_nonoise)
# Results for t > 0:
for ndt in range(1, n_dt + 1):
t1 = time.time()
# Initialise quantum circuit
circuit = qiskit.QuantumCircuit(L)
# Initialise time step quantum circuit
cc = trotter_evolution_H(L, V, Omega, dt)
# Initialise the Néel state |010101...>
for ii in range(1,L,2):
circuit.x(ii)
# Time evolve n_dt time steps
for _ in range(ndt):
circuit = circuit.compose(cc)
# Measure in computational basis
circuit.measure_all()
# Run circuit on the backends
NO_NOISE = False
# Unmitigated simulation
unmitigated, se_unmitigated = simulate_circuits(circuit)
# Mitigated simulation with the custom linear error bars
folded_circuits = [
zne.scaling.fold_gates_at_random(circuit, scale)
for scale in scale_factors
]
mitigated, se_mitigated = simulate_circuits(folded_circuits)
# Mitigated simulation with error bars from repeated simulations
LinFac = zne.inference.LinearFactory(scale_factors=[1.,1.25,1.5])
mit_list = []
for _ in range(n_runs):
mit = zne.execute_with_zne(circuit, ibmq_executor, factory=LinFac)
mit_list.append(mit)
mit_alt, se_mit_alt = [np.mean(mit_list),
np.std(mit_list, ddof=1)/np.sqrt(n_runs)]
NO_NOISE = True
# No noise simulation
nonoise, se_nonoise = simulate_circuits(circuit)
# Store results
unmitigated_list.append(unmitigated)
mitigated_list.append(mitigated)
mitigated_alt_list.append(mit_alt)
nonoise_list.append(nonoise)
se_unmitigated_list.append(se_unmitigated)
se_mitigated_list.append(se_mitigated)
se_mitigated_alt_list.append(se_mit_alt)
se_nonoise_list.append(se_nonoise)
print("Time elapsed: {:.3f}s".format(time.time() - T1))
```

```
Time elapsed: 320.083s
```

```
# Creation of the figure
time_list = np.arange(0, n_dt + 1, 1)
plt.figure()
plt.errorbar(time_list, unmitigated_list, yerr=se_unmitigated_list,
marker =".", c="r", linestyle =":", label = "Unmitigated")
plt.errorbar(time_list, mitigated_list, yerr=se_mitigated_list,
marker = ".", c="C1", linestyle ="--", label = "Mitigated (ZNE)")
plt.errorbar(time_list, mitigated_alt_list, yerr=se_mitigated_alt_list,
marker = ".", c="C0", linestyle ="--", label = "Mitigated (ZNE) alternative")
plt.errorbar(time_list, nonoise_list, yerr=se_nonoise_list,
marker = '.', c="k", label = "Ideal Trotter")
plt.xlim([time_list[0], time_list[-1]])
plt.yticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75],
['-1', '-0.75', '-0.5', '-0.25', '0', '0.25', '0.5', '0.75'])
plt.xlabel(r"$Vt$")
plt.ylabel(r"$<Z_\pi>/L$")
plt.legend()
plt.show()
```