diff --git a/docs/source/examples/quantum_simulation_1d_ising.md b/docs/source/examples/quantum_simulation_1d_ising.md
index d47ac33c51..093c9f60cc 100644
--- a/docs/source/examples/quantum_simulation_1d_ising.md
+++ b/docs/source/examples/quantum_simulation_1d_ising.md
@@ -4,18 +4,15 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.15.2
+    jupytext_version: 1.14.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
 
-
  # Using ZNE and learning-based methods to mitigate the 1D transverse-longitudinal Ising model
 
-
-
 In this tutorial, we employ ZNE, CDR, and VNCDR mitigation techniques to address errors in the simulation of the 1-D Transverse-Longitudinal Ising model using Mitiq. It is important to note that the results presented here are not original, but rather an attempt to reproduce some of the findings outlined in the paper available at {cite}`Sopena_2023_Quantum`.
 
 One of the primary applications of quantum computers is simulating dynamics in many-body systems. This is particularly significant because as the system size increases, the number of parameters grows exponentially. As a result, classical computers struggle to efficiently simulate such dynamics. However, we are currently in the Noisy Intermediate-Scale Quantum (NISQ) era, which means we lack the necessary resources for fault-tolerant quantum computing. Nevertheless, Quantum Error Mitigation techniques have been developed to address noise using minimal qubit resources. These techniques harness the power of classical computers to handle and mitigate quantum noise. In quantum simulation, our main interest is usually finding the average value of an observable. However, NISQ hardware can only provide us with noisy results. In mitigation techniques, we combine these noisy results with the computational power of classical computers to combat the noise. In this tutorial, we specifically utilize Zero Noise Extrapolation (ZNE), Corrected Dynamical Reduction (CDR), and Variational Noise-Corrected Dynamical Reduction (VNCDR) techniques to mitigate errors in the simulation of a 1-D Ising Hamiltonian.
@@ -25,9 +22,11 @@ The Hamiltonian for the quantum one-dimensional Ising model, with both transvers
 \begin{equation}
 H = H_{ZZ} + H_Z + H_X,
 \end{equation}
+
 \begin{equation}
 H = -J\bigg[\sum_{i=1}^{L-1}Z_iZ_{i+1} + h_Z \sum_{i=1}^{L}Z_i  + h_X\sum_{i=1}^L X_i\bigg].
 \end{equation}
+
 where J is an exchange coupling constant, which sets the microscopic energy scale and $h_X$ and $h_Z$ are the transverse and longitudinal relative field strengths, respectively. This model is integrable for $h_Z = 0$ while for $h_Z \neq 0$ it is only integrable in the continuum when $h_X = 1$.
 The dynamics of this model is governed by the Schrödinger equation
 \begin{equation}
@@ -43,93 +42,70 @@ U(\Delta t) \approx e^{-iH_{ZZ}\Delta t}e^{-iH_{Z}\Delta t}e^{-iH_{X}\Delta t},
 \end{equation}
 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.
 
-
-
 For the first step we import some packages.
 
-```{code-cell}
-
-
-from mitiq.observable import Observable, PauliString
-
-from mitiq import zne
-
-from cirq.ops.common_channels import DepolarizingChannel
-
-from mitiq.interface.mitiq_cirq import compute_density_matrix
-
-import cirq
-
-from mitiq import zne
-
-from mitiq import zne, cdr, Observable, PauliString
-
+```{code-cell} ipython3
 import time
 
+import cirq
+import matplotlib.pyplot as plt
 import numpy as np
+from cirq.ops.common_channels import DepolarizingChannel
 
-import matplotlib.pyplot as plt
+from mitiq import Observable, PauliString, cdr, zne
+from mitiq.interface.mitiq_cirq import compute_density_matrix
+from mitiq.observable import Observable, PauliString
 
 import warnings
-warnings.filterwarnings('ignore')
+warnings.filterwarnings("ignore")
 ```
 
-
-
 To start with coding, we define function "trotter_evolution_H" which is  $U(\Delta t)$
 
-```{code-cell}
-
-
-
-
-'''
+```{code-cell} ipython3
+"""
 in the paper the example trotter circuit in fig. 12 is exp(iH*dt) instead of exp(-iH*dt).
 To be consistent with their result I write the function for exp(iH*dt) too
-'''
+"""
+
 
-def trotter_evolution_H(L: int, h_z: float, h_x: float, Jdt: float) -> cirq.Circuit:
-    '''Return the circuit that performs a time step.
+def trotter_evolution_H(
+    L: int, h_z: float, h_x: float, Jdt: float
+) -> cirq.Circuit:
+    """Return the circuit that performs a time step.
 
     Args:
         L: Length of the Ising chain
         h_z: z interaction strength
         h_x: x interaction strength
         jdt: zz interaction strength time dt
-    '''
+    """
     # First define L qubits
     qubits = cirq.LineQubit.range(L)
     cq = cirq.Circuit()
 
     # Apply Rx gates:
     for ii in range(L):
-        cq.append(cirq.rx(-2*h_x*Jdt).on(qubits[ii]))
-
+        cq.append(cirq.rx(-2 * h_x * Jdt).on(qubits[ii]))
 
     # Apply Rz gates:
     for ii in range(L):
-        cq.append(cirq.rz(-2*h_z*Jdt).on(qubits[ii]))
-
+        cq.append(cirq.rz(-2 * h_z * Jdt).on(qubits[ii]))
 
     # We implement Rzz gate using two CNOT gates
 
-    for ii in range(1, L-1, 2):
-        cq.append(cirq.CNOT(qubits[ii], qubits[ii+1]))
-        cq.append(cirq.rz(-2*Jdt).on(qubits[ii + 1]))
-        cq.append(cirq.CNOT(qubits[ii], qubits[ii+1]))
+    for ii in range(1, L - 1, 2):
+        cq.append(cirq.CNOT(qubits[ii], qubits[ii + 1]))
+        cq.append(cirq.rz(-2 * Jdt).on(qubits[ii + 1]))
+        cq.append(cirq.CNOT(qubits[ii], qubits[ii + 1]))
 
-    for ii in range(0, L-1, 2):
-        cq.append(cirq.CNOT(qubits[ii], qubits[ii+1]))
-        cq.append(cirq.rz(-2*Jdt).on(qubits[ii + 1]))
-        cq.append(cirq.CNOT(qubits[ii], qubits[ii+1]))
+    for ii in range(0, L - 1, 2):
+        cq.append(cirq.CNOT(qubits[ii], qubits[ii + 1]))
+        cq.append(cirq.rz(-2 * Jdt).on(qubits[ii + 1]))
+        cq.append(cirq.CNOT(qubits[ii], qubits[ii + 1]))
     return cq
-
-
-
-
 ```
 
-
 Instead of real hardware, we use a (Mitiq-wrapped) Cirq simulator to obtain the simulation results, and add depolarizing noise.
 
 If you are interested, you can check out the tutorial here to find out how to do the same using IBM hardware
@@ -137,34 +113,28 @@ https://mitiq.readthedocs.io/en/stable/examples/cirq-ibmq-backends.html
 
 In addition, we define an exact simulator as we need to compare the exact and mitigated results.
 
-```{code-cell}
-
-
-# Define an exact simulator
+```{code-cell} ipython3
 def exact_simulator(circuit: cirq.Circuit) -> np.ndarray:
     return compute_density_matrix(circuit, noise_level=(0.0,))
-# We need a simulator for the noise evolution too, in principle this is the real quantum hardware the we use
-def noisy_simulator(circuit: cirq.Circuit) -> np.ndarray:
-    return compute_density_matrix(circuit, DepolarizingChannel ,noise_level=(0.007,))
-```
-
 
 
+def noisy_simulator(circuit: cirq.Circuit) -> np.ndarray:
+    return compute_density_matrix(
+        circuit, DepolarizingChannel, noise_level=(0.007,)
+    )
+```
 
 We set the Hamiltonian parameters and simulate the Hamiltonian using the trotterization technique explained earlier. Then, we plot the average value of different observables. We begin by calculating the average value of $Z_2$ when the initial state consists of all spins up.
 
-```{code-cell}
-
-
+```{code-cell} ipython3
 L = 5
 h_z = 0.9
 h_x = 0.5
 Jdt = 0.5
-n_dt = 6 #number of time steps
+n_dt = 6  # number of time steps
 ```
 
-```{code-cell}
-
+```{code-cell} ipython3
 # Create qubits
 qubits = cirq.LineQubit.range(L)
 
@@ -175,66 +145,69 @@ test_circuit = cirq.Circuit()
 obs = Observable(PauliString("IIZ"))
 
 # we create different lists variables for average value of Z_2 and its mitigated quantities
-unmitigated_measurement=[1]
-exact_measurement=[1]
-mitigated_measurement_cdr=[1]
-mitigated_measurement_vncdr=[1]
-mitigated_measurement_zne=[1]
+unmitigated_measurement = [1]
+exact_measurement = [1]
+mitigated_measurement_cdr = [1]
+mitigated_measurement_vncdr = [1]
+mitigated_measurement_zne = [1]
 
 # repeat the trotter evolution n_dt times and compute compute the unmitigated, exact and mitigated quantities for each step using Mitiq
 for ii in range(n_dt):
-  test_circuit+=trotter_evolution_H(L,h_z,h_x,Jdt)
-  unmitigated_measurement.append(obs.expectation(test_circuit, noisy_simulator).real)
-  exact_measurement.append(obs.expectation(test_circuit,exact_simulator))
-  mitigated_measurement_vncdr.append(cdr.execute_with_cdr(
-    test_circuit,
-    noisy_simulator,
-    observable=obs,
-    simulator=exact_simulator,
-    scale_factors=(1,3)
-).real)
-  mitigated_measurement_cdr.append(cdr.execute_with_cdr(
-    test_circuit,
-    noisy_simulator,
-    observable=obs,
-    simulator=exact_simulator,
-).real)
-  mitigated_measurement_zne.append(zne.execute_with_zne(
-    test_circuit,
-    noisy_simulator,
-    observable=obs,
-).real)
-
-  # this is to keep track of how fast the simulation is
-  print(ii)
-
-
-
-
+    test_circuit += trotter_evolution_H(L, h_z, h_x, Jdt)
+    unmitigated_measurement.append(
+        obs.expectation(test_circuit, noisy_simulator).real
+    )
+    exact_measurement.append(obs.expectation(test_circuit, exact_simulator))
+    mitigated_measurement_vncdr.append(
+        cdr.execute_with_cdr(
+            test_circuit,
+            noisy_simulator,
+            observable=obs,
+            simulator=exact_simulator,
+            scale_factors=(1, 3),
+        ).real
+    )
+    mitigated_measurement_cdr.append(
+        cdr.execute_with_cdr(
+            test_circuit,
+            noisy_simulator,
+            observable=obs,
+            simulator=exact_simulator,
+        ).real
+    )
+    mitigated_measurement_zne.append(
+        zne.execute_with_zne(
+            test_circuit,
+            noisy_simulator,
+            observable=obs,
+        ).real
+    )
+
+    # this is to keep track of how fast the simulation is
+    print(ii)
 ```
 
-
-
 We plot all the exact, unmitigated, and different mitigated measurements to compare the results. As suggested in the original paper {cite}`Sopena_2023_Quantum` , the evolution of the $⟨Z_i(t)⟩$ is as follows
 
 $
-⟨Z_i(t)⟩=a_1 e^{-a_2 t}cos(a_3t)+a_4t+a_5
+⟨Z_i(t)⟩=a_1 e^{-a_2 t}\cos(a_3t)+a_4t+a_5
 $
 
-we used the trotterized simulation results to estimate parametrs $a_1$,...,$a_5$.
-
-```{code-cell}
+we used the trotterized simulation results to estimate parametrs $a_1, \dots ,a_5$.
 
+```{code-cell} ipython3
+import matplotlib.pyplot as plt
 import numpy
 from scipy.optimize import curve_fit
-import matplotlib.pyplot as plt
+
 
 # Define the average value function
 def sigma_z_t(x, a_1, a_2, a_3, a_4, a_5):
     return a_1 * numpy.exp(-a_2 * x) * numpy.cos(a_3 * x) + a_4 * x + a_5
 
+
 # Generate some sample data
-x_data = np.linspace(0, n_dt*Jdt, n_dt+1)
+x_data = np.linspace(0, n_dt * Jdt, n_dt + 1)
 
 # Fit the cosine function to the data
 popt, pcov = curve_fit(sigma_z_t, x_data, exact_measurement)
@@ -243,36 +216,31 @@ popt, pcov = curve_fit(sigma_z_t, x_data, exact_measurement)
 a_1_opt, a_2_opt, a_3_opt, a_4_opt, a_5_opt = popt
 
 # Generate the fitted curve
-x_fit = np.linspace(0, n_dt*Jdt, 300)
+x_fit = np.linspace(0, n_dt * Jdt, 300)
 y_fit = sigma_z_t(x_fit, *popt)
 
 # Plot the  data sets and the fitted curve
-plt.scatter(x_data, unmitigated_measurement, label='Unmitigated')
-plt.scatter(x_data, exact_measurement, label='Exact_Trotterization')
-plt.scatter(x_data, mitigated_measurement_cdr, label='Mitigated_cdr')
-plt.scatter(x_data, mitigated_measurement_vncdr, label='Mitigated_vncdr')
-plt.scatter(x_data, mitigated_measurement_zne, label='Mitigated_zne')
-plt.plot(x_fit, y_fit, 'r-', label='Fitted Curve')
-plt.xlabel('n_dt*Jdt')
-plt.ylabel('$<Z_2>$')
+plt.scatter(x_data, unmitigated_measurement, label="Unmitigated")
+plt.scatter(x_data, exact_measurement, label="Exact_Trotterization")
+plt.scatter(x_data, mitigated_measurement_cdr, label="Mitigated_cdr")
+plt.scatter(x_data, mitigated_measurement_vncdr, label="Mitigated_vncdr")
+plt.scatter(x_data, mitigated_measurement_zne, label="Mitigated_zne")
+plt.plot(x_fit, y_fit, "r-", label="Fitted Curve")
+plt.xlabel(r"$n_dt*Jdt$")
+plt.ylabel(r"$\langle Z_2 \rangle$")
 plt.legend()
 plt.show()
 ```
 
-
-
 We repeat the same thing using a different initial state. this time the initial state is $\vert 00011⟩$
 
-```{code-cell}
-
+```{code-cell} ipython3
 # Create qubits
 qubits = cirq.LineQubit.range(L)
 
 # Create a circuit
 test_circuit = cirq.Circuit()
 
-
-
 # initialize the input state in|00011>
 test_circuit.append(cirq.I(qubits[0]))
 test_circuit.append(cirq.I(qubits[1]))
@@ -284,51 +252,59 @@ test_circuit.append(cirq.X(qubits[4]))
 obs = Observable(PauliString("IIZII"))
 
 # repeat the trotter evolution n_dt times and compute compute the unmitigated, exact and mitigated quantities for each step
-unmitigated_measurement=[1]
-exact_measurement=[1]
-mitigated_measurement_cdr=[1]
-mitigated_measurement_vncdr=[1]
-mitigated_measurement_zne=[1]
+unmitigated_measurement = [1]
+exact_measurement = [1]
+mitigated_measurement_cdr = [1]
+mitigated_measurement_vncdr = [1]
+mitigated_measurement_zne = [1]
 for ii in range(n_dt):
-  test_circuit+=trotter_evolution_H(L,h_z,h_x,Jdt)
-  unmitigated_measurement.append(obs.expectation(test_circuit, noisy_simulator).real)
-  exact_measurement.append(obs.expectation(test_circuit,exact_simulator))
-  mitigated_measurement_cdr.append(cdr.execute_with_cdr(
-    test_circuit,
-    noisy_simulator,
-    observable=obs,
-    simulator=exact_simulator,
-).real)
-  mitigated_measurement_vncdr.append(cdr.execute_with_cdr(
-    test_circuit,
-    noisy_simulator,
-    observable=obs,
-    simulator=exact_simulator,
-    scale_factors=(1,3),
-).real)
-  mitigated_measurement_zne.append(zne.execute_with_zne(
-    test_circuit,
-    noisy_simulator,
-    observable=obs,
-).real)
-
-  # this is to keep track of how fast the simulation is
-  print(ii)
-
+    test_circuit += trotter_evolution_H(L, h_z, h_x, Jdt)
+    unmitigated_measurement.append(
+        obs.expectation(test_circuit, noisy_simulator).real
+    )
+    exact_measurement.append(obs.expectation(test_circuit, exact_simulator))
+    mitigated_measurement_cdr.append(
+        cdr.execute_with_cdr(
+            test_circuit,
+            noisy_simulator,
+            observable=obs,
+            simulator=exact_simulator,
+        ).real
+    )
+    mitigated_measurement_vncdr.append(
+        cdr.execute_with_cdr(
+            test_circuit,
+            noisy_simulator,
+            observable=obs,
+            simulator=exact_simulator,
+            scale_factors=(1, 3),
+        ).real
+    )
+    mitigated_measurement_zne.append(
+        zne.execute_with_zne(
+            test_circuit,
+            noisy_simulator,
+            observable=obs,
+        ).real
+    )
+
+    # this is to keep track of how fast the simulation is
+    print(ii)
 ```
 
-```{code-cell}
-
+```{code-cell} ipython3
 import numpy
 from scipy.optimize import curve_fit
 import matplotlib.pyplot as plt
 
+
 # Define the cosine function
 def func(x, a_1, a_2, a_3, a_4, a_5):
     return a_1 * numpy.exp(-a_2 * x) * numpy.cos(a_3 * x) + a_4 * x + a_5
 
+
 # Generate some sample data
-x_data = np.linspace(0, n_dt*Jdt, n_dt+1)
+x_data = np.linspace(0, n_dt * Jdt, n_dt + 1)
 
 # Fit the cosine function to the data
 popt, pcov = curve_fit(func, x_data, exact_measurement)
@@ -337,43 +313,38 @@ popt, pcov = curve_fit(func, x_data, exact_measurement)
 a_1_opt, a_2_opt, a_3_opt, a_4_opt, a_5_opt = popt
 
 # Generate the fitted curve
-x_fit = np.linspace(0, n_dt*Jdt, 300)
+x_fit = np.linspace(0, n_dt * Jdt, 300)
 y_fit = func(x_fit, *popt)
 
 # Plot the original data and the fitted curve
-plt.scatter(x_data, unmitigated_measurement, label='Unmitigated')
-plt.scatter(x_data, exact_measurement, label='Exact_Trotterization')
-plt.scatter(x_data, mitigated_measurement_cdr, label='Mitigated_cdr')
-plt.scatter(x_data, mitigated_measurement_vncdr, label='Mitigated_vncdr')
-plt.scatter(x_data, mitigated_measurement_zne, label='Mitigated_zne')
-plt.plot(x_fit, y_fit, 'r-', label='Fitted Curve')
-plt.xlabel('n_dt * Jdt')
-plt.ylabel('$<Z_2>$')
+plt.scatter(x_data, unmitigated_measurement, label="Unmitigated")
+plt.scatter(x_data, exact_measurement, label="Exact_Trotterization")
+plt.scatter(x_data, mitigated_measurement_cdr, label="Mitigated_cdr")
+plt.scatter(x_data, mitigated_measurement_vncdr, label="Mitigated_vncdr")
+plt.scatter(x_data, mitigated_measurement_zne, label="Mitigated_zne")
+plt.plot(x_fit, y_fit, "r-", label="Fitted Curve")
+plt.xlabel("n_dt * Jdt")
+plt.ylabel("$<Z_2>$")
 plt.legend()
 plt.show()
 ```
 
-
-
 One can see that for both different initial states the mitigated results perform much better than the unmitigated results.
 
-
 In order to show the temporal evolution of the position of fermions and mesons in this 1-D model, one should measure the probability distribution of kinks {cite}`Sopena_2023_Quantum`
 \begin{equation}
 \Delta_i=\frac{1}{2} (I-Z_i Z_{i+1}).
 \end{equation}
 We assess the effectiveness of various mitigation techniques in the simulation of operatore $Δ_i$
 
-```{code-cell}
-
-
-delta=[1]*4
+```{code-cell} ipython3
+delta = [1] * 4
 
 # define the \delta_i^zz observables
-delta[0] = Observable(PauliString("ZZIII",-1/2),PauliString("I",1/2))
-delta[1] = Observable(PauliString("IZZII",-1/2),PauliString("I",1/2))
-delta[2] = Observable(PauliString("IIZZI",-1/2),PauliString("I",1/2))
-delta[3] = Observable(PauliString("IIIZZ",-1/2),PauliString("I",1/2))
+delta[0] = Observable(PauliString("ZZIII", -1 / 2), PauliString("I", 1 / 2))
+delta[1] = Observable(PauliString("IZZII", -1 / 2), PauliString("I", 1 / 2))
+delta[2] = Observable(PauliString("IIZZI", -1 / 2), PauliString("I", 1 / 2))
+delta[3] = Observable(PauliString("IIIZZ", -1 / 2), PauliString("I", 1 / 2))
 
 # Create qubits
 qubits = cirq.LineQubit.range(L)
@@ -397,52 +368,51 @@ mitigated_measurement_zne = [[1] * (n_dt) for _ in range(4)]
 
 # repeat the trotter evolution n_dt times and compute compute the unmitigated, exact and mitigated quantities for each step
 for ii in range(n_dt):
-  for jj in range(4):
-    unmitigated_measurement[jj][ii]=(delta[jj].expectation(test_circuit, noisy_simulator).real)
-    exact_measurement[jj][ii]=(delta[jj].expectation(test_circuit,exact_simulator))
-    mitigated_measurement_cdr[jj][ii]=(cdr.execute_with_cdr(
-      test_circuit,
-      noisy_simulator,
-      observable=delta[jj],
-      simulator=exact_simulator
-  ).real)
-    mitigated_measurement_vncdr[jj][ii]=(cdr.execute_with_cdr(
-      test_circuit,
-      noisy_simulator,
-      observable=delta[jj],
-      simulator=exact_simulator,
-      scale_factors=(1,3)
-  ).real)
-    mitigated_measurement_zne[jj][ii]=(zne.execute_with_zne(
-      test_circuit,
-      noisy_simulator,
-      observable=delta[jj],
-  ).real)
-  test_circuit+=trotter_evolution_H(L,h_z,h_x,Jdt)
-
-  print(ii)
-
-
+    for jj in range(4):
+        unmitigated_measurement[jj][ii] = (
+            delta[jj].expectation(test_circuit, noisy_simulator).real
+        )
+        exact_measurement[jj][ii] = delta[jj].expectation(
+            test_circuit, exact_simulator
+        )
+        mitigated_measurement_cdr[jj][ii] = cdr.execute_with_cdr(
+            test_circuit,
+            noisy_simulator,
+            observable=delta[jj],
+            simulator=exact_simulator,
+        ).real
+        mitigated_measurement_vncdr[jj][ii] = cdr.execute_with_cdr(
+            test_circuit,
+            noisy_simulator,
+            observable=delta[jj],
+            simulator=exact_simulator,
+            scale_factors=(1, 3),
+        ).real
+        mitigated_measurement_zne[jj][ii] = zne.execute_with_zne(
+            test_circuit,
+            noisy_simulator,
+            observable=delta[jj],
+        ).real
+    test_circuit += trotter_evolution_H(L, h_z, h_x, Jdt)
+
+    print(ii)
 ```
 
-
-
 Then, we plot different simulations using color plots
 
-```{code-cell}
-
+```{code-cell} ipython3
 # Plot exact_measurement using colors
-plt.imshow(np.transpose(exact_measurement), cmap='viridis')
+plt.imshow(np.transpose(exact_measurement), cmap="viridis")
 
 # Invert the y-axis
 plt.gca().invert_yaxis()
 
 # Add a title
-plt.title('Exact Trotterization')
+plt.title("Exact Trotterization")
 
 # Add x and y labels
-plt.xlabel(r'$\Delta_i$')
-plt.ylabel('n_dt')
+plt.xlabel(r"$\Delta_i$")
+plt.ylabel("n_dt")
 
 # Add a colorbar for reference
 plt.colorbar()
@@ -451,17 +421,17 @@ plt.colorbar()
 plt.show()
 
 # Plot unmitigated_measurement using colors
-plt.imshow(np.transpose(unmitigated_measurement), cmap='viridis')
+plt.imshow(np.transpose(unmitigated_measurement), cmap="viridis")
 
 # Invert the y-axis
 plt.gca().invert_yaxis()
 
 # Add a title
-plt.title('Unmitigated')
+plt.title("Unmitigated")
 
 # Add x and y labels
-plt.xlabel(r'$\Delta_i$')
-plt.ylabel('n_dt')
+plt.xlabel(r"$\Delta_i$")
+plt.ylabel("n_dt")
 
 # Add a colorbar for reference
 plt.colorbar()
@@ -471,17 +441,17 @@ plt.show()
 
 
 # Plot mitigated_measurement_zne using colors
-plt.imshow(np.transpose(mitigated_measurement_zne), cmap='viridis')
+plt.imshow(np.transpose(mitigated_measurement_zne), cmap="viridis")
 
 # Invert the y-axis
 plt.gca().invert_yaxis()
 
 # Add a title
-plt.title('Mitigated ZNE')
+plt.title("Mitigated ZNE")
 
 # Add x and y labels
-plt.xlabel(r'$\Delta_i$')
-plt.ylabel('n_dt')
+plt.xlabel(r"$\Delta_i$")
+plt.ylabel("n_dt")
 
 # Add a colorbar for reference
 plt.colorbar()
@@ -491,17 +461,17 @@ plt.show()
 
 
 # Plot mitigated_measurement_cdr using colors
-plt.imshow(np.transpose(mitigated_measurement_cdr), cmap='viridis')
+plt.imshow(np.transpose(mitigated_measurement_cdr), cmap="viridis")
 
 # Invert the y-axis
 plt.gca().invert_yaxis()
 
 # Add a title
-plt.title('Mitigated CDR')
+plt.title("Mitigated CDR")
 
 # Add x and y labels
-plt.xlabel(r'$\Delta_i$')
-plt.ylabel('n_dt')
+plt.xlabel(r"$\Delta_i$")
+plt.ylabel("n_dt")
 
 # Add a colorbar for reference
 plt.colorbar()
@@ -510,17 +480,17 @@ plt.colorbar()
 plt.show()
 
 # Plot mitigated_measurement_vncdr using colors
-plt.imshow(np.transpose(mitigated_measurement_vncdr), cmap='viridis')
+plt.imshow(np.transpose(mitigated_measurement_vncdr), cmap="viridis")
 
 # Invert the y-axis
 plt.gca().invert_yaxis()
 
 # Add a title
-plt.title('Mitigated VNCDR')
+plt.title("Mitigated VNCDR")
 
 # Add x and y labels
-plt.xlabel(r'$\Delta_i$')
-plt.ylabel('n_dt')
+plt.xlabel(r"$\Delta_i$")
+plt.ylabel("n_dt")
 
 # Add a colorbar for reference
 plt.colorbar()
@@ -529,26 +499,28 @@ plt.colorbar()
 plt.show()
 ```
 
-
-
 The plots are consist of different cells, and the color of each cell represents the value of $\Delta_i$ at that particular time.
 
 To enhance the visualization of the performance of various mitigation techniques, the absolute value of the difference between the exact measurement and different mitigation techniques is presented.
 
-```{code-cell}
-
+```{code-cell} ipython3
 # Plot the difference between exact and unmitigated_measurement using colors
-plt.imshow(np.abs(np.transpose(exact_measurement)-np.transpose(unmitigated_measurement)), cmap='viridis')
+plt.imshow(
+    np.abs(
+        np.transpose(exact_measurement) - np.transpose(unmitigated_measurement)
+    ),
+    cmap="viridis",
+)
 
 # Invert the y-axis
 plt.gca().invert_yaxis()
 
 # Add a title
-plt.title('Difference between exact Trotterization and unmitigated')
+plt.title("Difference between exact Trotterization and unmitigated")
 
 # Add x and y labels
-plt.xlabel(r'$\Delta_i$')
-plt.ylabel('n_dt')
+plt.xlabel(r"$\Delta_i$")
+plt.ylabel("n_dt")
 
 # Add a colorbar for reference
 plt.colorbar()
@@ -556,19 +528,24 @@ plt.colorbar()
 # Show the plot
 plt.show()
 
-
 # Plot the difference between exact and mitigated_measurement_zne using colors
-plt.imshow(np.abs(np.transpose(exact_measurement)-np.transpose(mitigated_measurement_zne)), cmap='viridis')
+plt.imshow(
+    np.abs(
+        np.transpose(exact_measurement)
+        - np.transpose(mitigated_measurement_zne)
+    ),
+    cmap="viridis",
+)
 
 # Invert the y-axis
 plt.gca().invert_yaxis()
 
 # Add a title
-plt.title('Difference between exact Trotterization and mitigated ZNE')
+plt.title("Difference between exact Trotterization and mitigated ZNE")
 
 # Add x and y labels
-plt.xlabel(r'$\Delta_i$')
-plt.ylabel('n_dt')
+plt.xlabel(r"$\Delta_i$")
+plt.ylabel("n_dt")
 
 # Add a colorbar for reference
 plt.colorbar()
@@ -576,19 +553,24 @@ plt.colorbar()
 # Show the plot
 plt.show()
 
-
 # Plot the difference between exact and mitigated_measurement_cdr using colors
-plt.imshow(np.abs(np.transpose(exact_measurement)-np.transpose(mitigated_measurement_cdr)), cmap='viridis')
+plt.imshow(
+    np.abs(
+        np.transpose(exact_measurement)
+        - np.transpose(mitigated_measurement_cdr)
+    ),
+    cmap="viridis",
+)
 
 # Invert the y-axis
 plt.gca().invert_yaxis()
 
 # Add a title
-plt.title('Difference between exact Trotterization and mitigated CDR')
+plt.title("Difference between exact Trotterization and mitigated CDR")
 
 # Add x and y labels
-plt.xlabel(r'$\Delta_i$')
-plt.ylabel('n_dt')
+plt.xlabel(r"$\Delta_i$")
+plt.ylabel("n_dt")
 
 # Add a colorbar for reference
 plt.colorbar()
@@ -597,17 +579,23 @@ plt.colorbar()
 plt.show()
 
 # Plot the difference between exact and mitigated_measurement_vncdr using colors
-plt.imshow(np.abs(np.transpose(exact_measurement)-np.transpose(mitigated_measurement_vncdr)), cmap='viridis')
+plt.imshow(
+    np.abs(
+        np.transpose(exact_measurement)
+        - np.transpose(mitigated_measurement_vncdr)
+    ),
+    cmap="viridis",
+)
 
 # Invert the y-axis
 plt.gca().invert_yaxis()
 
 # Add a title
-plt.title('Difference between exact Trotterization and mitigated VNCDR')
+plt.title("Difference between exact Trotterization and mitigated VNCDR")
 
 # Add x and y labels
-plt.xlabel(r'$\Delta_i$')
-plt.ylabel('n_dt')
+plt.xlabel(r"$\Delta_i$")
+plt.ylabel("n_dt")
 
 # Add a colorbar for reference
 plt.colorbar()
@@ -616,5 +604,4 @@ plt.colorbar()
 plt.show()
 ```
 
-
 As one can see, the VNCDR method out performs the other methods and is much better than the unmitigated result.
diff --git a/docs/source/examples/rshadows_tutorial.md b/docs/source/examples/rshadows_tutorial.md
index 6873b39ecd..a0f4470f58 100644
--- a/docs/source/examples/rshadows_tutorial.md
+++ b/docs/source/examples/rshadows_tutorial.md
@@ -10,6 +10,7 @@ kernelspec:
   language: python
   name: python3
 ---
+
 # Robust Shadow Estimation with Mitiq
 
 **Corresponding to:** Min Li (minl2@illinois.edu)
@@ -29,6 +30,7 @@ from mitiq import MeasurementResult
 from mitiq.interface.mitiq_cirq.cirq_utils import (
     sample_bitstrings as cirq_sample_bitstrings,
 )
+
 # set random seed
 np.random.seed(666)
 ```
@@ -47,8 +49,10 @@ file_directory = "./resources"
 if not run_quantum_processing:
     saved_data_name = "saved_data-rshadows"
 
-    with zipfile.ZipFile(f"{file_directory}/rshadows-tutorial-{saved_data_name}.zip") as zf:
-        saved_data = pickle.load(zf.open(f"{saved_data_name}.pkl"))    
+    with zipfile.ZipFile(
+        f"{file_directory}/rshadows-tutorial-{saved_data_name}.zip"
+    ) as zf:
+        saved_data = pickle.load(zf.open(f"{saved_data_name}.pkl"))
 ```
 
 The *robust shadow estimation*{cite}`chen2021robust` approach put forth in *Predicting Many Properties of a Quantum System from Very Few Measurements* {cite}`huang2020predicting` exhibits noise resilience. The inherent randomization in the protocol simplifies the noise, transforming it into a Pauli noise channel that can be characterized relatively straightforwardly. Once the noisy channel $\widehat{\mathcal{M}}$ is characterized $\widehat{\mathcal{M}}$, it is incorporated into the channel inversion $\widehat{\mathcal{M}}^{-1}$, resulting in an unbiased state estimator. The sampling error in the determination of the Pauli channel contributes to the variance of this estimator. 
@@ -71,7 +75,8 @@ qubits: List[cirq.Qid] = cirq.LineQubit.range(num_qubits)
 
 if download_ising_circuits:
     with open(f"{file_directory}/rshadows-tutorial-1D_Ising_g=1_{num_qubits}qubits.pkl", "rb") as file:
-        circuit = pickle.load(file)
+        old_cirq_circuit = pickle.load(file)
+        circuit = cirq.Circuit(old_cirq_circuit.all_operations())
     g = 1
 
 # or user can import from tensorflow_quantum
@@ -102,18 +107,20 @@ def cirq_executor(
     sampler=cirq.Simulator(),
 ) -> MeasurementResult:
     """
-    This function returns the measurement outcomes of a circuit with noisy channel added right before measurement.
+    This function returns the measurement outcomes of a circuit with noisy
+    channel added right before measurement.
     Args:
         circuit: The circuit to execute.
     Returns:
-        A one shot MeasurementResult object containing the measurement outcomes.
+        A one shot MeasurementResult object containing the measurement
+        outcomes.
     """
-    
-    circuit = circuit.copy()
-    qubits = sorted(list(circuit.all_qubits()))
+
+    tmp_circuit = circuit.copy()
+    qubits = sorted(list(tmp_circuit.all_qubits()))
     if noise_level[0] > 0:
         noisy_circuit = cirq.Circuit()
-        operations = list(circuit)
+        operations = list(tmp_circuit)
         n_ops = len(operations)
         for i, op in enumerate(operations):
             if i == n_ops - 1:
@@ -123,16 +130,16 @@ def cirq_executor(
                     )
                 )
             noisy_circuit.append(op)
-        circuit = noisy_circuit
+        tmp_circuit = noisy_circuit
     # circuit.append(cirq.Moment(*noise_model_function(*noise_level).on_each(*qubits)))
     executor = cirq_sample_bitstrings(
-        circuit,
+        tmp_circuit,
         noise_model_function=None,
         noise_level=(0,),
         shots=1,
         sampler=sampler,
     )
-    
+
     return executor
 ```
 
@@ -147,7 +154,6 @@ from mitiq.utils import operator_ptm_vector_rep
 operator_ptm_vector_rep(cirq.I._unitary_() / np.sqrt(2))
 ```
 
-
 ### 2.2 Pauli Twirling of Quantum Channel and Pauli Fidelity:
 The classical shadow estimation involves Pauli twirling of a quantum channel represented by $\mathcal{G} \subset U(d)$, with PTM representation $\mathcal{U}$. This twirling allows direct computation of $\widehat{\mathcal{M}}$ for the noisy channel $\Lambda$:
 \begin{equation}
@@ -242,9 +248,7 @@ plt.xticks(
 plt.ylabel("Pauli fidelity")
 plt.legend()
 ```
-<!--     
-![png](../img/rshadows_compare_channels.png)
-     -->
+
 
 ## 3. Calibration of the operator expectation value estimation
 The expectation value for a series of operators, denoted as $\{O_\iota\}_{\iota\leq M}$, has a snapshot version of expectation value estimation by random Pauli measurement $\widetilde{\mathcal{M}}=\bigotimes_{i}\widetilde{\mathcal{M}}_{P_i}$ and Pauli-twirling calibration $\widehat{\mathcal{M}}^{-1}=\sum_{b\in\{0,1\}^n}f_b^{-1}\bigotimes_{i}\Pi_{b_i\in b}$, which is given by
@@ -263,14 +267,12 @@ Therefore, the final expression of the expectation value estimation can be simpl
 \end{align}
  Additionally, when $P_i =X_i$ (r.e.s.p. $Y_i,\;Z_i$), $U_i^{(2)}$ must correspond to $X$ (r.e.s.p. $Y,\;Z$)-basis measurement to yield a non-zero value, which is easy to check considered that the $P$-basis measurement channel has a PTM rep: $\widetilde{\mathcal{M}}_{P}=\frac{1}{2}(|I\rangle\!\rangle\langle\!\langle I|+|P\rangle\!\rangle\langle\!\langle P|)$, observiously the only measurement that didn't vanished by the operator's $i$-th qubit component in PTM rep: $P\rightarrow \langle\!\langle P|$, is the local $P$-basis measurement.  
 
-
 Next steps are same with classical protocol,with the statistical method of taking an average called "median of means" to achieve an acceptable failure probability of estimation, which need $R_2=N_2K_2$ snapshots, where we use the subscript "2" to denote the index of classical shadow protocol. Acctually,
 \begin{align}
 \hat{o}_\iota(N_2,K_2):=\mathrm{median}\{\hat{o}_\iota^{(1)},\cdots,\hat{o}_\iota^{(K_2)}\}~~\mathrm{where}~~\hat{o}_\iota^{(j)}=N_2^{-1}\sum_{k=N_2(j-1)+1}^{N_2j}\mathrm{Tr}(O_\iota\hat{\rho}_k),\qquad \forall~1\leq j\leq K_2,
 \end{align}
 where we have $K_2$ estimators each of which is the average of $N_2$ single-round estimators $\hat{o}_i^{(j)}$, and take the median of these $K_2$ estimators as our final estimator $\hat{o}_\iota(N_2,K_2)$. We can calculate the median of means of each irreps with projection $\Pi_b=\bigotimes_{i=1}^n\Pi_{b_i}$, 
 
-
 ### 3.1 Ground State Energy Estimation of Ising model with the RSHADOWS algorithm
 
 In this section, we will use the robust shadows estimation algorithm to estimate the ground state energy of the Ising model. We will use the `compare_shadow_methods` function to compare the performance of the robust shadows estimation algorithm and the classical shadows estimation algorithm. The `compare_shadow_methods` function takes the following parameters:
@@ -464,11 +466,6 @@ df_hamiltonian = df_hamiltonian.reset_index()
 noise_model = "bit_flip"
 ```
 
-    /var/folders/fj/9qx1s04j6n713s58ml36xbbm0000gn/T/ipykernel_23400/2684437156.py:1: FutureWarning: The default value of numeric_only in DataFrameGroupBy.sum is deprecated. In a future version, numeric_only will default to False. Either specify numeric_only or select only columns which should be valid for the function.
-      df_hamiltonian = df_energy.groupby(["noise_level", "method"]).sum()
-
-
-
 ```{code-cell} ipython3
 # Define a color palette
 palette = {"exact": "black", "robust": "red", "standard": "green"}
@@ -483,24 +480,17 @@ sns.lineplot(
     markers=True,
     style="method",
     dashes=False,
-    ci=95,
+    errorbar=("ci", 95),
 )
 plt.title(f"Hamiltonian Estimation for {noise_model} Noise")
 plt.xlabel("Noise Level")
-plt.ylabel("Energy Value")
+plt.ylabel("Energy Value");
 ```
-<!-- 
-    
-![png](../img/rshadows_energy.png)
-     -->
-
-
 
 
 ### 3.2 Two Point Correlation Function Estimation with RShadows
 Let's estimate two point correlation fuction: $\langle Z_0 Z_i\rangle$ of a 16-spin 1D Ising model with transverse field on critical point ground state.
 
-
 Import groud state of 1-D Ising model with periodic boundary condition
 
 
@@ -509,7 +499,8 @@ num_qubits = 16
 qubits = cirq.LineQubit.range(num_qubits)
 if download_ising_circuits:
     with open(f"{file_directory}/rshadows-tutorial-1D_Ising_g=1_{num_qubits}qubits.pkl", "rb") as file:
-        circuit = pickle.load(file)
+        old_cirq_circuit = pickle.load(file)
+        circuit = cirq.Circuit(old_cirq_circuit.all_operations())
     g = 1
 else:
     qbs = cirq.GridQubit.rect(num_qubits, 1)
@@ -622,13 +613,9 @@ sns.lineplot(
     markers=True,
     style="method",
     dashes=False,
-    ci=95,
+    errorbar=("ci", 95),
 )
-plt.title(f"Correlation Function Estimation w/ 0.3 Depolarization Noise")
+plt.title("Correlation Function Estimation w/ 0.3 Depolarization Noise")
 plt.xlabel(r"Correlation Function $\langle Z_0Z_i \rangle$")
-plt.ylabel("Correlation")
+plt.ylabel("Correlation");
 ```
-<!-- 
-    
-![png](../img/rshadows_2pt_func.png)
-     -->
diff --git a/mitiq/cdr/clifford_training_data.py b/mitiq/cdr/clifford_training_data.py
index 6be540f3c2..a900e0111d 100644
--- a/mitiq/cdr/clifford_training_data.py
+++ b/mitiq/cdr/clifford_training_data.py
@@ -60,22 +60,18 @@ def generate_training_circuits(
         random_state = np.random.RandomState(random_state)
 
     # Find the non-Clifford operations in the circuit.
-    operations = np.array(list(circuit.all_operations()))
-    non_clifford_indices_and_ops = np.array(
-        [
-            [i, op]
-            for i, op in enumerate(operations)
-            if not cirq.has_stabilizer_effect(op)
-        ]
-    )
+    operations = list(circuit.all_operations())
+    non_clifford_indices_and_ops = [
+        (i, op)
+        for i, op in enumerate(operations)
+        if not cirq.has_stabilizer_effect(op)
+    ]
 
     if len(non_clifford_indices_and_ops) == 0:
         return [circuit] * num_training_circuits
 
-    non_clifford_indices = np.int32(non_clifford_indices_and_ops[:, 0])
-    non_clifford_ops = cast(
-        List[cirq.ops.Operation], non_clifford_indices_and_ops[:, 1]
-    )
+    non_clifford_indices = [i for i, _ in non_clifford_indices_and_ops]
+    non_clifford_ops = [op for _, op in non_clifford_indices_and_ops]
 
     # Replace (some of) the non-Clifford operations.
     near_clifford_circuits = []
@@ -88,7 +84,9 @@ def generate_training_circuits(
             random_state,
             **kwargs,
         )
-        operations[non_clifford_indices] = new_ops
+        for non_clifford_index, new_op in zip(non_clifford_indices, new_ops):
+            operations[non_clifford_index] = new_op
+
         near_clifford_circuits.append(Circuit(operations))
 
     return near_clifford_circuits
diff --git a/requirements.txt b/requirements.txt
index 88aa708379..164376f4ba 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,4 @@
 numpy>=1.22.0
 scipy>=1.5.0,<=1.11.4
-cirq>=1.0.0,<1.3.0
+cirq>=1.0.0,<1.4.0
 tabulate