Simulează împrăștierea neutronilor în materiale cuantice cu circuite cuantice

Estimare de utilizare: 13 minute pe un procesor Heron r2 (NOTĂ: Aceasta este doar o estimare. Timpul tău de rulare poate varia.)

Care tutorial ar trebui să folosesc?

Folosește acest tutorial pentru a învăța implementarea pas cu pas, cu calculul clasic rulând local pe laptopul tău și execuția pe hardware pe un QPU IBM Quantum®. Exemplul la scară largă necesită memorie substanțială, iar compilarea cuantică aproximativă (AQC) poate dura câteva ore pe resurse tipice de laptop. Pentru a delega compresia către resursele de calcul și memorie Qiskit Serverless, folosește tutorialul Serverless însoțitor.

Rezultate ale învățării

• Cum se conectează spectrele de împrăștiere inelastică a neutronilor (INS) la factorii de structură dinamică (DSF) ai modelelor de spin cuantic.

• Cum să pregătești o stare fundamentală, să aplici o perturbație locală și să efectuezi evoluție temporală Trotter pe un circuit cuantic.

• Cum să folosești compilarea cuantică aproximativă (AQC) cu qiskit-addon-aqc-tensor pentru a comprima circuite Trotter adânci pentru execuția pe hardware.

• Cum să extragi funcția Green retardată (RGF) din valorile de așteptare ale qubiților și să o transformi Fourier într-un DSF.

Cerințe preliminare

• Bazele informației cuantice

• Proiectarea algoritmilor variaționali

• Introducere în primitivele Qiskit (Estimator și Sampler)

Context

În acest tutorial, reproducem rezultatele din Lee et al., arXiv:2603.15608.

Împrăștierea neutronilor ajută cercetătorii să caracterizeze materiale relevante pentru sustenabilitate, inclusiv electrozii bateriilor. Acest tutorial explorează modul în care circuitele cuantice pot simula excitații magnetice și pot conecta modelele teoretice la măsurătorile de împrăștiere. Exemplul său cu magnet cuantic dezvoltă metode pentru studierea comportamentului materialelor, contribuind la setul mai amplu de instrumente de cercetare pentru descoperirea de materiale.

Împrăștierea inelastică a neutronilor și factorul de structură dinamică

Împrăștierea inelastică a neutronilor (INS) este una dintre cele mai puternice sonde experimentale ale excitațiilor magnetice în materialele cuantice. Când un fascicul de neutroni termici sau reci lovește un cristal, neutronii individuali schimbă atât impuls $\mathbf{q}$, cât și energie $\omega$ cu subsistemul magnetic. Intensitatea de împrăștiere măsurată este proporțională cu factorul de structură dinamică (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,$

care codifică corelațiile complete spațiu-timp ale gradelor de libertate de spin.

KCuF$_3$: un magnet Luttinger-liquid canonic

Fluorura de cupru-potasiu (KCuF$_3$) este un antiferromagnet cvasi-unidimensional în care lanțuri de ioni Cu$^{2+}$ de spin-$\frac{1}{2}$ interacționează printr-un schimb Heisenberg $J$ cu vecinii cei mai apropiați, în timp ce cuplajul interlanț este de doar $\sim 2.7\%$ din $J$. La $T = 6\;\mathrm{K}$, unde sunt disponibile date INS, spectrul este dominat de excitații de spinon fracționalizate, caracteristice unui lichid Tomonaga-Luttinger. Deoarece dinamica intralanț este bine descrisă de hamiltonianul XXZ unidimensional de spin-$\frac{1}{2}$ la punctul izotrop ($\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$ servește ca un benchmark ideal pentru simularea cuantică: hamiltonianul este suficient de simplu pentru a fi implementat pe un procesor cuantic, dar starea fundamentală este puternic încâlcită (entangled), iar spectrul de excitații prezintă un continuum larg de doi spinoni.

Notă: acest tutorial stabilește $J = 1$ ca unitate de energie și adoptă normalizarea $H = J\sum_i[\ldots]$, care corespunde hamiltonianului local $H_\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ folosit în implementarea circuitului din lucrare (Fig. S3 din materialul suplimentar). Hamiltonianul complet din lucrare (Ec. 3) conține un factor global suplimentar de 2, astfel încât $J$ din lucrare este dublul lui $J$ folosit aici.

Ce simulăm și măsurăm

Mărimea fizică pe care o calculăm este funcția Green retardată (RGF), definită ca funcția de corelație spin-spin dependentă de timp

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

unde $j_c$ este un sit de referință (centrul lanțului), iar $S_j^\alpha(t) = e^{iHt}S_j^\alpha e^{-iHt}$ este operatorul de spin în imaginea Heisenberg. În acest tutorial ne concentrăm pe componenta $zz$ ($\alpha = \beta = z$). Pe un calculator cuantic, RGF este accesată prin pregătirea stării fundamentale, aplicarea unei perturbații locale la $j_c$, evoluția în timp a stării perturbate și măsurarea valorii de așteptare a unui singur qubit $\langle\sigma_j^z\rangle$ la fiecare sit $j$ pentru fiecare pas de timp. Ideea cheie este că fiecare $\langle\sigma_j^z\rangle$ dă diferența față de magnetizarea stării fundamentale; deoarece antiferromagnetul Heisenberg izotrop are magnetizare netă pe sit zero ($\langle\sigma_j^z\rangle_\mathrm{GS} = 0$), valoarea măsurată brută dă direct $G^R(j, j_c, t)$ fără nicio scădere explicită.

Prin colectarea $G^R(j, j_c, t)$ pentru toate siturile și pașii de timp, asamblăm un set de date bidimensional care este apoi transformat Fourier atât în spațiu, cât și în timp pentru a produce factorul de structură dinamică $S(q,\omega)$. DSF este mărimea măsurată direct într-un experiment INS: ne spune ce excitații magnetice există la fiecare impuls $q$ și energie $\omega$. Pentru lanțul Heisenberg izotrop, spectrul exact de excitații este un continuum de doi spinoni, o bandă largă de intensitate de împrăștiere a cărei formă servește ca un benchmark riguros de la un capăt la altul pentru simularea cuantică: validează pregătirea stării fundamentale, perturbația, evoluția temporală Trotter și protocolul de măsurare, toate deodată.

Fluxul de lucru al simulării cuantice

Fluxul de lucru reflectă fizica unui eveniment INS. (1) Pregătim starea fundamentală cu mai multe corpuri $|\psi_{\mathrm{GS}}\rangle$ pe $n$ qubiți, (2) aplicăm o perturbație locală de răsturnare a spinului $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ în centrul lanțului pentru a imita transferul de spin al neutronului, (3) evoluăm sub $H$ pentru pași de timp discreți folosind trotterizarea de ordinul doi și (4) măsurăm $\langle\sigma_i^z\rangle$ pe fiecare qubit la fiecare pas pentru a obține RGF. O transformată Fourier discretă bidimensională produce apoi DSF $S(q,\omega)$.

Observabila pe care o măsurăm la fiecare pas de timp este $\sigma_i^z$ pe fiecare qubit $i$. În Qiskit, aceasta este reprezentată ca o listă de operatori SparsePauliOp: un singur $Z$ pe un qubit, încorporat în șirul de identitate de $n$ qubiți pentru fiecare sit. Aceste observabile sunt construite o singură dată pentru fiecare instanță de problemă, în timpul pasului de mapare a problemei (Pasul 1), și reutilizate pentru fiecare circuit la acea scară.

Compilarea cuantică aproximativă (AQC)

Circuitele Trotter adânci pot fi comprimate cu compilarea cuantică aproximativă (AQC), care înlocuiește primele câteva straturi Trotter cu un ansatz parametrizat mai scurt, ai cărui parametri sunt optimizați clasic pentru a maximiza fidelitatea la nivel MPS cu circuitul original adânc. Pașii Trotter rămași sunt adăugați exact, producând un circuit "mixt" AQC + Trotter cu semnificativ mai puține porți pe doi qubiți.

Simulare MPS

Pentru un sistem unidimensional, metodele de tip matrix-product-state (MPS) pot simula eficient atât pregătirea stării fundamentale (cu grupul de renormalizare a matricei de densitate, sau DMRG), cât și evoluția temporală la nivel de circuit. Prin controlul dimensiunii de legătură $\chi$, echilibrăm precizia cu costul computațional. În acest tutorial folosim simularea MPS în cadrul qiskit-addon-aqc-tensor pentru a calcula ansatz-uri AQC de înaltă fidelitate care comprimă circuitele Trotter adânci pentru execuția pe hardware.

Cerințe

Înainte de a începe acest tutorial, asigură-te că ai instalate următoarele:

• Qiskit SDK cu suport pentru vizualizare

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

• qiskit-addon-aqc-tensor cu extrasele quimb și JAX (pip install 'qiskit-addon-aqc-tensor[quimb-jax]')

Configurare

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

Exemplu de simulator la scară mică

Demonstrăm mai întâi întregul flux de lucru pe 10 qubiți, optimizând ansatz-ul HVA pentru starea fundamentală prin folosirea simulării MPS și utilizând simulatorul Qiskit de tip vector de stare (statevector) pentru evoluția temporală. DMRG oferă o energie de referință a stării fundamentale și un MPS de referință. Acest exemplu la scară mică ne permite să validăm fiecare pas înainte de a extinde scara.

Pasul 1: Maparea intrărilor clasice la o problemă cuantică

Începem prin definirea modelului fizic și construirea circuitelor cuantice.

Hamiltonian. KCuF$_3$ este modelat de hamiltonianul XXZ 1D la punctul izotrop ($J = 1$, $\epsilon = 1$, stabilind $J = 1$ ca unitate de energie).

Starea fundamentală. Folosim un circuit de tip ansatz variațional hamiltonian (HVA), build_ground_state_ansatz, ca circuit de pregătire a stării fundamentale. Parametrii HVA sunt optimizați clasic prin maximizarea fidelității stării $|\langle\psi_{\mathrm{HVA}}(\theta)|\psi_{\mathrm{DMRG}}\rangle|^2$ cu MPS-ul de stare fundamentală DMRG, unde starea HVA este evaluată prin simularea circuitului ca o stare de tip produs de matrice (MPS) folosind CircuitMPS din quimb. Folosim scipy.optimize.minimize pentru optimizare. Energia $\langle H \rangle$ a ansatz-ului optimizat este de asemenea calculată, ca referință.

Porți Trotter. Fiecare termen de interacțiune între vecinii cei mai apropiați $e^{-i\Delta t\, H_{\mathrm{pair}}}$ cu $H_{\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ este construit cu PauliEvolutionGate(H_pair, time=time_step). Qiskit sintetizează aceasta în descompunerea optimă cu trei porți CNOT în timpul transpilării.

Perturbație. O poartă $R_z(\pi/2)$ pe qubitul central implementează $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$, imitând răsturnarea locală a spinului produsă de un neutron împrăștiat.

Observabile. Construim o observabilă $\sigma^z$ pentru fiecare sit de qubit. Aceste obiecte SparsePauliOp sunt transmise primitivei Estimator în Pasul 3 pentru a extrage $\langle\sigma_i^z\rangle$ la fiecare pas de timp.

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

Pasul 2: Optimizarea problemei pentru execuția pe hardware cuantic

Pentru hardware real, circuitele Trotter de mai sus ar fi prea adânci. Compilarea cuantică aproximativă (AQC) rezolvă această problemă prin înlocuirea primelor $k$ straturi Trotter (inclusiv circuitul stării fundamentale) cu un ansatz parametrizat mai scurt, optimizat pentru a maximiza fidelitatea la nivel MPS cu circuitul original adânc. Pașii Trotter rămași sunt adăugați exact, producând un circuit "AQC + Trotter" mai puțin adânc.

Pentru a echilibra expresivitatea cu adâncimea circuitului, sunt folosite două ansatz-uri: un ansatz cu un strat (generat dintr-un singur pas Trotter) comprimă primii pași de timp, iar un ansatz cu două straturi mai adânc (generat din doi pași Trotter) comprimă următorii câțiva pași unde este necesară o fidelitate mai mare.

Fluxul de lucru AQC are patru sub-pași:

1. Construirea circuitelor țintă — primele $k_1 + k_2$ circuite din Pasul 1 servesc direct drept ținte AQC.

2. Calcularea MPS-ului țintă — simulează fiecare circuit țintă ca o stare de tip produs de matrice folosind quimb.tensor.CircuitMPS.

3. Generarea și optimizarea ansatz-urilor — generate_ansatz_from_circuit creează un ansatz cu un strat și un ansatz cu două straturi, parametrizate; parametrii sunt optimizați folosind L-BFGS-B cu gradienți accelerați de JAX pentru a minimiza $1 - |\langle\psi_{\mathrm{ansatz}}|\psi_{\mathrm{target}}\rangle|^2$. Primii $k_1$ pași folosesc ansatz-ul cu un strat, iar următorii $k_2$ pași folosesc ansatz-ul cu două straturi; în cadrul fiecărei etape, fiecare pas pornește cu parametrii optimizați ai pasului anterior, iar parametrii sunt resetați la valorile implicite ale etapei la granița dintre etape.

4. Asamblarea circuitelor mixte — pentru pașii de timp de dincolo de punctul de control AQC, adaugă straturi Trotter exacte la circuitul AQC optimizat cu două straturi, folosind 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

Pasul 3: Executarea folosind primitivele Qiskit

Simulăm fiecare circuit compilat AQC folosind primitiva StatevectorEstimator.

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

Pasul 4: Post-procesarea și returnarea rezultatului în formatul clasic dorit

Extragem acum valoarea de așteptare $\langle\sigma_i^z\rangle$ pentru fiecare qubit $i$ la fiecare pas de timp $t$. Aceste valori formează matricea funcției Green retardate $G^R(j, j_c, t)$. Apoi, transformăm Fourier RGF în factorul de structură dinamică $S(q, \omega)$ și reprezentăm grafic atât RGF, cât și DSF. get_dsf aplică intern simetria de oglindire și taie valorile negative.

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

Execuție hardware la scară mare

Extindem acum scara la 50 de qubiți. La această scară, optimizările ating fidelități mai mici decât în exemplul la scară mică: fidelitatea ansatz-ului stării fundamentale scade la aproximativ 0.65, iar fidelitățile AQC scad la aproximativ 0.7 la punctele de control ulterioare. Acest lucru este de așteptat și poți îmbunătăți aceste fidelități prin creșterea numărului de straturi ale ansatz-ului stării fundamentale (gs_n_layers) sau a numărului de iterații de optimizare (maxiter), cu un cost clasic suplimentar. Reține de asemenea că fidelitățile AQC cu două straturi ies mai mici decât cele cu un strat. Aceasta nu este o regresie: pașii de timp ulteriori generează mai multă încâlcire (entanglement) și sunt pur și simplu mai greu de comprimat, motiv pentru care ansatz-ul mai expresiv cu două straturi este folosit pentru ei.

Optimizarea AQC poate, de asemenea, să dureze mai multe ore de timp de calcul clasic (aproximativ șase ore în rularea prezentată aici, cea mai mare parte petrecută la punctele de control cu două straturi din Pasul 2c). Pentru a reduce timpul real de execuție, ia în considerare rularea acestui notebook pe hardware clasic mai puternic, cum ar fi un sistem de calcul de înaltă performanță (HPC). Alternativ, poți reduce scara la o instanță de problemă mai mică (de exemplu, mai puțini qubiți sau mai puțini pași de timp); rezultatele tale vor diferi atunci de cele prezentate aici.

Codul de mai jos urmează aceeași structură în patru pași ca exemplul la scară mică. Parametrii stării fundamentale HVA sunt din nou optimizați folosind simularea MPS. Pe QPU, activăm dynamical decoupling (DD), Pauli twirling și twirled readout error extinction (TREX) pentru suprimarea și atenuarea erorilor. Tabelul următor rezumă modul în care experimentul la scară mare diferă de cel la scară mică:

	Scară mică	Scară mare
Qubiți	10	50
Pași de timp	10	20
Puncte de control AQC (1 strat + 2 straturi)	3 + 2 = 5	6 + 4 = 10
Straturi ale ansatz-ului stării fundamentale	3	5
Dimensiunea maximă de legătură MPS	32	128
Estimator	StatevectorEstimator	QPU cu DD, Pauli twirling și TREX

Mesaje de compilare lentă XLA

În timpul optimizării AQC (Pasul 2c de mai jos) este posibil să vezi mesaje stderr precum:

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

Acestea sunt diagnostice benigne de la XLA, compilatorul care stă la baza autodiff-ului JAX al qiskit-addon-aqc-tensor. La 50 de qubiți cu dimensiunea de legătură MPS de 128, XLA petrece câteva minute compilând funcția de gradient prima dată când este trasată. Compilarea reușește, iar rezultatele optimizării nu sunt afectate.

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

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

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

gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

result = job.result()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Rezultatele hardware reproduc trăsăturile cheie ale continuumului de doi spinoni: intensitatea de împrăștiere este concentrată aproape de vectorul de undă antiferromagnetic $q = \pi$ la energie joasă și este mărginită inferior de dispersia sinusoidală a spinonilor, cu un continuum larg de greutate spectrală deasupra ei, mai degrabă decât un singur mod ascuțit. Aceasta este aceeași structură măsurată prin împrăștiere inelastică a neutronilor pe KCuF$_3$ și validează întregul flux de lucru de pregătire a stării fundamentale, perturbație, evoluție Trotter comprimată AQC și măsurare cu atenuare a erorilor la scara de 50 de qubiți.

Pași următori

Dacă ai găsit această lucrare interesantă, s-ar putea să fii interesat de următoarele materiale:

Recomandări

• Simulează împrăștierea neutronilor cu un flux de lucru Serverless de dinamică AQC + Trotter — un tutorial complementar care rulează același experiment printr-un șablon de funcție Qiskit Serverless implementat

• Lee et al., "Benchmarking quantum simulation with neutron-scattering experiments" (arXiv:2603.15608) — lucrarea de referință pe care se bazează acest tutorial

• Tehnici de atenuare și suprimare a erorilor — DD, Pauli twirling și TREX folosite în experimentele hardware

• Compilarea cuantică aproximativă pentru circuite de evoluție temporală — tutorial despre AQC-Tensor

• Documentația AQC-Tensor