Simulasikan hamburan neutron pada material kuantum dengan sirkuit kuantum

Perkiraan penggunaan: 13 menit pada prosesor Heron r2 (CATATAN: Ini hanya perkiraan. Runtime Anda bisa berbeda.)

Tutorial mana yang harus saya gunakan?

Gunakan tutorial ini untuk mempelajari implementasinya langkah demi langkah, dengan komputasi klasik berjalan secara lokal di laptop Anda dan eksekusi hardware pada QPU IBM Quantum®. Contoh skala besar memerlukan memori yang cukup besar, dan approximate quantum compiling (AQC) dapat memakan waktu beberapa jam pada sumber daya laptop biasa. Untuk mengalihkan kompresi ke sumber daya komputasi dan memori Qiskit Serverless, gunakan tutorial Serverless pendampingnya.

Hasil pembelajaran

• Bagaimana spektrum hamburan neutron inelastis (INS) berhubungan dengan faktor struktur dinamis (DSF) dari model spin kuantum.

• Bagaimana menyiapkan keadaan dasar, menerapkan gangguan lokal, dan melakukan evolusi waktu Trotter pada sirkuit kuantum.

• Bagaimana menggunakan approximate quantum compiling (AQC) dengan qiskit-addon-aqc-tensor untuk mengompresi sirkuit Trotter yang dalam untuk eksekusi perangkat keras.

• Bagaimana mengekstrak fungsi Green tertunda (RGF) dari nilai ekspektasi qubit dan mentransformasikan Fourier menjadi DSF.

Prerequisites

• Dasar-dasar informasi kuantum

• Desain algoritma variasional

• Pengantar primitives Qiskit (Estimator dan Sampler)

Background

Dalam tutorial ini, kita mereproduksi hasil dari Lee et al., arXiv:2603.15608.

Hamburan neutron membantu peneliti mengkarakterisasi material yang relevan dengan keberlanjutan, termasuk elektroda baterai. Tutorial ini mengeksplorasi bagaimana circuit kuantum dapat mensimulasikan eksitasi magnetik dan menghubungkan model teoretis dengan pengukuran hamburan. Contoh magnet kuantum dalam tutorial ini mengembangkan metode untuk mempelajari perilaku material, berkontribusi pada toolkit penelitian yang lebih luas untuk penemuan material.

Hamburan neutron inelastik dan faktor struktur dinamis

Hamburan neutron inelastis (INS) adalah salah satu probe eksperimental paling kuat untuk eksitasi magnetik dalam material kuantum. Ketika berkas neutron termal atau dingin menabrak kristal, neutron individu bertukar momentum $\mathbf{q}$ dan energi $\omega$ dengan subsistem magnetik. Intensitas hamburan yang terukur berbanding lurus dengan faktor struktur dinamis (DSF),

$S^{\alpha\beta}(q,\omega) = \sum_{j} e^{-iq\,j}\int_{-\infty}^{\infty} dt\; e^{i\omega t}\, \langle S_0^{\alpha}(0)\, S_j^{\beta}(t)\rangle,$

yang mengkodekan korelasi ruang-waktu penuh dari derajat kebebasan spin.

KCuF$_3$: magnet cairan-Luttinger kanonik

Kalium tembaga fluorida (KCuF$_3$) adalah antiferomagnet kuasi-satu-dimensi di mana rantai ion Cu$^{2+}$ spin-$\frac{1}{2}$ berinteraksi melalui pertukaran Heisenberg tetangga terdekat $J$, sementara kopling antar-rantai hanya $\sim 2.7\%$ dari $J$. Pada $T = 6\;\mathrm{K}$ di mana data INS tersedia, spektrum didominasi oleh eksitasi spinon terfraksionasi yang khas dari cairan Tomonaga-Luttinger. Karena dinamika intrachain digambarkan dengan baik oleh Hamiltonian XXZ spin-$\frac{1}{2}$ satu dimensi pada titik isotropik ($\epsilon = 1$),

$H = J\sum_{i}\left[S_i^Z S_{i+1}^Z + \epsilon\left(S_i^X S_{i+1}^X + S_i^Y S_{i+1}^Y\right)\right],$

KCuF$_3$ berfungsi sebagai benchmark ideal untuk simulasi kuantum: Hamiltonian-nya cukup sederhana untuk diimplementasikan pada prosesor kuantum, namun keadaan dasarnya sangat terbelit dan spektrum eksitasinya menampilkan kontinum two-spinon yang lebar.

Catatan: tutorial ini menetapkan $J = 1$ sebagai satuan energi dan mengadopsi normalisasi $H = J\sum_i[\ldots]$, yang sesuai dengan Hamiltonian lokal $H_\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ yang digunakan dalam implementasi sirkuit paper (Gambar S3 pada suplementari). Hamiltonian penuh paper (Persamaan 3) memiliki faktor keseluruhan tambahan sebesar 2, sehingga $J$ pada paper adalah dua kali $J$ yang digunakan di sini.

Apa yang kita simulasikan dan ukur

Kuantitas fisik yang kita hitung adalah fungsi Green tertunda (RGF), didefinisikan sebagai fungsi korelasi spin-spin yang bergantung waktu

$G^R_{\alpha,\beta}(j, j_c, t) = -\frac{i}{2}\,\langle\psi_{\mathrm{GS}}|\,S_j^\alpha(t)\,S_{j_c}^\beta(0) - S_{j_c}^\beta(0)\,S_j^\alpha(t)\,|\psi_{\mathrm{GS}}\rangle,$

di mana $j_c$ adalah situs referensi (pusat rantai) dan $S_j^\alpha(t) = e^{iHt}S_j^\alpha e^{-iHt}$ adalah operator spin dalam gambaran Heisenberg. Dalam tutorial ini kita berfokus pada komponen $zz$ ($\alpha = \beta = z$). Pada komputer kuantum, RGF diakses dengan menyiapkan keadaan dasar, menerapkan gangguan lokal pada $j_c$, mengevolusikan waktu keadaan yang terganggu, dan mengukur nilai ekspektasi qubit tunggal $\langle\sigma_j^z\rangle$ pada setiap situs $j$ untuk setiap langkah waktu. Wawasan kuncinya adalah setiap $\langle\sigma_j^z\rangle$ memberikan selisih dari magnetisasi keadaan dasar; karena antiferomagnet Heisenberg isotropik memiliki magnetisasi per-situs bersih nol ($\langle\sigma_j^z\rangle_\mathrm{GS} = 0$), nilai yang terukur secara langsung memberikan $G^R(j, j_c, t)$ tanpa pengurangan eksplisit apa pun.

Dengan mengumpulkan $G^R(j, j_c, t)$ di semua situs dan langkah waktu, kita menyusun dataset dua dimensi yang kemudian ditransformasikan Fourier baik dalam ruang maupun waktu untuk menghasilkan faktor struktur dinamis $S(q,\omega)$. DSF adalah kuantitas yang secara langsung diukur dalam eksperimen INS: ia memberi tahu kita eksitasi magnetik apa yang ada pada setiap momentum $q$ dan energi $\omega$. Untuk rantai Heisenberg isotropik, spektrum eksitasi eksaknya adalah kontinum two-spinon, sebuah pita lebar intensitas hamburan yang bentuknya berfungsi sebagai benchmark end-to-end yang ketat untuk simulasi kuantum: ia memvalidasi penyiapan keadaan dasar, gangguan, evolusi waktu Trotter, dan protokol pengukuran sekaligus.

Alur kerja simulasi kuantum

Alur kerja ini mencerminkan fisika dari peristiwa INS. Kita (1) menyiapkan keadaan dasar banyak-benda $|\psi_{\mathrm{GS}}\rangle$ pada $n$ qubit, (2) menerapkan gangguan spin-flip lokal $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ pada pusat rantai untuk meniru transfer spin neutron, (3) mengevolusikan di bawah $H$ untuk langkah waktu diskret menggunakan Trotterisasi orde-kedua, dan (4) mengukur $\langle\sigma_i^z\rangle$ pada setiap qubit di setiap langkah untuk mendapatkan RGF. Transformasi Fourier diskret dua dimensi kemudian menghasilkan DSF $S(q,\omega)$.

Observable yang kita ukur pada setiap langkah waktu adalah $\sigma_i^z$ pada setiap qubit $i$. Dalam Qiskit ini direpresentasikan sebagai daftar operator SparsePauliOp: satu $Z$ qubit-tunggal yang disematkan dalam string identitas $n$-qubit untuk setiap situs. Observable ini dibangun sekali per instans masalah selama langkah pemetaan masalah (Langkah 1) dan digunakan kembali untuk setiap sirkuit pada skala tersebut.

Approximate quantum compiling (AQC)

Sirkuit Trotter yang dalam dapat dikompresi dengan approximate quantum compiling (AQC), yang menggantikan beberapa lapisan Trotter pertama dengan ansatz terparameterisasi yang lebih pendek yang parameternya dioptimalkan secara klasik untuk memaksimalkan fidelitas level-MPS dengan sirkuit dalam yang asli. Langkah Trotter yang tersisa ditambahkan secara eksak, menghasilkan sirkuit "campuran" AQC + Trotter dengan jauh lebih sedikit gerbang dua-qubit.

Simulasi MPS

Untuk sistem satu dimensi, metode matrix-product-state (MPS) dapat secara efisien mensimulasikan baik persiapan ground-state (dengan density matrix renormalization group, atau DMRG) maupun evolusi waktu tingkat circuit. Dengan mengontrol bond dimension $\chi$, kita menukar (trade-off) akurasi dengan biaya komputasi. Dalam tutorial ini kita menggunakan simulasi MPS dalam qiskit-addon-aqc-tensor untuk menghitung ansatz AQC berfidelitas tinggi yang mengompresi circuit Trotter yang dalam untuk eksekusi hardware.

Persyaratan

Sebelum memulai tutorial ini, pastikan Anda telah menginstal hal berikut:

• Qiskit SDK dengan dukungan visualisasi

• Qiskit Runtime (pip install qiskit-ibm-runtime)

• qiskit-addon-aqc-tensor dengan ekstra quimb dan JAX (pip install 'qiskit-addon-aqc-tensor[quimb-jax]')

Persiapan

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_addon_aqc_tensor": "qiskit-addon-aqc-tensor", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "qiskit_quimb": "qiskit-quimb", "quimb": "quimb", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import timeit
import warnings
from collections.abc import Iterator, Sequence
from functools import partial

import matplotlib.pyplot as plt
import numpy as np
import quimb.tensor as qtn
import scipy.optimize
from numpy.typing import NDArray
from qiskit import QuantumCircuit
from qiskit.circuit import CircuitInstruction, ParameterVector, Qubit
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.primitives import StatevectorEstimator
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_addon_aqc_tensor import generate_ansatz_from_circuit
from qiskit_addon_aqc_tensor.objective import MaximizeStateFidelity
from qiskit_addon_aqc_tensor.simulation import tensornetwork_from_circuit
from qiskit_addon_aqc_tensor.simulation.quimb import QuimbSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_quimb import quimb_circuit
from scipy.sparse import SparseEfficiencyWarning

# scipy.sparse.linalg.expm (reached via PauliEvolutionGate.to_matrix) prefers
# CSC and warns when handed the CSR matrix from SparsePauliOp.to_matrix; the
# result is unaffected, so silence the cosmetic warning.
warnings.filterwarnings("ignore", category=SparseEfficiencyWarning)

def xxz_hamiltonian_mpo(
    n_qubits: int, interaction: float = 1.0, anisotropy: float = 1.0
) -> qtn.MatrixProductOperator:
    """1D XXZ Hamiltonian as a quimb MPO.

    Builds the Hamiltonian using ``qtn.SpinHam1D``.

    Args:
        n_qubits: Number of sites.
        interaction: Overall interaction strength (J in the paper).
        anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the
            isotropic Heisenberg point; ``ε = 0`` is the Ising limit.

    Returns:
        The Hamiltonian as a matrix product operator.
    """
    builder = qtn.SpinHam1D(S=1 / 2)
    builder += interaction * anisotropy * 0.5, "+", "-"
    builder += interaction * anisotropy * 0.5, "-", "+"
    builder += interaction, "Z", "Z"
    return builder.build_mpo(L=n_qubits)

def build_ground_state_ansatz(n_qubits: int, n_layers: int) -> QuantumCircuit:
    """Hamiltonian variational ansatz (HVA) circuit for ground-state preparation.

    Starts from a product of singlet pairs and applies alternating
    odd/even layers of parameterized XXZ pair-evolution gates.

    The returned circuit is parameterized: it carries a ``ParameterVector``
    named ``"theta"`` of length ``2 * n_layers`` whose values must be
    assigned (e.g. via ``circuit.assign_parameters``) before simulation.
    Element ``2 * r`` is the odd-layer angle and element ``2 * r + 1`` is the
    even-layer angle of layer ``r``.

    Args:
        n_qubits: Number of qubits (must be even).
        n_layers: Number of HVA layers.

    Returns:
        The parameterized HVA preparation circuit.
    """
    theta = ParameterVector("theta", 2 * n_layers)
    circuit = QuantumCircuit(n_qubits)
    # Initial singlet product state
    for i in range(n_qubits // 2):
        circuit.x(2 * i)
        circuit.x(2 * i + 1)
        circuit.h(2 * i + 1)
        circuit.cx(2 * i + 1, 2 * i)
    # Variational layers: each pair gate is exp(-i theta (XX+YY+ZZ)/2)
    pair_ham = SparsePauliOp(
        ["XX", "YY", "ZZ"], coeffs=[0.5, 0.5, 0.5]
    )  # H_pair (HVA form)
    for r in range(n_layers):
        for i in range(1, (n_qubits + 1) // 2):  # odd layer
            circuit.append(
                PauliEvolutionGate(pair_ham, time=theta[2 * r]),
                [2 * i - 1, 2 * i],
            )
        for i in range(n_qubits // 2):  # even layer
            circuit.append(
                PauliEvolutionGate(pair_ham, time=theta[2 * r + 1]),
                [2 * i, 2 * i + 1],
            )
    return circuit

def optimize_ground_state_ansatz(
    ansatz: QuantumCircuit,
    x0: NDArray[np.floating],
    target_mps: qtn.MatrixProductState,
    *,
    max_bond: int | None = None,
    cutoff: float = 1e-10,
    method: str = "COBYQA",
    options: dict | None = None,
) -> scipy.optimize.OptimizeResult:
    """Optimize HVA parameters by maximizing fidelity with a target MPS.

    The HVA circuit is simulated as a matrix product state with the given
    bond-dimension truncation, and the parameters are optimized to maximize
    the state fidelity ``|<psi_HVA | target_mps>|**2`` with the DMRG ground
    state. Both states are normalized, so the minimized objective is the
    infidelity ``1 - |<psi_HVA | target_mps>|**2``.

    Args:
        ansatz: Parameterized HVA circuit from ``build_ground_state_ansatz``.
            The length of ``x0`` must equal ``ansatz.num_parameters``.
        x0: Initial parameters.
        target_mps: Target MPS (DMRG ground state) to maximize fidelity with.
        max_bond: Maximum MPS bond dimension during gate application.
        cutoff: Singular-value cutoff during gate application.
        method: ``scipy.optimize.minimize`` method.
        options: Options dict forwarded to ``scipy.optimize.minimize``.

    Returns:
        The Scipy OptimizeResult.
    """

    def infidelity(params: NDArray[np.floating]) -> float:
        circuit = ansatz.assign_parameters(params)
        circuit_mps = quimb_circuit(
            circuit.decompose(["PauliEvolution"]),
            quimb_circuit_class=qtn.CircuitMPS,
            max_bond=max_bond,
            cutoff=cutoff,
        )
        return 1 - abs(circuit_mps.psi.H @ target_mps) ** 2

    return scipy.optimize.minimize(
        infidelity, np.asarray(x0), method=method, options=options
    )

def trotter_evolution(
    qubits: Sequence[Qubit],
    interaction: float,
    anisotropy: float,
    time_step: float,
    n_steps: int,
) -> Iterator[CircuitInstruction]:
    """Second-order Trotter steps of the XXZ pair Hamiltonian.

    While the paper used a hand-optimized circuit for the Trotter steps, we use
    PauliEvolutionGate here for simplicity and generality. The final two-qubit gate
    count and gate depth are equivalent when transpiled with ``optimization_level=3``.

    Args:
        qubits: Qubits to act on (length ``n_qubits``).
        interaction: Overall interaction strength (J in the paper).
        anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the
            isotropic Heisenberg point; ``ε = 0`` is the Ising limit.
        time_step: Per-step Trotter time.
        n_steps: Number of Trotter steps.

    Yields:
        ``CircuitInstruction``s implementing the Trotter steps.
    """
    if n_steps == 0:
        return
    n_qubits = len(qubits)
    pair_ham = SparsePauliOp(
        ["XX", "YY", "ZZ"],
        coeffs=[
            0.25 * interaction * anisotropy,
            0.25 * interaction * anisotropy,
            0.25 * interaction,
        ],
    )
    half_evo = PauliEvolutionGate(pair_ham, time=time_step / 2)
    full_evo = PauliEvolutionGate(pair_ham, time=time_step)
    for i in range(n_qubits // 2):  # half even layer
        yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))
    for i in range(n_qubits // 2 - 1):  # full odd layer
        yield CircuitInstruction(
            full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])
        )
    for _ in range(n_steps - 1):  # interior steps
        for i in range(n_qubits // 2):
            yield CircuitInstruction(
                full_evo, (qubits[2 * i], qubits[2 * i + 1])
            )
        for i in range(n_qubits // 2 - 1):
            yield CircuitInstruction(
                full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])
            )
    for i in range(n_qubits // 2):  # half even layer
        yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))

def get_dsf(
    n_qubits: int,
    rgf_mat: NDArray[np.floating],
    time_step: float,
    n_steps: int,
    n_points_momentum: int,
    n_points_frequency: int,
) -> NDArray[np.floating]:
    """Compute the dynamical structure factor from the retarded Green's function.

    Uses the center-site approximation and a discrete Fourier transform.
    The result is symmetrized about the momentum axis and clipped to
    non-negative values, ready for plotting.

    Args:
        n_qubits: Number of qubits (sites).
        rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.
        time_step: Trotter time-step size.
        n_steps: Number of time steps.
        n_points_momentum: Number of momentum points.
        n_points_frequency: Number of frequency points.

    Returns:
        DSF array of shape ``(n_points_frequency, n_points_momentum)``,
        symmetrized about the momentum axis and clipped to non-negative values.
    """
    max_frequency = np.pi / time_step
    momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)
    frequency_range = np.linspace(0, max_frequency, n_points_frequency)
    result = np.zeros((frequency_range.shape[0], momentum_range.shape[0]))
    center = n_qubits // 2 - 1
    for iw, w in enumerate(frequency_range):
        exponent = np.exp(1j * w * time_step * np.arange(1, n_steps + 1))
        # S = sigma/2, so two spin operators contribute a factor of (1/2)(1/2) = 1/4.
        rgf_omega = (
            np.dot(rgf_mat.T, exponent) * time_step / 4
        )  # S(omega): time Fourier slice of the Green's function
        for iq, q in enumerate(momentum_range):
            momentum_phases = np.exp(
                -1j * q * np.arange(-center, center + 2, 1)
            )
            result[iw, iq] = np.imag(np.dot(rgf_omega, momentum_phases))
    result = -(result + result[:, ::-1]) / 2
    result = np.clip(result, a_min=0, a_max=None)
    return result

def plot_dsf(
    dsf: NDArray[np.floating],
    time_step: float,
    n_points_momentum: int,
    n_points_frequency: int,
    title: str | None = None,
) -> None:
    """Heat-map of the dynamical structure factor.

    Args:
        dsf: DSF array of shape ``(n_points_frequency, n_points_momentum)``.
        time_step: Trotter time-step size.
        n_points_momentum: Number of momentum points.
        n_points_frequency: Number of frequency points.
        title: Optional plot title.
    """
    max_frequency = np.pi / time_step
    momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)
    frequency_range = np.linspace(0, max_frequency, n_points_frequency)
    x, y = np.meshgrid(momentum_range, frequency_range)
    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")
    ax.set_ylim(0, 3.6)
    ax.set_xlim(0, 2 * np.pi)
    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 title:
        ax.set_title(title, fontsize=14)
    plt.tight_layout()
    plt.show()

def plot_rgf(
    n_qubits: int,
    rgf_mat: NDArray[np.floating],
    time_step: float,
    n_steps: int,
    title: str | None = None,
) -> None:
    """Heat-map of the retarded Green's function in real space and time.

    Args:
        n_qubits: Number of qubits (sites).
        rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.
        time_step: Trotter time-step size.
        n_steps: Number of time steps.
        title: Optional plot title.
    """
    fig, ax = plt.subplots(figsize=(8, 6))
    qubit_axis = np.arange(n_qubits)
    t_axis = np.arange(1, n_steps + 1) * time_step
    x, y = np.meshgrid(qubit_axis, t_axis)
    c = ax.pcolormesh(
        x,
        y,
        np.real(rgf_mat),
        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("Qubit", fontsize=16)
    ax.xaxis.set_major_locator(
        plt.matplotlib.ticker.MaxNLocator(integer=True)
    )
    ax.set_ylabel(r"Time", fontsize=16)
    if title:
        ax.set_title(title, fontsize=14)
    plt.tight_layout()
    plt.show()

def uniform_2q_depth(circuit: QuantumCircuit) -> int:
    """Two-qubit gate depth in a standardized basis."""
    pass_manager = generate_preset_pass_manager(
        optimization_level=0, basis_gates=["cz", "id", "rz", "sx", "x"]
    )
    return pass_manager.run(circuit).depth(
        lambda inst: inst.operation.num_qubits == 2
    )

Contoh simulator skala kecil

Kita pertama-tama mendemonstrasikan alur kerja penuh pada 10 qubit, mengoptimalkan ansatz ground-state HVA dengan menggunakan simulasi MPS dan menggunakan simulator statevector Qiskit untuk evolusi waktu. DMRG memberikan energi ground-state referensi dan MPS referensi. Contoh skala kecil ini memungkinkan kita memvalidasi setiap langkah sebelum memperbesar skala.

Langkah 1: Memetakan input klasik ke masalah kuantum

Kita mulai dengan mendefinisikan model fisik dan membangun circuit kuantum.

Hamiltonian. KCuF$_3$ dimodelkan oleh Hamiltonian XXZ 1D pada titik isotropik ($J = 1$, $\epsilon = 1$, menetapkan $J = 1$ sebagai satuan energi).

Ground state. Kita menggunakan circuit Hamiltonian variational ansatz (HVA), build_ground_state_ansatz, sebagai circuit persiapan ground-state. Parameter HVA dioptimalkan secara klasik dengan memaksimalkan fidelitas state $|\langle\psi_{\mathrm{HVA}}(\theta)|\psi_{\mathrm{DMRG}}\rangle|^2$ dengan MPS ground-state DMRG, di mana state HVA dievaluasi dengan mensimulasikan circuit sebagai matrix product state menggunakan CircuitMPS milik quimb. Kita menggunakan scipy.optimize.minimize untuk optimasi. Energi $\langle H \rangle$ dari ansatz yang telah dioptimalkan juga dihitung sebagai referensi.

Trotter gates. Setiap suku interaksi nearest-neighbor $e^{-i\Delta t\, H_{\mathrm{pair}}}$ dengan $H_{\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ dibangun dengan PauliEvolutionGate(H_pair, time=time_step). Qiskit mensintesis ini menjadi dekomposisi tiga-CNOT yang optimal selama transpilasi.

Perturbasi. Sebuah gate $R_z(\pi/2)$ pada qubit tengah mengimplementasikan $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$, meniru spin-flip lokal yang dihasilkan oleh neutron yang terhambur.

Observabel. Kita membangun observabel $\sigma^z$ untuk setiap situs qubit. Objek SparsePauliOp ini diteruskan ke primitive Estimator pada Langkah 3 untuk mengekstrak $\langle\sigma_i^z\rangle$ pada setiap langkah waktu.

# -- Physical parameters --
n_qubits = 10
interaction = 1.0  # J
anisotropy = 1.0  # ε (isotropic point)
time_step = 0.6
n_steps = 10
mps_max_bond = 32
mps_cutoff = 1e-8
center = n_qubits // 2 - 1

# -- Hamiltonian MPO --
ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)

# -- Reference ground-state energy via DMRG --
dmrg = qtn.DMRG2(ham_mpo)
dmrg.solve(tol=1e-8)
print(f"Ground-state energy (DMRG): {dmrg.energy:.6f}")

# -- Build ground state ansatz circuit --
gs_n_layers = 3
gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)

# -- Optimize ground state ansatz parameters to maximize fidelity with the DMRG MPS --
# Initialize odd-layer angles near 0 (where the inter-pair gate is the
# identity) and even-layer angles near pi/2 (where the intra-pair gate
# is a SWAP, since 0.5*(XX + YY + ZZ) = SWAP - I/2).
rng = np.random.default_rng(12345)
x0 = np.tile([0, np.pi / 2], gs_n_layers) + rng.normal(
    scale=0.1, size=2 * gs_n_layers
)
print("Optimizing ground state ansatz...")
t0 = timeit.default_timer()
result = optimize_ground_state_ansatz(
    gs_ansatz,
    x0,
    dmrg.state,
    max_bond=mps_max_bond,
    cutoff=mps_cutoff,
    options=dict(maxiter=100),
)
t1 = timeit.default_timer()
print(f"Finished optimizing ground state ansatz in {t1 - t0} seconds.")
print(f"Ground state ansatz fidelity: {1 - result.fun:.6f}")

gs_circuit = gs_ansatz.assign_parameters(result.x)
gs_circuit_mps = quimb_circuit(
    gs_circuit.decompose(["PauliEvolution"]),
    quimb_circuit_class=qtn.CircuitMPS,
    max_bond=mps_max_bond,
    cutoff=mps_cutoff,
)
gs_ansatz_energy = qtn.expec_TN_1D(
    gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi
)
print(f"Ground state ansatz energy: {gs_ansatz_energy:.6f}")

# -- Build circuits for each time step --
perturbed = gs_circuit.copy()
perturbed.rz(np.pi / 2, center)

circuits = []
for t in range(1, n_steps + 1):
    circuit = perturbed.copy()
    for instr in trotter_evolution(
        circuit.qubits, interaction, anisotropy, time_step, t
    ):
        circuit.append(instr)
    circuits.append(circuit)
print(
    f"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = "
    f"{uniform_2q_depth(circuits[-1])}"
)

# -- Observables: Z on each qubit site --
observables = [
    SparsePauliOp.from_sparse_list([("Z", [i], 1)], num_qubits=n_qubits)
    for i in range(n_qubits)
]
print(f"Defined {len(observables)} Z observables.")

Ground-state energy (DMRG): -4.258035
Optimizing ground state ansatz...
Finished optimizing ground state ansatz in 8.850687621976249 seconds.
Ground state ansatz fidelity: 0.984277
Ground state ansatz energy: -4.232565
Built 10 circuits, deepest 2q depth (uniform basis) = 163
Defined 10 Z observables.

Langkah 2: Mengoptimalkan masalah untuk eksekusi hardware kuantum

Untuk hardware nyata, circuit Trotter di atas akan terlalu dalam. Approximate quantum compiling (AQC) mengatasi ini dengan mengganti $k$ lapisan Trotter pertama (termasuk circuit ground-state) dengan ansatz terparameterisasi yang lebih pendek yang dioptimalkan untuk memaksimalkan fidelitas tingkat MPS dengan circuit dalam aslinya. Langkah-langkah Trotter yang tersisa ditambahkan secara tepat, menghasilkan circuit "AQC + Trotter" yang lebih dangkal.

Untuk menyeimbangkan ekspresivitas terhadap kedalaman circuit, dua ansatz digunakan: ansatz satu lapisan (dihasilkan dari satu langkah Trotter) mengompresi langkah waktu paling awal, dan ansatz dua lapisan yang lebih dalam (dihasilkan dari dua langkah Trotter) mengompresi beberapa langkah berikutnya di mana fidelitas yang lebih tinggi diperlukan.

Alur kerja AQC memiliki empat sub-langkah:

1. Membangun circuit target — $k_1 + k_2$ circuit pertama dari Langkah 1 berfungsi langsung sebagai target AQC.

2. Menghitung MPS target — mensimulasikan setiap circuit target sebagai matrix-product state menggunakan quimb.tensor.CircuitMPS.

3. Menghasilkan dan mengoptimalkan ansatz — generate_ansatz_from_circuit membuat ansatz terparameterisasi satu lapisan dan dua lapisan; parameter dioptimalkan dengan menggunakan L-BFGS-B dengan gradien yang dipercepat JAX untuk meminimalkan $1 - |\langle\psi_{\mathrm{ansatz}}|\psi_{\mathrm{target}}\rangle|^2$. $k_1$ langkah pertama menggunakan ansatz satu lapisan dan $k_2$ langkah berikutnya menggunakan ansatz dua lapisan; dalam setiap tahap, setiap langkah melakukan warm-start dari parameter yang dioptimalkan pada langkah sebelumnya, dan parameter direset ke default tahap pada batas tahap.

4. Merakit circuit campuran — untuk langkah waktu di luar checkpoint AQC, tambahkan lapisan Trotter eksak ke circuit AQC dua lapisan yang telah dioptimalkan menggunakan trotter_evolution.

# Number of time steps to compress into an AQC ansatz with one layer
aqc_n_steps_1 = 3
# Number of time steps to compress into an AQC ansatz with two layers
aqc_n_steps_2 = 2
aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2

# ── Step 2a: Target circuits (first aqc_n_steps_total circuits from Step 1) ──
target_circuits = {
    k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)
}
# ── Step 2b: Compute target MPS ──
# Decompose PauliEvolutionGate to RXX/RYY/RZZ before passing to the AQC MPS
# backend, which does not understand PauliEvolutionGate natively.
aqc_sim = QuimbSimulator(
    quimb_circuit_factory=partial(
        qtn.CircuitMPS,
        gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),
    ),
    autodiff_backend="jax",
)
print("\nStep 2b — target MPS:")
target_mps = {}
for k in range(1, aqc_n_steps_total + 1):
    target_mps[k] = tensornetwork_from_circuit(
        target_circuits[k].decompose(["PauliEvolution"]), aqc_sim
    )
    print(f"  k={k}: max bond = {target_mps[k].psi.max_bond()}")

# ── Step 2c: Generate ansätze (one-layer and two-layer) and optimize parameters ──
ansatz_1, initial_params_1 = generate_ansatz_from_circuit(
    target_circuits[1].decompose(["PauliEvolution"]),
    qubits_initially_zero=True,
)
initial_params_1 = np.array(initial_params_1)
ansatz_2, initial_params_2 = generate_ansatz_from_circuit(
    target_circuits[2].decompose(["PauliEvolution"]),
    qubits_initially_zero=True,
)
initial_params_2 = np.array(initial_params_2)
print(
    f"\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, "
    f"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}"
)
print(
    f"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, "
    f"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}"
)

aqc_circuits = {}
aqc_params = {}
for k in range(1, aqc_n_steps_total + 1):
    if k <= aqc_n_steps_1:
        ansatz, base_params = ansatz_1, initial_params_1
    else:
        ansatz, base_params = ansatz_2, initial_params_2
    # Warm-start from the previous step only within the same stage
    same_stage = (k - 1 >= 1) and (
        (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)
    )
    x0 = aqc_params[k - 1] if same_stage else base_params
    obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)
    t0 = timeit.default_timer()
    result = scipy.optimize.minimize(
        obj.loss_function,
        x0,
        method="L-BFGS-B",
        jac=True,
        options=dict(maxiter=100),
    )
    elapsed = timeit.default_timer() - t0
    aqc_params[k] = result.x
    aqc_circuits[k] = ansatz.assign_parameters(result.x)
    print(
        f"  k={k}: fidelity = {1 - result.fun:.4f}, "
        f"2q depth (uniform basis) = "
        f"{uniform_2q_depth(aqc_circuits[k])}, "
        f"{elapsed:.1f}s"
    )

# ── Step 2d: Assemble full circuit set (AQC + Trotter) ──
all_circuits = []
for k in range(1, aqc_n_steps_total + 1):
    all_circuits.append(aqc_circuits[k])
base = aqc_circuits[aqc_n_steps_total]
for k in range(1, n_steps - aqc_n_steps_total + 1):
    circuit = base.copy()
    for instr in trotter_evolution(
        circuit.qubits, interaction, anisotropy, time_step, k
    ):
        circuit.append(instr)
    all_circuits.append(circuit)

full_depths = [
    uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)
]
aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]
full_2q = [
    circuits[k - 1].decompose(["PauliEvolution"]).num_nonlocal_gates()
    for k in range(1, n_steps + 1)
]
aqc_2q = [
    circuit.decompose(["PauliEvolution"]).num_nonlocal_gates()
    for circuit in all_circuits
]
print(f"\nStep 2d — assembled {len(all_circuits)} circuits")
print(
    f"  At step {n_steps} (uniform basis): "
    f"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; "
    f"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}"
)

steps_axis = np.arange(1, n_steps + 1)
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(steps_axis, full_depths, "-o", color="black", label="Full Trotter")
ax.plot(
    steps_axis, aqc_depths, "-o", color="cadetblue", label="AQC + Trotter"
)
ax.set_xlabel("Trotter steps", fontsize=13)
ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))
ax.set_ylabel("2q gate depth (uniform basis)", fontsize=13)
ax.set_title(f"AQC circuit-depth reduction ({n_qubits} qubits)", fontsize=13)
ax.legend(fontsize=11)
plt.tight_layout()
plt.show()

Step 2b — target MPS:
  k=1: max bond = 22
  k=2: max bond = 22
  k=3: max bond = 26
  k=4: max bond = 27
  k=5: max bond = 30

Step 2c — one-layer ansatz: 389 parameters, 2q depth (uniform basis) = 27
Step 2c — two-layer ansatz: 470 parameters, 2q depth (uniform basis) = 33
  k=1: fidelity = 1.0000, 2q depth (uniform basis) = 27, 11.8s
  k=2: fidelity = 0.9980, 2q depth (uniform basis) = 27, 13.2s
  k=3: fidelity = 0.9890, 2q depth (uniform basis) = 27, 13.7s
  k=4: fidelity = 0.9983, 2q depth (uniform basis) = 33, 16.4s
  k=5: fidelity = 0.9957, 2q depth (uniform basis) = 33, 16.4s

Step 2d — assembled 10 circuits
  At step 10 (uniform basis): Trotter 2q depth = 163, AQC+Trotter 2q depth = 99; 2q gates: Trotter = 371, AQC+Trotter = 200

Langkah 3: Mengeksekusi menggunakan primitive Qiskit

Kita mensimulasikan setiap circuit yang dikompilasi AQC menggunakan primitive StatevectorEstimator.

estimator = StatevectorEstimator()
pubs = [(circuit, observables) for circuit in all_circuits]
job = estimator.run(pubs)
result = job.result()

Langkah 4: Pasca-pemrosesan dan mengembalikan hasil dalam format klasik yang diinginkan

Sekarang kita mengekstrak nilai ekspektasi $\langle\sigma_i^z\rangle$ untuk setiap qubit $i$ pada setiap langkah waktu $t$. Nilai-nilai ini membentuk matriks fungsi Green tertunda $G^R(j, j_c, t)$. Kemudian, kita melakukan transformasi Fourier RGF menjadi faktor struktur dinamis $S(q, \omega)$ dan memplot baik RGF maupun DSF. get_dsf menerapkan simetri cermin dan memotong nilai negatif secara internal.

rgf_mat = np.stack([pub_result.data.evs for pub_result in result])

# -- Compute DSF --
n_points_momentum, n_points_frequency = 100, 100
spectrum = get_dsf(
    n_qubits,
    rgf_mat,
    time_step,
    n_steps,
    n_points_momentum,
    n_points_frequency,
)

# -- Plot retarded Green's function --
plot_rgf(
    n_qubits,
    rgf_mat,
    time_step,
    n_steps,
    title=f"Retarded Green's function — {n_qubits} qubits (simulation)",
)

# -- Plot DSF --
plot_dsf(
    spectrum,
    time_step,
    n_points_momentum,
    n_points_frequency,
    title=f"Dynamical structure factor — {n_qubits} qubits (simulation)",
)

Eksekusi hardware skala besar

Sekarang kita memperbesar skala ke 50 qubit. Pada skala ini, optimasi mencapai fidelitas yang lebih rendah dibandingkan pada contoh skala kecil: fidelitas ansatz ground-state turun menjadi sekitar 0,65, dan fidelitas AQC menurun menjadi sekitar 0,7 pada checkpoint yang lebih akhir. Hal ini memang wajar, dan Anda dapat meningkatkan fidelitas ini dengan menambah jumlah lapisan ansatz ground-state (gs_n_layers) atau jumlah iterasi optimasi (maxiter) dengan biaya klasik tambahan. Perhatikan juga bahwa fidelitas AQC dua lapisan lebih rendah dibandingkan yang satu lapisan. Ini bukan regresi: langkah waktu yang lebih akhir menghasilkan lebih banyak entanglement dan memang lebih sulit dikompresi, itulah sebabnya ansatz dua lapisan yang lebih ekspresif digunakan untuk langkah-langkah tersebut.

Optimasi AQC juga dapat memakan waktu komputasi klasik beberapa jam (sekitar enam jam pada run yang ditampilkan di sini, sebagian besar dihabiskan pada checkpoint dua lapisan di Langkah 2c). Untuk mengurangi waktu wall-clock, pertimbangkan untuk menjalankan notebook ini pada hardware klasik yang lebih kuat, seperti sistem high-performance computing (HPC). Alternatifnya, Anda dapat memperkecil skala ke instance masalah yang lebih kecil (misalnya, lebih sedikit qubit atau lebih sedikit langkah waktu); hasil Anda kemudian akan berbeda dari yang ditampilkan di sini.

Kode di bawah ini mengikuti struktur empat langkah yang sama seperti contoh skala kecil. Parameter ground-state HVA kembali dioptimalkan dengan menggunakan simulasi MPS. Pada QPU, kita mengaktifkan dynamical decoupling (DD), Pauli twirling, dan twirled readout error extinction (TREX) untuk penekanan dan mitigasi error. Tabel berikut merangkum bagaimana eksperimen skala besar berbeda dari yang skala kecil:

	Skala kecil	Skala besar
Qubit	10	50
Langkah waktu	10	20
Checkpoint AQC (1-lapisan + 2-lapisan)	3 + 2 = 5	6 + 4 = 10
Lapisan ansatz ground state	3	5
Bond dimension maksimum MPS	32	128
Estimator	StatevectorEstimator	QPU dengan DD, Pauli twirling, dan TREX

Pesan kompilasi lambat XLA

Selama optimasi AQC (Langkah 2c di bawah) Anda mungkin melihat pesan stderr seperti:

[Compiling module jit_MakeArrayFn for CPU] Very slow compile? ...
The operation took 2m14s

Ini adalah diagnostik yang tidak berbahaya dari XLA, compiler yang mendukung autodiff JAX milik qiskit-addon-aqc-tensor. Pada 50 qubit dengan bond dimension MPS 128, XLA menghabiskan beberapa menit mengompilasi fungsi gradien saat pertama kali dilacak. Kompilasi berhasil, dan hasil optimasi tidak terpengaruh.

# ── Parameters ──────────────────────────────────────────────────────────────
n_qubits = 50  # 10 → 50
interaction = 1.0  # J
anisotropy = 1.0  # ε (isotropic point)
time_step = 0.6
n_steps = 20  # 10 → 20
aqc_n_steps_1 = 6  # 3 → 6
aqc_n_steps_2 = 4  # 2 → 4
aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2
gs_n_layers = 5  # 3 → 5
mps_max_bond = 128  # 32 -> 128
mps_cutoff = 1e-8
center = n_qubits // 2 - 1

# ── Step 1: Map ──────────────────────────────────────────────────────────────
ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)

dmrg = qtn.DMRG2(ham_mpo)
dmrg.solve(tol=1e-8)
print(f"Ground-state energy (DMRG): {dmrg.energy:.6f}")

gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)

rng = np.random.default_rng(12345)
x0 = np.tile([0.0, np.pi / 2], gs_n_layers) + rng.normal(
    scale=0.1, size=2 * gs_n_layers
)
print("Optimizing ground state ansatz...")
t0 = timeit.default_timer()
result = optimize_ground_state_ansatz(
    gs_ansatz,
    x0,
    dmrg.state,
    max_bond=mps_max_bond,
    cutoff=mps_cutoff,
    options=dict(maxiter=100),
)
t1 = timeit.default_timer()
print(f"Finished optimizing ground state ansatz in {t1 - t0} seconds.")
print(f"Ground state ansatz fidelity: {1 - result.fun:.6f}")

gs_circuit = gs_ansatz.assign_parameters(result.x)
gs_circuit_mps = quimb_circuit(
    gs_circuit.decompose(["PauliEvolution"]),
    quimb_circuit_class=qtn.CircuitMPS,
    max_bond=mps_max_bond,
    cutoff=mps_cutoff,
)
gs_ansatz_energy = qtn.expec_TN_1D(
    gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi
)
print(f"Ground state ansatz energy: {gs_ansatz_energy:.6f}")

perturbed = gs_circuit.copy()
perturbed.rz(np.pi / 2, center)

circuits = []
for t in range(1, n_steps + 1):
    circuit = perturbed.copy()
    for instr in trotter_evolution(
        circuit.qubits, interaction, anisotropy, time_step, t
    ):
        circuit.append(instr)
    circuits.append(circuit)
print(
    f"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = "
    f"{uniform_2q_depth(circuits[-1])}"
)

observables = [
    SparsePauliOp.from_sparse_list([("Z", [i], 1)], num_qubits=n_qubits)
    for i in range(n_qubits)
]
print(f"Defined {len(observables)} Z observables.")

# ── Step 2: AQC ──────────────────────────────────────────────────────────────
target_circuits = {
    k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)
}

aqc_sim = QuimbSimulator(
    quimb_circuit_factory=partial(
        qtn.CircuitMPS,
        gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),
    ),
    autodiff_backend="jax",
)
print("\nStep 2b — target MPS:")
target_mps = {}
for k in range(1, aqc_n_steps_total + 1):
    target_mps[k] = tensornetwork_from_circuit(
        target_circuits[k].decompose(["PauliEvolution"]), aqc_sim
    )
    print(f"  k={k}: max bond = {target_mps[k].psi.max_bond()}")

ansatz_1, initial_params_1 = generate_ansatz_from_circuit(
    target_circuits[1].decompose(["PauliEvolution"]),
    qubits_initially_zero=True,
)
initial_params_1 = np.array(initial_params_1)
ansatz_2, initial_params_2 = generate_ansatz_from_circuit(
    target_circuits[2].decompose(["PauliEvolution"]),
    qubits_initially_zero=True,
)
initial_params_2 = np.array(initial_params_2)
print(
    f"\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, "
    f"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}"
)
print(
    f"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, "
    f"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}"
)

aqc_circuits = {}
aqc_params = {}
for k in range(1, aqc_n_steps_total + 1):
    if k <= aqc_n_steps_1:
        ansatz, base_params = ansatz_1, initial_params_1
    else:
        ansatz, base_params = ansatz_2, initial_params_2
    # Warm-start from the previous step only within the same stage
    same_stage = (k - 1 >= 1) and (
        (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)
    )
    x0 = aqc_params[k - 1] if same_stage else base_params
    obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)
    t0 = timeit.default_timer()
    result = scipy.optimize.minimize(
        obj.loss_function,
        x0,
        method="L-BFGS-B",
        jac=True,
        options=dict(maxiter=100),
    )
    elapsed = timeit.default_timer() - t0
    aqc_params[k] = result.x
    aqc_circuits[k] = ansatz.assign_parameters(result.x)
    print(
        f"  k={k}: fidelity = {1 - result.fun:.4f}, "
        f"2q depth (uniform basis) = "
        f"{uniform_2q_depth(aqc_circuits[k])}, "
        f"{elapsed:.1f}s"
    )

all_circuits = []
for k in range(1, aqc_n_steps_total + 1):
    all_circuits.append(aqc_circuits[k])
base = aqc_circuits[aqc_n_steps_total]
for k in range(1, n_steps - aqc_n_steps_total + 1):
    circuit = base.copy()
    for instr in trotter_evolution(
        circuit.qubits, interaction, anisotropy, time_step, k
    ):
        circuit.append(instr)
    all_circuits.append(circuit)

full_depths = [
    uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)
]
aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]
full_2q = [
    circuits[k - 1].decompose(["PauliEvolution"]).num_nonlocal_gates()
    for k in range(1, n_steps + 1)
]
aqc_2q = [
    circuit.decompose(["PauliEvolution"]).num_nonlocal_gates()
    for circuit in all_circuits
]
print(f"\nStep 2d — assembled {len(all_circuits)} circuits")
print(
    f"  At step {n_steps} (uniform basis): "
    f"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; "
    f"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}"
)

steps_axis = np.arange(1, n_steps + 1)
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(steps_axis, full_depths, "-o", color="black", label="Full Trotter")
ax.plot(
    steps_axis, aqc_depths, "-o", color="cadetblue", label="AQC + Trotter"
)
ax.set_xlabel("Trotter steps", fontsize=13)
ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))
ax.set_ylabel("2q gate depth (uniform basis)", fontsize=13)
ax.set_title(f"AQC circuit-depth reduction ({n_qubits} qubits)", fontsize=13)
ax.legend(fontsize=11)
plt.tight_layout()
plt.show()

# ── Step 3: Execute on IBM Quantum hardware ───────────────────────────────────
# (replaces StatevectorEstimator)
service = QiskitRuntimeService()
backend = service.least_busy(
    min_num_qubits=n_qubits,
    operational=True,
    simulator=False,
    filters=lambda x: x.configuration().processor_type["family"] == "Heron",
)
print(f"Backend: {backend.name}")

pm = generate_preset_pass_manager(optimization_level=3, backend=backend)
isa_circuits = pm.run(all_circuits, num_processes=1)
isa_2q_depths = [
    isa_circuit.depth(lambda inst: inst.operation.num_qubits == 2)
    for isa_circuit in isa_circuits
]
print(
    f"Transpiled 2q depth (deepest, ISA on {backend.name}): "
    f"{max(isa_2q_depths)} "
    f"(2q gates: {max(isa_circuit.num_nonlocal_gates() for isa_circuit in isa_circuits)})"
)

estimator = Estimator(backend)
estimator.options.environment.job_tags = ["TUT_SNS"]
estimator.options.dynamical_decoupling.enable = True
estimator.options.dynamical_decoupling.sequence_type = "XY4"
estimator.options.twirling.enable_gates = True
estimator.options.twirling.num_randomizations = 1000
estimator.options.twirling.shots_per_randomization = 128
estimator.options.resilience.measure_mitigation = True
estimator.options.resilience.measure_noise_learning.num_randomizations = 32
estimator.options.resilience.measure_noise_learning.shots_per_randomization = 100

pubs = [
    (
        isa_circuit,
        [obs.apply_layout(isa_circuit.layout) for obs in observables],
    )
    for isa_circuit in isa_circuits
]
job = estimator.run(pubs)
print(f"Job ID: {job.job_id()}")

result = job.result()

# ── Step 4: Post-process ──────────────────────────────────────────────────────
rgf_mat = np.stack([pub_result.data.evs for pub_result in result])

n_points_momentum, n_points_frequency = 100, 100
spectrum = get_dsf(
    n_qubits,
    rgf_mat,
    time_step,
    n_steps,
    n_points_momentum,
    n_points_frequency,
)

plot_rgf(
    n_qubits,
    rgf_mat,
    time_step,
    n_steps,
    title=f"Retarded Green's function — {n_qubits} qubits (QPU)",
)
plot_dsf(
    spectrum,
    time_step,
    n_points_momentum,
    n_points_frequency,
    title=rf"KCuF$_3$ DSF — {n_qubits} qubits (QPU)"
    "\n(AQC + DD + Pauli twirling + TREX)",
)

Ground-state energy (DMRG): -21.972109
Optimizing ground state ansatz...
Finished optimizing ground state ansatz in 132.2076231740648 seconds.
Ground state ansatz fidelity: 0.645956
Ground state ansatz energy: -21.616744
Built 20 circuits, deepest 2q depth (uniform basis) = 307
Defined 50 Z observables.

Step 2b — target MPS:
  k=1: max bond = 44
  k=2: max bond = 46
  k=3: max bond = 53
  k=4: max bond = 62
  k=5: max bond = 75
  k=6: max bond = 96
  k=7: max bond = 118
  k=8: max bond = 128
  k=9: max bond = 128
  k=10: max bond = 128

Step 2c — one-layer ansatz: 2971 parameters, 2q depth (uniform basis) = 39
Step 2c — two-layer ansatz: 3412 parameters, 2q depth (uniform basis) = 45
  k=1: fidelity = 1.0000, 2q depth (uniform basis) = 39, 215.6s
  k=2: fidelity = 0.9676, 2q depth (uniform basis) = 39, 269.8s
  k=3: fidelity = 0.9014, 2q depth (uniform basis) = 39, 293.5s
  k=4: fidelity = 0.8755, 2q depth (uniform basis) = 39, 279.6s
  k=5: fidelity = 0.8450, 2q depth (uniform basis) = 39, 266.4s
  k=6: fidelity = 0.8146, 2q depth (uniform basis) = 39, 371.0s

E0626 02:59:07.229070 1508496 slow_operation_alarm.cc:73] 
********************************
[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.
********************************
E0626 02:59:21.713796 1508477 slow_operation_alarm.cc:140] The operation took 2m14.484932074s

********************************
[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.
********************************

k=7: fidelity = 0.7312, 2q depth (uniform basis) = 45, 10865.3s
  k=8: fidelity = 0.7720, 2q depth (uniform basis) = 45, 1450.3s
  k=9: fidelity = 0.7316, 2q depth (uniform basis) = 45, 3388.6s

E0626 07:21:05.148724 1508496 slow_operation_alarm.cc:73] 
********************************
[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.
********************************
E0626 07:21:24.286095 1508477 slow_operation_alarm.cc:140] The operation took 2m19.137566498s

********************************
[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.
********************************

k=10: fidelity = 0.6726, 2q depth (uniform basis) = 45, 3396.1s

Step 2d — assembled 20 circuits
  At step 20 (uniform basis): Trotter 2q depth = 307, AQC+Trotter 2q depth = 171; 2q gates: Trotter = 3775, AQC+Trotter = 1913

Backend: ibm_fez
Transpiled 2q depth (deepest, ISA on ibm_fez): 105 (2q gates: 2574)
Job ID: d8v39vhropqc738biotg

Hasil hardware mereproduksi fitur-fitur kunci dari kontinum two-spinon: intensitas hamburan terkonsentrasi di dekat vektor gelombang antiferromagnetik $q = \pi$ pada energi rendah dan dibatasi dari bawah oleh dispersi spinon sinusoidal, dengan kontinum bobot spektral yang luas di atasnya alih-alih satu mode yang tajam. Ini adalah struktur yang sama yang diukur oleh hamburan neutron inelastis pada KCuF$_3$, dan ini memvalidasi seluruh alur kerja persiapan ground-state, perturbasi, evolusi Trotter yang dikompresi AQC, dan pengukuran yang dimitigasi error pada skala 50-qubit.

Langkah selanjutnya

Jika Anda merasa pekerjaan ini menarik, Anda mungkin tertarik dengan materi berikut:

Rekomendasi

• Mensimulasikan hamburan neutron dengan alur kerja Serverless dinamika AQC + Trotter — tutorial pendamping yang menjalankan eksperimen yang sama ini melalui template fungsi Qiskit Serverless yang di-deploy

• Lee et al., "Benchmarking quantum simulation with neutron-scattering experiments" (arXiv:2603.15608) — makalah referensi tempat tutorial ini didasarkan

• Teknik mitigasi dan penekanan error — DD, Pauli twirling, dan TREX yang digunakan dalam eksperimen hardware

• Kompilasi kuantum aproksimatif untuk circuit evolusi waktu — tutorial tentang AQC-Tensor

• Dokumentasi AQC-Tensor