{ "cells": [ { "cell_type": "markdown", "id": "8e82ead7", "metadata": {}, "source": [ "---\n", "title: \"量子回路を用いて量子材料における中性子散乱をシミュレートする\"\n", "description: \"トロッター回路とMPSシミュレーションを用いて、量子磁石の動的構造因子を計算する。\"\n", "---\n", "\n", "{/* cspell:ignore DSFs spinon DMRG viridis fontsize vmax vmin cadetblue antiferromagnet COBYQA */}\n", "\n", "\n", "\n", "# 量子回路を用いて量子材料における中性子散乱をシミュレートする\n", "\n", "*推定実行時間：Heron r2 プロセッサで 13 分（注：これはあくまで推定値です。 （実行時間は状況によって異なる場合があります。）*\n", "\n" ] }, { "cell_type": "markdown", "id": "11033666", "metadata": {}, "source": [ "\n", "\n", "## 学習成果\n", "\n", "このチュートリアルを修了すると、以下の内容を理解できるようになります：\n", "\n", "* 非弾性中性子散乱（INS）スペクトルと量子スピンモデルの動的構造因子（DSF）との関連性。\n", "* 量子回路において、基底状態を準備し、局所的な摂動を与え、トロッターの時間発展を行う方法。\n", "* 近似量子コンパイル（AQC）を用いて `qiskit-addon-aqc-tensor` 、ハードウェア実行向けのディープ・トロッター回路を圧縮する方法。\n", "* 量子ビットの期待値から遅延グリーン関数（RGF）を抽出し、フーリエ変換を行ってDSFに変換する方法。\n", "\n", "\n", "\n", "## 前提条件\n", "\n", "以下のトピックについて、あらかじめ理解しておくことをお勧めします：\n", "\n", "* [量子情報の基礎](/learning/courses/basics-of-quantum-information)\n", "* [変分アルゴリズム設計](/learning/courses/variational-algorithm-design)\n", "* [Qiskit primitives （EstimatorおよびSampler）の概要](/docs/guides/qiskit-runtime-primitives)\n", "\n", "\n", "\n", "## 背景\n", "\n", "このチュートリアルでは、以下の結果を再現します [Leeら、 arXiv:2603.15608](https://arxiv.org/abs/2603.15608).\n", "\n", "\n", "\n", "### 非弾性中性子散乱と動的構造因子\n", "\n", "非弾性中性子散乱（INS）は、量子材料における磁気励起を調べる上で最も強力な実験的手法の一つである。 熱中性子または冷中性子のビームが結晶に衝突すると、個々の中性子は、磁気サブシステムと運動量 $\\mathbf{q}$ およびエネルギー $\\omega$ の両方を交換する。 測定された散乱強度は、動的構造因子（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", "これは、スピン自由度の時空間全体の相関を記述するものである。\n", "\n", "\n", "\n", "### KCuF$_3$：標準的なルティンガー液体の磁性体\n", "\n", "フッ化銅カリウム（ KCuF$_3$ ）は、準一次元反強磁性体であり、スピン $\\frac{1}{2}$ Cu $^{2+}$ イオンからなる鎖が、最近接ハイゼンベルグ交換 $J$ を介して相互作用する一方、鎖間結合は $J$ の $\\sim 2.7\\%$ に過ぎない。 $T = 6\\;\\mathrm{K}$ において、INSデータが得られているが、そのスペクトルは、トモナガ・ルティンガー液体に特徴的な分画化されたスピノン励起によって支配されている。 等方点において、鎖内ダイナミクスは一次元スピン- $\\frac{1}{2}$ XXZハミルトニアンによってよく記述されるため（ $\\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$は、量子シミュレーションにとって理想的なベンチマークとなっている。そのハミルトニアンは量子プロセッサ上で実装できるほど単純である一方、基底状態は強くもつれ合っており、励起スペクトルには広い2スピノン連続体が現れる。\n", "\n", "注：このチュートリアルでは、エネルギー単位として $J = 1$ を設定し、 $H = J\\sum_i[\\ldots]$ という正規化を採用しています。これは、論文の回路実装（図）で使用されている局所ハミルトニアン $H_\\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ に対応しています。 S3 （補足資料より）。 本論文の完全ハミルトニアン（式 3) にはさらに2という総合係数が適用されるため、この論文の $J$ は、ここで使用されている $J$ の2倍となる。\n", "\n", "\n", "\n", "### シミュレーションおよび測定対象\n", "\n", "ここで計算する物理量は、 **遅延グリーン関数** （RGF）であり、これは時間依存のスピン間相関関数として定義される\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", "ここで、 $j_c$ は基準点（鎖の中心）であり、 $S_j^\\alpha(t) = e^{iHt}S_j^\\alpha e^{-iHt}$ はハイゼンベルク像のスピン演算子である。 このチュートリアルでは、 $zz$ コンポーネント（ $\\alpha = \\beta = z$ ）に焦点を当てます。量子コンピュータ上でRGFにアクセスするには、基底状態を準備し、 $j_c$ で局所的な摂動を与え、摂動を受けた状態を時間発展させ、各時間ステップごとに、すべてのサイト $j$ における単一量子ビットの期待値 $\\langle\\sigma_j^z\\rangle$ を測定します。 重要な点は、各 $\\langle\\sigma_j^z\\rangle$ が基底状態の磁化との差を表しているということである。等方性ハイゼンベルク反強磁性体は、サイトあたりの正味の磁化がゼロである（ $\\langle\\sigma_j^z\\rangle_\\mathrm{GS} = 0$ ）ため、測定された生データから、明示的な差し引きを行うことなく、直接 $G^R(j, j_c, t)$ が得られる。\n", "\n", "すべてのサイトおよび時間ステップにわたって $G^R(j, j_c, t)$ を収集することで、2次元データセットを構築し、これを空間および時間の両方でフーリエ変換することで、 **動的構造因子**$S(q,\\omega)$ を導出する。 DSFはINS実験で直接測定される量であり、各運動量 $q$ およびエネルギー $\\omega$ においてどのような磁気励起が存在するかを示す。等方性ハイゼンベルグ鎖の場合、正確な励起スペクトルは **2スピノン連続体**であり、その形状は量子シミュレーションに対する厳格なエンドツーエンドのベンチマークとして機能する。これにより、基底状態の準備、摂動、トロッター時間発展、および測定プロトコルのすべてを一度に検証することができる。\n", "\n", "\n", "\n", "### 量子シミュレーションのワークフロー\n", "\n", "このワークフローは、INSイベントの物理的挙動を反映しています。 我々は、（1） $n$ 個の量子ビット上で多体基底状態 $|\\psi_{\\mathrm{GS}}\\rangle$ を準備し、（2）中性子のスピン転移を模倣するために、鎖の中心で局所的なスピン反転摂動 $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ を加え、（3）2次トロッター化を用いて $H$ の下で離散的な時間ステップごとに系を進化させ、（4）各ステップで全量子ビットについて $\\langle\\sigma_i^z\\rangle$ を測定し、RGF を求める。 その後、2次元離散フーリエ変換を行うと、DSF $S(q,\\omega)$ が得られる。\n", "\n", "各時間ステップで測定する観測量は、各量子ビット $i$ における $\\sigma_i^z$ です。Qiskitでは、これは演算子の `SparsePauliOp` リストとして表現されます。具体的には、各サイトごとに、 $n$ のn量子ビットの恒等演算子列の中に、単一量子ビットの $Z$ が埋め込まれた形となります。 これらの観測量は、問題マッピング段階（ステップ1）において、問題インスタンスごとに1回構築され、そのスケールにおけるすべての回路で再利用される。\n", "\n", "\n", "\n", "### 近似量子コンパイル（AQC）\n", "\n", "ディープ・トロッター回路は、 **近似量子コンパイル（AQC）** によって圧縮することができる。AQCでは、最初の数層のトロッター層を、より短いパラメータ化されたアンザッツに置き換える。このアンザッツのパラメータは、元のディープ回路とのMPSレベルの忠実度を最大化するよう、古典的に最適化されている。 残りのトロッター手順をそのまま付加することで、2量子ビットゲートの数が大幅に少ない「混合型」のAQC＋トロッター回路が生成される。\n", "\n", "\n", "\n", "### MPSシミュレーション\n", "\n", "1次元系の場合、行列積状態（MPS）法を用いることで、基底状態の準備（密度行列再正規化群、DMRGを用いた）と回路レベルの時間発展の両方を効率的にシミュレーションすることができる。 結合の次元 $\\chi$ を制御することで、精度と計算コストのバランスを調整します。 このチュートリアルでは、MPSシミュレーションを用いて `qiskit-addon-aqc-tensor` 、ハードウェア実行のためにディープ・トロッター回路を圧縮する高精度なAQCアンザッツを計算します。\n", "\n" ] }, { "cell_type": "markdown", "id": "de1ad5a2", "metadata": {}, "source": [ "\n", "\n", "## 要件\n", "\n", "このチュートリアルを始める前に、以下のものがインストールされていることを確認してください：\n", "\n", "* Qiskit SDK [可視化](/docs/api/qiskit/visualization)機能を搭載\n", "* Qiskit Runtime (`pip install qiskit-ibm-runtime`)\n", "* `qiskit-addon-aqc-tensor` with `quimb` および JAX extras (`pip install 'qiskit-addon-aqc-tensor[quimb-jax]'`)\n", "\n" ] }, { "cell_type": "markdown", "id": "9938e4bd", "metadata": {}, "source": [ "\n", "\n", "## セットアップ\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", "## 小規模シミュレータの例\n", "\n", "まず、 **10キュービット**を用いてワークフロー全体を実証し、MPSシミュレーションを用いてHVA基底状態アンザッツを最適化し、Qiskitの状態ベクトルシミュレータを用いて時間発展を計算する。 DMRGは、基準となる基底状態エネルギーと基準となるMPSを算出する。 この小規模な事例により、規模を拡大する前に各段階を検証することができます。\n", "\n" ] }, { "cell_type": "markdown", "id": "3818ed8f", "metadata": {}, "source": [ "\n", "\n", "### ステップ1：古典的な入力を量子問題に写像する\n", "\n", "まず、物理モデルを定義し、量子回路を構築することから始めます。\n", "\n", "**ハミルトニアン。** KCuF$_3$ は、等方点（ $J = 1$、 $\\epsilon = 1$、エネルギー単位を $J = 1$ と設定）における 1D XXZ ハミルトニアンによってモデル化される。\n", "\n", "**基底状態。** `build_ground_state_ansatz`基底状態の準備回路として、ハミルトニアン変分アンザッツ（HVA）回路 を使用する。 `CircuitMPS`HVAパラメータは、DMRG基底状態MPSを用いて状態忠実度 $|\\langle\\psi_{\\mathrm{HVA}}(\\theta)|\\psi_{\\mathrm{DMRG}}\\rangle|^2$ を最大化する古典的な手法によって最適化される。この際、HVA状態は、quimbを用いて回路を行列積状態としてシミュレーションすることで評価される。 最適化には を使用 `scipy.optimize.minimize` します。 参考として、最適化されたアンザッツのエネルギー $\\langle H \\rangle$ も計算した。\n", "\n", "**トロッターゲート。** `PauliEvolutionGate(H_pair, time=time_step)`各最近傍相互作用項 $e^{-i\\Delta t\\, H_{\\mathrm{pair}}}$ （ $H_{\\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ ）は、次のように構成される。 Qiskitは、トランスパイラ処理の過程で、これを最適な3つのCNOTへの分解に合成します。\n", "\n", "**摂動。** 中央の量子ビットに $R_z(\\pi/2)$ ゲートを適用することで、 $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ が実現され、散乱中性子によって生じる局所的なスピン反転が再現される。\n", "\n", "**オブザーバブル。** 各量子ビットサイトについて、 $\\sigma^z$ の観測量を定義する。 これらの `SparsePauliOp` オブジェクトは、ステップ3でEstimatorプリミティブに渡され、各時間ステップごとに $\\langle\\sigma_i^z\\rangle$ が抽出されます。\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", "### ステップ 2：量子ハードウェアでの実行に向けて問題を最適化する\n", "\n", "実際のハードウェアの場合、上記のトロッター回路は深くなりすぎるだろう。 **近似量子コンパイル（AQC）** は、最初の $k$ トロッター層（基底状態回路を含む）を、元の深層回路とのMPSレベルでのフィデリティを最大化するよう最適化された、より短いパラメータ化アンザッツに置き換えることで、この問題に対処する。 残りのトロッターのステップをそのまま追加することで、より浅い「AQC + トロッター」回路が生成される。\n", "\n", "表現力と回路の深さのバランスをとるために、2つのアンザッツが用いられる。 **1層のアンザッツ** （1回のトロッターステップから生成される）は、最も初期の時間ステップを圧縮し、 **より**深い2層のアンザッツ（2回のトロッターステップから生成される）は、より高い精度が求められるその後の数ステップを圧縮する。\n", "\n", "AQCのワークフローには、4つのサブステップがあります：\n", "\n", "1. **ターゲット回路の構築** — ステップ1で作成した最初の $k_1 + k_2$ 回路は、そのままAQCのターゲットとして機能します。\n", "2. `quimb.tensor.CircuitMPS`**ターゲットMPSの計算** — 各ターゲット回路を、行列積状態としてシミュレートする。\n", "3. **アンザッツの生成と最適化** — `generate_ansatz_from_circuit` 1層および2層のパラメータ化されたアンザッツを作成します。パラメータは、 $1 - |\\langle\\psi_{\\mathrm{ansatz}}|\\psi_{\\mathrm{target}}\\rangle|^2$ を最小化するために、JAXによる勾配計算の高速化を適用したL-BFGS-Bを用いて最適化されます。 最初の $k_1$ ステップでは1層のアンザッツが使用され、次の $k_2$ ステップでは2層のアンザッツが使用されます。各ステージ内では、各ステップが前のステップで最適化されたパラメータからウォームスタートし、ステージの境界ではパラメータがステージのデフォルト値にリセットされます。\n", "4. `trotter_evolution`**混合回路を組み立てる** — AQCチェックポイント以降の時間ステップについては、. を使用して、最適化された2層のAQC回路に正確なトロッター層を追加する。\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", "### ステップ 3: `Qiskit primitives` を使用して実行する\n", "\n", "`StatevectorEstimator` AQCでコンパイルされた各回路を、プリミティブを用いてシミュレーションします。\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", "### ステップ4：後処理を行い、希望する従来の形式で結果を返す\n", "\n", "ここで、各時間ステップ $t$ におけるすべての量子ビット $i$ について、期待値 $\\langle\\sigma_i^z\\rangle$ を算出する。これらの値は、遅延グリーン関数行列 $G^R(j, j_c, t)$ を構成する。次に、この遅延グリーン関数（RGF）をフーリエ変換して動的構造因子 $S(q, \\omega)$ に変換し、RGF と DSF の両方をプロットする。 `get_dsf` ミラー対称性を適用し、負の値は内部で切り捨てます。\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", "## 大規模なハードウェア実行\n", "\n", "ここで\\*\\*、量子ビット数を\\*\\* 50まで増やします。 このスケールでは、最適化によるフィデリティは小規模な例よりも低くなります。基底状態のアンザッツのフィデリティはおよそ 0.65 まで低下し、AQCのフィデリティも後半のチェックポイントではおよそ 0.7 まで低下します。 これは予想される現象であり、追加の古典的計算コストを支払うことで、基底状態のアンザッツ層の数（`gs_n_layers`）や最適化反復回数（`maxiter`）を増やすことにより、これらの精度を向上させることができます。 また、2層のAQCの忠実度は、1層のものよりも低くなる点にも留意してください。 これは退行ではありません。後の時間ステップほどエンタングルメントが増加し、圧縮が単純に難しくなるため、それらには表現力の高い2層アンザッツが用いられているのです。\n", "\n", "AQCの最適化には、従来の計算時間でも数時間かかる場合があります（ここで示した実行例では約6時間かかり、その大部分はステップ 2c の2層のチェックポイントに費やされました）。 実処理時間を短縮するには、このノートブックを、高性能計算（HPC）システムなど、より高性能な従来のハードウェア上で実行することを検討してください。 あるいは、問題の規模を縮小（例えば、量子ビット数を減らしたり、時間ステップ数を減らしたり）することも可能です。その場合、得られる結果はここに示したものとは異なります。\n", "\n", "以下のコードは、小規模な例と同じ4段階の構成に従っています。 HVAの基底状態パラメータは、MPSシミュレーションを用いて再度最適化された。 QPU上では、エラーの抑制および軽減のために、動的デカップリング（DD）、パウリ・トゥワーリング、およびトゥワーリング付き読み出しエラー消去（TREX）を実現しています。 以下の表は、大規模実験と小規模実験の違いをまとめたものです：\n", "\n", "|   | 小規模 | 大規模な |\n", "| ------------------ | ---------------------- | ----------------------- |\n", "| 量子ビット | 10 | 50 |\n", "| 時間ステップ | 10 | 20 |\n", "| AQCチェックポイント（1層＋2層） | 3 ＋ 2 ＝ 5 | 6 ＋ 4 ＝ 10 |\n", "| 基底状態のアンザッツ層 | 3 | 5 |\n", "| MPSの最大結合寸法 | 32 | 128 |\n", "| 推定法 | `StatevectorEstimator` | DD搭載のQPU、パウリの旋回、そしてTREX |\n", "\n" ] }, { "cell_type": "markdown", "id": "9e85fbd1", "metadata": {}, "source": [ "\n", " AQCの最適化中（以下のステップ 2c ）に、次のようなメッセージが表示 `stderr` される場合があります：\n", "\n", " ```\n", " [Compiling module jit_MakeArrayFn for CPU] Very slow compile? ...\n", " The operation took 2m14s\n", " ```\n", "\n", " `qiskit-addon-aqc-tensor`これらは、JAXのautodiffを支えるコンパイラであるXLAによる、問題のない診断結果です。 50キュービット、MPSの結合次元が128の場合、XLAは勾配関数を初めてトレースする際に、コンパイルに数分かかる。 コンパイルは正常に完了し、最適化の結果にも影響はありません。\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": [ "ハードウェアによる結果からは、2スピノン連続体の主要な特徴が再現されている。すなわち、散乱強度は低エネルギー域において反強磁性波数 $q = \\pi$ 付近に集中しており、その下限は正弦波状のスピノン分散によって制限されている。また、その上方には単一の鋭いモードではなく、スペクトル重みの広い連続体が存在する。 これは、 KCuF$_3$ における非弾性中性子散乱によって測定されたものと同じ構造であり、50量子ビット規模における基底状態の準備、摂動、AQC圧縮トロッター進化、および誤差低減測定からなる一連のワークフローの有効性を裏付けるものである。\n", "\n" ] }, { "cell_type": "markdown", "id": "ls-nextsteps", "metadata": {}, "source": [ "\n", "\n", "## 次のステップ\n", "\n", "この作品に興味を持たれた方は、以下の資料もご参考になるかもしれません：\n", "\n", "\n", " * [Leeら、 「中性子散乱実験を用いた量子シミュレーションのベンチマーク」（ arXiv:2603.15608 ）](https://arxiv.org/abs/2603.15608) — 本チュートリアルが基にしている参考論文\n", " * [エラーの軽減および抑制手法](/docs/guides/error-mitigation-and-suppression-techniques) — ハードウェア実験で用いられたDD、パウリ・トゥワーリング、およびTREX\n", " * [時間発展回路のための近似量子コンパイル](/docs/tutorials/approximate-quantum-compilation-for-time-evolution) — AQC-Tensor チュートリアル\n", " * [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 }