{
"cells": [
{
"cell_type": "markdown",
"id": "title",
"metadata": {},
"source": [
"# Trotter dynamics with Hamming-weight phasing\n",
"\n",
"**Download Notebook** - {nb-download}`trotter_hamming_weight_phasing.ipynb`\n",
"\n",
"This notebook compares two implementations of the same first-order Trotter step for the open-chain Ising Hamiltonian\n",
"\n",
"$$H = J \\sum_{i=0}^{n-2} Z_i Z_{i+1}. $$\n",
"\n",
"The baseline uses `trotter_first_order`, which constructs one ordinary Pauli exponential per Zixy term. The second implementation groups the nearest-neighbor terms into even and odd brick-wall layers. Within each layer it computes all `ZZ` parities in parallel, replaces the identical `Rz` rotations with Hamming-weight phasing, and uncomputes the parities.\n",
"\n",
"Because all terms in this simple Hamiltonian commute, first-order Trotterization is exact. This lets us isolate the effect of changing the circuit implementation."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "imports",
"metadata": {
"execution": {
"iopub.execute_input": "2026-08-14T08:09:17.812711Z",
"iopub.status.busy": "2026-08-14T08:09:17.812627Z",
"iopub.status.idle": "2026-08-14T08:09:19.900413Z",
"shell.execute_reply": "2026-08-14T08:09:19.899955Z"
}
},
"outputs": [],
"source": [
"from typing import no_type_check\n",
"import warnings\n",
"\n",
"warnings.filterwarnings(\"ignore\", category=SyntaxWarning)\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import zixy.qubit.pauli as zqp\n",
"from guppylang import guppy\n",
"from guppylang.defs import GuppyFunctionDefinition\n",
"from guppylang.std.angles import angle\n",
"from guppylang.std.builtins import array, comptime\n",
"from guppylang.std.debug import state_output\n",
"from guppylang.std.quantum import cx, discard_array, h, qubit\n",
"from selene_sim import Quest\n",
"\n",
"from guppyalgos.primitives.arithmetic.hamming_weight import num_hamming_weight_bits\n",
"from guppyalgos.primitives.rotations.hamming_weight_phasing import hamming_weight_phase\n",
"from guppyalgos.algorithms.time_evolution.trotter import ham_sim_trotter, trotter_first_order\n",
"from guppyalgos.utils import qarray, transversal"
]
},
{
"cell_type": "markdown",
"id": "hamiltonian-intro",
"metadata": {},
"source": [
"## Build the nearest-neighbor Hamiltonian with Zixy\n",
"\n",
"We use eight state qubits for the dynamics comparison. The Guppy angle convention implemented by `pauli_exp` gives $\\exp[-i(\\pi/2)\\, J\\, \\Delta t\\, Z_iZ_{i+1}]$ for each term."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "hamiltonian",
"metadata": {
"execution": {
"iopub.execute_input": "2026-08-14T08:09:19.901643Z",
"iopub.status.busy": "2026-08-14T08:09:19.901506Z",
"iopub.status.idle": "2026-08-14T08:09:19.905687Z",
"shell.execute_reply": "2026-08-14T08:09:19.905291Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(0.7, Z0 Z1), (0.7, Z1 Z2), (0.7, Z2 Z3), (0.7, Z3 Z4), (0.7, Z4 Z5), (0.7, Z5 Z6), (0.7, Z6 Z7)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N_STATE_QUBITS = 8\n",
"COUPLING = 0.7\n",
"TIME_STEP = 0.12\n",
"MAX_STEPS = 12\n",
"\n",
"\n",
"def nearest_neighbor_ising_hamiltonian(\n",
" n_qubits: int, coupling: float\n",
") -> zqp.RealTermSum:\n",
" if n_qubits < 3:\n",
" raise ValueError(\"This brick-wall example needs at least three qubits\")\n",
" source = \", \".join(\n",
" f\"({coupling}, Z{left} Z{left + 1})\"\n",
" for left in range(n_qubits - 1)\n",
" )\n",
" return zqp.RealTermSum.from_str(source, n_qubits)\n",
"\n",
"\n",
"hamiltonian = nearest_neighbor_ising_hamiltonian(\n",
" N_STATE_QUBITS, COUPLING\n",
")\n",
"hamiltonian"
]
},
{
"cell_type": "markdown",
"id": "brick-wall-explanation",
"metadata": {},
"source": [
"## Convert a Pauli-exponential step to brick-wall Hamming-weight phasing\n",
"\n",
"For one edge, the ordinary `ZZ` Pauli exponential is\n",
"\n",
"$$\\operatorname{CX}_{i,i+1}\\; R_z(\\theta)_{i+1}\\; \\operatorname{CX}_{i,i+1}. $$\n",
"\n",
"Edges in one brick-wall layer are disjoint, so their CX gates and parity targets can be handled together. If a layer contains $m$ identical rotations, Hamming-weight phasing replaces those $m$ arbitrary rotations with `num_hamming_weight_bits(m)` rotations, at the cost of reversible adders and clean ancillas."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "brick-wall-parser",
"metadata": {
"execution": {
"iopub.execute_input": "2026-08-14T08:09:19.906992Z",
"iopub.status.busy": "2026-08-14T08:09:19.906920Z",
"iopub.status.idle": "2026-08-14T08:09:19.909830Z",
"shell.execute_reply": "2026-08-14T08:09:19.909440Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Even layer: ((0, 1), (2, 3), (4, 5), (6, 7))\n",
"Odd layer: ((1, 2), (3, 4), (5, 6))\n"
]
}
],
"source": [
"def brick_wall_layers(\n",
" hamiltonian: zqp.RealTermSum,\n",
") -> tuple[tuple[tuple[int, int], ...], tuple[tuple[int, int], ...], float]:\n",
" terms: list[zqp.RealTerm] = list(hamiltonian.to_terms())\n",
" if not terms:\n",
" raise ValueError(\"The Hamiltonian must contain at least one ZZ term\")\n",
"\n",
" coupling = float(terms[0].coeff)\n",
" even_edges: list[tuple[int, int]] = []\n",
" odd_edges: list[tuple[int, int]] = []\n",
"\n",
" for term in terms:\n",
" paulis = term.string.get_dict()\n",
" support = tuple(sorted(paulis))\n",
" if (\n",
" len(support) != 2\n",
" or support[1] != support[0] + 1\n",
" or any(pauli != zqp.Z for pauli in paulis.values())\n",
" ):\n",
" raise ValueError(\"Expected only nearest-neighbor ZZ terms\")\n",
" if not np.isclose(float(term.coeff), coupling):\n",
" raise ValueError(\"Hamming-weight phasing requires uniform coupling\")\n",
" (even_edges if support[0] % 2 == 0 else odd_edges).append(support)\n",
"\n",
" return tuple(even_edges), tuple(odd_edges), coupling\n",
"\n",
"\n",
"even_edges, odd_edges, _ = brick_wall_layers(hamiltonian)\n",
"print(\"Even layer:\", even_edges)\n",
"print(\"Odd layer: \", odd_edges)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "hwp-step",
"metadata": {
"execution": {
"iopub.execute_input": "2026-08-14T08:09:19.910683Z",
"iopub.status.busy": "2026-08-14T08:09:19.910616Z",
"iopub.status.idle": "2026-08-14T08:09:19.918310Z",
"shell.execute_reply": "2026-08-14T08:09:19.917951Z"
}
},
"outputs": [],
"source": [
"def hamming_weight_zz_layer(\n",
" edges: tuple[tuple[int, int], ...],\n",
" n_state_qubits: int,\n",
") -> GuppyFunctionDefinition:\n",
" n_edges = len(edges)\n",
" phase_targets = hamming_weight_phase(n_edges)\n",
"\n",
" @guppy\n",
" @no_type_check\n",
" def layer(\n",
" state_qreg: array[qubit, n_state_qubits], theta: angle\n",
" ) -> None:\n",
" lefts = comptime(array(left for left, _ in edges))\n",
" rights = comptime(array(right for _, right in edges))\n",
" controls = array(\n",
" state_qreg.take(lefts[i]) for i in range(n_edges)\n",
" )\n",
" targets = array(\n",
" state_qreg.take(rights[i]) for i in range(n_edges)\n",
" )\n",
"\n",
" transversal(cx, controls, targets)\n",
" phase_targets(targets, theta)\n",
" transversal(cx, controls, targets)\n",
"\n",
" for i in range(n_edges):\n",
" state_qreg.put(controls.take(i), lefts[i])\n",
" state_qreg.put(targets.take(i), rights[i])\n",
" controls.discard_all_taken()\n",
" targets.discard_all_taken()\n",
"\n",
" return layer\n",
"\n",
"\n",
"def hamming_weight_ising_trotter_step(\n",
" hamiltonian: zqp.RealTermSum,\n",
" n_state_qubits: int,\n",
") -> GuppyFunctionDefinition:\n",
" even_edges, odd_edges, coupling = brick_wall_layers(hamiltonian)\n",
" if not even_edges or not odd_edges:\n",
" raise ValueError(\"Both brick-wall layers must be non-empty\")\n",
"\n",
" even_layer = hamming_weight_zz_layer(even_edges, n_state_qubits)\n",
" odd_layer = hamming_weight_zz_layer(odd_edges, n_state_qubits)\n",
"\n",
" @guppy\n",
" @no_type_check\n",
" def trotter_step(\n",
" state_qreg: array[qubit, n_state_qubits], time_step: float\n",
" ) -> None:\n",
" theta = angle(comptime(coupling) * time_step)\n",
" even_layer(state_qreg, theta)\n",
" odd_layer(state_qreg, theta)\n",
"\n",
" return trotter_step\n",
"\n",
"\n",
"normal_step = trotter_first_order(hamiltonian, N_STATE_QUBITS)\n",
"hwp_step = hamming_weight_ising_trotter_step(\n",
" hamiltonian, N_STATE_QUBITS\n",
")"
]
},
{
"cell_type": "markdown",
"id": "dynamics-explanation",
"metadata": {},
"source": [
"## Build the full dynamics with `ham_sim_trotter`\n",
"\n",
"We prepare $|+\\rangle^{\\otimes n}$ and track the mean $X$ magnetization\n",
"\n",
"$$\\langle \\bar X \\rangle = \\frac{1}{n}\\sum_i \\langle X_i \\rangle. $$\n",
"\n",
"For every time point, `ham_sim_trotter` composes the requested number of steps into a full Hamiltonian-simulation function. We use exactly the same helper for the ordinary and Hamming-weight-phased steps. The latter emulator is given enough capacity for the largest brick-wall layer's clean ancillas."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "compile-dynamics",
"metadata": {
"execution": {
"iopub.execute_input": "2026-08-14T08:09:19.919700Z",
"iopub.status.busy": "2026-08-14T08:09:19.919631Z",
"iopub.status.idle": "2026-08-14T08:09:19.921742Z",
"shell.execute_reply": "2026-08-14T08:09:19.921407Z"
}
},
"outputs": [],
"source": [
"def output_state(\n",
" hamiltonian_simulation: GuppyFunctionDefinition,\n",
" simulator_qubits: int,\n",
") -> np.ndarray:\n",
" @guppy\n",
" @no_type_check\n",
" def main() -> None:\n",
" state_qreg = qarray(N_STATE_QUBITS)\n",
" transversal(h, state_qreg)\n",
" hamiltonian_simulation(state_qreg)\n",
" state_output(\"state\", state_qreg)\n",
" discard_array(state_qreg)\n",
"\n",
" result = main.emulator(simulator_qubits).run()\n",
" states = Quest.extract_states_dict(result.results[0].entries)\n",
" return states[\"state\"].get_single_state()\n",
"\n",
"\n",
"NORMAL_SIMULATOR_QUBITS = N_STATE_QUBITS\n",
"HWP_SIMULATOR_QUBITS = (\n",
" N_STATE_QUBITS + max(len(even_edges), len(odd_edges))\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "run-dynamics",
"metadata": {
"execution": {
"iopub.execute_input": "2026-08-14T08:09:19.922828Z",
"iopub.status.busy": "2026-08-14T08:09:19.922685Z",
"iopub.status.idle": "2026-08-14T08:12:47.622348Z",
"shell.execute_reply": "2026-08-14T08:12:47.622040Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Maximum statevector error: 6.835e-15\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" step | \n",
" time | \n",
" Normal Pauli exponentials | \n",
" Hamming-weight phasing | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 0 | \n",
" 0.00 | \n",
" 1.000000 | \n",
" 1.000000 | \n",
"
\n",
" \n",
" | 1 | \n",
" 1 | \n",
" 0.12 | \n",
" 0.940317 | \n",
" 0.940317 | \n",
"
\n",
" \n",
" | 2 | \n",
" 2 | \n",
" 0.24 | \n",
" 0.775754 | \n",
" 0.775754 | \n",
"
\n",
" \n",
" | 3 | \n",
" 3 | \n",
" 0.36 | \n",
" 0.545950 | \n",
" 0.545950 | \n",
"
\n",
" \n",
" | 4 | \n",
" 4 | \n",
" 0.48 | \n",
" 0.305267 | \n",
" 0.305267 | \n",
"
\n",
" \n",
" | 5 | \n",
" 5 | \n",
" 0.60 | \n",
" 0.108557 | \n",
" 0.108557 | \n",
"
\n",
" \n",
" | 6 | \n",
" 6 | \n",
" 0.72 | \n",
" -0.003023 | \n",
" -0.003023 | \n",
"
\n",
" \n",
" | 7 | \n",
" 7 | \n",
" 0.84 | \n",
" -0.012361 | \n",
" -0.012361 | \n",
"
\n",
" \n",
" | 8 | \n",
" 8 | \n",
" 0.96 | \n",
" 0.069876 | \n",
" 0.069876 | \n",
"
\n",
" \n",
" | 9 | \n",
" 9 | \n",
" 1.08 | \n",
" 0.209057 | \n",
" 0.209057 | \n",
"
\n",
" \n",
" | 10 | \n",
" 10 | \n",
" 1.20 | \n",
" 0.356858 | \n",
" 0.356858 | \n",
"
\n",
" \n",
" | 11 | \n",
" 11 | \n",
" 1.32 | \n",
" 0.465143 | \n",
" 0.465143 | \n",
"
\n",
" \n",
" | 12 | \n",
" 12 | \n",
" 1.44 | \n",
" 0.499605 | \n",
" 0.499605 | \n",
"
\n",
" \n",
"
\n",
"