Zum Hauptinhalt springen

Neutronenstreuung mit einem AQC + Trotter-Dynamics-Serverless-Workflow simulieren

Nutzungsschätzung: 18 Minuten auf einem Heron-r3-Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine Laufzeit kann variieren.)

Lernziele

  • Wie ein inelastisches Neutronenstreuungsspektrum auf den dynamischen Strukturfaktor S(q,ω)S(q, \omega) eines eindimensionalen Quantenmagneten abgebildet wird.

  • Wie man den Grundzustand von KCuF3_3 (isotropes Heisenberg-Modell) mit der Dichtematrix-Renormierungsgruppe (DMRG) und der Fidelity-Maximierung mit Matrixproduktzuständen (MPS) vorbereitet.

  • Wie man Trotter-Zeitentwicklung, Circuit-Kompression durch Approximate Quantum Compilation (AQC) und fehlergemilderte Ausführung als einzelnen Funktionsaufruf durchführt.

  • Wie man die ortsaufgelöste σz(t)\langle \sigma_z \rangle(t)-Zeitreihe zu S(q,ω)S(q, \omega) nachbearbeitet und das Zwei-Spinon-Kontinuum identifiziert.

Voraussetzungen

Hintergrund

Inelastische Neutronenstreuung misst den dynamischen Strukturfaktor S(q,ω)S(q, \omega), die Raum-Zeit-Fourier-Transformierte der Spin-Spin-Korrelationsfunktion, sodass die Reproduktion von S(q,ω)S(q, \omega) aus einem mikroskopischen Spinmodell einen direkten, falsifizierbaren Test einer Quantensimulation darstellt. Dieses Tutorial untersucht KCuF3_3, eine antiferromagnetische Heisenberg-Kette mit Spin-12\frac{1}{2}, deren Anregungen keine einzelnen Spin-Flips, sondern Paare fraktionalisierter Spinonen sind: Statt einer scharfen Magnon-Dispersion zeigt S(q,ω)S(q, \omega) ein breites Zwei-Spinon-Kontinuum, das unten durch π2sinq\tfrac{\pi}{2}|\sin q| und oben durch πsin(q/2)\pi|\sin(q/2)| begrenzt ist. Das sind die gestrichelten Kurven in den folgenden Plots. Die vollständige Physik sowie der Vergleich mit gemessenen Neutronendaten werden im ursprünglichen Tutorial und in Lee et al., arXiv:2603.15608, behandelt.

Der Quanten-Workflow spiegelt das Streuungsexperiment wider:

  1. Bereite den Grundzustand ψ0|\psi_0\rangle der Kette vor.

  2. Versetze ihm einen Kick durch eine lokale Störung an der zentralen Stelle, eine π/2\pi/2-ZZ-Rotation, die den Impuls- und Energieübertrag durch das Neutron nachahmt.

  3. Entwickle ihn zeitlich unter dem Heisenberg-Hamiltonian, eiHte^{-iHt}, mit einer Trotter-Produktformel.

  4. Miss die ortsaufgelöste Magnetisierung σzj(t)\langle \sigma_z^j \rangle(t). Als Funktion der Stelle jj und der Zeit tt ist dies genau die retardierte Greensche Funktion GR(j,jc,t)G^R(j, j_c, t), sodass vor der Fourier-Transformation in Schritt 5 keine Umrechnung nötig ist.

  5. Fourier-transformiere GRG^R zu S(q,ω)S(q, \omega).

In Schritt 3 können Probleme auftreten, wenn exakte Trotter-Circuits für lange Zeitentwicklungen zu tief für die Hardware werden. AQC mit Tensornetzwerken löst dies, indem ein Block von Trotter-Schritten in einen festen, flachen parametrisierten Ansatz komprimiert wird, dessen Zustands-Fidelity zur exakten Zeitentwicklung klassisch mit einem MPS-Simulator maximiert wird (arXiv:2301.08609). Die AQC-Dynamics-Vorlage verpackt diesen gesamten Quantenkern (Trotter-Synthese, AQC-Kompression und fehlergemilderte Ausführung) hinter einem einzigen Aufruf:

PRE (dieses Notebook)FUNCTION (aqc-dynamics-function)POST (dieses Notebook)
Grundzustand aus DMRG plus MPS-Fidelity-Maximierung, mit dem Neutronen-Kick, der in denselben Circuit eingebaut istTrotter-Synthese → AQC-Kompression → Ausführung auf statevector, fake oder runtime, gibt σzj(t)\langle \sigma_z^j \rangle(t) pro Stelle zurückS(q,ω)S(q, \omega), der dynamische Strukturfaktor

Die experimentspezifische Arbeit bleibt hier im Notebook: die Grundzustandsvorbereitung (PRE) und die S(q,ω)S(q, \omega)-Nachbearbeitung (POST). Die beiden quantenintensiven Schritte, Kompression und Ausführung, laufen innerhalb der Funktion.

Dieses Tutorial ist ein Begleit-Tutorial zu Neutronenstreuung in Quantenmaterialien mit Quanten-Circuits simulieren, das dasselbe Experiment inline aufbaut: dasselbe KCuF3_3-Modell, Grundzustandsvorbereitung, Neutronen-Kick und Nachbearbeitung, wobei die Trotter-Synthese, AQC-Kompression und fehlergemilderte Ausführung Schritt für Schritt ausgeschrieben sind. Lies jenes Tutorial, um zu erfahren, wie die AQC-Kompression funktioniert. Lies dieses hier, um dasselbe Experiment über eine bereitgestellte Funktionsvorlage auszuführen: Der Quantenkern wird zu einem einzigen Funktionsaufruf, und die mehrstündige AQC-Kompression läuft innerhalb des Serverless-Workers statt auf deinem eigenen Rechner, sodass du während der Ausführung weder ein HPC-System noch einen offenen Kernel benötigst. Derselbe Aufruf steuert auch andere 1D-Dynamics-Experimente.

Voraussetzungen

Stelle vor Beginn dieses Tutorials sicher, dass du Folgendes hast:

  • Die Funktion, die in deinem Qiskit-Serverless-Konto bereitgestellt ist. Führe zuerst die zugehörige Funktionsvorlage aus: AQC + Trotter-Dynamics-Funktionsvorlage bereitstellen und ausführen. Dieser Leitfaden führt dich durch das Beschaffen der Quelldateien und das Hochladen der Funktion in dein Konto. Dieses Tutorial ruft nur die bereitgestellte Funktion auf.

  • Für QiskitServerless gespeicherte IBM Quantum®-Anmeldedaten (siehe die Funktionsvorlage). Beide Beispiele in diesem Tutorial rufen die bereitgestellte Funktion auf, daher benötigen beide sie.

  • Qiskit SDK v2.0 oder höher (pip install qiskit).

  • Der Qiskit IBM Catalog-Client (pip install qiskit-ibm-catalog).

  • NumPy, SciPy und Matplotlib (pip install numpy scipy matplotlib). SciPy 1.14 oder höher wird für den COBYQA-Optimierer benötigt, der bei der Grundzustandsvorbereitung verwendet wird.

  • Den AQC-Tensornetzwerk-Stack, da die Grundzustandsvorbereitung in Schritt 1 lokal in diesem Notebook läuft: pip install 'qiskit-addon-aqc-tensor[quimb-jax]==0.3.1'.

Der erste Aufruf einer neu bereitgestellten Funktion wartet, während der Serverless-Worker seine Abhängigkeiten installiert, daher ist bei diesem Durchlauf mit zusätzlicher Latenz zu rechnen.

Setup

Importiere die Bibliotheken und definiere die später verwendeten experimentspezifischen Hilfsfunktionen: build_gs_ansatz (der Hamiltonian Variational Ansatz, oder HVA, für die Grundzustandsvorbereitung), prepare_ground_state (DMRG plus MPS-Fidelity-Maximierung) sowie get_spectrum, plot_green und plot_spectrum (die S(q,ω)S(q, \omega)-Nachbearbeitung). Diese sind aus dem ursprünglichen Neutronenstreuungs-Tutorial übernommen.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-addon-aqc-tensor qiskit-ibm-catalog quimb scipy
from functools import partial

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize

import quimb.tensor as qtn
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit_addon_aqc_tensor.simulation import tensornetwork_from_circuit
from qiskit_addon_aqc_tensor.simulation.quimb import QuimbSimulator
from qiskit_ibm_catalog import QiskitServerless
# Dynamical structure factor via discrete Fourier transform

def get_spectrum(n, Gjjc, dt, time_steps, q_steps, w_steps):
"""Compute the dynamical structure factor from the retarded Green's function.

Uses the center-site approximation and a discrete Fourier transform.
"""
green = Gjjc / 4 # sigma -> S=1/2
omega_max = np.pi / dt
qpoints = np.arange(0, 2 * np.pi, 2 * np.pi / q_steps)
omegas = np.arange(0, omega_max, omega_max / w_steps)
green_map = np.zeros((omegas.shape[0], qpoints.shape[0]))
center = n // 2 - 1
for iw, w in enumerate(omegas):
exponent = np.exp(1j * w * dt * np.arange(1, time_steps + 1))
S_w = np.dot(green.T, exponent) * dt
for iq, q in enumerate(qpoints):
q_matrix = np.exp(-1j * q * np.arange(-center, center + 2, 1))
green_map[iw, iq] = np.imag(np.dot(S_w, q_matrix))
return green_map

# Plotting helpers

def plot_spectrum(
dsf,
dt,
q_steps,
w_steps,
lower_bound=False,
upper_bound=False,
title=None,
):
"""Heat-map of the dynamical structure factor."""
omega_max = np.pi / dt
qpoints = np.arange(0, 2 * np.pi, 2 * np.pi / q_steps)
omegas = np.arange(0, omega_max, omega_max / w_steps)
x, y = np.meshgrid(qpoints, omegas)
fig, ax = plt.subplots(figsize=(8, 5))
c = ax.pcolormesh(x, y, dsf / np.max(dsf), cmap="viridis", shading="auto")
fig.colorbar(c, ax=ax, label="Normalized intensity")
if lower_bound:
ax.plot(
qpoints,
np.pi * np.abs(np.sin(qpoints)) / 2,
"--",
color="white",
lw=1.5,
label="Lower bound",
)
if upper_bound:
ax.plot(
qpoints,
np.pi * np.abs(np.sin(qpoints / 2)),
"--",
color="red",
lw=1.5,
label="Upper bound",
)
ax.set_ylim(0, 3.6)
ax.set_xlim(0, 2 * np.pi - 2 * np.pi / q_steps)
ax.set_xlabel(r"$q$", fontsize=16)
ax.set_ylabel(r"$\tilde{\omega} = \omega / J$", fontsize=16)
ax.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])
ax.set_xticklabels(["0", r"$\pi/2$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"])
if lower_bound or upper_bound:
ax.legend(loc="upper right", fontsize=11)
if title:
ax.set_title(title, fontsize=14)
plt.tight_layout()
plt.show()

def plot_green(n, Gjjc, time_steps, dt, title=None):
"""Heat-map of the retarded Green's function in real space and time."""
fig, ax = plt.subplots(figsize=(8, 6))
t_axis = np.arange(1, time_steps + 1) * dt
site_axis = np.arange(n)
x, y = np.meshgrid(t_axis, site_axis)
c = ax.pcolormesh(
x,
y,
np.real(Gjjc).T,
cmap="RdBu",
vmax=0.5,
vmin=-0.5,
shading="auto",
)
fig.colorbar(c, ax=ax, label=r"Re $G^R(j, j_c, t)$")
ax.set_xlabel(r"Time ($t / J^{-1}$)", fontsize=16)
ax.set_ylabel("Site index $j$", fontsize=16)
if title:
ax.set_title(title, fontsize=14)
plt.tight_layout()
plt.show()

# Variational ground-state ansatz (HVA)

def _apply_xxz_pair_gate(qc, q0, q1, theta):
"""Apply the parameterized XXZ-type two-qubit gate used in the HVA."""
qc.cx(q0, q1)
qc.rz(theta, q1)
qc.h(q0)
qc.rz(theta + np.pi / 2, q0)
qc.cx(q0, q1)
qc.rz(-theta, q1)
qc.h(q1)
qc.cx(q1, q0)
qc.rz(np.pi / 2, q1)
qc.rz(-np.pi / 2, q0)
qc.h(q1)
qc.h(q0)

def build_gs_ansatz(n, params, layers):
"""Build the Hamiltonian variational ansatz (HVA) circuit for
ground-state preparation of the 1D Heisenberg model.

Starts from a product of singlet pairs and applies alternating
odd/even layers of parameterized XXZ gates. For layer r,
params[2 * r] is the odd-layer (inter-pair) angle and
params[2 * r + 1] is the even-layer (intra-pair) angle.
"""
qc = QuantumCircuit(n)
# Initial singlet product state
for i in range(n // 2):
qc.x(2 * i)
qc.x(2 * i + 1)
qc.h(2 * i + 1)
qc.cx(2 * i + 1, 2 * i)
# Variational layers
for r in range(layers):
for i in range(1, (n + 1) // 2): # odd layer
_apply_xxz_pair_gate(qc, 2 * i - 1, 2 * i, params[2 * r])
for i in range(n // 2): # even layer
_apply_xxz_pair_gate(qc, 2 * i, 2 * i + 1, params[2 * r + 1])
return qc

def prepare_ground_state(n, gs_layers=5, max_bond=128, cutoff=1e-8):
"""Prepare the KCuF3 (isotropic Heisenberg) ground state as a QuantumCircuit.

Runs DMRG (quimb MPO + DMRG2) to get the chain's ground state, then optimizes
the HVA angles to maximize the MPS overlap |<psi_ansatz|psi_DMRG>|^2. No exact
diagonalization, so it scales to larger n.
"""
J = Jz = 1.0
builder = qtn.SpinHam1D(S=1 / 2)
builder += J * 0.5, "+", "-"
builder += J * 0.5, "-", "+"
builder += Jz, "Z", "Z"
H_mpo = builder.build_mpo(L=n)
dmrg = qtn.DMRG2(H_mpo)
dmrg.solve(tol=1e-8, verbosity=0)

gs_sim = QuimbSimulator(
quimb_circuit_factory=partial(
qtn.CircuitMPS, gate_opts=dict(cutoff=cutoff, max_bond=max_bond)
),
autodiff_backend="jax",
)

def gs_infidelity(params):
psi = tensornetwork_from_circuit(
build_gs_ansatz(n, params, gs_layers), gs_sim
).psi
return 1 - abs(psi.H @ dmrg.state) ** 2

# Seed and optimizer match the original tutorial. Each layer starts at
# [0, pi/2]: an odd-layer angle of 0 makes the inter-pair gate the identity,
# and an even-layer angle of pi/2 makes the intra-pair gate a SWAP (since
# 0.5 * (XX + YY + ZZ) = SWAP - I/2). That puts the seed at the singlet-pair
# product limit, which is already a decent approximation to the Heisenberg
# ground state, so the optimizer only has to refine it. The small jitter
# (fixed RNG seed, so runs are reproducible) breaks the exact symmetry
# between layers; COBYQA then runs for up to 100 iterations.
rng = np.random.default_rng(12345)
x0 = np.tile([0.0, np.pi / 2], gs_layers) + rng.normal(
scale=0.1, size=2 * gs_layers
)
result_gs = scipy.optimize.minimize(
gs_infidelity, x0, method="COBYQA", options={"maxiter": 100}
)
print(f"DMRG ground-state energy: {dmrg.energy:.6f}")
print(f"GS fidelity: {1 - result_gs.fun:.4f}")
return build_gs_ansatz(n, result_gs.x, gs_layers)

print("Setup complete - helpers defined.")
Setup complete - helpers defined.

Die Funktionsvorlage laden

Verbinde dich mit Qiskit Serverless und lade die bereitgestellte aqc-dynamics-function. Beide Beispiele in diesem Tutorial rufen dasselbe fn-Handle auf, daher wird die Funktion hier nur einmal geladen.

# Credentials are read from the account saved once via QiskitServerless.save_account(...)
serverless = QiskitServerless()
fn = serverless.load("aqc-dynamics-function")

Kleinmaßstäbliches Simulatorbeispiel

Wir führen zunächst den vollständigen Workflow auf einer kleinen Kette mit 10 Stellen unter Verwendung des exakten statevector-Backends aus. Dies validiert die PRE → FUNCTION → POST-Pipeline, bevor QPU-Zeit aufgewendet wird.

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Baue den KCuF3_3-Hamiltonian als SparsePauliOp auf (isotropes Heisenberg-Modell: XX+YY+ZZXX + YY + ZZ mit Kopplung 14\tfrac14 auf jeder Bindung zwischen nächsten Nachbarn; die Strings sind Pauli-Operatoren, sodass 14\tfrac14 die Spin-12\frac{1}{2}-Kopplung ergibt). Bereite den Grundzustand mit DMRG plus MPS-Fidelity-Maximierung vor und baue dann den Neutronen-Kick ein: eine π/2\pi/2-ZZ-Rotation an der zentralen Stelle. Der vorbereitete Circuit ist das, was wir der Funktion als initial_state übergeben. Wir belassen observables bei seinem Standardwert (ortsaufgelöstes ZZ), was genau das σzj(t)\langle \sigma_z^j \rangle(t)-Readout ist, das der Neutronen-Workflow benötigt.

n = 10
dt = 0.6 # physical time per Trotter step (also the omega-axis unit in POST)
time_steps = 10
center = n // 2 - 1

# MPS-simulator settings, shared by the ground-state prep here and the AQC
# compression inside the function (matches the original tutorial).
mps_max_bond = 32
mps_cutoff = 1e-8

# 1D isotropic Heisenberg (KCuF3) Hamiltonian on n qubits
H = SparsePauliOp.from_sparse_list(
[(p, [i, i + 1], 0.25) for i in range(n - 1) for p in ("XX", "YY", "ZZ")],
num_qubits=n,
)

# Ground state (DMRG + fidelity max) + neutron kick baked into the same circuit
gs_circuit = prepare_ground_state(
n, gs_layers=3, max_bond=mps_max_bond, cutoff=mps_cutoff
)
gs_circuit.rz(
np.pi / 2, center
) # exp(-i (pi/2)/2 Z_center): the neutron perturbation
print(
f"Prepared {n}-qubit ground state with the neutron kick at site {center}."
)
DMRG ground-state energy: -4.258035
GS fidelity: 0.9841
Prepared 10-qubit ground state with the neutron kick at site 4.

Schritte 2 und 3: Komprimieren und Ausführen mit der Funktionsvorlage

In einem handgeschriebenen Workflow sind dies zwei getrennte Phasen: die Circuits für die Hardware optimieren (Schritt 2) und sie ausführen (Schritt 3). Die Funktionsvorlage fasst beides in einem einzigen Aufruf zusammen. Sie führt Trotter-Synthese, AQC-Kompression und Hardware-Transpilation durch und führt dann die Circuits aus (hier auf dem exakten Simulator, später mit eingebauter Fehlerminderung auf Hardware). Die beiden Feinabstimmungsparameter sind aqc_segments (der Kompressionsplan) und aqc_options (die MPS- und Optimierereinstellungen). Jedes Segment {"n_steps": k, "ansatz_steps": m} komprimiert k aufeinanderfolgende Trotter-Schritte in einen Ansatz, der aus einem m-Schritt-Trotter-Ziel aufgebaut ist, und alle Schritte über sum(n_steps) hinaus laufen als reines Trotter. Frühe Schritte mit geringer Verschränkung lassen sich gut in einen flachen Ansatz (ansatz_steps=1) komprimieren, daher komprimieren wir hier die ersten drei Schritte in einen einlagigen Ansatz und die nächsten zwei in einen tieferen zweilagigen Ansatz; die restlichen fünf der 10 Trotter-Schritte laufen als reines Trotter. Für aqc_options orientieren wir uns am ursprünglichen Tutorial: MPS-Bindungsdimension max_bond=32, cutoff=1e-8, und ein L-BFGS-B-Optimierer, begrenzt auf 100 Iterationen.

Rufe die im Setup geladene Funktion auf. backend="statevector" führt den exakten Referenzpfad aus: keine QPU-Zeit, wobei die Circuits auf einem exakten Statevector-Simulator innerhalb des Serverless-Workers laufen (zum Aufrufen wird weiterhin ein gespeichertes Qiskit-Serverless-Konto benötigt). initial_state trägt den vorbereiteten Grundzustand (einschließlich des Kicks); observables wird weggelassen, sodass die Funktion das standardmäßige ortsaufgelöste ZZ misst.

job = fn.run(
t_steps=time_steps,
aqc_segments=[
{
"n_steps": 3,
"ansatz_steps": 1,
}, # early steps -> shallow 1-layer ansatz
{
"n_steps": 2,
"ansatz_steps": 2,
}, # later steps -> deeper 2-layer ansatz
],
aqc_options={
"max_bond": mps_max_bond, # MPS bond dimension for AQC compression
"cutoff": mps_cutoff,
"optimizer_settings": {
"method": "L-BFGS-B",
"jac": True,
"options": {"maxiter": 100},
},
},
dt=dt,
hamiltonian=H,
initial_state=gs_circuit, # prepared ground state including the neutron kick
# observables omitted -> default per-site Z (the neutron sigma_z readout)
backend="statevector",
)
print(job.status()) # rerun this cell until status says DONE
DONE
# The per-site <sigma_z>(t) the function returns is the retarded Green's function
# G(j, j_c, t). The workflow samples t = 1..time_steps, so drop the t = 0 row (the
# prepared+kicked state before any evolution) before post-processing.
result = job.result()
print(
"AQC fidelities:",
{k: round(v, 4) for k, v in result["metadata"]["aqc_fidelities"].items()},
)

ev = np.array(result["expectation_values"])
Gjjc = ev[1:] # shape (time_steps, n)
print("Green's function shape:", Gjjc.shape)
AQC fidelities: {'1': 1.0, '2': 0.9999, '3': 0.9992, '4': 0.9998, '5': 0.9995}
Green's function shape: (10, 10)

Schritt 4: Nachbearbeiten und Ergebnis im gewünschten klassischen Format zurückgeben

Fourier-transformiere die Greensche Funktion zu S(q,ω)S(q, \omega), symmetriere sie durch Spiegelung und schneide negative Werte ab: die übliche Neutronen-Nachbearbeitung. Die Spiegelung ist exakt, da S(q,ω)=S(q,ω)S(q, \omega) = S(-q, \omega) für dieses Modell gilt, und die verbleibenden negativen Werte sind Artefakte der Fourier-Transformation einer endlichen, diskret abgetasteten Zeitreihe, weshalb sie auf null gesetzt werden. Bei diesem kleinen exakten Durchlauf wird das Zwei-Spinon-Kontinuum nur grob aufgelöst, aber die Mechanik ist identisch mit dem folgenden Hardware-Durchlauf.

q_res, w_res = 100, 100
spectrum = get_spectrum(n, Gjjc, dt, time_steps, q_res, w_res)
spectrum = -(spectrum + spectrum[:, ::-1]) / 2 # mirror symmetry
spectrum = np.clip(spectrum, a_min=0, a_max=None) # clip negatives

plot_green(
n,
Gjjc,
time_steps,
dt,
title=f"Retarded Green's function - {n} qubits (AQC, statevector)",
)
plot_spectrum(
spectrum,
dt,
q_res,
w_res,
lower_bound=True,
upper_bound=True,
title=f"Dynamical structure factor - {n} qubits (AQC, statevector)",
)

Output of the previous code cell

Output of the previous code cell

Großmaßstäbliches Hardwarebeispiel

Derselbe Workflow skaliert ohne Änderung des wissenschaftlichen Codes hoch: eine Kette mit 30 Stellen, doppelte Trotter-Tiefe (20 Schritte), ein Kompressionsplan, der die Ansatz-Tiefe variiert (ein tieferer Ansatz für die späteren, stärker verschränkten Schritte), sowie die Ausführung auf einem IBM Quantum-Prozessor mit der eingebauten Fehlerminderung der Funktion (Dynamical Decoupling, Pauli-Twirling und Twirled Readout Error eXtinction (TREX)). Wir durchlaufen dieselben vier Schritte wie im Simulatorbeispiel und verwenden das fn-Handle aus dem Setup weiter.

Kleiner MaßstabGroßer Maßstab
Qubits1030
Trotter-Schritte1020
AQC-komprimierte Schritte (1-lagig + 2-lagig)3 + 2 = 56 + 4 = 10
Grundzustands-Ansatz-Schichten35
Maximale MPS-Bindungsdimension32128
BackendstatevectorQPU mit DD, Pauli-Twirling und TREX

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Baue denselben KCuF3_3-Heisenberg-SparsePauliOp auf und bereite den Grundzustand vor, nun mit einem tieferen Ansatz gs_layers=5 für die längere Kette, und baue dann den π/2\pi/2-ZZ-Neutronen-Kick an der zentralen Stelle ein. Dies ist identisch mit der kleinmaßstäblichen Abbildung, jedoch bei n=30n = 30.

Erwarte eine niedrigere Grundzustands-Fidelity als beim Durchlauf mit 10 Stellen: hier etwa 0,82 gegenüber 0,98 bei der kleineren Kette, da fünf HVA-Schichten einen Grundzustand mit 30 Stellen nicht vollständig erfassen können. Das ist eher zu erwarten als ein Fehlschlag, und das ursprüngliche Tutorial akzeptiert aus demselben Grund etwa 0,65 bei 50 Stellen. Eine Erhöhung von gs_layers oder der COBYQA-Iterationsobergrenze verbessert dies, allerdings zu zusätzlichen klassischen Kosten.

n = 30
dt = 0.6
time_steps = 20
center = n // 2 - 1

# Same MPS settings as the original large-scale run: a larger bond for the
# longer, more-entangled chain (shared by GS prep and AQC compression).
mps_max_bond = 128
mps_cutoff = 1e-8

# Same KCuF3 Hamiltonian and ground-state prep, on a larger chain
H = SparsePauliOp.from_sparse_list(
[(p, [i, i + 1], 0.25) for i in range(n - 1) for p in ("XX", "YY", "ZZ")],
num_qubits=n,
)
gs_circuit = prepare_ground_state(
n, gs_layers=5, max_bond=mps_max_bond, cutoff=mps_cutoff
)
gs_circuit.rz(np.pi / 2, center) # neutron kick at the center site
print(
f"Prepared {n}-qubit ground state with the neutron kick at site {center}."
)
DMRG ground-state energy: -13.111355
GS fidelity: 0.8201
Prepared 30-qubit ground state with the neutron kick at site 14.

Schritte 2 und 3: Komprimieren und Ausführen mit der Funktionsvorlage

Derselbe einzelne Aufruf wie im Simulator-Beispiel, nun mit backend_name, der auf einen IBM Quantum-Prozessor zeigt, sodass die Funktion dort transpiliert und ausgeführt wird. Der Kompressionsplan variiert die Ansatz-Tiefe: Die ersten sechs (schwach verschränkten) Trotter-Schritte werden zu einem flachen einlagigen Ansatz komprimiert, die nächsten vier zu einem tieferen zweilagigen Ansatz, und die verbleibenden 10 der 20 Schritte laufen als reines Trotter. aqc_options erhöht die MPS-Bindungsdimension auf max_bond=128 für die längere, stärker verschränkte Kette (passend zum Original) und behält denselben L-BFGS-B-Optimierer bei, begrenzt auf 100 Iterationen. Die estimator_options aktivieren die integrierte Fehlerminderung: Dynamical Decoupling (XY4), Gate-Twirling und TREX-Messfehlerminderung. Die Standardwerte der Funktion stimmen bereits mit dem Original-Tutorial für all diese Punkte überein, außer beim TREX-Lernbudget (measure_noise_learning). Der gesamte Block wird trotzdem ausgeschrieben, weil vom Aufrufer bereitgestellte estimator_options die Standardwerte der Funktion vollständig ersetzen, statt mit ihnen zusammengeführt zu werden. Würde man also einen Schlüssel weglassen, griffe stattdessen der Standard von IBM Quantum Compute statt der Funktionsstandard.

# Steps 2 + 3: the function compresses (varied ansatz) and executes on hardware.
job = fn.run(
t_steps=time_steps,
aqc_segments=[
{
"n_steps": 6,
"ansatz_steps": 1,
}, # early steps -> shallow 1-layer ansatz
{
"n_steps": 4,
"ansatz_steps": 2,
}, # later steps -> deeper 2-layer ansatz
],
aqc_options={
"max_bond": mps_max_bond, # 128 for the longer chain
"cutoff": mps_cutoff,
"optimizer_settings": {
"method": "L-BFGS-B",
"jac": True,
"options": {"maxiter": 100},
},
},
dt=dt,
hamiltonian=H,
initial_state=gs_circuit,
backend_name="ibm_pittsburgh",
# Mitigation settings from the original tutorial. Only the two
# measure_noise_learning values differ from the function's defaults; the rest
# restates them, because a caller-supplied estimator_options dict replaces the
# function's defaults wholesale rather than merging into them.
estimator_options={
"environment": {"job_tags": ["TUT-SNS"]},
"dynamical_decoupling": {"enable": True, "sequence_type": "XY4"},
"twirling": {
"enable_gates": True,
"num_randomizations": 1000,
"shots_per_randomization": 128,
},
"resilience": {
"measure_mitigation": True,
"measure_noise_learning": {
"num_randomizations": 32,
"shots_per_randomization": 100,
},
},
},
)
print("job ID (save this to reconnect later):", job.job_id)
job ID (save this to reconnect later): 43ed8d07-6d7d-4f33-b70a-7f31b765b310
Wiederverbindung mit einem lange laufenden Job

Der großskalige Lauf ist nicht schnell, und der Großteil der Zeit ist klassisch statt auf der QPU. Die AQC-Kompression läuft innerhalb der Funktion, bevor überhaupt etwas die QPU erreicht: Bei 30 Stellen mit max_bond=128 dauerte das in unserem Lauf fast vier Stunden, gegenüber den rund 18 Minuten QPU-Zeit, die in der Nutzungsschätzung am Anfang dieses Tutorials angegeben sind. Die Warteschlangenzeit kommt bei beidem noch hinzu. Du musst dieses Notebook oder den Kernel während der Ausführung nicht geöffnet lassen.

Kopiere die von der vorherigen Zelle ausgegebene Job-ID und speichere sie. Mit den nächsten drei Zellen kannst du den Lauf später wieder aufnehmen:

  1. Wiederverbinden, nur in einer neuen Kernel-Sitzung nötig: Führe die Setup-Zellen erneut aus, um serverless neu zu erstellen, und baue dann das job-Handle anhand der gespeicherten ID neu auf. Überspringe diese Zelle, wenn du dich noch in der Sitzung befindest, in der du den Auftrag übermittelt hast, da das Handle bereits aktiv ist.

  2. Status prüfen: Führe die Zelle erneut aus, bis DONE gemeldet wird.

  3. Ergebnis abrufen: Führe dies erst aus, sobald der Status DONE ist.

Die folgende Wiederverbindungszelle enthält einen Platzhalter. Ersetze ihn durch deine eigene job_id:

# Reconnect to a previously submitted job by its ID. Only needed in a NEW kernel
# session; if you are still in the session where you submitted, the `job` handle
# from the preceding cell is already live, so skip this cell. Replace the ID that follows with your own.
job = serverless.get_job_by_id("<your job ID>")
# Check where the job is. Re-run this until it reports DONE before fetching the
# result in the following cell: QUEUED -> INITIALIZING -> RUNNING: OPTIMIZING_FOR_HARDWARE ->
# RUNNING: WAITING_FOR_QPU -> RUNNING: EXECUTING_QPU -> RUNNING: POST_PROCESSING
# -> DONE.
print(job.status())
DONE
# Run this only once the preceding status cell reports DONE. result() blocks until
# the job finishes, so calling it earlier just waits (possibly for hours).
result = job.result()
print(
"AQC fidelities:",
{k: round(v, 4) for k, v in result["metadata"]["aqc_fidelities"].items()},
)

ev = np.array(result["expectation_values"])
Gjjc = ev[1:] # drop the t = 0 row -> shape (time_steps, n)
AQC fidelities: {'1': 1.0, '2': 0.9994, '3': 0.9944, '4': 0.9853, '5': 0.9747, '6': 0.959, '7': 0.9495, '8': 0.9542, '9': 0.9533, '10': 0.9451}

Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

Identische Nachbearbeitung wie beim Simulatorlauf: Die Green'sche Funktion wird per Fourier-Transformation in S(q,ω)S(q, \omega) umgewandelt, spiegelsymmetrisiert und negative Werte werden abgeschnitten. Durch die längere Kette und die längere Zeitentwicklung ist das Zwei-Spinon-Kontinuum deutlich besser aufgelöst. Es sollte das Band zwischen den gestrichelten Grenzen ausfüllen, am hellsten in der Nähe von q=πq = \pi.

n = result["metadata"]["n"]
q_res, w_res = 100, 100
spectrum = get_spectrum(n, Gjjc, dt, time_steps, q_res, w_res)
spectrum = -(spectrum + spectrum[:, ::-1]) / 2 # mirror symmetry
spectrum = np.clip(spectrum, a_min=0, a_max=None) # clip negatives

plot_green(
n,
Gjjc,
time_steps,
dt,
title=f"Retarded Green's function - {n} qubits (AQC, hardware)",
)
plot_spectrum(
spectrum,
dt,
q_res,
w_res,
lower_bound=True,
upper_bound=True,
title=f"Dynamical structure factor - {n} qubits (AQC, hardware)",
)

Ausgabe der vorherigen Codezelle

Ausgabe der vorherigen Codezelle

Anhang

Das vorangegangene Hardware-Beispiel läuft mit einer einzigen Kettenlänge. Die drei folgenden Spektren stammen aus früheren Hardware-Läufen desselben Workflows auf ibm_pittsburgh mit 10, 20 und 30 Stellen, wobei alle anderen Eingaben unverändert blieben: 20 Trotter-Schritte mit dt = 0.6, der Kompressionsplan aus sechs einlagigen plus vier zweilagigen AQC-komprimierten Schritten und max_bond = 128. Dies sind aufgezeichnete Ergebnisse, keine Ausgabe der vorangegangenen Zellen.

Bei allen drei Größen werden dieselben Einstellungen verwendet, sodass die Spektren direkt vergleichbar sind. Eine Anpassung pro Kettenlänge, zum Beispiel mit mehr Grundzustands-Ansatz-Schichten oder einem größeren max_bond, kann bessere Ergebnisse liefern als alle hier gezeigten.

Dynamischer Strukturfaktor bei 10 Stellen, ein einzelner scharfer heller Peak bei q = pi nahe der unteren Grenze

Dynamischer Strukturfaktor bei 20 Stellen, spektrales Gewicht füllt das Band zwischen den beiden gestrichelten Zwei-Spinon-Grenzen

Dynamischer Strukturfaktor bei 30 Stellen, das Kontinuum feiner aufgelöst mit schwächerem Kontrast und etwas Gewicht außerhalb der Grenzen

Nächste Schritte

Empfehlungen