{ "cells": [ { "cell_type": "markdown", "id": "8e82ead7", "metadata": {}, "source": [ "---\n", "title: \"Simuler la diffusion des neutrons dans les matériaux quantiques à l'aide de circuits quantiques\"\n", "description: \"Calculer le facteur de structure dynamique d'un aimant quantique à l'aide de circuits de Trotter et d'une simulation MPS.\"\n", "---\n", "\n", "{/* cspell:ignore DSFs spinon DMRG viridis fontsize vmax vmin cadetblue antiferromagnet COBYQA */}\n", "\n", "\n", "\n", "# Simuler la diffusion des neutrons dans les matériaux quantiques à l'aide de circuits quantiques\n", "\n", "*Estimation de la durée d'exécution : 13 minutes sur un processeur Heron r2 (REMARQUE : il s'agit uniquement d'une estimation. (Votre temps d'exécution peut varier.)*\n", "\n" ] }, { "cell_type": "markdown", "id": "11033666", "metadata": {}, "source": [ "\n", "\n", "## Acquis d'apprentissage\n", "\n", "À l'issue de ce tutoriel, vous devriez être en mesure de comprendre les éléments suivants :\n", "\n", "* Comment les spectres de diffusion inélastique des neutrons (INS) sont liés aux facteurs de structure dynamiques (DSF) des modèles de spin quantiques.\n", "* Comment préparer un état fondamental, appliquer une perturbation locale et effectuer une évolution temporelle de Trotter sur un circuit quantique.\n", "* Comment utiliser la compilation quantique approximative (AQC) avec `qiskit-addon-aqc-tensor` pour compresser des circuits de Trotter profonds en vue de leur exécution sur matériel.\n", "* Comment extraire la fonction de Green retardée (RGF) à partir des valeurs attendues des qubits et la transformer par transformée de Fourier en une fonction de Green de Schrödinger (DSF).\n", "\n", "\n", "\n", "## Prérequis\n", "\n", "Nous vous recommandons de vous familiariser avec les sujets suivants :\n", "\n", "* [Notions de base sur l’information quantique](/learning/courses/basics-of-quantum-information)\n", "* [Conception d'algorithmes variationnels](/learning/courses/variational-algorithm-design)\n", "* [Introduction à « Qiskit primitives » (Estimateur et échantillonneur)](/docs/guides/qiskit-runtime-primitives)\n", "\n", "\n", "\n", "## Arrière-plan\n", "\n", "Dans ce tutoriel, nous reproduisons les résultats de [Lee et al., arXiv:2603.15608](https://arxiv.org/abs/2603.15608).\n", "\n", "\n", "\n", "### Diffusion inélastique des neutrons et facteur de structure dynamique\n", "\n", "La diffusion inélastique des neutrons (INS) est l'un des outils expérimentaux les plus puissants pour l'étude des excitations magnétiques dans les matériaux quantiques. Lorsqu'un faisceau de neutrons thermiques ou froids frappe un cristal, chaque neutron échange à la fois de la quantité de mouvement $\\mathbf{q}$ et de l'énergie $\\omega$ avec le sous-système magnétique. L'intensité de diffusion mesurée est proportionnelle au facteur de structure dynamique (DSF),\n", "\n", "$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,$\n", "\n", "qui code l'ensemble des corrélations spatio-temporelles des degrés de liberté de spin.\n", "\n", "\n", "\n", "### KCuF$_3$ : un aimant canonique de type liquide de Luttinger\n", "\n", "Le fluorure de cuivre et de potassium ( KCuF$_3$) est un antiferromagnétique quasi-unidimensionnel dans lequel des chaînes d’ions spin- $\\frac{1}{2}$ Cu $^{2+}$ interagissent via une $J$ d’échange de Heisenberg entre voisins immédiats, tandis que le couplage interchaînes n’est que $\\sim 2.7\\%$ e de $J$. À $T = 6\\;\\mathrm{K}$, où des données INS sont disponibles, le spectre est dominé par des excitations de spinons fractionnées caractéristiques d’un liquide de Tomonaga-Luttinger. Étant donné que la dynamique intra-chaîne est bien décrite par l'hamiltonien XXZ de type spin- $\\frac{1}{2}$ e unidimensionnel au point isotrope ( $\\epsilon = 1$ ),\n", "\n", "$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],$\n", "\n", "KCuF$_3$ constitue un modèle de référence idéal pour la simulation quantique : son hamiltonien est suffisamment simple pour être mis en œuvre sur un processeur quantique, tandis que son état fondamental présente un fort enchevêtrement et que son spectre d'excitation comporte un large continuum à deux spinons.\n", "\n", "Remarque : ce tutoriel utilise l’ $J = 1$ e comme unité d’énergie et adopte la normalisation $H = J\\sum_i[\\ldots]$, qui correspond à l’hamiltonien local $H_\\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ utilisé dans la mise en œuvre du circuit présentée dans l’article (fig. S3 (de l'annexe). L'hamiltonien complet de l'article (équation 3) comporte un facteur global supplémentaire de 2; l' $J$ e de l'article est donc le double de l' $J$ e utilisée ici.\n", "\n", "\n", "\n", "### Ce que nous simulons et mesurons\n", "\n", "La grandeur physique que nous calculons est la **fonction de Green retardée** (RGF), définie comme la fonction de corrélation spin-spin dépendante du temps\n", "\n", "$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,$\n", "\n", "où $j_c$ est un site de référence (le centre de la chaîne) et $S_j^\\alpha(t) = e^{iHt}S_j^\\alpha e^{-iHt}$ est l'opérateur de spin selon la représentation de Heisenberg. Dans ce tutoriel, nous nous intéressons au composant « $zz$ » ( $\\alpha = \\beta = z$ ). Sur un ordinateur quantique, on accède au RGF en préparant l’état fondamental, en appliquant une perturbation locale à l’adresse $j_c$, en faisant évoluer l’état perturbé dans le temps, puis en mesurant la valeur attendue d’un qubit unique $\\langle\\sigma_j^z\\rangle$ à chaque site $j$ pour chaque pas de temps. L'idée essentielle est que chaque mesure de l' $\\langle\\sigma_j^z\\rangle$ e correspond à la différence par rapport à la de l'état fondamental; comme l'antiferromagnétique isotropique de Heisenberg présente une nette nulle par site ( $\\langle\\sigma_j^z\\rangle_\\mathrm{GS} = 0$ ), la valeur brute mesurée donne directement l' $G^R(j, j_c, t)$ sans qu'il soit nécessaire de procéder à une soustraction explicite.\n", "\n", "En collectant l’ $G^R(j, j_c, t)$ s sur l’ensemble des sites et des pas de temps, nous constituons un ensemble de données bidimensionnel qui est ensuite soumis à une transformée de Fourier à la fois dans l’espace et dans le temps afin de produire le **facteur de structure dynamique** $S(q,\\omega)$. Le DSF est la grandeur directement mesurée lors d’une expérience INS : il nous indique quelles excitations magnétiques existent à chaque $q$ d’impulsion et à chaque $\\omega$ d’énergie. Pour la chaîne de Heisenberg isotrope, le spectre d’excitation exact est un **continuum à deux spinons**, une large bande d’intensité de diffusion dont la forme sert de référence rigoureuse de bout en bout pour la simulation quantique : elle valide à la fois la préparation de l’état fondamental, la perturbation, l’évolution temporelle de Trotter et le protocole de mesure.\n", "\n", "\n", "\n", "### Processus de simulation quantique\n", "\n", "Le déroulement des opérations reflète les mécanismes physiques d'un événement INS. Nous (1) préparons l’ $|\\psi_{\\mathrm{GS}}\\rangle$ de l’état fondamental à N corps sur $n$ qubits, (2) appliquons une perturbation locale de renversement de spin $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ au centre de la chaîne afin de reproduire le transfert de spin du neutron, (3) faisons évoluer le système selon l’équation d’ $H$ e par pas de temps discrets à l’aide d’une trotterisation du second ordre, et (4) mesurons $\\langle\\sigma_i^z\\rangle$ sur chaque qubit à chaque pas afin d’obtenir la fonction de transfert de spin (RGF). Une transformée de Fourier discrète bidimensionnelle permet alors d'obtenir l' $S(q,\\omega)$ DSF.\n", "\n", "La grandeur observable que nous mesurons à chaque pas de temps est l' $\\sigma_i^z$ sur chaque qubit $i$. Dans Qiskit, cela est représenté par une liste d'opérateurs `SparsePauliOp` : un opérateur à un seul qubit $Z$ intégré dans la chaîne d'identité à $n$ qubits pour chaque site. Ces observables sont construits une seule fois par instance du problème lors de l'étape de cartographie du problème (étape 1) et réutilisés pour chaque circuit à cette échelle.\n", "\n", "\n", "\n", "### Compilation quantique approximative (AQC)\n", "\n", "Les circuits Deep Trotter peuvent être compressés à l'aide de **la compilation quantique approximative (AQC)**, qui remplace les premières couches de Trotter par un ansatz paramétré plus court, dont les paramètres sont optimisés de manière classique afin de maximiser la fidélité au niveau MPS par rapport au circuit profond d'origine. Les étapes restantes de l'algorithme de Trotter sont ajoutées à l'identique, ce qui donne un circuit « mixte » AQC + Trotter comportant nettement moins de portes à deux qubits.\n", "\n", "\n", "\n", "### Simulation MPS\n", "\n", "Pour un système unidimensionnel, les méthodes MPS (Matrix-Product-State) permettent de simuler efficacement à la fois la préparation de l'état fondamental (à l'aide du groupe de renormalisation de la matrice de densité, ou DMRG) et l'évolution temporelle au niveau du circuit. En ajustant la dimension de la liaison $\\chi$, on trouve un compromis entre la précision et le coût de calcul. Dans ce tutoriel, nous utilisons la simulation MPS pour `qiskit-addon-aqc-tensor` calculer des ansätze AQC haute fidélité qui permettent de compresser les circuits de Trotter profonds en vue de leur exécution matérielle.\n", "\n" ] }, { "cell_type": "markdown", "id": "de1ad5a2", "metadata": {}, "source": [ "\n", "\n", "## Exigences\n", "\n", "Avant de commencer ce tutoriel, assurez-vous d'avoir installé les éléments suivants :\n", "\n", "* Qiskit SDK avec prise en charge de [la visualisation](/docs/api/qiskit/visualization)\n", "* Qiskit Runtime (`pip install qiskit-ibm-runtime`)\n", "* `qiskit-addon-aqc-tensor` avec `quimb` et les fonctionnalités supplémentaires de JAX (`pip install 'qiskit-addon-aqc-tensor[quimb-jax]'`)\n", "\n" ] }, { "cell_type": "markdown", "id": "9938e4bd", "metadata": {}, "source": [ "\n", "\n", "## Configuration\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "71835ff6", "metadata": {}, "outputs": [], "source": [ "import timeit\n", "import warnings\n", "from collections.abc import Iterator, Sequence\n", "from functools import partial\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import quimb.tensor as qtn\n", "import scipy.optimize\n", "from numpy.typing import NDArray\n", "from qiskit import QuantumCircuit\n", "from qiskit.circuit import CircuitInstruction, ParameterVector, Qubit\n", "from qiskit.circuit.library import PauliEvolutionGate\n", "from qiskit.primitives import StatevectorEstimator\n", "from qiskit.quantum_info import SparsePauliOp\n", "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", "from qiskit_addon_aqc_tensor import generate_ansatz_from_circuit\n", "from qiskit_addon_aqc_tensor.objective import MaximizeStateFidelity\n", "from qiskit_addon_aqc_tensor.simulation import tensornetwork_from_circuit\n", "from qiskit_addon_aqc_tensor.simulation.quimb import QuimbSimulator\n", "from qiskit_ibm_runtime import EstimatorV2 as Estimator\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "from qiskit_quimb import quimb_circuit\n", "from scipy.sparse import SparseEfficiencyWarning\n", "\n", "# scipy.sparse.linalg.expm (reached via PauliEvolutionGate.to_matrix) prefers\n", "# CSC and warns when handed the CSR matrix from SparsePauliOp.to_matrix; the\n", "# result is unaffected, so silence the cosmetic warning.\n", "warnings.filterwarnings(\"ignore\", category=SparseEfficiencyWarning)\n", "\n", "\n", "def xxz_hamiltonian_mpo(\n", " n_qubits: int, interaction: float = 1.0, anisotropy: float = 1.0\n", ") -> qtn.MatrixProductOperator:\n", " \"\"\"1D XXZ Hamiltonian as a quimb MPO.\n", "\n", " Builds the Hamiltonian using ``qtn.SpinHam1D``.\n", "\n", " Args:\n", " n_qubits: Number of sites.\n", " interaction: Overall interaction strength (J in the paper).\n", " anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the\n", " isotropic Heisenberg point; ``ε = 0`` is the Ising limit.\n", "\n", " Returns:\n", " The Hamiltonian as a matrix product operator.\n", " \"\"\"\n", " builder = qtn.SpinHam1D(S=1 / 2)\n", " builder += interaction * anisotropy * 0.5, \"+\", \"-\"\n", " builder += interaction * anisotropy * 0.5, \"-\", \"+\"\n", " builder += interaction, \"Z\", \"Z\"\n", " return builder.build_mpo(L=n_qubits)\n", "\n", "\n", "def build_ground_state_ansatz(n_qubits: int, n_layers: int) -> QuantumCircuit:\n", " \"\"\"Hamiltonian variational ansatz (HVA) circuit for ground-state preparation.\n", "\n", " Starts from a product of singlet pairs and applies alternating\n", " odd/even layers of parameterized XXZ pair-evolution gates.\n", "\n", " The returned circuit is parameterized: it carries a ``ParameterVector``\n", " named ``\"theta\"`` of length ``2 * n_layers`` whose values must be\n", " assigned (e.g. via ``circuit.assign_parameters``) before simulation.\n", " Element ``2 * r`` is the odd-layer angle and element ``2 * r + 1`` is the\n", " even-layer angle of layer ``r``.\n", "\n", " Args:\n", " n_qubits: Number of qubits (must be even).\n", " n_layers: Number of HVA layers.\n", "\n", " Returns:\n", " The parameterized HVA preparation circuit.\n", " \"\"\"\n", " theta = ParameterVector(\"theta\", 2 * n_layers)\n", " circuit = QuantumCircuit(n_qubits)\n", " # Initial singlet product state\n", " for i in range(n_qubits // 2):\n", " circuit.x(2 * i)\n", " circuit.x(2 * i + 1)\n", " circuit.h(2 * i + 1)\n", " circuit.cx(2 * i + 1, 2 * i)\n", " # Variational layers: each pair gate is exp(-i theta (XX+YY+ZZ)/2)\n", " pair_ham = SparsePauliOp(\n", " [\"XX\", \"YY\", \"ZZ\"], coeffs=[0.5, 0.5, 0.5]\n", " ) # H_pair (HVA form)\n", " for r in range(n_layers):\n", " for i in range(1, (n_qubits + 1) // 2): # odd layer\n", " circuit.append(\n", " PauliEvolutionGate(pair_ham, time=theta[2 * r]),\n", " [2 * i - 1, 2 * i],\n", " )\n", " for i in range(n_qubits // 2): # even layer\n", " circuit.append(\n", " PauliEvolutionGate(pair_ham, time=theta[2 * r + 1]),\n", " [2 * i, 2 * i + 1],\n", " )\n", " return circuit\n", "\n", "\n", "def optimize_ground_state_ansatz(\n", " ansatz: QuantumCircuit,\n", " x0: NDArray[np.floating],\n", " target_mps: qtn.MatrixProductState,\n", " *,\n", " max_bond: int | None = None,\n", " cutoff: float = 1e-10,\n", " method: str = \"COBYQA\",\n", " options: dict | None = None,\n", ") -> scipy.optimize.OptimizeResult:\n", " \"\"\"Optimize HVA parameters by maximizing fidelity with a target MPS.\n", "\n", " The HVA circuit is simulated as a matrix product state with the given\n", " bond-dimension truncation, and the parameters are optimized to maximize\n", " the state fidelity ``||**2`` with the DMRG ground\n", " state. Both states are normalized, so the minimized objective is the\n", " infidelity ``1 - ||**2``.\n", "\n", " Args:\n", " ansatz: Parameterized HVA circuit from ``build_ground_state_ansatz``.\n", " The length of ``x0`` must equal ``ansatz.num_parameters``.\n", " x0: Initial parameters.\n", " target_mps: Target MPS (DMRG ground state) to maximize fidelity with.\n", " max_bond: Maximum MPS bond dimension during gate application.\n", " cutoff: Singular-value cutoff during gate application.\n", " method: ``scipy.optimize.minimize`` method.\n", " options: Options dict forwarded to ``scipy.optimize.minimize``.\n", "\n", " Returns:\n", " The Scipy OptimizeResult.\n", " \"\"\"\n", "\n", " def infidelity(params: NDArray[np.floating]) -> float:\n", " circuit = ansatz.assign_parameters(params)\n", " circuit_mps = quimb_circuit(\n", " circuit.decompose([\"PauliEvolution\"]),\n", " quimb_circuit_class=qtn.CircuitMPS,\n", " max_bond=max_bond,\n", " cutoff=cutoff,\n", " )\n", " return 1 - abs(circuit_mps.psi.H @ target_mps) ** 2\n", "\n", " return scipy.optimize.minimize(\n", " infidelity, np.asarray(x0), method=method, options=options\n", " )\n", "\n", "\n", "def trotter_evolution(\n", " qubits: Sequence[Qubit],\n", " interaction: float,\n", " anisotropy: float,\n", " time_step: float,\n", " n_steps: int,\n", ") -> Iterator[CircuitInstruction]:\n", " \"\"\"Second-order Trotter steps of the XXZ pair Hamiltonian.\n", "\n", " While the paper used a hand-optimized circuit for the Trotter steps, we use\n", " PauliEvolutionGate here for simplicity and generality. The final two-qubit gate\n", " count and gate depth are equivalent when transpiled with ``optimization_level=3``.\n", "\n", " Args:\n", " qubits: Qubits to act on (length ``n_qubits``).\n", " interaction: Overall interaction strength (J in the paper).\n", " anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the\n", " isotropic Heisenberg point; ``ε = 0`` is the Ising limit.\n", " time_step: Per-step Trotter time.\n", " n_steps: Number of Trotter steps.\n", "\n", " Yields:\n", " ``CircuitInstruction``s implementing the Trotter steps.\n", " \"\"\"\n", " if n_steps == 0:\n", " return\n", " n_qubits = len(qubits)\n", " pair_ham = SparsePauliOp(\n", " [\"XX\", \"YY\", \"ZZ\"],\n", " coeffs=[\n", " 0.25 * interaction * anisotropy,\n", " 0.25 * interaction * anisotropy,\n", " 0.25 * interaction,\n", " ],\n", " )\n", " half_evo = PauliEvolutionGate(pair_ham, time=time_step / 2)\n", " full_evo = PauliEvolutionGate(pair_ham, time=time_step)\n", " for i in range(n_qubits // 2): # half even layer\n", " yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))\n", " for i in range(n_qubits // 2 - 1): # full odd layer\n", " yield CircuitInstruction(\n", " full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])\n", " )\n", " for _ in range(n_steps - 1): # interior steps\n", " for i in range(n_qubits // 2):\n", " yield CircuitInstruction(\n", " full_evo, (qubits[2 * i], qubits[2 * i + 1])\n", " )\n", " for i in range(n_qubits // 2 - 1):\n", " yield CircuitInstruction(\n", " full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])\n", " )\n", " for i in range(n_qubits // 2): # half even layer\n", " yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))\n", "\n", "\n", "def get_dsf(\n", " n_qubits: int,\n", " rgf_mat: NDArray[np.floating],\n", " time_step: float,\n", " n_steps: int,\n", " n_points_momentum: int,\n", " n_points_frequency: int,\n", ") -> NDArray[np.floating]:\n", " \"\"\"Compute the dynamical structure factor from the retarded Green's function.\n", "\n", " Uses the center-site approximation and a discrete Fourier transform.\n", " The result is symmetrized about the momentum axis and clipped to\n", " non-negative values, ready for plotting.\n", "\n", " Args:\n", " n_qubits: Number of qubits (sites).\n", " rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.\n", " time_step: Trotter time-step size.\n", " n_steps: Number of time steps.\n", " n_points_momentum: Number of momentum points.\n", " n_points_frequency: Number of frequency points.\n", "\n", " Returns:\n", " DSF array of shape ``(n_points_frequency, n_points_momentum)``,\n", " symmetrized about the momentum axis and clipped to non-negative values.\n", " \"\"\"\n", " max_frequency = np.pi / time_step\n", " momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)\n", " frequency_range = np.linspace(0, max_frequency, n_points_frequency)\n", " result = np.zeros((frequency_range.shape[0], momentum_range.shape[0]))\n", " center = n_qubits // 2 - 1\n", " for iw, w in enumerate(frequency_range):\n", " exponent = np.exp(1j * w * time_step * np.arange(1, n_steps + 1))\n", " # S = sigma/2, so two spin operators contribute a factor of (1/2)(1/2) = 1/4.\n", " rgf_omega = (\n", " np.dot(rgf_mat.T, exponent) * time_step / 4\n", " ) # S(omega): time Fourier slice of the Green's function\n", " for iq, q in enumerate(momentum_range):\n", " momentum_phases = np.exp(\n", " -1j * q * np.arange(-center, center + 2, 1)\n", " )\n", " result[iw, iq] = np.imag(np.dot(rgf_omega, momentum_phases))\n", " result = -(result + result[:, ::-1]) / 2\n", " result = np.clip(result, a_min=0, a_max=None)\n", " return result\n", "\n", "\n", "def plot_dsf(\n", " dsf: NDArray[np.floating],\n", " time_step: float,\n", " n_points_momentum: int,\n", " n_points_frequency: int,\n", " title: str | None = None,\n", ") -> None:\n", " \"\"\"Heat-map of the dynamical structure factor.\n", "\n", " Args:\n", " dsf: DSF array of shape ``(n_points_frequency, n_points_momentum)``.\n", " time_step: Trotter time-step size.\n", " n_points_momentum: Number of momentum points.\n", " n_points_frequency: Number of frequency points.\n", " title: Optional plot title.\n", " \"\"\"\n", " max_frequency = np.pi / time_step\n", " momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)\n", " frequency_range = np.linspace(0, max_frequency, n_points_frequency)\n", " x, y = np.meshgrid(momentum_range, frequency_range)\n", " fig, ax = plt.subplots(figsize=(8, 5))\n", " c = ax.pcolormesh(x, y, dsf / np.max(dsf), cmap=\"viridis\", shading=\"auto\")\n", " fig.colorbar(c, ax=ax, label=\"Normalized intensity\")\n", " ax.set_ylim(0, 3.6)\n", " ax.set_xlim(0, 2 * np.pi)\n", " ax.set_xlabel(r\"$q$\", fontsize=16)\n", " ax.set_ylabel(r\"$\\tilde{\\omega} = \\omega / J$\", fontsize=16)\n", " ax.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])\n", " ax.set_xticklabels([\"0\", r\"$\\pi/2$\", r\"$\\pi$\", r\"$3\\pi/2$\", r\"$2\\pi$\"])\n", " if title:\n", " ax.set_title(title, fontsize=14)\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "\n", "def plot_rgf(\n", " n_qubits: int,\n", " rgf_mat: NDArray[np.floating],\n", " time_step: float,\n", " n_steps: int,\n", " title: str | None = None,\n", ") -> None:\n", " \"\"\"Heat-map of the retarded Green's function in real space and time.\n", "\n", " Args:\n", " n_qubits: Number of qubits (sites).\n", " rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.\n", " time_step: Trotter time-step size.\n", " n_steps: Number of time steps.\n", " title: Optional plot title.\n", " \"\"\"\n", " fig, ax = plt.subplots(figsize=(8, 6))\n", " qubit_axis = np.arange(n_qubits)\n", " t_axis = np.arange(1, n_steps + 1) * time_step\n", " x, y = np.meshgrid(qubit_axis, t_axis)\n", " c = ax.pcolormesh(\n", " x,\n", " y,\n", " np.real(rgf_mat),\n", " cmap=\"RdBu\",\n", " vmax=0.5,\n", " vmin=-0.5,\n", " shading=\"auto\",\n", " )\n", " fig.colorbar(c, ax=ax, label=r\"Re $G^R(j, j_c, t)$\")\n", " ax.set_xlabel(\"Qubit\", fontsize=16)\n", " ax.xaxis.set_major_locator(\n", " plt.matplotlib.ticker.MaxNLocator(integer=True)\n", " )\n", " ax.set_ylabel(r\"Time\", fontsize=16)\n", " if title:\n", " ax.set_title(title, fontsize=14)\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "\n", "def uniform_2q_depth(circuit: QuantumCircuit) -> int:\n", " \"\"\"Two-qubit gate depth in a standardized basis.\"\"\"\n", " pass_manager = generate_preset_pass_manager(\n", " optimization_level=0, basis_gates=[\"cz\", \"id\", \"rz\", \"sx\", \"x\"]\n", " )\n", " return pass_manager.run(circuit).depth(\n", " lambda inst: inst.operation.num_qubits == 2\n", " )" ] }, { "cell_type": "markdown", "id": "d28d01a6", "metadata": {}, "source": [ "\n", "\n", "## Exemple de simulateur à petite échelle\n", "\n", "Nous présentons tout d'abord le flux de travail complet sur **10 qubits**, en optimisant l'ansatz de l'état fondamental HVA à l'aide d'une simulation MPS et en utilisant le simulateur de vecteurs d'état de Qiskit pour l'évolution temporelle. Le DMRG fournit une énergie de référence de l'état fondamental et un MPS de référence. Cet exemple à petite échelle nous permet de valider chaque étape avant de passer à une plus grande échelle.\n", "\n" ] }, { "cell_type": "markdown", "id": "3818ed8f", "metadata": {}, "source": [ "\n", "\n", "### Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique\n", "\n", "Nous commençons par définir le modèle physique et par construire les circuits quantiques.\n", "\n", "**hamiltonien.** KCuF$_3$ est modélisé par l'hamiltonien XXZ de type « 1D » au point isotrope ( $J = 1$, $\\epsilon = 1$, en prenant $J = 1$ comme unité d'énergie).\n", "\n", "**État fondamental.** Nous utilisons un circuit d'approche variationnelle `build_ground_state_ansatz` hamiltonienne (HVA) comme circuit de préparation de l'état fondamental. `CircuitMPS`Les paramètres HVA sont optimisés de manière classique en maximisant la fidélité d'état $|\\langle\\psi_{\\mathrm{HVA}}(\\theta)|\\psi_{\\mathrm{DMRG}}\\rangle|^2$ à l'aide de la méthode MPS de l'état fondamental DMRG, où l'état HVA est évalué en simulant le circuit sous la forme d'un état de produit matriciel à l'aide de quimb. Nous utilisons `scipy.optimize.minimize` pour l'optimisation. L' $\\langle H \\rangle$ e énergétique de l'ansatz optimisé est également calculée à titre de référence.\n", "\n", "**Portails pour trotteurs.** Chaque terme d'interaction « plus proche voisin » $e^{-i\\Delta t\\, H_{\\mathrm{pair}}}$, où $H_{\\mathrm{pair}} = (J/4)(XX + YY + ZZ)$, est construit à l'aide de `PauliEvolutionGate(H_pair, time=time_step)`. Qiskit synthétise cela en une décomposition optimale en trois CNOT lors de la transpilation.\n", "\n", "**Perturbation.** Une porte « $R_z(\\pi/2)$ » appliquée au qubit central met en œuvre une opération « $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ », imitant le renversement de spin local provoqué par un neutron diffusé.\n", "\n", "**Observables.** Nous construisons une observable de type « $\\sigma^z$ » pour chaque site de qubit. Ces `SparsePauliOp` objets sont transmis à la primitive « Estimator » à l'étape 3 afin d'extraire l' $\\langle\\sigma_i^z\\rangle$ s à chaque pas de temps.\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "58305fbc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ground-state energy (DMRG): -4.258035\n", "Optimizing ground state ansatz...\n", "Finished optimizing ground state ansatz in 8.850687621976249 seconds.\n", "Ground state ansatz fidelity: 0.984277\n", "Ground state ansatz energy: -4.232565\n", "Built 10 circuits, deepest 2q depth (uniform basis) = 163\n", "Defined 10 Z observables.\n" ] } ], "source": [ "# -- Physical parameters --\n", "n_qubits = 10\n", "interaction = 1.0 # J\n", "anisotropy = 1.0 # ε (isotropic point)\n", "time_step = 0.6\n", "n_steps = 10\n", "mps_max_bond = 32\n", "mps_cutoff = 1e-8\n", "center = n_qubits // 2 - 1\n", "\n", "# -- Hamiltonian MPO --\n", "ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)\n", "\n", "# -- Reference ground-state energy via DMRG --\n", "dmrg = qtn.DMRG2(ham_mpo)\n", "dmrg.solve(tol=1e-8)\n", "print(f\"Ground-state energy (DMRG): {dmrg.energy:.6f}\")\n", "\n", "# -- Build ground state ansatz circuit --\n", "gs_n_layers = 3\n", "gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)\n", "\n", "# -- Optimize ground state ansatz parameters to maximize fidelity with the DMRG MPS --\n", "# Initialize odd-layer angles near 0 (where the inter-pair gate is the\n", "# identity) and even-layer angles near pi/2 (where the intra-pair gate\n", "# is a SWAP, since 0.5*(XX + YY + ZZ) = SWAP - I/2).\n", "rng = np.random.default_rng(12345)\n", "x0 = np.tile([0, np.pi / 2], gs_n_layers) + rng.normal(\n", " scale=0.1, size=2 * gs_n_layers\n", ")\n", "print(\"Optimizing ground state ansatz...\")\n", "t0 = timeit.default_timer()\n", "result = optimize_ground_state_ansatz(\n", " gs_ansatz,\n", " x0,\n", " dmrg.state,\n", " max_bond=mps_max_bond,\n", " cutoff=mps_cutoff,\n", " options=dict(maxiter=100),\n", ")\n", "t1 = timeit.default_timer()\n", "print(f\"Finished optimizing ground state ansatz in {t1 - t0} seconds.\")\n", "print(f\"Ground state ansatz fidelity: {1 - result.fun:.6f}\")\n", "\n", "gs_circuit = gs_ansatz.assign_parameters(result.x)\n", "gs_circuit_mps = quimb_circuit(\n", " gs_circuit.decompose([\"PauliEvolution\"]),\n", " quimb_circuit_class=qtn.CircuitMPS,\n", " max_bond=mps_max_bond,\n", " cutoff=mps_cutoff,\n", ")\n", "gs_ansatz_energy = qtn.expec_TN_1D(\n", " gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi\n", ")\n", "print(f\"Ground state ansatz energy: {gs_ansatz_energy:.6f}\")\n", "\n", "# -- Build circuits for each time step --\n", "perturbed = gs_circuit.copy()\n", "perturbed.rz(np.pi / 2, center)\n", "\n", "circuits = []\n", "for t in range(1, n_steps + 1):\n", " circuit = perturbed.copy()\n", " for instr in trotter_evolution(\n", " circuit.qubits, interaction, anisotropy, time_step, t\n", " ):\n", " circuit.append(instr)\n", " circuits.append(circuit)\n", "print(\n", " f\"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = \"\n", " f\"{uniform_2q_depth(circuits[-1])}\"\n", ")\n", "\n", "# -- Observables: Z on each qubit site --\n", "observables = [\n", " SparsePauliOp.from_sparse_list([(\"Z\", [i], 1)], num_qubits=n_qubits)\n", " for i in range(n_qubits)\n", "]\n", "print(f\"Defined {len(observables)} Z observables.\")" ] }, { "cell_type": "markdown", "id": "39191552", "metadata": {}, "source": [ "\n", "\n", "### Étape 2 : Optimiser le problème en vue de son exécution sur un matériel quantique\n", "\n", "Pour du matériel réel, les circuits de Trotter présentés ci-dessus seraient trop complexes. **La compilation quantique approximative (AQC)** résout ce problème en remplaçant les premières couches de Trotter de l' $k$ e (y compris le circuit à l'état fondamental) par un ansatz paramétré plus court, optimisé pour maximiser la fidélité au niveau MPS par rapport au circuit profond d'origine. Les étapes restantes de l'algorithme de Trotter sont ajoutées à l'identique, ce qui donne un circuit « AQC + Trotter » moins profond.\n", "\n", "Afin de trouver un équilibre entre l'expressivité et la profondeur du circuit, deux approches sont utilisées : une **approche à une couche** (générée à partir d'une seule étape de Trotter) compresse les premières étapes temporelles, tandis qu'une **approche** plus profonde à deux couches (générée à partir de deux étapes de Trotter) compresse les étapes suivantes, pour lesquelles une plus grande fidélité est requise.\n", "\n", "Le processus AQC comprend quatre sous-étapes :\n", "\n", "1. **Construire des circuits cibles** — les premiers circuits d’ $k_1 + k_2$ s de l’étape 1 servent directement de cibles pour l’AQC.\n", "2. **Calculer le MPS cible** — simuler chaque circuit cible sous la forme d'un état de produit matriciel à l'aide de `quimb.tensor.CircuitMPS`.\n", "3. **Générer et optimiser des ansätze** — `generate_ansatz_from_circuit` crée un ansatz paramétré à une couche et un autre à deux couches; les paramètres sont optimisés à l’aide de l’algorithme L-BFGS-B avec des gradients accélérés par JAX afin de minimiser l’ $1 - |\\langle\\psi_{\\mathrm{ansatz}}|\\psi_{\\mathrm{target}}\\rangle|^2$. Les premières étapes de l’ $k_1$ utilisent l’ansatz à une couche et les étapes suivantes $k_2$ utilisent l’ansatz à deux couches; au sein de chaque étape, chaque itération démarre à chaud à partir des paramètres optimisés de l’itération précédente, et les paramètres sont réinitialisés aux valeurs par défaut de l’étape à la limite de celle-ci.\n", "4. **Assembler des circuits mixtes** : pour les pas de temps situés au-delà du point de contrôle AQC, ajouter des couches de Trotter exactes au circuit AQC optimisé à deux couches à l'aide de `trotter_evolution`.\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "e78f8d9a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Step 2b — target MPS:\n", " k=1: max bond = 22\n", " k=2: max bond = 22\n", " k=3: max bond = 26\n", " k=4: max bond = 27\n", " k=5: max bond = 30\n", "\n", "Step 2c — one-layer ansatz: 389 parameters, 2q depth (uniform basis) = 27\n", "Step 2c — two-layer ansatz: 470 parameters, 2q depth (uniform basis) = 33\n", " k=1: fidelity = 1.0000, 2q depth (uniform basis) = 27, 11.8s\n", " k=2: fidelity = 0.9980, 2q depth (uniform basis) = 27, 13.2s\n", " k=3: fidelity = 0.9890, 2q depth (uniform basis) = 27, 13.7s\n", " k=4: fidelity = 0.9983, 2q depth (uniform basis) = 33, 16.4s\n", " k=5: fidelity = 0.9957, 2q depth (uniform basis) = 33, 16.4s\n", "\n", "Step 2d — assembled 10 circuits\n", " At step 10 (uniform basis): Trotter 2q depth = 163, AQC+Trotter 2q depth = 99; 2q gates: Trotter = 371, AQC+Trotter = 200\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Number of time steps to compress into an AQC ansatz with one layer\n", "aqc_n_steps_1 = 3\n", "# Number of time steps to compress into an AQC ansatz with two layers\n", "aqc_n_steps_2 = 2\n", "aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2\n", "\n", "# ── Step 2a: Target circuits (first aqc_n_steps_total circuits from Step 1) ──\n", "target_circuits = {\n", " k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)\n", "}\n", "# ── Step 2b: Compute target MPS ──\n", "# Decompose PauliEvolutionGate to RXX/RYY/RZZ before passing to the AQC MPS\n", "# backend, which does not understand PauliEvolutionGate natively.\n", "aqc_sim = QuimbSimulator(\n", " quimb_circuit_factory=partial(\n", " qtn.CircuitMPS,\n", " gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),\n", " ),\n", " autodiff_backend=\"jax\",\n", ")\n", "print(\"\\nStep 2b — target MPS:\")\n", "target_mps = {}\n", "for k in range(1, aqc_n_steps_total + 1):\n", " target_mps[k] = tensornetwork_from_circuit(\n", " target_circuits[k].decompose([\"PauliEvolution\"]), aqc_sim\n", " )\n", " print(f\" k={k}: max bond = {target_mps[k].psi.max_bond()}\")\n", "\n", "# ── Step 2c: Generate ansätze (one-layer and two-layer) and optimize parameters ──\n", "ansatz_1, initial_params_1 = generate_ansatz_from_circuit(\n", " target_circuits[1].decompose([\"PauliEvolution\"]),\n", " qubits_initially_zero=True,\n", ")\n", "initial_params_1 = np.array(initial_params_1)\n", "ansatz_2, initial_params_2 = generate_ansatz_from_circuit(\n", " target_circuits[2].decompose([\"PauliEvolution\"]),\n", " qubits_initially_zero=True,\n", ")\n", "initial_params_2 = np.array(initial_params_2)\n", "print(\n", " f\"\\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, \"\n", " f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}\"\n", ")\n", "print(\n", " f\"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, \"\n", " f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}\"\n", ")\n", "\n", "aqc_circuits = {}\n", "aqc_params = {}\n", "for k in range(1, aqc_n_steps_total + 1):\n", " if k <= aqc_n_steps_1:\n", " ansatz, base_params = ansatz_1, initial_params_1\n", " else:\n", " ansatz, base_params = ansatz_2, initial_params_2\n", " # Warm-start from the previous step only within the same stage\n", " same_stage = (k - 1 >= 1) and (\n", " (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)\n", " )\n", " x0 = aqc_params[k - 1] if same_stage else base_params\n", " obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)\n", " t0 = timeit.default_timer()\n", " result = scipy.optimize.minimize(\n", " obj.loss_function,\n", " x0,\n", " method=\"L-BFGS-B\",\n", " jac=True,\n", " options=dict(maxiter=100),\n", " )\n", " elapsed = timeit.default_timer() - t0\n", " aqc_params[k] = result.x\n", " aqc_circuits[k] = ansatz.assign_parameters(result.x)\n", " print(\n", " f\" k={k}: fidelity = {1 - result.fun:.4f}, \"\n", " f\"2q depth (uniform basis) = \"\n", " f\"{uniform_2q_depth(aqc_circuits[k])}, \"\n", " f\"{elapsed:.1f}s\"\n", " )\n", "\n", "# ── Step 2d: Assemble full circuit set (AQC + Trotter) ──\n", "all_circuits = []\n", "for k in range(1, aqc_n_steps_total + 1):\n", " all_circuits.append(aqc_circuits[k])\n", "base = aqc_circuits[aqc_n_steps_total]\n", "for k in range(1, n_steps - aqc_n_steps_total + 1):\n", " circuit = base.copy()\n", " for instr in trotter_evolution(\n", " circuit.qubits, interaction, anisotropy, time_step, k\n", " ):\n", " circuit.append(instr)\n", " all_circuits.append(circuit)\n", "\n", "full_depths = [\n", " uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)\n", "]\n", "aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]\n", "full_2q = [\n", " circuits[k - 1].decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n", " for k in range(1, n_steps + 1)\n", "]\n", "aqc_2q = [\n", " circuit.decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n", " for circuit in all_circuits\n", "]\n", "print(f\"\\nStep 2d — assembled {len(all_circuits)} circuits\")\n", "print(\n", " f\" At step {n_steps} (uniform basis): \"\n", " f\"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; \"\n", " f\"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}\"\n", ")\n", "\n", "steps_axis = np.arange(1, n_steps + 1)\n", "fig, ax = plt.subplots(figsize=(7, 4))\n", "ax.plot(steps_axis, full_depths, \"-o\", color=\"black\", label=\"Full Trotter\")\n", "ax.plot(\n", " steps_axis, aqc_depths, \"-o\", color=\"cadetblue\", label=\"AQC + Trotter\"\n", ")\n", "ax.set_xlabel(\"Trotter steps\", fontsize=13)\n", "ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))\n", "ax.set_ylabel(\"2q gate depth (uniform basis)\", fontsize=13)\n", "ax.set_title(f\"AQC circuit-depth reduction ({n_qubits} qubits)\", fontsize=13)\n", "ax.legend(fontsize=11)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "43ad6048", "metadata": {}, "source": [ "\n", "\n", "### Étape 3 : Exécutez la commande à l'aide d' Qiskit primitives\n", "\n", "Nous simulons chaque circuit compilé par AQC à l'aide de la `StatevectorEstimator` primitive.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "b8f9e41c", "metadata": {}, "outputs": [], "source": [ "estimator = StatevectorEstimator()\n", "pubs = [(circuit, observables) for circuit in all_circuits]\n", "job = estimator.run(pubs)\n", "result = job.result()" ] }, { "cell_type": "markdown", "id": "cb12662c", "metadata": {}, "source": [ "\n", "\n", "### Étape 4 : Traitement ultérieur et restitution du résultat dans le format classique souhaité\n", "\n", "Nous calculons ensuite la valeur attendue $\\langle\\sigma_i^z\\rangle$ pour chaque qubit $i$ à chaque pas de temps $t$. Ces valeurs forment la matrice de la fonction de Green retardée $G^R(j, j_c, t)$. Nous appliquons ensuite la transformée de Fourier à la fonction de Green retardée pour obtenir le facteur de structure dynamique $S(q, \\omega)$, puis nous représentons graphiquement à la fois la fonction de Green retardée et le facteur de structure dynamique. `get_dsf` applique la symétrie miroir et limite les valeurs négatives en interne.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "4585d139", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rgf_mat = np.stack([pub_result.data.evs for pub_result in result])\n", "\n", "# -- Compute DSF --\n", "n_points_momentum, n_points_frequency = 100, 100\n", "spectrum = get_dsf(\n", " n_qubits,\n", " rgf_mat,\n", " time_step,\n", " n_steps,\n", " n_points_momentum,\n", " n_points_frequency,\n", ")\n", "\n", "# -- Plot retarded Green's function --\n", "plot_rgf(\n", " n_qubits,\n", " rgf_mat,\n", " time_step,\n", " n_steps,\n", " title=f\"Retarded Green's function — {n_qubits} qubits (simulation)\",\n", ")\n", "\n", "# -- Plot DSF --\n", "plot_dsf(\n", " spectrum,\n", " time_step,\n", " n_points_momentum,\n", " n_points_frequency,\n", " title=f\"Dynamical structure factor — {n_qubits} qubits (simulation)\",\n", ")" ] }, { "cell_type": "markdown", "id": "ls-intro", "metadata": {}, "source": [ "\n", "\n", "## Exécution matérielle à grande échelle\n", "\n", "Nous passons désormais à **50 qubits**. À cette échelle, les optimisations atteignent des fidélités inférieures à celles de l’exemple à petite échelle : la fidélité de l’ansatz de l’état fondamental chute à environ 0.65, et les fidélités AQC diminuent jusqu’à environ 0.7 aux derniers points de contrôle. C'est tout à fait normal, et vous pouvez améliorer ces précisions en augmentant le nombre de couches d'ansatz de l'état fondamental (`gs_n_layers`) ou le nombre d'itérations d'optimisation (`maxiter`) au prix d'un surcoût classique supplémentaire. Il convient également de noter que les fidélités de l'AQC à deux couches s'avèrent inférieures à celles de l'AQC à une seule couche. Il ne s'agit pas d'une régression : les pas de temps ultérieurs génèrent davantage d'intrication et sont tout simplement plus difficiles à compresser, ce qui explique pourquoi on utilise pour eux l'ansatz à deux couches, plus expressif.\n", "\n", "L'optimisation AQC peut également nécessiter plusieurs heures de calcul classique (environ six heures dans l'exécution présentée ici, dont la majeure partie est consacrée aux points de contrôle à deux couches de l'étape 2c ). Pour réduire le temps réel, envisagez d'exécuter ce notebook sur du matériel classique plus puissant, tel qu'un système de calcul haute performance (HPC). Vous pouvez également réduire la portée du problème (par exemple, en diminuant le nombre de qubits ou le nombre d'étapes temporelles); vos résultats seront alors différents de ceux présentés ici.\n", "\n", "Le code ci-dessous suit la même structure en quatre étapes que l'exemple à petite échelle. Les paramètres de l'état fondamental de la HVA sont à nouveau optimisés à l'aide d'une simulation MPS. Sur le QPU, nous activons le découplage dynamique (DD), le « Pauli twirling » et l'extinction des erreurs de lecture par rotation (TREX) afin de supprimer et d'atténuer les erreurs. Le tableau suivant résume les différences entre l'expérience à grande échelle et celle à petite échelle :\n", "\n", "|   | À petite échelle | À grande échelle |\n", "| --------------------------------------------- | ---------------------- | ----------------------------------- |\n", "| Qubits | 10 | 50 |\n", "| Pas de temps | 10 | 20 |\n", "| Points de contrôle AQC (1 couche + 2 couches) | 3 + 2 = 5 | 6 + 4 = 10 |\n", "| Couches d'ansatz de l'état fondamental | 3 | 5 |\n", "| Dimension maximale des liaisons MPS | 32 | 128 |\n", "| Estimateur | `StatevectorEstimator` | QPU avec DD, effet de Pauli et TREX |\n", "\n" ] }, { "cell_type": "markdown", "id": "9e85fbd1", "metadata": {}, "source": [ "\n", " Au cours de l'optimisation AQC (étape « 2c » ci-dessous), vous pourriez voir `stderr` s'afficher des messages tels que :\n", "\n", " ```\n", " [Compiling module jit_MakeArrayFn for CPU] Very slow compile? ...\n", " The operation took 2m14s\n", " ```\n", "\n", " Il s'agit de messages de diagnostic sans gravité provenant de XLA, le compilateur qui sous-tend `qiskit-addon-aqc-tensor`la fonctionnalité « autodiff » de JAX. Avec 50 qubits et une dimension de liaison MPS de 128, XLA met quelques minutes à compiler la fonction de gradient lors de son premier calcul. La compilation aboutit, et les résultats de l'optimisation ne sont pas affectés.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "fbe57e1a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ground-state energy (DMRG): -21.972109\n", "Optimizing ground state ansatz...\n", "Finished optimizing ground state ansatz in 132.2076231740648 seconds.\n", "Ground state ansatz fidelity: 0.645956\n", "Ground state ansatz energy: -21.616744\n", "Built 20 circuits, deepest 2q depth (uniform basis) = 307\n", "Defined 50 Z observables.\n", "\n", "Step 2b — target MPS:\n", " k=1: max bond = 44\n", " k=2: max bond = 46\n", " k=3: max bond = 53\n", " k=4: max bond = 62\n", " k=5: max bond = 75\n", " k=6: max bond = 96\n", " k=7: max bond = 118\n", " k=8: max bond = 128\n", " k=9: max bond = 128\n", " k=10: max bond = 128\n", "\n", "Step 2c — one-layer ansatz: 2971 parameters, 2q depth (uniform basis) = 39\n", "Step 2c — two-layer ansatz: 3412 parameters, 2q depth (uniform basis) = 45\n", " k=1: fidelity = 1.0000, 2q depth (uniform basis) = 39, 215.6s\n", " k=2: fidelity = 0.9676, 2q depth (uniform basis) = 39, 269.8s\n", " k=3: fidelity = 0.9014, 2q depth (uniform basis) = 39, 293.5s\n", " k=4: fidelity = 0.8755, 2q depth (uniform basis) = 39, 279.6s\n", " k=5: fidelity = 0.8450, 2q depth (uniform basis) = 39, 266.4s\n", " k=6: fidelity = 0.8146, 2q depth (uniform basis) = 39, 371.0s\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "E0626 02:59:07.229070 1508496 slow_operation_alarm.cc:73] \n", "********************************\n", "[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.\n", "********************************\n", "E0626 02:59:21.713796 1508477 slow_operation_alarm.cc:140] The operation took 2m14.484932074s\n", "\n", "********************************\n", "[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.\n", "********************************\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " k=7: fidelity = 0.7312, 2q depth (uniform basis) = 45, 10865.3s\n", " k=8: fidelity = 0.7720, 2q depth (uniform basis) = 45, 1450.3s\n", " k=9: fidelity = 0.7316, 2q depth (uniform basis) = 45, 3388.6s\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "E0626 07:21:05.148724 1508496 slow_operation_alarm.cc:73] \n", "********************************\n", "[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.\n", "********************************\n", "E0626 07:21:24.286095 1508477 slow_operation_alarm.cc:140] The operation took 2m19.137566498s\n", "\n", "********************************\n", "[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.\n", "********************************\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " k=10: fidelity = 0.6726, 2q depth (uniform basis) = 45, 3396.1s\n", "\n", "Step 2d — assembled 20 circuits\n", " At step 20 (uniform basis): Trotter 2q depth = 307, AQC+Trotter 2q depth = 171; 2q gates: Trotter = 3775, AQC+Trotter = 1913\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Backend: ibm_fez\n", "Transpiled 2q depth (deepest, ISA on ibm_fez): 105 (2q gates: 2574)\n", "Job ID: d8v39vhropqc738biotg\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# ── Parameters ──────────────────────────────────────────────────────────────\n", "n_qubits = 50 # 10 → 50\n", "interaction = 1.0 # J\n", "anisotropy = 1.0 # ε (isotropic point)\n", "time_step = 0.6\n", "n_steps = 20 # 10 → 20\n", "aqc_n_steps_1 = 6 # 3 → 6\n", "aqc_n_steps_2 = 4 # 2 → 4\n", "aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2\n", "gs_n_layers = 5 # 3 → 5\n", "mps_max_bond = 128 # 32 -> 128\n", "mps_cutoff = 1e-8\n", "center = n_qubits // 2 - 1\n", "\n", "# ── Step 1: Map ──────────────────────────────────────────────────────────────\n", "ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)\n", "\n", "dmrg = qtn.DMRG2(ham_mpo)\n", "dmrg.solve(tol=1e-8)\n", "print(f\"Ground-state energy (DMRG): {dmrg.energy:.6f}\")\n", "\n", "gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)\n", "\n", "rng = np.random.default_rng(12345)\n", "x0 = np.tile([0.0, np.pi / 2], gs_n_layers) + rng.normal(\n", " scale=0.1, size=2 * gs_n_layers\n", ")\n", "print(\"Optimizing ground state ansatz...\")\n", "t0 = timeit.default_timer()\n", "result = optimize_ground_state_ansatz(\n", " gs_ansatz,\n", " x0,\n", " dmrg.state,\n", " max_bond=mps_max_bond,\n", " cutoff=mps_cutoff,\n", " options=dict(maxiter=100),\n", ")\n", "t1 = timeit.default_timer()\n", "print(f\"Finished optimizing ground state ansatz in {t1 - t0} seconds.\")\n", "print(f\"Ground state ansatz fidelity: {1 - result.fun:.6f}\")\n", "\n", "gs_circuit = gs_ansatz.assign_parameters(result.x)\n", "gs_circuit_mps = quimb_circuit(\n", " gs_circuit.decompose([\"PauliEvolution\"]),\n", " quimb_circuit_class=qtn.CircuitMPS,\n", " max_bond=mps_max_bond,\n", " cutoff=mps_cutoff,\n", ")\n", "gs_ansatz_energy = qtn.expec_TN_1D(\n", " gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi\n", ")\n", "print(f\"Ground state ansatz energy: {gs_ansatz_energy:.6f}\")\n", "\n", "perturbed = gs_circuit.copy()\n", "perturbed.rz(np.pi / 2, center)\n", "\n", "circuits = []\n", "for t in range(1, n_steps + 1):\n", " circuit = perturbed.copy()\n", " for instr in trotter_evolution(\n", " circuit.qubits, interaction, anisotropy, time_step, t\n", " ):\n", " circuit.append(instr)\n", " circuits.append(circuit)\n", "print(\n", " f\"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = \"\n", " f\"{uniform_2q_depth(circuits[-1])}\"\n", ")\n", "\n", "observables = [\n", " SparsePauliOp.from_sparse_list([(\"Z\", [i], 1)], num_qubits=n_qubits)\n", " for i in range(n_qubits)\n", "]\n", "print(f\"Defined {len(observables)} Z observables.\")\n", "\n", "# ── Step 2: AQC ──────────────────────────────────────────────────────────────\n", "target_circuits = {\n", " k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)\n", "}\n", "\n", "aqc_sim = QuimbSimulator(\n", " quimb_circuit_factory=partial(\n", " qtn.CircuitMPS,\n", " gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),\n", " ),\n", " autodiff_backend=\"jax\",\n", ")\n", "print(\"\\nStep 2b — target MPS:\")\n", "target_mps = {}\n", "for k in range(1, aqc_n_steps_total + 1):\n", " target_mps[k] = tensornetwork_from_circuit(\n", " target_circuits[k].decompose([\"PauliEvolution\"]), aqc_sim\n", " )\n", " print(f\" k={k}: max bond = {target_mps[k].psi.max_bond()}\")\n", "\n", "ansatz_1, initial_params_1 = generate_ansatz_from_circuit(\n", " target_circuits[1].decompose([\"PauliEvolution\"]),\n", " qubits_initially_zero=True,\n", ")\n", "initial_params_1 = np.array(initial_params_1)\n", "ansatz_2, initial_params_2 = generate_ansatz_from_circuit(\n", " target_circuits[2].decompose([\"PauliEvolution\"]),\n", " qubits_initially_zero=True,\n", ")\n", "initial_params_2 = np.array(initial_params_2)\n", "print(\n", " f\"\\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, \"\n", " f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}\"\n", ")\n", "print(\n", " f\"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, \"\n", " f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}\"\n", ")\n", "\n", "aqc_circuits = {}\n", "aqc_params = {}\n", "for k in range(1, aqc_n_steps_total + 1):\n", " if k <= aqc_n_steps_1:\n", " ansatz, base_params = ansatz_1, initial_params_1\n", " else:\n", " ansatz, base_params = ansatz_2, initial_params_2\n", " # Warm-start from the previous step only within the same stage\n", " same_stage = (k - 1 >= 1) and (\n", " (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)\n", " )\n", " x0 = aqc_params[k - 1] if same_stage else base_params\n", " obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)\n", " t0 = timeit.default_timer()\n", " result = scipy.optimize.minimize(\n", " obj.loss_function,\n", " x0,\n", " method=\"L-BFGS-B\",\n", " jac=True,\n", " options=dict(maxiter=100),\n", " )\n", " elapsed = timeit.default_timer() - t0\n", " aqc_params[k] = result.x\n", " aqc_circuits[k] = ansatz.assign_parameters(result.x)\n", " print(\n", " f\" k={k}: fidelity = {1 - result.fun:.4f}, \"\n", " f\"2q depth (uniform basis) = \"\n", " f\"{uniform_2q_depth(aqc_circuits[k])}, \"\n", " f\"{elapsed:.1f}s\"\n", " )\n", "\n", "all_circuits = []\n", "for k in range(1, aqc_n_steps_total + 1):\n", " all_circuits.append(aqc_circuits[k])\n", "base = aqc_circuits[aqc_n_steps_total]\n", "for k in range(1, n_steps - aqc_n_steps_total + 1):\n", " circuit = base.copy()\n", " for instr in trotter_evolution(\n", " circuit.qubits, interaction, anisotropy, time_step, k\n", " ):\n", " circuit.append(instr)\n", " all_circuits.append(circuit)\n", "\n", "full_depths = [\n", " uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)\n", "]\n", "aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]\n", "full_2q = [\n", " circuits[k - 1].decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n", " for k in range(1, n_steps + 1)\n", "]\n", "aqc_2q = [\n", " circuit.decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n", " for circuit in all_circuits\n", "]\n", "print(f\"\\nStep 2d — assembled {len(all_circuits)} circuits\")\n", "print(\n", " f\" At step {n_steps} (uniform basis): \"\n", " f\"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; \"\n", " f\"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}\"\n", ")\n", "\n", "steps_axis = np.arange(1, n_steps + 1)\n", "fig, ax = plt.subplots(figsize=(7, 4))\n", "ax.plot(steps_axis, full_depths, \"-o\", color=\"black\", label=\"Full Trotter\")\n", "ax.plot(\n", " steps_axis, aqc_depths, \"-o\", color=\"cadetblue\", label=\"AQC + Trotter\"\n", ")\n", "ax.set_xlabel(\"Trotter steps\", fontsize=13)\n", "ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))\n", "ax.set_ylabel(\"2q gate depth (uniform basis)\", fontsize=13)\n", "ax.set_title(f\"AQC circuit-depth reduction ({n_qubits} qubits)\", fontsize=13)\n", "ax.legend(fontsize=11)\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# ── Step 3: Execute on IBM Quantum hardware ───────────────────────────────────\n", "# (replaces StatevectorEstimator)\n", "service = QiskitRuntimeService()\n", "backend = service.least_busy(\n", " min_num_qubits=n_qubits,\n", " operational=True,\n", " simulator=False,\n", " filters=lambda x: x.configuration().processor_type[\"family\"] == \"Heron\",\n", ")\n", "print(f\"Backend: {backend.name}\")\n", "\n", "pm = generate_preset_pass_manager(optimization_level=3, backend=backend)\n", "isa_circuits = pm.run(all_circuits, num_processes=1)\n", "isa_2q_depths = [\n", " isa_circuit.depth(lambda inst: inst.operation.num_qubits == 2)\n", " for isa_circuit in isa_circuits\n", "]\n", "print(\n", " f\"Transpiled 2q depth (deepest, ISA on {backend.name}): \"\n", " f\"{max(isa_2q_depths)} \"\n", " f\"(2q gates: {max(isa_circuit.num_nonlocal_gates() for isa_circuit in isa_circuits)})\"\n", ")\n", "\n", "estimator = Estimator(backend)\n", "estimator.options.environment.job_tags = [\"TUT_SNS\"]\n", "estimator.options.dynamical_decoupling.enable = True\n", "estimator.options.dynamical_decoupling.sequence_type = \"XY4\"\n", "estimator.options.twirling.enable_gates = True\n", "estimator.options.twirling.num_randomizations = 1000\n", "estimator.options.twirling.shots_per_randomization = 128\n", "estimator.options.resilience.measure_mitigation = True\n", "estimator.options.resilience.measure_noise_learning.num_randomizations = 32\n", "estimator.options.resilience.measure_noise_learning.shots_per_randomization = 100\n", "\n", "pubs = [\n", " (\n", " isa_circuit,\n", " [obs.apply_layout(isa_circuit.layout) for obs in observables],\n", " )\n", " for isa_circuit in isa_circuits\n", "]\n", "job = estimator.run(pubs)\n", "print(f\"Job ID: {job.job_id()}\")\n", "\n", "result = job.result()\n", "\n", "# ── Step 4: Post-process ──────────────────────────────────────────────────────\n", "rgf_mat = np.stack([pub_result.data.evs for pub_result in result])\n", "\n", "n_points_momentum, n_points_frequency = 100, 100\n", "spectrum = get_dsf(\n", " n_qubits,\n", " rgf_mat,\n", " time_step,\n", " n_steps,\n", " n_points_momentum,\n", " n_points_frequency,\n", ")\n", "\n", "plot_rgf(\n", " n_qubits,\n", " rgf_mat,\n", " time_step,\n", " n_steps,\n", " title=f\"Retarded Green's function — {n_qubits} qubits (QPU)\",\n", ")\n", "plot_dsf(\n", " spectrum,\n", " time_step,\n", " n_points_momentum,\n", " n_points_frequency,\n", " title=rf\"KCuF$_3$ DSF — {n_qubits} qubits (QPU)\"\n", " \"\\n(AQC + DD + Pauli twirling + TREX)\",\n", ")" ] }, { "cell_type": "markdown", "id": "ls-discussion", "metadata": {}, "source": [ "Les résultats obtenus sur le matériel reproduisent les caractéristiques principales du continuum à deux spinons : l'intensité de diffusion est concentrée près de l' $q = \\pi$ du vecteur d'onde antiferromagnétique à basse énergie et est limitée par le bas par la dispersion sinusoïdale des spinons, avec au-dessus un large continuum de poids spectral plutôt qu'un mode unique et bien défini. Il s’agit de la même structure que celle mesurée par diffusion inélastique de neutrons sur KCuF$_3$, ce qui valide l’ensemble du processus, depuis la préparation de l’état fondamental jusqu’à la perturbation, en passant par l’évolution de Trotter compressée par AQC et la mesure à erreur atténuée, à l’échelle de 50 qubits.\n", "\n" ] }, { "cell_type": "markdown", "id": "ls-nextsteps", "metadata": {}, "source": [ "\n", "\n", "## Etapes suivantes\n", "\n", "Si ce travail vous a paru intéressant, les documents suivants pourraient vous intéresser :\n", "\n", "\n", " * [Lee et al., « Évaluation comparative de la simulation quantique à l'aide d'expériences de diffusion des neutrons » ( arXiv:2603.15608 )](https://arxiv.org/abs/2603.15608) — l'article de référence sur lequel s'appuie ce tutoriel\n", " * [Techniques d'atténuation et de suppression des erreurs](/docs/guides/error-mitigation-and-suppression-techniques) — DD, « Pauli twirling » et TREX utilisées dans les expériences matérielles\n", " * [Compilation quantique approximative pour les circuits d'évolution temporelle](/docs/tutorials/approximate-quantum-compilation-for-time-evolution) — tutoriel sur AQC-Tensor\n", " * [Documentation AQC-Tensor](https://qiskit.github.io/qiskit-addon-aqc-tensor/)\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" }, "hours": 1.5, "qpuSeconds": 780 }, "nbformat": 4, "nbformat_minor": 5 }