Simulasikan serakan neutron dalam bahan kuantum dengan litar kuantum

Anggaran penggunaan: 13 minit pada pemproses Heron r2 (NOTA: Ini hanyalah anggaran. Masa jalan anda mungkin berbeza.)

Tutorial mana yang patut saya gunakan?

Gunakan tutorial ini untuk mempelajari pelaksanaannya langkah demi langkah, dengan pengkomputeran klasik berjalan secara tempatan pada komputer riba anda dan pelaksanaan perkakasan pada QPU IBM Quantum®. Contoh berskala besar memerlukan memori yang besar, dan penyusunan kuantum anggaran (AQC) boleh mengambil masa beberapa jam pada sumber komputer riba biasa. Untuk memindahkan pemampatan kepada sumber pengkomputeran dan memori Qiskit Serverless, gunakan tutorial Serverless yang berkaitan.

Hasil pembelajaran

Bagaimana spektrum serakan neutron tak elastik (INS) berkait dengan faktor struktur dinamik (DSF) bagi model spin kuantum.
Cara menyediakan keadaan asas, menerapkan gangguan setempat, dan melaksanakan evolusi masa Trotter pada litar kuantum.
Cara menggunakan penyusunan kuantum anggaran (AQC) dengan qiskit-addon-aqc-tensor untuk memampatkan litar Trotter yang mendalam bagi pelaksanaan perkakasan.
Cara mengekstrak fungsi Green tertangguh (RGF) daripada nilai jangkaan qubit dan menukarkannya kepada DSF melalui jelmaan Fourier.

Prasyarat

Latar belakang

Dalam tutorial ini, kita menghasilkan semula keputusan Lee et al., arXiv:2603.15608.

Serakan neutron membantu penyelidik mencirikan bahan yang relevan dengan kemampanan, termasuk elektrod bateri. Tutorial ini meneroka cara litar kuantum boleh mensimulasikan pengujaan magnet dan menghubungkan model teori dengan ukuran serakan. Contoh magnet-kuantumnya membangunkan kaedah untuk mengkaji gelagat bahan, menyumbang kepada set alat penyelidikan yang lebih luas untuk penemuan bahan.

Serakan neutron tak elastik dan faktor struktur dinamik

Serakan neutron tak elastik (INS) adalah salah satu penyiasat eksperimen paling berkuasa bagi pengujaan magnet dalam bahan kuantum. Apabila alur neutron termal atau sejuk melanggar hablur, neutron individu bertukar kedua-dua momentum $\mathbf{q}$ dan tenaga $\omega$ dengan subsistem magnet. Keamatan serakan yang diukur adalah berkadar dengan faktor struktur dinamik (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 mengekod korelasi ruang-masa penuh bagi darjah kebebasan spin.

KCuF $_3$ : magnet cecair-Luttinger kanonik

Kalium kuprum florida (KCuF $_3$ ) adalah antiferomagnet kuasi-satu-dimensi di mana rantai ion Cu $^{2+}$ spin- $\frac{1}{2}$ berinteraksi melalui pertukaran Heisenberg jiran-terdekat $J$ , manakala gandingan antara-rantai hanyalah $\sim 2.7\%$ daripada $J$ . Pada $T = 6\;\mathrm{K}$ di mana data INS tersedia, spektrum didominasi oleh pengujaan spinon terpecah yang menjadi ciri cecair Tomonaga-Luttinger. Kerana dinamik intra-rantai dijelaskan 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 penanda aras ideal untuk simulasi kuantum: Hamiltonian cukup ringkas untuk dilaksanakan pada pemproses kuantum, namun keadaan asas berkekusutan kuat dan spektrum pengujaan mempamerkan kontinum dwi-spinon yang luas.

Nota: tutorial ini menetapkan $J = 1$ sebagai unit tenaga dan menggunapakai penormalan $H = J\sum_i[\ldots]$ , yang sepadan dengan Hamiltonian setempat $H_\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ yang digunakan dalam pelaksanaan litar kertas kerja tersebut (Rajah S3 dalam suplemen). Hamiltonian penuh kertas kerja (Pers. 3) membawa faktor keseluruhan tambahan sebanyak 2, jadi $J$ kertas kerja adalah dua kali ganda $J$ yang digunakan di sini.

Apa yang kita simulasikan dan ukur

Kuantiti fizikal yang kita kira ialah fungsi Green tertangguh (RGF), ditakrifkan sebagai fungsi korelasi spin-spin bergantung-masa

$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 tapak rujukan (pusat rantai) dan $S_j^\alpha(t) = e^{iHt}S_j^\alpha e^{-iHt}$ adalah operator spin gambaran-Heisenberg. Dalam tutorial ini kita menumpukan pada komponen $zz$ ( $\alpha = \beta = z$ ). Pada komputer kuantum, RGF diakses dengan menyediakan keadaan asas, menerapkan gangguan setempat pada $j_c$ , mengevolusi-masa keadaan yang terganggu, dan mengukur nilai jangkaan qubit-tunggal $\langle\sigma_j^z\rangle$ pada setiap tapak $j$ bagi setiap langkah masa. Wawasan utamanya ialah setiap $\langle\sigma_j^z\rangle$ memberikan perbezaan daripada kemagnetan keadaan-asas; kerana antiferomagnet Heisenberg isotropik mempunyai kemagnetan per-tapak bersih sifar ( $\langle\sigma_j^z\rangle_\mathrm{GS} = 0$ ), nilai mentah yang diukur secara langsung menghasilkan $G^R(j, j_c, t)$ tanpa sebarang penolakan eksplisit.

Dengan mengumpulkan $G^R(j, j_c, t)$ merentasi semua tapak dan langkah masa, kita menyusun set data dua dimensi yang kemudian dijelmakan-Fourier dalam kedua-dua ruang dan masa untuk menghasilkan faktor struktur dinamik $S(q,\omega)$ . DSF adalah kuantiti yang diukur secara langsung dalam eksperimen INS: ia memberitahu kita pengujaan magnet apa yang wujud pada setiap momentum $q$ dan tenaga $\omega$ . Bagi rantai Heisenberg isotropik, spektrum pengujaan tepat ialah kontinum dwi-spinon, jalur keamatan serakan yang luas yang bentuknya berfungsi sebagai penanda aras hujung-ke-hujung yang ketat bagi simulasi kuantum: ia mengesahkan penyediaan keadaan-asas, gangguan, evolusi masa Trotter, dan protokol pengukuran sekali gus.

Aliran kerja simulasi kuantum

Aliran kerja ini mencerminkan fizik peristiwa INS. Kita (1) menyediakan keadaan asas jasad-banyak $|\psi_{\mathrm{GS}}\rangle$ pada $n$ qubit, (2) menerapkan gangguan songsang-spin setempat $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ pada pusat rantai untuk meniru pemindahan spin neutron, (3) mengevolusikan di bawah $H$ untuk langkah masa diskret menggunakan Trotterisasi peringkat-kedua, dan (4) mengukur $\langle\sigma_i^z\rangle$ pada setiap qubit pada setiap langkah untuk mendapatkan RGF. Jelmaan Fourier diskret dua dimensi kemudiannya menghasilkan DSF $S(q,\omega)$ .

Observable yang kita ukur pada setiap langkah masa ialah $\sigma_i^z$ pada setiap qubit $i$ . Dalam Qiskit ini diwakili sebagai senarai operator SparsePauliOp: satu $Z$ qubit-tunggal terbenam dalam rentetan identiti $n$ -qubit bagi setiap tapak. Observables ini dibina sekali sahaja bagi setiap kejadian masalah semasa langkah pemetaan-masalah (Langkah 1) dan digunakan semula untuk setiap litar pada skala tersebut.

Penyusunan kuantum anggaran (AQC)

Litar Trotter yang mendalam boleh dimampatkan dengan penyusunan kuantum anggaran (AQC), yang menggantikan beberapa lapisan Trotter pertama dengan ansatz berparameter yang lebih pendek dan parameternya dioptimumkan secara klasik untuk memaksimumkan fideliti peringkat-MPS dengan litar mendalam asal. Baki langkah Trotter dilampirkan secara tepat, menghasilkan litar AQC + Trotter "campuran" dengan get dua-qubit yang jauh lebih sedikit.

Simulasi MPS

Bagi sistem satu dimensi, kaedah keadaan-produk-matriks (MPS) boleh mensimulasikan dengan cekap kedua-dua penyediaan keadaan-asas (dengan kumpulan renormalisasi matriks ketumpatan, atau DMRG) dan evolusi masa peringkat-litar. Dengan mengawal dimensi ikatan $\chi$ , kita mengimbangi ketepatan dengan kos pengiraan. Dalam tutorial ini kita menggunakan simulasi MPS dalam qiskit-addon-aqc-tensor untuk mengira ansätze AQC berfideliti tinggi yang memampatkan litar Trotter mendalam bagi pelaksanaan perkakasan.

Keperluan

Sebelum memulakan tutorial ini, pastikan anda mempunyai perisian berikut dipasang:

Qiskit SDK dengan sokongan visualisasi
Qiskit Runtime (pip install qiskit-ibm-runtime)
qiskit-addon-aqc-tensor dengan tambahan quimb dan JAX (pip install 'qiskit-addon-aqc-tensor[quimb-jax]')

Persediaan

# 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 berskala kecil

Kita mula-mula menunjukkan aliran kerja penuh pada 10 qubit, mengoptimumkan ansatz keadaan-asas HVA dengan menggunakan simulasi MPS dan menggunakan simulator statevector Qiskit bagi evolusi masa. DMRG memberikan tenaga keadaan-asas rujukan dan MPS rujukan. Contoh berskala kecil ini membolehkan kita mengesahkan setiap langkah sebelum menskala naik.

Langkah 1: Memetakan input klasik kepada masalah kuantum

Kita mula dengan mentakrifkan model fizikal dan membina litar kuantum.

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

Keadaan asas. Kita menggunakan litar ansatz variasi Hamiltonian (HVA), build_ground_state_ansatz, sebagai litar penyediaan keadaan-asas. Parameter HVA dioptimumkan secara klasik dengan memaksimumkan fideliti keadaan $|\langle\psi_{\mathrm{HVA}}(\theta)|\psi_{\mathrm{DMRG}}\rangle|^2$ dengan MPS keadaan-asas DMRG, di mana keadaan HVA dinilai dengan mensimulasikan litar sebagai keadaan produk matriks menggunakan CircuitMPS quimb. Kita menggunakan scipy.optimize.minimize untuk pengoptimuman. Tenaga $\langle H \rangle$ bagi ansatz yang dioptimumkan juga dikira sebagai rujukan.

Get Trotter. Setiap sebutan interaksi jiran-terdekat $e^{-i\Delta t\, H_{\mathrm{pair}}}$ dengan $H_{\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ dibina dengan PauliEvolutionGate(H_pair, time=time_step). Qiskit mensintesis ini menjadi penguraian tiga-CNOT optimum semasa transpilasi.

Gangguan. Get $R_z(\pi/2)$ pada qubit tengah melaksanakan $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ , meniru songsang-spin setempat yang dihasilkan oleh neutron yang diserakkan.

Observables. Kita membina observable $\sigma^z$ bagi setiap tapak qubit. Objek SparsePauliOp ini diserahkan kepada primitif Estimator dalam Langkah 3 untuk mengekstrak $\langle\sigma_i^z\rangle$ pada setiap langkah masa.

# -- 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: Mengoptimumkan masalah untuk pelaksanaan perkakasan kuantum

Bagi perkakasan sebenar, litar Trotter di atas akan menjadi terlalu mendalam. Penyusunan kuantum anggaran (AQC) menangani ini dengan menggantikan $k$ lapisan Trotter pertama (termasuk litar keadaan-asas) dengan ansatz berparameter yang lebih pendek yang dioptimumkan untuk memaksimumkan fideliti peringkat-MPS dengan litar mendalam asal. Baki langkah Trotter dilampirkan secara tepat, menghasilkan litar "AQC + Trotter" yang lebih cetek.

Untuk mengimbangi keupayaan ekspresi dengan kedalaman litar, dua ansätze digunakan: ansatz satu-lapisan (dijana daripada satu langkah Trotter) memampatkan langkah masa paling awal, dan ansatz dua-lapisan yang lebih dalam (dijana daripada dua langkah Trotter) memampatkan beberapa langkah seterusnya yang memerlukan fideliti lebih tinggi.

Aliran kerja AQC mempunyai empat sub-langkah:

Bina litar sasaran — $k_1 + k_2$ litar pertama daripada Langkah 1 berfungsi secara langsung sebagai sasaran AQC.
Kira MPS sasaran — simulasikan setiap litar sasaran sebagai keadaan produk-matriks menggunakan quimb.tensor.CircuitMPS.
Jana dan optimumkan ansätze — generate_ansatz_from_circuit mencipta ansatz satu-lapisan dan dua-lapisan berparameter; parameter dioptimumkan dengan menggunakan L-BFGS-B dengan kecerunan yang dipercepatkan JAX untuk meminimumkan $1 - |\langle\psi_{\mathrm{ansatz}}|\psi_{\mathrm{target}}\rangle|^2$ . $k_1$ langkah pertama menggunakan ansatz satu-lapisan dan $k_2$ langkah seterusnya menggunakan ansatz dua-lapisan; dalam setiap peringkat, setiap langkah dimulakan dengan parameter teroptimum daripada langkah sebelumnya, dan parameter ditetapkan semula kepada lalai peringkat pada sempadan peringkat.
Susun litar campuran — bagi langkah masa melebihi titik semak AQC, lampirkan lapisan Trotter tepat kepada litar AQC dua-lapisan yang dioptimumkan 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

Output of the previous code cell

Langkah 3: Laksanakan menggunakan primitif Qiskit

Kita mensimulasikan setiap litar yang disusun-AQC menggunakan primitif StatevectorEstimator.

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

Langkah 4: Pasca-proses dan kembalikan keputusan dalam format klasik yang dikehendaki

Kita kini mengekstrak nilai jangkaan $\langle\sigma_i^z\rangle$ bagi setiap qubit $i$ pada setiap langkah masa $t$ . Nilai-nilai ini membentuk matriks fungsi Green tertangguh $G^R(j, j_c, t)$ . Kemudian, kita menjelmakan-Fourier RGF menjadi faktor struktur dinamik $S(q, \omega)$ dan memplot kedua-dua RGF dan DSF. get_dsf menerapkan simetri cermin dan memotong nilai negatif secara dalaman.

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)",
)

Output of the previous code cell

Pelaksanaan perkakasan berskala besar

Kita kini menskala naik kepada 50 qubit. Pada skala ini, pengoptimuman mencapai fideliti yang lebih rendah berbanding contoh berskala kecil: fideliti ansatz keadaan-asas menurun kepada kira-kira 0.65, dan fideliti AQC menurun kepada kira-kira 0.7 pada titik semak kemudian. Ini dijangka, dan anda boleh meningkatkan fideliti ini dengan menambah bilangan lapisan ansatz keadaan-asas (gs_n_layers) atau bilangan lelaran pengoptimuman (maxiter) dengan kos klasik tambahan. Perhatikan juga bahawa fideliti AQC dua-lapisan didapati lebih rendah daripada yang satu-lapisan. Ini bukan kemerosotan: langkah masa kemudian menjana lebih banyak kekusutan dan hanya lebih sukar untuk dimampatkan, itulah sebabnya ansatz dua-lapisan yang lebih berekspresi digunakan untuknya.

Pengoptimuman AQC juga boleh mengambil masa beberapa jam pengiraan klasik (kira-kira enam jam dalam larian yang ditunjukkan di sini, kebanyakannya dibelanjakan pada titik semak dua-lapisan dalam Langkah 2c). Untuk mengurangkan masa jam-dinding, pertimbangkan untuk menjalankan notebook ini pada perkakasan klasik yang lebih berkuasa, seperti sistem pengkomputeran berprestasi-tinggi (HPC). Sebagai alternatif, anda boleh menskala turun kepada kejadian masalah yang lebih kecil (contohnya, lebih sedikit qubit atau langkah masa); keputusan anda kemudiannya akan berbeza daripada yang ditunjukkan di sini.

Kod di bawah mengikut struktur empat-langkah yang sama seperti contoh berskala kecil. Parameter keadaan-asas HVA sekali lagi dioptimumkan dengan menggunakan simulasi MPS. Pada QPU, kita mendayakan penyahgandingan dinamik (DD), twirling Pauli, dan pemadaman ralat bacaan berpintal (TREX) untuk penindasan dan mitigasi ralat. Jadual berikut meringkaskan bagaimana eksperimen berskala besar berbeza daripada yang berskala kecil:

	Skala kecil	Skala besar
Qubit	10	50
Langkah masa	10	20
Titik semak AQC (1-lapisan + 2-lapisan)	3 + 2 = 5	6 + 4 = 10
Lapisan ansatz keadaan asas	3	5
Dimensi ikatan maksimum MPS	32	128
Estimator	`StatevectorEstimator`	QPU dengan DD, twirling Pauli, dan TREX

Mesej XLA slow-compile

Semasa pengoptimuman AQC (Langkah 2c di bawah) anda mungkin melihat mesej stderr seperti:

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

Ini adalah diagnostik tidak berbahaya daripada XLA, pengkompil yang menyokong autodiff JAX bagi qiskit-addon-aqc-tensor. Pada 50 qubit dengan dimensi ikatan MPS 128, XLA mengambil masa beberapa minit untuk mengkompil fungsi kecerunan pada kali pertama ia dijejaki. Kompilasi berjaya, dan keputusan pengoptimuman tidak terjejas.

# ── 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

Output of the previous code cell

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

Output of the previous code cell

Keputusan perkakasan menghasilkan semula ciri utama kontinum dwi-spinon: keamatan serakan tertumpu berhampiran vektor gelombang antiferomagnet $q = \pi$ pada tenaga rendah dan dibatasi dari bawah oleh dispersi spinon sinusoidal, dengan kontinum berat spektrum yang luas di atasnya berbanding satu mod tajam tunggal. Ini adalah struktur yang sama diukur oleh serakan neutron tak elastik pada KCuF $_3$ , dan ia mengesahkan aliran kerja penuh penyediaan keadaan-asas, gangguan, evolusi Trotter yang dimampat-AQC, dan pengukuran yang dimitigasi-ralat pada skala 50-qubit.

Langkah seterusnya

Jika anda mendapati kerja ini menarik, anda mungkin berminat dengan bahan berikut:

Cadangan

Simulasikan serakan neutron dengan aliran kerja Serverless dinamik AQC + Trotter — tutorial rakan yang menjalankan eksperimen yang sama ini melalui templat fungsi Qiskit Serverless yang telah dilancarkan
Lee et al., "Benchmarking quantum simulation with neutron-scattering experiments" (arXiv:2603.15608) — kertas kerja rujukan yang menjadi asas tutorial ini
Teknik mitigasi dan penindasan ralat — DD, twirling Pauli, dan TREX yang digunakan dalam eksperimen perkakasan
Penyusunan kuantum anggaran untuk litar evolusi masa — tutorial mengenai AQC-Tensor
Dokumentasi AQC-Tensor

Hasil pembelajaran​

Prasyarat​

Latar belakang​

Serakan neutron tak elastik dan faktor struktur dinamik​

KCuF3_33​: magnet cecair-Luttinger kanonik​

Apa yang kita simulasikan dan ukur​

Aliran kerja simulasi kuantum​

Penyusunan kuantum anggaran (AQC)​

Simulasi MPS​

Keperluan​

Persediaan​

Contoh simulator berskala kecil​

Langkah 1: Memetakan input klasik kepada masalah kuantum​

Langkah 2: Mengoptimumkan masalah untuk pelaksanaan perkakasan kuantum​

Langkah 3: Laksanakan menggunakan primitif Qiskit​

Langkah 4: Pasca-proses dan kembalikan keputusan dalam format klasik yang dikehendaki​

Pelaksanaan perkakasan berskala besar​

Langkah seterusnya​

Hasil pembelajaran

Prasyarat

Latar belakang

Serakan neutron tak elastik dan faktor struktur dinamik

KCuF $_3$ : magnet cecair-Luttinger kanonik

Apa yang kita simulasikan dan ukur

Aliran kerja simulasi kuantum

Penyusunan kuantum anggaran (AQC)

Simulasi MPS

Keperluan

Persediaan

Contoh simulator berskala kecil

Langkah 1: Memetakan input klasik kepada masalah kuantum

Langkah 2: Mengoptimumkan masalah untuk pelaksanaan perkakasan kuantum

Langkah 3: Laksanakan menggunakan primitif Qiskit

Langkah 4: Pasca-proses dan kembalikan keputusan dalam format klasik yang dikehendaki

Pelaksanaan perkakasan berskala besar

Langkah seterusnya