Simuluj neutronový rozptyl v kvantových materiálech pomocí kvantových obvodů

Odhad využití: 13 minut na procesoru Heron r2 (POZNÁMKA: Toto je pouze odhad. Tvůj skutečný čas běhu se může lišit.)

Který tutoriál mám použít?

Použij tento tutoriál k seznámení s implementací krok za krokem, kdy klasický výpočet běží lokálně na tvém notebooku a hardwarové vykonání na QPU IBM Quantum®. Rozsáhlý příklad vyžaduje značnou paměť a přibližné kvantové kompilování (AQC) může na typickém laptopu trvat několik hodin. Pro přesun komprese na výpočetní a paměťové prostředky Qiskit Serverless použij doprovodný tutoriál Serverless.

Výstupy z učení

Jak spektra nepružného neutronového rozptylu (INS) souvisejí s dynamickými strukturními faktory (DSF) kvantových spinových modelů.
Jak připravit základní stav, aplikovat lokální poruchu a provést Trotterovu časovou evoluci na kvantovém obvodu.
Jak použít přibližnou kvantovou kompilaci (AQC) s balíčkem qiskit-addon-aqc-tensor ke kompresi hlubokých Trotterových obvodů pro hardwarové spuštění.
Jak extrahovat retardovanou Greenovu funkci (RGF) z očekávaných hodnot qubitů a fourierovsky ji transformovat na DSF.

Předpoklady

Kontext

V tomto tutoriálu reprodukujeme výsledky z Lee et al., arXiv:2603.15608.

Neutronový rozptyl pomáhá výzkumníkům charakterizovat materiály relevantní pro udržitelnost, včetně elektrod baterií. Tento tutoriál zkoumá, jak mohou kvantové obvody simulovat magnetické excitace a propojit teoretické modely s měřeními rozptylu. Jeho příklad kvantového magnetu rozvíjí metody pro studium chování materiálů a přispívá k širšímu výzkumnému souboru nástrojů pro objevování materiálů.

Nepružný neutronový rozptyl a dynamický strukturní faktor

Nepružný neutronový rozptyl (INS) je jednou z nejsilnějších experimentálních sond magnetických excitací v kvantových materiálech. Když svazek tepelných nebo studených neutronů dopadne na krystal, jednotlivé neutrony si s magnetickým subsystémem vyměňují jak hybnost $\mathbf{q}$ , tak energii $\omega$ . Naměřená rozptylová intenzita je úměrná dynamickému strukturnímu faktoru (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,$

který zaznamenává úplné prostoročasové korelace spinových stupňů volnosti.

KCuF $_3$ : kanonický magnet typu Luttingerovy kapaliny

Fluorid draselno-měďnatý (KCuF $_3$ ) je kvazi-jednorozměrný antiferomagnet, v němž řetězce spin- $\frac{1}{2}$ iontů Cu $^{2+}$ interagují prostřednictvím Heisenbergovy výměnné vazby $J$ mezi nejbližšími sousedy, zatímco mezřetězcová vazba činí pouze $\sim 2,7\%$ z $J$ . Při $T = 6\;\mathrm{K}$ , kde jsou dostupná data INS, spektru dominují frakcionalizované spinonové excitace charakteristické pro Tomonaga-Luttingerovu kapalinu. Protože je dynamika uvnitř řetězce dobře popsána jednorozměrným spin- $\frac{1}{2}$ XXZ Hamiltoniánem v izotropním bodě ( $\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],$

slouží KCuF $_3$ jako ideální benchmark pro kvantovou simulaci: Hamiltonián je dost jednoduchý na implementaci na kvantovém procesoru, přesto je základní stav silně provázaný a excitační spektrum vykazuje široké dvouspinonové kontinuum.

Poznámka: tento tutoriál nastavuje $J = 1$ jako jednotku energie a používá normalizaci $H = J\sum_i[\ldots]$ , což odpovídá lokálnímu Hamiltoniánu $H_\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ použitému v implementaci obvodu v článku (obr. S3 doplňkových materiálů). Úplný Hamiltonián článku (rov. 3) nese navíc celkový faktor 2, takže $J$ v článku je dvojnásobkem $J$ použitého zde.

Co simulujeme a měříme

Fyzikální veličina, kterou počítáme, je retardovaná Greenova funkce (RGF), definovaná jako časově závislá spin-spinová korelační funkce

$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,$

kde $j_c$ je referenční místo (střed řetězce) a $S_j^\alpha(t) = e^{iHt}S_j^\alpha e^{-iHt}$ je spinový operátor v Heisenbergově obrazu. V tomto tutoriálu se zaměřujeme na $zz$ složku ( $\alpha = \beta = z$ ). Na kvantovém počítači se k RGF přistupuje přípravou základního stavu, aplikací lokální poruchy v $j_c$ , časovou evolucí narušeného stavu a měřením jednoqubitové očekávané hodnoty $\langle\sigma_j^z\rangle$ na každém místě $j$ pro každý časový krok. Klíčový poznatek je, že každé $\langle\sigma_j^z\rangle$ dává rozdíl oproti magnetizaci základního stavu; protože izotropní Heisenbergův antiferomagnet má nulovou čistou magnetizaci na jedno místo ( $\langle\sigma_j^z\rangle_\mathrm{GS} = 0$ ), dává surová naměřená hodnota přímo $G^R(j, j_c, t)$ bez jakéhokoli explicitního odečítání.

Sesbíráním $G^R(j, j_c, t)$ přes všechna místa a časové kroky sestavíme dvourozměrnou datovou sadu, která je poté fourierovsky transformována v prostoru i v čase, aby vznikl dynamický strukturní faktor $S(q,\omega)$ . DSF je veličina přímo měřená v experimentu INS: říká nám, jaké magnetické excitace existují při jaké hybnosti $q$ a energii $\omega$ . Pro izotropní Heisenbergův řetězec je přesné excitační spektrum dvouspinonové kontinuum, široké pásmo rozptylové intenzity, jehož tvar slouží jako přísný end-to-end benchmark pro kvantovou simulaci: ověřuje najednou přípravu základního stavu, poruchu, Trotterovu časovou evoluci i protokol měření.

Workflow kvantové simulace

Workflow zrcadlí fyziku události INS. (1) Připravíme mnohočásticový základní stav $|\psi_{\mathrm{GS}}\rangle$ na $n$ qubitech, (2) aplikujeme lokální poruchu spinového překlopení $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ ve středu řetězce, čímž napodobíme přenos spinu neutronem, (3) necháme stav vyvíjet s hamiltoniánem $H$ v diskrétních časových krocích pomocí Trotterizace druhého řádu a (4) měříme $\langle\sigma_i^z\rangle$ na každém qubitu v každém kroku, abychom získali RGF. Dvourozměrná diskrétní Fourierova transformace poté dává DSF $S(q,\omega)$ .

Pozorovatelná, kterou měříme v každém časovém kroku, je $\sigma_i^z$ na každém qubitu $i$ . V Qiskitu je to reprezentováno jako seznam operátorů SparsePauliOp: jeden jednoqubitový $Z$ vložený do $n$ -qubitového identitního řetězce pro každé místo. Tyto pozorovatelné se sestrojí jednou pro danou instanci problému během kroku mapování problému (krok 1) a znovu se použijí pro každý obvod v tomto měřítku.

Přibližná kvantová kompilace (AQC)

Hluboké Trotterovy obvody lze komprimovat pomocí přibližné kvantové kompilace (AQC), která nahrazuje prvních několik Trotterových vrstev kratším parametrizovaným ansatzem, jehož parametry jsou klasicky optimalizovány tak, aby maximalizovaly fidelitu na úrovni MPS s původním hlubokým obvodem. Zbývající Trotterovy kroky jsou přidány přesně, čímž vznikne „smíšený“ obvod AQC + Trotter s podstatně méně dvouqubitovými hradly.

MPS simulace

Pro jednorozměrný systém dokážou metody matrix-product-state (MPS) efektivně simulovat jak přípravu základního stavu (pomocí density matrix renormalization group, neboli DMRG), tak časovou evoluci na úrovni obvodu. Řízením bond dimension $\chi$ vyvažujeme přesnost proti výpočetní náročnosti. V tomto tutoriálu používáme MPS simulaci v rámci qiskit-addon-aqc-tensor k výpočtu vysoce věrných AQC ansätze, které komprimují hluboké Trotterovy obvody pro hardwarové spuštění.

Požadavky

Než začneš s tímto tutoriálem, ujisti se, že máš nainstalováno následující:

Qiskit SDK s podporou vizualizace
Qiskit Runtime (pip install qiskit-ibm-runtime)
qiskit-addon-aqc-tensor s extras quimb a JAX (pip install 'qiskit-addon-aqc-tensor[quimb-jax]')

Nastavení

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

Příklad malého rozsahu na simulátoru

Nejprve demonstrujeme celý workflow na 10 qubitech, optimalizací HVA ansatzu základního stavu pomocí MPS simulace a použitím Qiskit statevector simulátoru pro časovou evoluci. DMRG dává referenční energii základního stavu a referenční MPS. Tento příklad malého rozsahu nám umožňuje ověřit každý krok před přechodem na větší systém.

Krok 1: Namapuj klasické vstupy na kvantový problém

Začneme definicí fyzikálního modelu a sestavením kvantových obvodů.

Hamiltonián. KCuF $_3$ je modelován 1D XXZ Hamiltoniánem v izotropním bodě ( $J = 1$ , $\epsilon = 1$ , s nastavením $J = 1$ jako jednotky energie).

Základní stav. Jako obvod pro přípravu základního stavu používáme obvod Hamiltonian variational ansatz (HVA), build_ground_state_ansatz. Parametry HVA jsou optimalizovány klasicky maximalizací fidelity stavu $|\langle\psi_{\mathrm{HVA}}(\theta)|\psi_{\mathrm{DMRG}}\rangle|^2$ s DMRG základním stavem MPS, přičemž se stav HVA vyhodnocuje simulací obvodu jako matrix product state pomocí quimb CircuitMPS. Pro optimalizaci používáme scipy.optimize.minimize. Pro referenci je také vypočtena energie $\langle H \rangle$ optimalizovaného ansatzu.

Trotterova hradla. Každý interakční člen mezi nejbližšími sousedy $e^{-i\Delta t\, H_{\mathrm{pair}}}$ s $H_{\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ je sestrojen pomocí PauliEvolutionGate(H_pair, time=time_step). Qiskit to při transpilaci syntetizuje do optimálního rozkladu se třemi CNOT hradly.

Porucha. Hradlo $R_z(\pi/2)$ na centrálním qubitu implementuje $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ , čímž napodobuje lokální překlopení spinu způsobené rozptýleným neutronem.

Pozorovatelné. Sestrojíme pozorovatelnou $\sigma^z$ pro každé místo qubitu. Tyto objekty SparsePauliOp jsou v kroku 3 předány primitivu Estimator k extrakci $\langle\sigma_i^z\rangle$ v každém časovém kroku.

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

Krok 2: Optimalizuj problém pro spuštění na kvantovém hardwaru

Pro skutečný hardware by výše uvedené Trotterovy obvody byly příliš hluboké. Přibližná kvantová kompilace (AQC) to řeší nahrazením prvních $k$ Trotterových vrstev (včetně obvodu základního stavu) kratším parametrizovaným ansatzem optimalizovaným tak, aby maximalizoval fidelitu na úrovni MPS s původním hlubokým obvodem. Zbývající Trotterovy kroky jsou přidány přesně, čímž vznikne mělčí obvod „AQC + Trotter“.

Aby se vyvážila expresivita proti hloubce obvodu, používají se dva ansätze: jednovrstvý ansatz (generovaný z jednoho Trotterova kroku) komprimuje nejranější časové kroky a hlubší dvouvrstvý ansatz (generovaný ze dvou Trotterových kroků) komprimuje několik dalších kroků, kde je potřeba vyšší fidelita.

Workflow AQC má čtyři dílčí kroky:

Sestav cílové obvody — prvních $k_1 + k_2$ obvodů z kroku 1 slouží přímo jako cíle AQC.
Vypočti cílový MPS — simuluj každý cílový obvod jako matrix-product state pomocí quimb.tensor.CircuitMPS.
Vygeneruj a optimalizuj ansätze — generate_ansatz_from_circuit vytvoří jednovrstvý a dvouvrstvý parametrizovaný ansatz; parametry jsou optimalizovány pomocí L-BFGS-B s gradienty urychlenými pomocí JAX tak, aby se minimalizovalo $1 - |\langle\psi_{\mathrm{ansatz}}|\psi_{\mathrm{target}}\rangle|^2$ . Prvních $k_1$ kroků používá jednovrstvý ansatz a dalších $k_2$ kroků používá dvouvrstvý ansatz; v rámci každé fáze každý krok vychází (warm-start) z optimalizovaných parametrů předchozího kroku a na hranici fáze se parametry resetují na výchozí hodnoty fáze.
Sestav smíšené obvody — pro časové kroky nad rámec kontrolního bodu AQC přidej přesné Trotterovy vrstvy k optimalizovanému dvouvrstvému obvodu AQC pomocí 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

Krok 3: Spusť pomocí primitiv Qiskitu

Simulujeme každý AQC-zkompilovaný obvod pomocí primitiva StatevectorEstimator.

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

Krok 4: Zpracuj výsledek a vrať ho v požadovaném klasickém formátu

Nyní extrahujeme očekávanou hodnotu $\langle\sigma_i^z\rangle$ pro každý qubit $i$ v každém časovém kroku $t$ . Tyto hodnoty tvoří matici retardované Greenovy funkce $G^R(j, j_c, t)$ . Poté fourierovsky transformujeme RGF na dynamický strukturní faktor $S(q, \omega)$ a vykreslíme jak RGF, tak DSF. get_dsf interně aplikuje zrcadlovou symetrii a ořezává záporné hodnoty.

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

Rozsáhlé hardwarové spuštění

Nyní škálujeme nahoru na 50 qubitů. V tomto měřítku dosahují optimalizace nižší fidelity než v příkladu s malým měřítkem: fidelita ansatzu základního stavu klesá zhruba na 0,65 a fidelity AQC klesají zhruba na 0,7 na pozdějších kontrolních bodech. To je očekávané a tyto fidelity lze zlepšit zvýšením počtu vrstev ansatzu základního stavu (gs_n_layers) nebo počtu optimalizačních iterací (maxiter) za cenu dalších klasických výpočtů. Všimni si také, že dvouvrstvé fidelity AQC vycházejí nižší než jednovrstvé. Nejde o regresi: pozdější časové kroky generují více provázanosti a jsou prostě obtížnější na kompresi, a proto se pro ně používá expresivnější dvouvrstvý ansatz.

Optimalizace AQC může také zabrat několik hodin klasického výpočetního času (v tomto ukázaném běhu zhruba šest hodin, přičemž většina času byla strávena na dvouvrstvých kontrolních bodech v kroku 2c). Pro zkrácení skutečné doby běhu zvaž spuštění tohoto notebooku na výkonnějším klasickém hardwaru, například na systému pro vysoce výkonné výpočty (HPC). Alternativně můžeš zmenšit instanci problému (například méně qubitů nebo méně časových kroků); tvé výsledky se pak budou lišit od těch zde ukázaných.

Kód níže sleduje stejnou čtyřkrokovou strukturu jako příklad s malým měřítkem. Parametry HVA ansatzu základního stavu jsou opět optimalizovány pomocí MPS simulace. Na QPU zapínáme dynamické odsdružení (DD), Pauliho twirling a twirled readout error extinction (TREX) pro potlačení a mitigaci chyb. Následující tabulka shrnuje, čím se experiment s velkým měřítkem liší od experimentu s malým měřítkem:

	Malé měřítko	Velké měřítko
Qubity	10	50
Časové kroky	10	20
Kontrolní body AQC (1vrstvý + 2vrstvý)	3 + 2 = 5	6 + 4 = 10
Vrstvy ansatzu základního stavu	3	5
Maximální bond dimension MPS	32	128
Estimator	`StatevectorEstimator`	QPU s DD, Pauliho twirlingem a TREX

Zprávy o pomalé kompilaci XLA

Během optimalizace AQC (krok 2c níže) se ti mohou zobrazit zprávy stderr, například:

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

Jde o neškodné diagnostické zprávy z XLA, kompilátoru, na kterém stojí JAX autodiff v qiskit-addon-aqc-tensor. Při 50 qubitech s MPS bond dimension 128 stráví XLA pár minut kompilací gradientní funkce při jejím prvním trasování. Kompilace se zdaří a výsledky optimalizace nejsou ovlivněny.

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

Hardwarové výsledky reprodukují klíčové rysy dvouspinonového kontinua: rozptylová intenzita je soustředěna poblíž antiferomagnetického vlnového vektoru $q = \pi$ při nízké energii a je zdola omezena sinusoidální spinonovou disperzí, přičemž nad ní se rozkládá široké kontinuum spektrální váhy místo jediného ostrého módu. Jde o stejnou strukturu, jaká je měřena nepružným neutronovým rozptylem na KCuF $_3$ , a validuje celý workflow přípravy základního stavu, poruchy, AQC-komprimované Trotterovy evoluce a mitigovaného měření v měřítku 50 qubitů.

Další kroky

Pokud tě tato práce zaujala, mohl by tě zajímat následující materiál:

Doporučení

Simuluj neutronový rozptyl pomocí serverless workflow s AQC + Trotter dynamikou — doprovodný tutoriál, který spouští stejný experiment prostřednictvím nasazeného function template Qiskit Serverless
Lee et al., „Benchmarking quantum simulation with neutron-scattering experiments“ (arXiv:2603.15608) — referenční článek, na kterém je tento tutoriál založen
Techniky mitigace a potlačení chyb — DD, Pauliho twirling a TREX použité v hardwarových experimentech
Přibližná kvantová kompilace pro obvody časové evoluce — tutoriál o AQC-Tensor
Dokumentace AQC-Tensor

Výstupy z učení​

Předpoklady​

Kontext​

Nepružný neutronový rozptyl a dynamický strukturní faktor​

KCuF3_33​: kanonický magnet typu Luttingerovy kapaliny​

Co simulujeme a měříme​

Workflow kvantové simulace​

Přibližná kvantová kompilace (AQC)​

MPS simulace​

Požadavky​

Nastavení​

Příklad malého rozsahu na simulátoru​

Krok 1: Namapuj klasické vstupy na kvantový problém​

Krok 2: Optimalizuj problém pro spuštění na kvantovém hardwaru​

Krok 3: Spusť pomocí primitiv Qiskitu​

Krok 4: Zpracuj výsledek a vrať ho v požadovaném klasickém formátu​

Rozsáhlé hardwarové spuštění​

Další kroky​

Výstupy z učení

Předpoklady

Kontext

Nepružný neutronový rozptyl a dynamický strukturní faktor

KCuF $_3$ : kanonický magnet typu Luttingerovy kapaliny

Co simulujeme a měříme

Workflow kvantové simulace

Přibližná kvantová kompilace (AQC)

MPS simulace

Požadavky

Nastavení

Příklad malého rozsahu na simulátoru

Krok 1: Namapuj klasické vstupy na kvantový problém

Krok 2: Optimalizuj problém pro spuštění na kvantovém hardwaru

Krok 3: Spusť pomocí primitiv Qiskitu

Krok 4: Zpracuj výsledek a vrať ho v požadovaném klasickém formátu

Rozsáhlé hardwarové spuštění

Další kroky