# High Energy Physics: Implementing \(SU(2)\) Lattice Gauge Theory Trotter Circuits¶

Quantum computers have the ability to tackle problems beyond the reach of classical systems. One such application lies in the realm of quantum field theories, where lattice gauge theories provide a rigorous framework to derive non-perturbative results. Before diving in, you may find it helpful to review the Definitions at the end of this article if you’re unfamiliar with high energy physics.

## Introduction¶

Lattice gauge theory (LGT) is a vital computational tool in quantum field theory, enabling the acquisition of rigorous results without relying on perturbation theory. It’s particularly significant in the context of quantum chromodynamics (QCD), where lattice computations have been instrumental in uncovering new properties of known hadrons, evaluating the likelihood of undiscovered hadrons, and providing critical insights for the search for new physics beyond the standard model.

Typically, LGT studies are conducted in Euclidean spacetime, where time is treated as an imaginary component. This approach is effective for many observables, but it limits the ability to analyze real-time dynamics. However, a Hamiltonian formulation in real-time could overcome this limitation, given adequate computing resources. This is where the potential of future quantum computers becomes particularly intriguing. Quantum computers could handle and evolve quantum states far beyond the capabilities of classical computers, thus offering a promising avenue for exploring real-time dynamics within LGT.

The paper Self-mitigating Trotter circuits for SU(2) lattice gauge theory on a quantum computer by Rahman, Lewis, Mendicelli, and Powell delves into the implementation of SU(2) lattice gauge theory in Minkowski spacetime using quantum computers, showcasing the real-time evolution of excitations across a lattice. While the theoretical foundations laid out in the paper are robust, practical implementation on contemporary quantum hardware requires specialized techniques and considerations.

This knowledge article, complete with detailed code, explores the implementation of methods described in this paper on the Quantinuum H-Series devices. Utilizing TKET, we detail the process of defining the SU(2) 2-Plaquette (2-qubit) Model Hamiltonian, setting up quantum circuits, and employing Trotterization to simulate the quantum dynamics of the system. For a single Trotter step, we simulate the ideal execution of a quantum circuit using a statevector simulator. The two Trotter step case is executed on the Quantinuum H1-2 emulator, and we also demonstrate an example of symbolic circuit construction for these circuits. The final section details the implementation of the SU(2) 5-Plaquette (five-qubit) Model Hamiltonian using the `Qiskit`

library.

## Setup¶

This notebook was generated using `pytket=1.23.0`

, `pytket-quantinuum=0.26.0`

, and `pytket-qiskit=0.47.0`

.

First, we import the libraries and modules we will use:

`Pauli`

utilities and`QubitPauliOperator`

are used for Hamiltonian representation.`Circuit`

utilities from`pytket`

are used for constructing and analyzing quantum circuits.Additional utilities are used for generating term sequences and gate statistics.

Standard libraries are used for numerical operations and visualization.

```
from pytket.pauli import Pauli, QubitPauliString
from pytket.utils import QubitPauliOperator
from pytket.circuit import Circuit, Qubit, CircBox
from pytket.utils import gen_term_sequence_circuit
from pytket.utils.stats import gate_counts
from pytket import OpType
import numpy as np
from scipy.linalg import expm
import matplotlib.pyplot as plt
```

```
# Uncomment to verify installation of pytket extensions: pytket-quantinuum and pytket-qiskit.
# %pip show pytket pytket-quantinuum pytket-qiskit
```

*Note:* Plotting circuits requires an internet connection. An alternative for displaying circuits is the `pytket-offline-display`

package for rendering circuits offline with no internet connection, for details see here.

## Hamiltonian Formulation¶

In quantum mechanics, the Hamiltonian operator plays a pivotal role as it represents the total energy of a system and governs how quantum states evolve with time. The time evolution of a quantum state, denoted as \(|ψ(t)⟩\), is described by the Schrödinger equation:

where \(i\) stands for the imaginary unit, \(\hbar\) denotes the reduced Planck constant, and \(\hat{H}\) signifies the Hamiltonian operator of the system. In the realm of quantum simulations, particularly those involving lattice gauge theories, the Hamiltonian is indispensable. It not only elucidates the system’s dynamics but also forms the foundation for creating quantum circuits that simulate this evolution.

## SU(2) 2-Plaquette Model Hamiltonian¶

The SU(2) 2-plaquette model offers a streamlined representation of lattice gauge theories. The associated Hamiltonian for this model encapsulates the core dynamics of the system, making it apt for simulation on quantum hardware. The Hamiltonian’s general expression is:

where \(J_{ij}\) represents interaction coefficients between lattice sites \(i\) and \(j\), \(\sigma_i^a\) and \(\sigma_j^a\) are the Pauli matrices for site \(i\) or \(j\), with \(a\) taking values \(x\), \(y\) or \(z\), and \(h_i\) are the external field terms for site \(i\).

Here we show the 2-qubit Hamiltonian for the SU(2) 2-plaquette model:

The matrix form represents the action of the Hamiltonian on a 2-qubit system.

The expanded form shows the Hamiltonian in terms of tensor products of Pauli matrices, providing insight into the system’s dynamics and interactions.

This Hamiltonian encapsulates the dynamics of the system and serves as the foundation for subsequent quantum simulations. Specifically, the terms associated with Pauli \(Z\) matrices, like \(Z_0\) and \(Z_1\), capture the on-site energies and interactions between the qubits. The terms involving the Pauli \(X\) matrices, such as \(X_0\) and \(X_1\), represent the interactions influenced by external fields parameterized by \(x\). This Hamiltonian provides a framework for analyzing and simulating the dynamics of the SU(2) 2-plaquette model, making it a cornerstone for the quantum simulations presented in the paper and this notebook.

## Implementation of the SU(2) 2-Plaquette Model Hamiltonian¶

The Hamiltonian for the SU(2) 2-plaquette model is implemented in `pytket`

in this notebook using the `QubitPauliOperator`

to write the operator on qubits and `QubitPauliString`

to represent the individual Pauli string terms. Here we use \(x = 2\).

```
hamiltonian = QubitPauliOperator(
{
QubitPauliString([], []): +7 * 3 / 8,
QubitPauliString([Qubit(0)], [Pauli.Z]): -3 * 3 / 8,
QubitPauliString([Qubit(0), Qubit(1)], [Pauli.Z, Pauli.Z]): -3 / 8,
QubitPauliString([Qubit(1)], [Pauli.Z]): -3 * 3 / 8,
QubitPauliString([Qubit(0)], [Pauli.X]): -3,
QubitPauliString([Qubit(0), Qubit(1)], [Pauli.X, Pauli.Z]): -1,
QubitPauliString([Qubit(1)], [Pauli.X]): -3,
QubitPauliString([Qubit(0), Qubit(1)], [Pauli.Z, Pauli.X]): -1,
}
)
```

We can display the `QubitPauliOperator`

representation of the Hamiltonian:

```
print(hamiltonian)
```

```
{(): 2.62500000000000, (Zq[0]): -1.12500000000000, (Zq[0], Zq[1]): -0.375000000000000, (Zq[1]): -1.12500000000000, (Xq[0]): -3, (Xq[0], Zq[1]): -1, (Xq[1]): -3, (Zq[0], Xq[1]): -1}
```

You can also visualize the non-zero elements of the Hamiltonian matrix using a `spy`

plot from the `matplotlib`

library , shedding light on its non-zero components and overall structure:

```
plt.spy(hamiltonian.to_sparse_matrix())
```

```
<matplotlib.lines.Line2D at 0x1424678d0>
```

### Matrix Formulation¶

We can also convert the Hamiltonian `QubitPauliOperator`

representation to a dense matrix and extract the Hamiltonian matrix elements. This transformation caters to computational tasks and subsequent quantum simulations.

```
hm = hamiltonian.to_sparse_matrix(2).todense()
hm.real
```

```
matrix([[ 0. , -4. , -4. , 0. ],
[-4. , 3. , 0. , -2. ],
[-4. , 0. , 3. , -2. ],
[ 0. , -2. , -2. , 4.5]])
```

### Time Evolution Operator¶

To simulate the time evolution of the quantum system, we compute the matrix exponential of the Hamiltonian matrix. This process represents how the quantum state evolves over a specific time interval. For instance, a factor like \(0.1\) can be interpreted as a discrete time step in the system’s evolution.

In our implementation, we use the `expm`

function to calculate the matrix exponential of the Hamiltonian for a time step of \(0.1\). The time evolution operator ( U(t) ) is expressed as:

For a time step of \(0.1\), this operator becomes:

This calculated operator, `exp_hm`

, will then be applied to simulate the quantum state’s evolution over the time step.

```
# Compute the matrix exponential of the Hamiltonian for the time step of 0.1
exp_hm = expm(-1j * hm * 0.1)
exp_hm
```

```
array([[ 0.84639558+0.0152074j , 0.05508436+0.36836211j,
0.05508436+0.36836211j, -0.07461882+0.01911594j],
[ 0.05508436+0.36836211j, 0.8615547 -0.27351321j,
-0.09378179+0.022007j , 0.06951527+0.17342834j],
[ 0.05508436+0.36836211j, -0.09378179+0.022007j ,
0.8615547 -0.27351321j, 0.06951527+0.17342834j],
[-0.07461882+0.01911594j, 0.06951527+0.17342834j,
0.06951527+0.17342834j, 0.8648741 -0.41980935j]])
```

### Circuit Construction & Visualization¶

After calculating the time evolution operator, the next logical step is to integrate this operator into a quantum circuit. This allows for direct simulations or runs on quantum hardware or simulator. The `pytket`

library provides tools to seamlessly achieve this. By wrapping our time evolution operator in a 2-qubit unitary box, we can then incorporate it into a quantum circuit, setting the stage for subsequent quantum simulations or analyses.

```
# Creating a 2-qubit unitary box that represents the time evolution operator:
from pytket.circuit import Unitary2qBox
u2box = Unitary2qBox(exp_hm)
# Initializing a 2-qubit circuit and embedding the time evolution operator into it:
circ = Circuit(2)
circ.add_unitary2qbox(u2box, 0, 1)
```

```
[Unitary2qBox q[0], q[1]; ]
```

Retrieving the unitary matrix representation of the quantum circuit:

```
circ.get_unitary()
```

```
array([[ 0.84639558+0.0152074j , 0.05508436+0.36836211j,
0.05508436+0.36836211j, -0.07461882+0.01911594j],
[ 0.05508436+0.36836211j, 0.8615547 -0.27351321j,
-0.09378179+0.022007j , 0.06951527+0.17342834j],
[ 0.05508436+0.36836211j, -0.09378179+0.022007j ,
0.8615547 -0.27351321j, 0.06951527+0.17342834j],
[-0.07461882+0.01911594j, 0.06951527+0.17342834j,
0.06951527+0.17342834j, 0.8648741 -0.41980935j]])
```

And obtaining the resulting quantum state after applying the circuit (from a default initial state):

```
circ.get_statevector()
```

```
array([ 0.84639558+0.0152074j , 0.05508436+0.36836211j,
0.05508436+0.36836211j, -0.07461882+0.01911594j])
```

### Decomposition¶

Recognizing the constraints of real quantum hardware, the constructed circuit is decomposed into a set of elementary gates, ensuring its adaptability and compatibility. Specifically, in quantum circuit design, ‘boxed’ operations encapsulate complex actions. To standardize our circuit for broader applications, we’ll replace these boxed operations with sequences of elementary gates using pytket’s `DecomposeBoxes`

.

```
# Importing the necessary tool for gate decomposition:
from pytket.passes import DecomposeBoxes
# Decomposing custom gates in the circuit to ensure compatibility with quantum hardware:
DecomposeBoxes().apply(circ)
```

```
True
```

We can also visualize the quantum circuit using `pytket`

’s rendering tool for Jupyter:

```
from pytket.circuit.display import render_circuit_jupyter
render_circuit_jupyter(circ)
```

### Gate Analysis¶

The quantum circuit, post-decomposition, is dissected to provide a breakdown of its comprising gate types. This analysis gives valuable insights into the circuit’s intricacy and underlying structure.

Insights into the composition and complexity of our circuit can be obtained with the following command:

```
gate_counts(circ)
```

```
Counter({<OpType.TK1: 40>: 4,
<OpType.X: 23>: 2,
<OpType.Z: 22>: 2,
<OpType.Y: 24>: 2,
<OpType.TK2: 41>: 1})
```

The breakdown of gate occurrences in our quantum circuit is as follows:

TK1 gates: 4

X gates: 2

TK2 gates: 1

Z gates: 2

Y gates: 2

This summary serve as a valuable reference when we delve into the topics of recompilation and optimization for execution on H-series devices.

*Note:* See pytket.circuit.OpType for the mathematical descriptions of TK1 and TK2 gates.

### Time evolution circuits¶

The given Hamiltonian, \(-i H t\), needs to be discretized for simulation on a quantum computer. This is where the Trotter-Suzuki approximation comes in handy, allowing us to approximate the time evolution operator \(e^{-i H t}\).

### First-Order Trotter Step¶

In the first-order Trotter step, the time evolution operator is approximated by $\( \left(\prod_j e^{-i H_j \delta t} \right)^{N_t} \)\( where \)\delta t = \frac{t}{N_t}\(. Here, \) N_t \( represents the number of Trotter steps, which is a parameter that determines the granularity of the time evolution approximation. A higher value of \) N_t $ leads to a finer approximation of the time evolution operator.

```
# Setting the maximum time for the simulation
t_max = 1
# Defining the number of Trotter steps for the approximation
n_trotter_steps = 10
# Calculating the time step for each Trotter step
time_step = t_max / n_trotter_steps
# Creating a time line from 0 to t_max with equal intervals for each Trotter step
time_line = np.linspace(0, t_max, n_trotter_steps + 1)
# Printing the calculated time step
print(f"Time step: {time_step}")
# Multiplying the Hamiltonian by the time step to get a slice of the Hamiltonian for each Trotter step
hamiltonian_slice = hamiltonian * time_step
```

```
Time step: 0.1
```

```
print(hamiltonian_slice)
```

```
{(): 0.262500000000000, (Zq[0]): -0.112500000000000, (Zq[0], Zq[1]): -0.0375000000000000, (Zq[1]): -0.112500000000000, (Xq[0]): -0.300000000000000, (Xq[0], Zq[1]): -0.100000000000000, (Xq[1]): -0.300000000000000, (Zq[0], Xq[1]): -0.100000000000000}
```

After slicing the Hamiltonian with the calculated `time_step`

, we convert the `hamiltonian_slice`

into a matrix format. This is done by first converting the Hamiltonian slice into a sparse matrix representation and then converting it to a dense matrix. The `.real`

part extracts the real component of the matrix, which is necessary for some quantum computing simulations and visualizations.

```
hamiltonian_slice.to_sparse_matrix().todense().real
```

```
matrix([[ 1.38777878e-17, -4.00000000e-01, -4.00000000e-01,
0.00000000e+00],
[-4.00000000e-01, 3.00000000e-01, 0.00000000e+00,
-2.00000000e-01],
[-4.00000000e-01, 0.00000000e+00, 3.00000000e-01,
-2.00000000e-01],
[ 0.00000000e+00, -2.00000000e-01, -2.00000000e-01,
4.50000000e-01]])
```

### Constructing the Quantum Circuit for a Single Trotter Step¶

After calculating the Hamiltonian slice for a single time step, we proceed to construct the quantum circuit that simulates this evolution. We use a utility `gen_term_sequence_circuit`

which builds a circuit for us starting from the Hamiltonian we want in the exponent of a single Trotter step, and the initial reference state we want to start from. The pytket API documentation for gen_term_sequence_circuit has more details.

```
# initializing a quantum circuit with two qubits
initial_state = Circuit(2)
# generate the circuit for the Trotter step
trotter_step_circ = gen_term_sequence_circuit(hamiltonian_slice, initial_state)
# visualize the circuit
render_circuit_jupyter(trotter_step_circ)
```

The circuit for a single Trotter step is composed of several `CircBox`

components. The `get_commands()`

method is used to explore the composition of the circuit.

```
trotter_step_circ.get_commands()
```

```
[CircBox q[0], q[1];, CircBox q[0], q[1];, CircBox q[0], q[1];]
```

Each `CircBox`

can be examined to reveal its internal structure, which typically consists of `PauliExpBox`

boxes. We can get the individual `PauliExpBox`

boxes from each `CircBox`

in the circuit

```
pauliexpbox0 = trotter_step_circ.get_commands()[0].op
box_circ = pauliexpbox0.get_circuit()
box_circ.get_commands()
```

```
[PauliExpBox q[0];, PauliExpBox q[0], q[1];, PauliExpBox q[1];]
```

In this nested `CircBox`

there are 3 `PauliExpBox`

, of which only one acts on 2 qubits

```
pauli_exp_box = box_circ.get_commands()[1].op
render_circuit_jupyter(pauli_exp_box.get_circuit())
```

For a more detailed analysis, we use the `DecomposeBoxes`

pass from `pytket.passes`

. This pass decomposes the circuit into fundamental quantum gates, allowing us to see the basic operations that constitute the Trotter step.

For one Trotter step:

```
test_circuit = trotter_step_circ.copy()
DecomposeBoxes().apply(test_circuit)
render_circuit_jupyter(test_circuit)
```

and we can count the total number of gates

```
gate_counts(test_circuit)
```

```
Counter({<OpType.H: 33>: 8, <OpType.Rz: 36>: 7, <OpType.CX: 42>: 6})
```

or look at specific statistics of the circuit

```
print(f"Number of gates = {test_circuit.n_gates}")
print(f"Circuit depth = {test_circuit.depth()}")
print(f"Number of CX gates = {test_circuit.n_gates_of_type(OpType.CX)}")
```

```
Number of gates = 21
Circuit depth = 19
Number of CX gates = 6
```

## Example 1: Noiseless emulation of a circuit for one Trotter step¶

We will first run a noiseless statevector emulation of the circuit on Quantinuum’s System Model H1 emulator. Quantinuum H-Series emulators model the physical and noise parameters of H-Series quantum computers. The emulators provide a way to turn the noise model off, while still modeling the physical properties of the device, such as ion transport.

```
from pytket.extensions.quantinuum import QuantinuumBackend
h_backend = QuantinuumBackend("H1-2E") # emulator
h_backend.login()
```

### One Trotter step¶

The Trotter step circuit can be enclosed in a `CircBox`

and added to any circuit (this way we can add many steps in a sequence).

```
trotter_step_box = CircBox(trotter_step_circ)
# define circuit for one Trotter step
circ = Circuit(2)
circ.add_gate(trotter_step_box, circ.qubits)
render_circuit_jupyter(circ)
```

When we compile the circuit for the backend, the `CircBox`

nested boxes get decomposed into native gates. Each backend in `pytket`

has a list of default *passes* that are applied to it. They can be applied by calling `backend.get_compiled_circuit()`

with different `optimisation_level`

arguments (the pytket-quantinuum docs list the default values). Here we use the default `optimisation_level`

for the `QuantinuumBackend`

, which is 2.

*Note:* The measurement operations from the `MeasurementSetup`

object are appended to the numerical circuit. Once this step is complete, the circuit is ready for submission.

```
circ_compiled = h_backend.get_compiled_circuit(circ.measure_all())
render_circuit_jupyter(circ_compiled)
```

Before submitting to the emulator, the total cost of running the set of circuits in H-System Quantum Credits (HQCs) can be checked beforehand.

```
n_shots = 100
cost = h_backend.cost(circ_compiled, n_shots=n_shots, syntax_checker="H1-2SC")
print("Cost of experiment in HQCs:", cost)
```

```
Cost of experiment in HQCs: 5.92
```

We submit the circuit to the Quantinuum H1-2 Emulator and make the emulation noiseless by passing the option `noisy_simulation=False`

.

```
handle = h_backend.process_circuit(
circ_compiled, n_shots=n_shots, noisy_simulation=False
)
```

```
circuit_status = h_backend.circuit_status(handle)
print(circuit_status)
```

```
StatusEnum.COMPLETED
```

```
result = h_backend.get_result(handle)
result.get_counts()
```

```
Counter({(0, 0): 39, (1, 0): 35, (0, 1): 24, (1, 1): 2})
```

### Example 2: Noisy emulation on the System Model H1 Emulator for two Trotter steps¶

Next we implement and run two Trotter steps on the H1-2E emulator with the emulator’s default noise model for two Trotter steps. The noise model mimicks the noise of the H1-2 quantum computer.

Let’s add two trotter steps to our circuit:

```
# two steps
circ = Circuit(2)
circ.add_gate(trotter_step_box, circ.qubits)
circ.add_gate(trotter_step_box, circ.qubits)
```

```
[CircBox q[0], q[1]; CircBox q[0], q[1]; ]
```

```
circ_compiled = h_backend.get_compiled_circuit(circ.measure_all())
```

```
print(
f"Our circuit has {circ_compiled.n_gates} quantum gates in total and an overall depth of {circ_compiled.depth()}."
)
print(
f"Of these gates {circ_compiled.n_gates_of_type(OpType.ZZPhase)} are two qubit gate, ZZPhase, counts."
)
```

```
Our circuit has 13 quantum gates in total and an overall depth of 8.
Of these gates 3 are two qubit gate, ZZPhase, counts.
```

```
render_circuit_jupyter(circ_compiled)
```

Before submitting to the emulator, the total cost of running the set of circuits in H-System Quantum Credits (HQCs) is checked beforehand.

```
n_shots = 100
cost = h_backend.cost(circ_compiled, n_shots=n_shots, syntax_checker="H1-2SC")
print("Cost of experiment in HQCs:", cost)
```

```
Cost of experiment in HQCs: 6.16
```

Now we run the circuits.

```
handle = h_backend.process_circuit(circ_compiled, n_shots=n_shots)
```

The status of the jobs can be checked with `ciruit_status`

method. This method requires the `ResultHandle`

to be passed as input. In this example, the job has completed and the results are reported as being ready to request.

```
circuit_status = h_backend.circuit_status(handle)
print(circuit_status.status)
```

```
StatusEnum.COMPLETED
```

Once a job’s status returns completed, results can be returned using the `get_result`

function in conjunction with the `get_counts`

call.

```
result = h_backend.get_result(handle)
result.get_counts()
```

```
Counter({(1, 0): 47, (1, 1): 31, (0, 1): 19, (0, 0): 3})
```

We note that when using `optimisation_level`

equal to 0, we compile to the native gates of the H-Series without optimization and the number of 1- and 2-qubit gates is significantly larger than compared to the default `optimisation_level`

of 2.

```
circ_compiled_0 = h_backend.get_compiled_circuit(circ, optimisation_level=0)
gate_counts(circ_compiled_0)
```

```
Counter({<OpType.Rz: 36>: 78,
<OpType.PhasedX: 68>: 40,
<OpType.ZZPhase: 73>: 12,
<OpType.Measure: 63>: 2})
```

## Example 3: Symbolic Manipulation in Quantum Circuit Simulation¶

In this section, we explore the use of symbolic manipulation for quantum circuit simulation, particularly focusing on the Trotter-Suzuki decomposition. We employ `sympy`

for defining symbolic variables, which represent key parameters like coupling strengths and time steps in our quantum system. These symbols allow us to dynamically adjust our circuit parameters.

The symbols are related to the coupling and the time step.

```
from sympy import symbols
syms = symbols("x dt")
```

We define various `QubitPauliString`

objects to represent different Pauli operators acting on our qubits.

```
zi = QubitPauliString([Qubit(0)], [Pauli.Z])
iz = QubitPauliString([Qubit(1)], [Pauli.Z])
xi = QubitPauliString([Qubit(0)], [Pauli.X])
ix = QubitPauliString([Qubit(1)], [Pauli.X])
zz = QubitPauliString([Qubit(0), Qubit(1)], [Pauli.Z, Pauli.Z])
xz = QubitPauliString([Qubit(0), Qubit(1)], [Pauli.X, Pauli.Z])
zx = QubitPauliString([Qubit(0), Qubit(1)], [Pauli.Z, Pauli.X])
```

Next, we define two dictionaries, `singles_syms`

and `doubles_sym`

, to represent parts of our Hamiltonian symbolically:

`singles_syms`

: Maps 1-qubit Pauli strings to their coefficients. These terms, involving`Z`

and`X`

Pauli operators on individual qubits, reflect the linear contributions in the Hamiltonian. For example, terms like`zi`

and`xi`

are mapped to symbolic expressions involving`dt`

and`x`

, representing time step and coupling constants, respectively.`doubles_sym`

: Represents 2-qubit interaction terms. Here, combinations of`Z`

and`X`

operators on pairs of qubits are mapped to their corresponding coefficients. These terms capture the inter-qubit interactions in the Hamiltonian, with coefficients also expressed as functions of`dt`

and`x`

.

```
singles_syms = {
zi: -21 / 8 * syms[1],
iz: -21 / 8 * syms[1],
xi: -3 / 2 * syms[0] * syms[1],
ix: -3 / 2 * syms[0] * syms[1],
}
```

```
doubles_sym = {
zz: -3 / 8 * syms[1],
xz: -1 / 2 * syms[0] * syms[1],
zx: -1 / 2 * syms[0] * syms[1],
}
```

We now combine the 1-qubit and 2-qubit interaction terms defined in `singles_syms`

and `doubles_sym`

to form the first-order Trotterized Hamiltonian:

```
hamiltonian_trotterized = QubitPauliOperator({**singles_syms, **doubles_sym})
```

```
hamiltonian_trotterized
```

```
{(Zq[0]): -2.625*dt, (Zq[1]): -2.625*dt, (Xq[0]): -1.5*dt*x, (Xq[1]): -1.5*dt*x, (Zq[0], Zq[1]): -0.375*dt, (Xq[0], Zq[1]): -0.5*dt*x, (Zq[0], Xq[1]): -0.5*dt*x}
```

We can generate the circuit by using the same function as before

```
circ = Circuit(2)
trotter_step_circ_symb = gen_term_sequence_circuit(hamiltonian_trotterized, circ)
```

```
render_circuit_jupyter(trotter_step_circ_symb)
```

Again, we can also decompose this symbolic circuit into basic gates and count these gates.

```
from pytket.passes import DecomposeBoxes
test_circuit = trotter_step_circ_symb.copy()
DecomposeBoxes().apply(test_circuit)
gate_counts(test_circuit)
```

```
Counter({<OpType.H: 33>: 8, <OpType.Rz: 36>: 7, <OpType.CX: 42>: 6})
```

```
render_circuit_jupyter(test_circuit)
```

We can also make this symbolic circuit into a box that we can repeat for multiple Trotter steps:

```
trotter_step_box_symb = CircBox(trotter_step_circ_symb)
```

This is a 2-step Trotter symbolic circuit:

```
# two steps
circ = Circuit(2)
circ.add_gate(trotter_step_box_symb, circ.qubits)
circ.add_gate(trotter_step_box_symb, circ.qubits)
```

```
[CircBox q[0], q[1]; CircBox q[0], q[1]; ]
```

In the final step of working with the symbolic circuit, we demonstrate how to substitute our symbolic parameters with actual numerical values:

**Defining Symbol Values**: We first define a list`x_dt_values`

containing the specific values for our symbols. In this example,`x`

is set to 2.0 and`dt`

to 0.01.**Creating the Symbol-Value Map**: Using`zip(syms, x_dt_values)`

, we pair each symbol with its corresponding value, and then convert this pairing into a dictionary`sym_map`

. This dictionary maps each symbolic variable to its numerical value.**Applying Substitutions**: We use`circ.symbol_substitution(sym_map)`

to apply these substitutions to our circuit. This method replaces each occurrence of the symbolic variables in`circ`

with their defined numerical values.

Lastly, we render the updated circuit using `render_circuit_jupyter(circ)`

. This visualization shows the circuit as it would be with the specific parameters applied.

```
x_dt_values = [2.0, 0.01]
sym_map = dict(zip(syms, x_dt_values))
```

```
sym_map
```

```
{x: 2.0, dt: 0.01}
```

```
circ.symbol_substitution(sym_map)
render_circuit_jupyter(circ)
```

## Implementation of the SU(2) 5-Plaquette Model Hamiltonian¶

Having explored the SU(2) 2-Plaquette Model, a foundational representation in lattice gauge theories, we now shift our focus to the more complex SU(2) 5-Plaquette Model. This transition marks a significant step towards understanding and simulating more intricate quantum systems. The 5-Plaquette Model, with its increased number of qubits, provides a richer and more detailed canvas to examine the dynamics of lattice gauge theories. This model’s Hamiltonian, characterized by its division into electric and magnetic components, captures the nuanced interactions within a lattice. Such a detailed representation is crucial for simulating real-world quantum phenomena and lays the groundwork for advanced quantum computing applications. By delving into the 5-Plaquette Model, we gain deeper insights into the behavior of complex quantum systems, paving the way for more accurate simulations and potentially groundbreaking discoveries in the realm of quantum physics. The following section presents the Hamiltonian of the SU(2) 5-Plaquette Model, implemented using `Qiskit`

, demonstrating how to interact with `Qiskit`

code.

For the SU(2) 5-Plaquette Model, the 5-qubit Hamilatonian is given as follows

This model’s Hamiltonian is explicitly broken down into electric \(h_E \) and magnetic \(h_B \) components, each comprising various terms involving Pauli \( \sigma_Z \) and \( \sigma_X \) operators.

In this section, we construct the Hamiltonian for our quantum system using Qiskit’s `SparsePauliOp`

build from a list of `Pauli`

strings and their associated coefficients. The Hamiltonian is divided into two main components: the electric component (`h_e`

) and the magnetic component (`h_b`

). Each component is expressed as list of Pauli terms acting on the qubits and their respective coefficients. We use `SparsePauliOp`

to write the Hamiltonian operator following the equations above quite literally (fix \(x=2\)):

The electric component

`h_e`

includes terms involving`Z`

Pauli operators, reflecting the electric interactions in the system.The magnetic component

`h_b`

comprises terms with`X`

and combined`XZ`

Pauli operators, representing the magnetic interactions.

The coefficients of each term are specific to our model, reflecting the strength of each interaction. We then combine these two components to form the complete Hamiltonian of the system, which is used for further analysis and simulation.

```
from qiskit.quantum_info import SparsePauliOp
h_e = [
("ZIIII", -9.0 / 8),
("IIIIZ", -9.0 / 8),
("IZIII", -3.0 / 4),
("IIZII", -3.0 / 4),
("IIIZI", -3.0 / 4),
("ZZIII", -3.0 / 8),
("IZZII", -3.0 / 8),
("IIZZI", -3.0 / 8),
("IIIZZ", -3.0 / 8),
]
h_b = [
("XIIII", -3.0),
("XZIII", -1.0),
("IIIIX", -3.0),
("IIIZX", -1.0),
("IXIII", -9.0 / 4),
("ZXIII", -3.0 / 4),
("IXZII", -3.0 / 4),
("ZXZII", -1.0 / 4),
("IIXII", -9.0 / 4),
("IZXII", -3.0 / 4),
("IIXZI", -3.0 / 4),
("IZXZI", -1.0 / 4),
("IIIXI", -9.0 / 4),
("IIZXI", -3.0 / 4),
("IIIXZ", -3.0 / 4),
("IIZXZ", -1.0 / 4),
]
hamiltonian = SparsePauliOp.from_list(h_e + h_b)
hamiltonian
```

```
SparsePauliOp(['ZIIII', 'IIIIZ', 'IZIII', 'IIZII', 'IIIZI', 'ZZIII', 'IZZII', 'IIZZI', 'IIIZZ', 'XIIII', 'XZIII', 'IIIIX', 'IIIZX', 'IXIII', 'ZXIII', 'IXZII', 'ZXZII', 'IIXII', 'IZXII', 'IIXZI', 'IZXZI', 'IIIXI', 'IIZXI', 'IIIXZ', 'IIZXZ'],
coeffs=[-1.125+0.j, -1.125+0.j, -0.75 +0.j, -0.75 +0.j, -0.75 +0.j, -0.375+0.j,
-0.375+0.j, -0.375+0.j, -0.375+0.j, -3. +0.j, -1. +0.j, -3. +0.j,
-1. +0.j, -2.25 +0.j, -0.75 +0.j, -0.75 +0.j, -0.25 +0.j, -2.25 +0.j,
-0.75 +0.j, -0.75 +0.j, -0.25 +0.j, -2.25 +0.j, -0.75 +0.j, -0.75 +0.j,
-0.25 +0.j])
```

A dense (or sparse) matrix representation of the operator can be accessed directly with the method `SparsePauliOp.to_matrix()`

:

```
hamiltonian.to_matrix(sparse=False)
```

```
array([[-6. +0.j, -4. +0.j, -4. +0.j, ..., -0. +0.j, -0. +0.j, -0. +0.j],
[-4. +0.j, -3. +0.j, -0. +0.j, ..., -0. +0.j, -0. +0.j, -0. +0.j],
[-4. +0.j, -0. +0.j, -3. +0.j, ..., -0. +0.j, -0. +0.j, -0. +0.j],
...,
[-0. +0.j, -0. +0.j, -0. +0.j, ..., 3. +0.j, -0. +0.j, -1. +0.j],
[-0. +0.j, -0. +0.j, -0. +0.j, ..., -0. +0.j, 1.5+0.j, -2. +0.j],
[-0. +0.j, -0. +0.j, -0. +0.j, ..., -1. +0.j, -2. +0.j, 3. +0.j]])
```

As a first step we can compute the ground state of the Hamiltonian by diagonalizing and extracting the lowest eigenvalue. There are many ways to do this but one can pass the (sparse) matrix obtained from `SparsePauliOp`

representing our Hamiltonian to a numerical eigensolver

```
from scipy.sparse.linalg import eigs
vals, vecs = eigs(hamiltonian.to_matrix(sparse=True), k=1)
print(f"Energy of the ground state {vals[0].real}")
```

```
Energy of the ground state -16.999769538608106
```

In this section, we build utility functions to translate `Qiskit`

’s Hamiltonian representation into the `QubitPauliOperator`

format used by `pytket`

.

`qps_from_sparsepauliop`

: This function converts a tensor of Pauli operators from`SparsePauliOp`

format to pytket’s`QubitPauliString`

. It iterates over qubit indices and corresponding Pauli operators, mapping them to the`pytket`

equivalent.`qpo_from_sparsepauliop`

: This function takes an`Qiskit`

`SparsePauliOp`

and converts it into a`QubitPauliOperator`

. It processes each term in the`Qiskit`

operator, using`qps_from_sparsepauliop`

for conversion, and associates it with the corresponding coefficient.

Finally, we apply `qpo_from_sparsepauliop`

to our previously defined Hamiltonian and print the result. This step confirms that the Hamiltonian is correctly converted into the format required by pytket for further quantum computations.

```
from pytket.pauli import Pauli, QubitPauliString
from pytket.utils import QubitPauliOperator
from pytket import Qubit
from collections import defaultdict
pauli_sym = {"I": Pauli.I, "X": Pauli.X, "Y": Pauli.Y, "Z": Pauli.Z}
def qps_from_sparsepauliop(paulis):
"""Convert SparsePauliOp tensor of Paulis to pytket QubitPauliString."""
qlist = []
plist = []
for q, p in enumerate(paulis):
if p != "I":
qlist.append(Qubit(q))
plist.append(pauli_sym[p])
return QubitPauliString(qlist, plist)
def qpo_from_sparsepauliop(sp_op):
"""Convert SparsePauliOp QubitOperator to pytket QubitPauliOperator."""
tk_op = defaultdict(complex)
for term, coeff in sp_op.to_list():
string = qps_from_sparsepauliop(term)
tk_op[string] += coeff
return QubitPauliOperator(tk_op)
hamiltonian_op = qpo_from_sparsepauliop(hamiltonian)
print(hamiltonian_op)
```

```
{(Zq[0]): -1.12500000000000, (Zq[4]): -1.12500000000000, (Zq[1]): -0.750000000000000, (Zq[2]): -0.750000000000000, (Zq[3]): -0.750000000000000, (Zq[0], Zq[1]): -0.375000000000000, (Zq[1], Zq[2]): -0.375000000000000, (Zq[2], Zq[3]): -0.375000000000000, (Zq[3], Zq[4]): -0.375000000000000, (Xq[0]): -3.00000000000000, (Xq[0], Zq[1]): -1.00000000000000, (Xq[4]): -3.00000000000000, (Zq[3], Xq[4]): -1.00000000000000, (Xq[1]): -2.25000000000000, (Zq[0], Xq[1]): -0.750000000000000, (Xq[1], Zq[2]): -0.750000000000000, (Zq[0], Xq[1], Zq[2]): -0.250000000000000, (Xq[2]): -2.25000000000000, (Zq[1], Xq[2]): -0.750000000000000, (Xq[2], Zq[3]): -0.750000000000000, (Zq[1], Xq[2], Zq[3]): -0.250000000000000, (Xq[3]): -2.25000000000000, (Zq[2], Xq[3]): -0.750000000000000, (Xq[3], Zq[4]): -0.750000000000000, (Zq[2], Xq[3], Zq[4]): -0.250000000000000}
```

The next line of code converts the `QubitPauliOperator`

representation of the Hamiltonian, obtained in the previous step, into a dense matrix format. The `to_sparse_matrix()`

method first transforms the Hamiltonian into a sparse matrix. We then use `todense()`

to convert this sparse matrix into a regular, dense matrix. This dense matrix representation is useful for certain types of numerical analyses and visual inspections of the Hamiltonian structure.

```
hamiltonian_op.to_sparse_matrix().todense()
```

```
matrix([[-6. +0.j, -4. +0.j, -4. +0.j, ..., 0. +0.j, 0. +0.j, 0. +0.j],
[-4. +0.j, -3. +0.j, 0. +0.j, ..., 0. +0.j, 0. +0.j, 0. +0.j],
[-4. +0.j, 0. +0.j, -3. +0.j, ..., 0. +0.j, 0. +0.j, 0. +0.j],
...,
[ 0. +0.j, 0. +0.j, 0. +0.j, ..., 3. +0.j, 0. +0.j, -1. +0.j],
[ 0. +0.j, 0. +0.j, 0. +0.j, ..., 0. +0.j, 1.5+0.j, -2. +0.j],
[ 0. +0.j, 0. +0.j, 0. +0.j, ..., -1. +0.j, -2. +0.j, 3. +0.j]])
```

```
np.allclose(
hamiltonian_op.to_sparse_matrix().todense(), hamiltonian.to_matrix(sparse=False)
)
```

```
True
```

We now proceed with the Trotter-Suzuki decomposition to approximate the time-evolution operator of the Hamiltonian in quantum simulations as we did above.

*Note:* `rescaling_factor_from_pytket = np.pi/2`

is defined for aligning simulation results with exact theoretical predictions. This factor compensates for an additional π/2 factor introduced by `gen_term_sequence_circuit`

in `pytket`

, which affects the time evolution steps. When comparing with exact results, divide the time step (\(\delta t\)) by this rescaling factor to account for this discrepancy.

```
# rescaling_factor_from_pytket = np.pi/2 ## use it to divide \deltat in case you want to compare with the exact results: gen_term_sequence adds Pauliexpbox with an extrat pi/2 factor
t_max = 1
n_trotter_steps = 10
time_step = t_max / n_trotter_steps
time_line = np.linspace(0, t_max, n_trotter_steps + 1)
print(f"Time step: {time_step}")
hamiltonian_slice = hamiltonian_op * time_step
```

```
Time step: 0.1
```

Next we construct the corresponding five-qubit circuit for two Trotter steps. A product of exponential of Pauli strings can be automatically generated by `pytket`

using `gen_term_sequence_circuits`

, including aggregating mutually commuting terms into groups.

```
initial_state = Circuit(5)
trotter_step_circ = gen_term_sequence_circuit(hamiltonian_slice, initial_state)
```

The Trotter step circuit can be enclosed in a `CircBox`

and added to any circuit (this way we can add many steps in a sequence)

```
trotter_step_box = CircBox(trotter_step_circ)
```

Now one Trotter step is a single box, and we can simply stack boxes to create circuits for more Trotter steps

```
# two steps
circ = Circuit(5)
circ.add_gate(trotter_step_box, circ.qubits)
circ.add_gate(trotter_step_box, circ.qubits)
```

```
[CircBox q[0], q[1], q[2], q[3], q[4]; CircBox q[0], q[1], q[2], q[3], q[4]; ]
```

```
from pytket.circuit.display import render_circuit_jupyter
render_circuit_jupyter(circ)
```

Now we sign into the Quantinuum Backend we will submit to, H1-2E, and compile the circuit for the backend.

```
from pytket.extensions.quantinuum import QuantinuumBackend
h_backend = QuantinuumBackend("H1-2E") # emulator
h_backend.login()
```

```
circ_compiled = h_backend.get_compiled_circuit(circ.measure_all())
render_circuit_jupyter(circ_compiled)
```

We can look at the statistics of the circuit before running it on the emulator.

```
gate_counts(circ_compiled)
```

```
Counter({<OpType.PhasedX: 68>: 59,
<OpType.ZZPhase: 73>: 34,
<OpType.Measure: 63>: 5})
```

Now we run the circuit on the Quantinuum H1-2 Emulator. Before submitting to the emulator, the total cost of running the circuit should be checked beforehand.

```
n_shots = 100
cost = h_backend.cost(circ_compiled, n_shots=n_shots, syntax_checker="H1-2SC")
print("Cost of experiment in HQCs:", cost)
```

```
Cost of experiment in HQCs: 13.98
```

Now we run the circuit with the noise model turned on

```
handle = h_backend.process_circuit(
circ_compiled, n_shots=n_shots, noisy_simulation=True
)
```

```
circuit_status = h_backend.circuit_status(handle)
print(circuit_status)
```

```
CircuitStatus(status=<StatusEnum.COMPLETED: 'Circuit has completed. Results are ready.'>, message='{"name": "job", "submit-date": "2024-01-11T06:24:50.716190", "result-date": "2024-01-11T06:25:15.622257", "queue-position": null, "cost": "13.98", "error": null}', error_detail=None, completed_time=None, queued_time=None, submitted_time=None, running_time=None, cancelled_time=None, error_time=None, queue_position=None)
```

```
result = h_backend.get_result(handle)
result.get_counts()
```

```
Counter({(1, 0, 1, 0, 1): 36,
(1, 0, 1, 1, 0): 6,
(1, 1, 1, 0, 1): 6,
(0, 0, 0, 0, 1): 5,
(0, 1, 1, 0, 1): 5,
(1, 0, 1, 1, 1): 5,
(0, 1, 1, 1, 1): 4,
(1, 1, 0, 1, 0): 4,
(1, 1, 0, 1, 1): 4,
(0, 0, 1, 0, 1): 3,
(0, 1, 0, 0, 1): 3,
(1, 0, 1, 0, 0): 3,
(1, 0, 0, 0, 0): 2,
(1, 0, 0, 1, 0): 2,
(1, 0, 0, 1, 1): 2,
(1, 1, 0, 0, 0): 2,
(0, 0, 0, 1, 0): 1,
(0, 0, 0, 1, 1): 1,
(0, 1, 0, 1, 1): 1,
(0, 1, 1, 0, 0): 1,
(1, 0, 0, 0, 1): 1,
(1, 1, 0, 0, 1): 1,
(1, 1, 1, 0, 0): 1,
(1, 1, 1, 1, 1): 1})
```

We can also plot the result using the following helper code:

```
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
def plot_counts(counts):
counts_record = [
{"State": str(state), "Count": count} for state, count in counts.items()
]
count_df = pd.DataFrame().from_records(counts_record)
sns.catplot(x="State", y="Count", kind="bar", data=count_df, aspect=12, height=2)
plt.show()
```

```
counts = result.get_counts()
```

```
plot_counts(counts)
```

And compare our results to the noiseless emulation, we see that our emulator results match the noiseless simulation closely.

```
noiseless_handle = h_backend.process_circuit(
circ_compiled, n_shots=n_shots, noisy_simulation=False
)
noiseless_result = h_backend.get_result(noiseless_handle)
noiseless_counts = noiseless_result.get_counts()
plot_counts(noiseless_counts)
```

One can also run the noiseless simulation using the Qiskit `Aer`

backend directly from `pytket`

using our `pytket-qiskit`

extension:

```
from pytket.extensions.qiskit import AerBackend
# define Aer backend
backend = AerBackend()
# compile circuit as usual, but now for the new backend
compiled_circ = backend.get_compiled_circuit(circ)
# submit and run the job: this runs locally and not in the cloud
handle = backend.process_circuit(compiled_circ, n_shots=100)
counts = backend.get_result(handle).get_counts()
plot_counts(counts)
# print(counts)
```

## Definitions¶

Definitions for the use case context in this article.

*Quantum Field Theory (QFT):*a mathematical and conceuptual framework in theoretical physics extending quantum mechanics, which deals with particles, over to fields, which are systems with infinite number of degrees of freedom. It is used in particle physics to construct physical models of subatomic particles and condensed matter physics to construct models of quasiparticles. Difficult due to the infinite number of degrees of freedom.*Lattice Gauge Theory (LGT):*a mathematical framework in QFT. Rather than having infinite degrees of freedom, space and time are discretized to a very large number of degrees of freedom, enabling its study on computers.*Perturbation Theory:*a method for finding an approximate solutions to a problem, by starting from the exact solution of a related, simpler problem*Quantum Chromodynamics:*theory that describes the action of the strong nuclear force, the fundamental interaction between subatomic particles of matter inside protons and neutrons.