Simulation de la diffusion neutronique à l'aide d'un workflow sans serveur combinant AQC et la dynamique de Trotter

Estimation de la durée d'exécution : 18 minutes sur un processeur Heron r3 (REMARQUE : il ne s'agit que d'une estimation. (Votre temps d'exécution peut varier.)

Acquis d'apprentissage

• Comment un spectre de diffusion inélastique des neutrons permet de déterminer le facteur de structure dynamique $S(q, \omega)$ d'un aimant quantique de type « 1D ».
• Comment préparer l' KCuF$_3$ (état de Heisenberg isotrope) à l'aide du groupe de renormalisation de la matrice de densité (DMRG) et de la maximisation de la fidélité des états de produit matriciel (MPS).
• Comment exécuter l'évolution temporelle de Trotter, la compression de circuits par compilation quantique approximative (AQC) et l'exécution atténuée en un seul appel de fonction.
• Comment traiter les séries chronologiques d' $\langle \sigma_z \rangle(t)$ s par site pour obtenir des données de type « $S(q, \omega)$ » et identifier le continuum à deux spinons.

Prérequis

• Connaissance des modèles Qiskit et SparsePauliOpde l'évolution temporelle de Trotter.
• Une connaissance de base des méthodes des réseaux de tenseurs (DMRG et MPS) est utile mais pas indispensable, tout comme la maîtrise de la bibliothèque qiskit-addon-aqc-tensor utilisée par la fonction pour compresser les circuits de Trotter.

Arrière-plan

La diffusion inélastique des neutrons permet de mesurer le facteur de structure dynamique $S(q, \omega)$, qui correspond à la transformée de Fourier spatio-temporelle de la fonction de corrélation spin-spin; par conséquent, la reproduction de l’ $S(q, \omega)$ issue d’un modèle de spin microscopique constitue un test direct et falsifiable d’une simulation quantique. Ce tutoriel porte sur KCuF$_3$, une chaîne de Heisenberg antiferromagnétique à spin $\frac{1}{2}$, dont les excitations ne sont pas des inversions de spin isolées, mais des paires de spinons fractionnés : au lieu d’une dispersion de magnons bien définie, $S(q, \omega)$ présente un large continuum à deux spinons, limité en bas par $\tfrac{\pi}{2}|\sin q|$ et en haut par $\pi|\sin(q/2)|$. Il s’agit des courbes en pointillés sur les graphiques ci-dessous. Les aspects physiques dans leur intégralité, ainsi que la comparaison avec les données mesurées sur les neutrons, sont abordés dans le tutoriel original et dans l'article de Lee et al., arXiv:2603.15608.

Le déroulement des opérations quantiques reflète l'expérience de diffusion :

1. Préparez l' $|\psi_0\rangle$ de l'état fondamental de la chaîne.
2. Commencez par une perturbation locale au niveau du site central, une rotation de type « $\pi/2$ » $Z$, imitant le transfert de quantité de mouvement et d’énergie provenant du neutron.
3. Évolution temporelle sous l'hamiltonien de Heisenberg, $e^{-iHt}$, à l'aide d'une formule de produit de Trotter.
4. Mesurer l’ $\langle \sigma_z^j \rangle(t)$ de magnétisation par site. En fonction de l’ $j$ du site et du temps $t$, il s’agit exactement de la fonction de Green retardée $G^R(j, j_c, t)$; aucune conversion n’est donc nécessaire avant la transformée de Fourier à l’étape 5.
5. La transformée de Fourier de $G^R$ donne $S(q, \omega)$.

Des problèmes peuvent survenir à l'étape 3, lorsque les circuits de Trotter exacts pour les longues évolutions deviennent trop profonds pour le matériel. L'AQC avec réseaux de tenseurs résout ce problème en compressant un bloc d'étapes de Trotter en un ansatz paramétré fixe et peu profond, dont la fidélité de l'état par rapport à l'évolution exacte est maximisée de manière classique à l'aide d'un simulateur MPS ( arXiv:2301.08609 ). Le modèle AQC Dynamics regroupe l'ensemble de ce noyau quantique (synthèse de Trotter, compression AQC et exécution atténuée) en un seul appel :

Les travaux spécifiques à l'expérience sont consignés ici, dans le cahier de laboratoire : la préparation à l'état fondamental (PRE) et le post-traitement de l' $S(q, \omega)$ e (POST). Les deux étapes les plus gourmandes en ressources quantiques, la compression et l'exécution, s'exécutent à l'intérieur de la fonction.

Ce tutoriel accompagne l'article « Simuler la diffusion neutronique dans les matériaux quantiques à l'aide de circuits quantiques », qui présente en détail la réalisation de cette même expérience : le même modèle KCuF$_3$, la préparation de l'état fondamental, la projection neutronique et le post-traitement, avec la synthèse de Trotter, la compression AQC et l'exécution atténuée décrites étape par étape. Lisez ce tutoriel pour découvrir le fonctionnement de la compression AQC. Lisez cet article pour réaliser la même expérience à l'aide d'un modèle de fonction déployé : le cœur quantique se résume alors à un simple appel de fonction, et la compression AQC, qui dure plusieurs heures, s'exécute au sein du travailleur « serverless » plutôt que sur votre machine; vous n'avez donc pas besoin d'un système HPC ni d'un noyau ouvert pendant son exécution. Cet appel de fonds finance également d'autres expériences sur la dynamique des « 1D s ».

Exigences

Avant de commencer ce tutoriel, assurez-vous de disposer des éléments suivants :

• La fonction déployée sur votre compte Qiskit Serverless. Exécutez d'abord le modèle de fonction associé : déployez et exécutez le modèle de fonction « AQC + dynamique de Trotter ». Ce guide vous explique comment récupérer les fichiers sources et mettre en ligne la fonction sur votre compte. Ce tutoriel se contente d'appeler la fonction déployée.

• IBM Quantum® identifiants enregistrés pour QiskitServerless (voir le modèle de fonction). Les deux exemples de ce tutoriel appellent la fonction déployée; ils en ont donc tous les deux besoin.

• Qiskit SDK v2.0 ou version ultérieure (pip install qiskit).

• Le client « Catalog » de Qiskit IBM (pip install qiskit-ibm-catalog).

• NumPy, SciPy, et Matplotlib (pip install numpy scipy matplotlib). SciPy 1.14 ou une version ultérieure est nécessaire pour l'optimiseur COBYQA utilisé dans la préparation de l'état fondamental.

• La pile de réseaux de tenseurs AQC, car la préparation de l'état fondamental à l'étape 1 s'exécute localement dans ce cahier : pip install 'qiskit-addon-aqc-tensor[quimb-jax]==0.3.1'.

Le premier appel à une fonction nouvellement déployée est mis en attente le temps que le worker Serverless installe ses dépendances; il faut donc s'attendre à une latence supplémentaire lors de cette exécution.

Configuration

Importez les bibliothèques et définissez les fonctions auxiliaires spécifiques à l'expérience qui seront utilisées par la suite : build_gs_ansatz (l'approche variationnelle hamiltonienne, ou HVA, pour la préparation de l'état fondamental), prepare_ground_state (DMRG associé à la maximisation de la fidélité MPS), ainsi que get_spectrum, plot_green, et plot_spectrum (le post-traitement « $S(q, \omega)$ »). Ces éléments sont tirés du tutoriel original sur la diffusion des neutrons.

from functools import partial

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

import quimb.tensor as qtn
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit_addon_aqc_tensor.simulation import tensornetwork_from_circuit
from qiskit_addon_aqc_tensor.simulation.quimb import QuimbSimulator
from qiskit_ibm_catalog import QiskitServerless

#  Dynamical structure factor via discrete Fourier transform


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

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


#  Plotting helpers


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


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


#  Variational ground-state ansatz (HVA)


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


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

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


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

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

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

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

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


print("Setup complete - helpers defined.")

Output:

Setup complete - helpers defined.

Charger le modèle de fonction

Connectez-vous à Qiskit Serverless et chargez le fichier déployé aqc-dynamics-function. Les deux exemples de ce tutoriel font appel au même fn handle; la fonction n'est donc chargée qu'une seule fois, ici.

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

Exemple de simulateur à petite échelle

Nous avons d'abord exécuté l'intégralité du workflow sur une petite chaîne de 10 sites en utilisant le statevector backend exact. Cela permet de valider le pipeline PRE → FUNCTION → POST avant de mobiliser des ressources du QPU.

Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique

Construisez l'hamiltonien KCuF$_3$ sous la forme SparsePauliOp d'un hamiltonien de Heisenberg isotrope : $XX + YY + ZZ$ avec un couplage de $\tfrac14$ sur chaque liaison entre voisins immédiats; les chaînes sont des opérateurs de Pauli, donc $\tfrac14$ donne le couplage de spin $\frac{1}{2}$. Préparez l'état fondamental à l'aide de la méthode DMRG combinée à la maximisation de la fidélité MPS, puis intégrez l'effet neutronique : une rotation de type « $\pi/2$ » $Z$ au niveau du site central. C'est ce circuit préparé que nous transmettons à la fonction sous la forme initial_state. Nous conservons la valeur par observables défaut ( $Z$ par site), qui correspond exactement à la valeur « $\langle \sigma_z^j \rangle(t)$ » requise par le flux de travail sur les neutrons.

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

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

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

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

Output:

DMRG ground-state energy: -4.258035
GS fidelity: 0.9841
Prepared 10-qubit ground state with the neutron kick at site 4.

Étapes 2 et 3 : Compresser et exécuter à l'aide du modèle de fonction

Dans un processus manuel, il s'agit de deux étapes distinctes : l'optimisation des circuits pour le matériel (étape 2) et leur exécution (étape 3). Le modèle de fonction regroupe les deux en un seul appel. Il effectue la synthèse de Trotter, la compression AQC et la transpilation matérielle, puis exécute les circuits (ici sur le simulateur exact, puis ultérieurement avec un système intégré d'atténuation des erreurs sur le matériel). Les deux paramètres de réglage sont aqc_segments (le plan de compression) et aqc_options (les paramètres du MPS et de l'optimiseur). Chaque segment {"n_steps": k, "ansatz_steps": m} compresse des étapes de Trotter k consécutives en un ansatz construit à partir d'une cible de Trotter à métapes, et toutes les étapes au-delà de sont sum(n_steps) exécutées selon la méthode de Trotter classique. Les premières étapes, caractérisées par un faible niveau d'intrication, se compressent bien en un ansatz peu profond (ansatz_steps=1), c'est pourquoi nous compressons ici les trois premières étapes en un ansatz à une seule couche et les deux suivantes en un ansatz plus profond à deux couches; les cinq étapes restantes parmi les dix étapes de Trotter sont exécutées selon la méthode de Trotter classique. En effet, nous aqc_options reprenons le tutoriel d'origine : la dimension de la liaison MPS max_bond=32, cutoff=1e-8, et un optimiseur L-BFGS-B plafonné à 100 itérations.

Appelez la fonction chargée dans Setup. backend="statevector" exécute le chemin de référence exact : aucun temps de traitement sur le QPU, les circuits s'exécutant sur un simulateur de vecteurs d'état exact au sein du worker « serverless » (un compte Qiskit Serverless enregistré reste toutefois nécessaire pour l'appeler). La fonction initial_state reprend l'état fondamental préparé (y compris le « kick »); observables est omis afin que la fonction mesure l' $Z$ par site par défaut.

job = fn.run(
    t_steps=time_steps,
    aqc_segments=[
        {
            "n_steps": 3,
            "ansatz_steps": 1,
        },  # early steps -> shallow 1-layer ansatz
        {
            "n_steps": 2,
            "ansatz_steps": 2,
        },  # later steps -> deeper 2-layer ansatz
    ],
    aqc_options={
        "max_bond": mps_max_bond,  # MPS bond dimension for AQC compression
        "cutoff": mps_cutoff,
        "optimizer_settings": {
            "method": "L-BFGS-B",
            "jac": True,
            "options": {"maxiter": 100},
        },
    },
    dt=dt,
    hamiltonian=H,
    initial_state=gs_circuit,  # prepared ground state including the neutron kick
    # observables omitted -> default per-site Z (the neutron sigma_z readout)
    backend="statevector",
)

print(job.status())  # rerun this cell until status says DONE

Output:

DONE

# The per-site <sigma_z>(t) the function returns is the retarded Green's function
# G(j, j_c, t). The workflow samples t = 1..time_steps, so drop the t = 0 row (the
# prepared+kicked state before any evolution) before post-processing.
result = job.result()
print(
    "AQC fidelities:",
    {k: round(v, 4) for k, v in result["metadata"]["aqc_fidelities"].items()},
)

ev = np.array(result["expectation_values"])
Gjjc = ev[1:]  # shape (time_steps, n)
print("Green's function shape:", Gjjc.shape)

Output:

AQC fidelities: {'1': 1.0, '2': 0.9999, '3': 0.9992, '4': 0.9998, '5': 0.9995}
Green's function shape: (10, 10)

Étape 4 : Traitement ultérieur et restitution du résultat dans le format classique souhaité

Appliquer la transformée de Fourier à la fonction de Green pour obtenir $S(q, \omega)$, effectuer une symétrie par réflexion et supprimer les valeurs négatives : c'est le post-traitement standard des données neutroniques. La symétrie est exacte car $S(q, \omega) = S(-q, \omega)$ pour ce modèle, et les valeurs négatives qui subsistent sont des artefacts issus de la transformation de Fourier d'une série temporelle finie et échantillonnée de manière discrète; elles sont donc tronquées à zéro. Dans cette petite simulation exacte, le continuum à deux spinons n'est que grossièrement résolu, mais le principe est identique à celui de la simulation sur matériel qui suit.

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

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

Output:

Exemple de matériel à grande échelle

Ce même flux de travail peut être étendu sans modifier le code scientifique : une chaîne de 30 sites, une profondeur de Trotter doublée (20 étapes), un plan de compression qui fait varier la profondeur de l’ansatz (un ansatz plus profond pour les étapes ultérieures, où l’intrication est plus forte), et une exécution sur un processeur IBM Quantum avec les mécanismes intégrés de mitigation des erreurs de la fonction (découplage dynamique, « Pauli twirling » et extinction des erreurs de lecture par rotation (TREX)). Nous suivons les quatre mêmes étapes que dans l'exemple du simulateur, en réutilisant le fn handle créé lors de la configuration.

Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique

Construisez le même KCuF$_3$ Heisenberg et SparsePauliOp préparez l'état fondamental, en utilisant cette fois un gs_layers=5 ansatz plus profond pour la chaîne plus longue, puis intégrez la impulsion neutronique $\pi/2$ $Z$ au site central. C'est identique à la cartographie à petite échelle, mais à l'adresse $n = 30$.

Il faut s'attendre à une fidélité de l'état fondamental inférieure à celle obtenue pour la simulation à 10 sites : environ 0.82 ici, contre 0.98 pour la chaîne plus courte, car cinq couches HVA ne permettent pas de reproduire intégralement l'état fondamental d'une chaîne de 30 sites. C'est tout à fait normal, et ce n'est pas un problème; d'ailleurs, le tutoriel d'origine admet environ 0.65 sur 50 sites pour la même raison. L'augmentation gs_layers de la limite d'itération COBYQA permet d'améliorer les résultats, mais entraîne un surcoût en termes de calcul classique.

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

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

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

Output:

DMRG ground-state energy: -13.111355
GS fidelity: 0.8201
Prepared 30-qubit ground state with the neutron kick at site 14.

Étapes 2 et 3 : Compresser et exécuter à l'aide du modèle de fonction

Il s'agit du même appel que dans l'exemple du simulateur, mais cette fois-ci avec backend_name pointant vers un processeur IBM Quantum, de sorte que la fonction soit compilée et exécutée à cet endroit. Le plan de compression fait varier la profondeur de l'ansatz : les six premières étapes de Trotter (à faible enchevêtrement) sont compressées en un ansatz simple de faible profondeur, les quatre suivantes en un ansatz à deux couches plus profond, et les 10 étapes restantes sur les 20 sont exécutées selon la méthode de Trotter standard. aqc_options augmente la dimension de la liaison MPS à max_bond=128 pour la chaîne plus longue et plus enchevêtrée (conformément à l'original), tout en conservant le même optimiseur L-BFGS-B plafonné à 100 itérations. estimator_options Activer les mécanismes intégrés d'atténuation des erreurs : le découplage dynamique ( XY4 ), le « gate twirling » et l'atténuation des mesures TREX. Les valeurs par défaut de la fonction correspondent déjà à celles du tutoriel d'origine pour tous ces éléments, à l'exception du budget d'apprentissage TREX (measure_noise_learning). L'ensemble du bloc est tout de même affiché, car un paramètre fourni par l'appelant estimator_options remplace intégralement les valeurs par défaut de la fonction au lieu de s'y fondre; ainsi, si une clé était omise, la valeur par défaut de IBM Quantum Compute serait utilisée à la place de celle de la fonction.

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

Output:

job ID (save this to reconnect later): 43ed8d07-6d7d-4f33-b70a-7f31b765b310

# Reconnect to a previously submitted job by its ID. Only needed in a NEW kernel
# session; if you are still in the session where you submitted, the `job` handle
# from the preceding cell is already live, so skip this cell. Replace the ID that follows with your own.
job = serverless.get_job_by_id("<your job ID>")

# Check where the job is. Re-run this until it reports DONE before fetching the
# result in the following cell: QUEUED -> INITIALIZING -> RUNNING: OPTIMIZING_FOR_HARDWARE ->
# RUNNING: WAITING_FOR_QPU -> RUNNING: EXECUTING_QPU -> RUNNING: POST_PROCESSING
# -> DONE.
print(job.status())

Output:

DONE

# Run this only once the preceding status cell reports DONE. result() blocks until
# the job finishes, so calling it earlier just waits (possibly for hours).
result = job.result()
print(
    "AQC fidelities:",
    {k: round(v, 4) for k, v in result["metadata"]["aqc_fidelities"].items()},
)

ev = np.array(result["expectation_values"])
Gjjc = ev[1:]  # drop the t = 0 row -> shape (time_steps, n)

Output:

AQC fidelities: {'1': 1.0, '2': 0.9994, '3': 0.9944, '4': 0.9853, '5': 0.9747, '6': 0.959, '7': 0.9495, '8': 0.9542, '9': 0.9533, '10': 0.9451}

Étape 4 : Traitement ultérieur et restitution du résultat dans le format classique souhaité

Traitement post-calcul identique à celui effectué lors de la simulation : appliquer la transformée de Fourier à la fonction de Green pour obtenir l $S(q, \omega)$ s, effectuer une symétrie par réflexion et écrêter les valeurs négatives. Grâce à la chaîne plus longue et à l'évolution, le continuum à deux spinons est bien mieux résolu. Elle devrait remplir la bande comprise entre les limites en pointillés, avec une luminosité maximale près de $q = \pi$.

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

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

Output:

Appendice

L'exemple de matériel ci-dessus utilise une seule longueur de chaîne. Les trois spectres ci-dessous proviennent de simulations antérieures de ce même flux de travail ibm_pittsburgh sur 10, 20 et 30 sites, tous les autres paramètres d'entrée étant maintenus fixes : 20 étapes de Trotter à dt = 0.6, le plan de compression comprenant six étapes à une couche et quatre étapes à deux couches compressées par AQC, et max_bond = 128. Il s'agit de résultats enregistrés, et non du résultat des cellules précédentes.

Les mêmes paramètres sont utilisés pour les trois tailles, ce qui permet de comparer directement les spectres. Les ajuster en fonction de la longueur de la chaîne, par exemple en augmentant le nombre de couches de l'ansatz de l'état fondamental ou en choisissant une valeur plus grande max_bond, peut donner de meilleurs résultats que ceux présentés ici.

Etapes suivantes