{
"cells": [
{
"cell_type": "markdown",
"id": "8e82ead7",
"metadata": {},
"source": [
"---\n",
"title: \"Simular a dispersão de nêutrons em materiais quânticos com circuitos quânticos\"\n",
"description: \"Calcular o fator de estrutura dinâmico de um ímã quântico utilizando circuitos de Trotter e simulação MPS.\"\n",
"---\n",
"\n",
"{/* cspell:ignore DSFs spinon DMRG viridis fontsize vmax vmin cadetblue antiferromagnet COBYQA */}\n",
"\n",
"\n",
"\n",
"# Simular a dispersão de nêutrons em materiais quânticos com circuitos quânticos\n",
"\n",
"*Estimativa de tempo de execução: 13 minutos em um processador Heron r2 (NOTA: Trata-se apenas de uma estimativa. (O tempo de execução pode variar.)*\n",
"\n",
"\n",
" Use este tutorial para aprender a implementação passo a passo, com o cálculo clássico sendo executado localmente no seu laptop e a execução em hardware em uma QPU da IBM Quantum®. O exemplo em grande escala requer uma quantidade considerável de memória, e a compilação quântica aproximada (AQC) pode levar várias horas em um laptop comum. Para transferir a compressão para os recursos de computação e memóri Qiskit Serverless, utilize o [tutorial](/docs/tutorials/simulate-neutron-scattering-with-a-serverless-workflow) complementar sobre Serverless.\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "11033666",
"metadata": {},
"source": [
"\n",
"\n",
"## Resultados do aprendizado\n",
"\n",
"* Como os espectros de espalhamento inelástico de nêutrons (INS) se relacionam com os fatores de estrutura dinâmicos (DSFs) dos modelos de spin quântico.\n",
"* Como preparar um estado fundamental, aplicar uma perturbação local e realizar a evolução temporal de Trotter em um circuito quântico.\n",
"* Como usar a compilação quântica aproximada (AQC) com `qiskit-addon-aqc-tensor` para compactar circuitos de Trotter profundos para execução em hardware.\n",
"* Como extrair a função de Green retardada (RGF) a partir dos valores esperados dos qubits e transformá-la por meio da transformada de Fourier em uma DSF.\n",
"\n",
"\n",
"\n",
"## Pré-requisitos\n",
"\n",
"* [Noções básicas de informação quântica](/learning/courses/basics-of-quantum-information)\n",
"* [Design de algoritmos variacionais](/learning/courses/variational-algorithm-design)\n",
"* [Introdução ao “ Qiskit primitives ” (Estimador e Amostrador)](/docs/guides/qiskit-runtime-primitives)\n",
"\n",
"\n",
"\n",
"## Segundo plano\n",
"\n",
"Neste tutorial, reproduzimos os resultados de [Lee et al., arXiv:2603.15608](https://arxiv.org/abs/2603.15608).\n",
"\n",
"A dispersão de nêutrons ajuda os pesquisadores a caracterizar materiais relevantes para a sustentabilidade, incluindo eletrodos de baterias. Este tutorial explora como os circuitos quânticos podem simular excitações magnéticas e relacionar modelos teóricos às medições de espalhamento. Seu exemplo de ímã quântico desenvolve métodos para o estudo do comportamento dos materiais, contribuindo para o conjunto mais amplo de ferramentas de pesquisa voltadas para a descoberta de novos materiais.\n",
"\n",
"\n",
"\n",
"### Espalhamento inelástico de nêutrons e o fator de estrutura dinâmico\n",
"\n",
"A dispersão inelástica de nêutrons (INS) é uma das técnicas experimentais mais poderosas para o estudo das excitações magnéticas em materiais quânticos. Quando um feixe de nêutrons térmicos ou frios incide sobre um cristal, os nêutrons individuais trocam tanto o momento $\\mathbf{q}$ quanto a energia $\\omega$ com o subsistema magnético. A intensidade de espalhamento medida é proporcional ao fator de estrutura dinâmica (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",
"que codifica todas as correlações espaço-temporais dos graus de liberdade de spin.\n",
"\n",
"\n",
"\n",
"### KCuF$_3$: um ímã canônico do tipo líquido de Luttinger\n",
"\n",
"O fluoreto de potássio e cobre ( KCuF$_3$ ) é um antiferromagneto quase unidimensional no qual cadeias de íons spin- $\\frac{1}{2}$ Cu $^{2+}$ interagem por meio de um $J$ de troca de Heisenberg entre vizinhos mais próximos, enquanto o acoplamento entre cadeias é apenas $\\sim 2.7\\%$ de $J$. Em $T = 6\\;\\mathrm{K}$, onde há dados de INS disponíveis, o espectro é dominado por excitações de spinons fracionadas, características de um líquido de Tomonaga-Luttinger. Como a dinâmica intracadeia é bem descrita pelo hamiltoniano unidimensional de spin- $\\frac{1}{2}$ XXZ no ponto isotrópico ( $\\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$ serve como um modelo de referência ideal para simulação quântica: o hamiltoniano é simples o suficiente para ser implementado em um processador quântico, mas o estado fundamental apresenta forte entrelaçamento e o espectro de excitação apresenta um amplo continuum de dois spinons.\n",
"\n",
"Observação: este tutorial define $J = 1$ como unidade de energia e adota a normalização $H = J\\sum_i[\\ldots]$, que corresponde ao hamiltoniano local $H_\\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ utilizado na implementação do circuito apresentada no artigo (Fig. S3 (do suplemento). O hamiltoniano completo do artigo (Eq. 3) apresenta um fator geral adicional de 2; portanto, o valor de “ $J$ ” do artigo é o dobro do valor de “ $J$ ” utilizado aqui.\n",
"\n",
"\n",
"\n",
"### O que simulamos e medimos\n",
"\n",
"A grandeza física que calculamos é a **função de Green retardada** (RGF), definida como a função de correlação spin-spin dependente do tempo\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",
"onde $j_c$ é um ponto de referência (o centro da cadeia) e $S_j^\\alpha(t) = e^{iHt}S_j^\\alpha e^{-iHt}$ é o operador de spin na representação de Heisenberg. Neste tutorial, vamos nos concentrar no componente “ $zz$ ” ( $\\alpha = \\beta = z$ ). Em um computador quântico, o RGF é acessado preparando-se o estado fundamental, aplicando-se uma perturbação local em $j_c$, evoluindo-se o estado perturbado no tempo e medindo-se o valor esperado do qubit único $\\langle\\sigma_j^z\\rangle$ em cada local $j$ para cada passo temporal. A ideia principal é que cada $\\langle\\sigma_j^z\\rangle$ representa a diferença em relação à magnetização do estado fundamental; como o antiferromagneto isotrópico de Heisenberg possui magnetização líquida por vazio igual a zero ( $\\langle\\sigma_j^z\\rangle_\\mathrm{GS} = 0$ ), o valor bruto medido resulta diretamente em $G^R(j, j_c, t)$ sem qualquer subtração explícita.\n",
"\n",
"Ao coletar o $G^R(j, j_c, t)$ em todos os locais e intervalos de tempo, montamos um conjunto de dados bidimensional que é então submetido à transformada de Fourier tanto no espaço quanto no tempo para produzir o **fator de estrutura dinâmica** $S(q,\\omega)$. O DSF é a grandeza medida diretamente em um experimento de INS: ele nos informa quais excitações magnéticas existem em cada momento $q$ e energia $\\omega$. Para a cadeia de Heisenberg isotrópica, o espectro de excitação exato é um **contínuo de dois spinons**, uma ampla faixa de intensidade de espalhamento cuja forma serve como um rigoroso parâmetro de referência de ponta a ponta para a simulação quântica: ela valida de uma só vez a preparação do estado fundamental, a perturbação, a evolução temporal de Trotter e o protocolo de medição.\n",
"\n",
"\n",
"\n",
"### Fluxo de trabalho de simulação quântica\n",
"\n",
"O fluxo de trabalho reflete a física de um evento INS. Nós (1) preparamos o estado fundamental de muitos corpos $|\\psi_{\\mathrm{GS}}\\rangle$ em $n$ qubits, (2) aplicamos uma perturbação local de inversão de spin $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ no centro da cadeia para simular a transferência de spin do nêutron, (3) evoluímos sob $H$ em passos de tempo discretos usando a trotterização de segunda ordem e (4) medimos $\\langle\\sigma_i^z\\rangle$ em cada qubit a cada passo para obter o RGF. A transformada discreta de Fourier bidimensional fornece, então, o DSF $S(q,\\omega)$.\n",
"\n",
"A observável que medimos a cada passo temporal é $\\sigma_i^z$ em cada qubit $i$. No Qiskit, isso é representado como uma lista de `SparsePauliOp` operadores: um operador de qubit único $Z$ incorporado na sequência de identidade de $n$ -qubit para cada local. Essas variáveis observáveis são construídas uma vez por instância do problema durante a etapa de mapeamento do problema (Etapa 1) e reutilizadas para todos os circuitos nessa escala.\n",
"\n",
"\n",
"\n",
"### Compilação quântica aproximada (AQC)\n",
"\n",
"Os circuitos Deep Trotter podem ser compactados por meio **da compilação quântica aproximada (AQC)**, que substitui as primeiras camadas do Trotter por um ansatz parametrizado mais curto, cujos parâmetros são otimizados classicamente para maximizar a fidelidade no nível do MPS em relação ao circuito profundo original. Os passos restantes do algoritmo de Trotter são anexados exatamente, resultando em um circuito “misto” de AQC + Trotter com um número substancialmente menor de portas de dois qubits.\n",
"\n",
"\n",
"\n",
"### Simulação MPS\n",
"\n",
"Para um sistema unidimensional, os métodos de estado de produto matricial (MPS) podem simular com eficiência tanto a preparação do estado fundamental (com o grupo de renormalização da matriz de densidade, ou DMRG) quanto a evolução temporal no nível do circuito. Ao controlar a dimensão da ligação $\\chi$, fazemos um equilíbrio entre precisão e custo computacional. Neste tutorial, utilizamos a simulação MPS para `qiskit-addon-aqc-tensor` calcular ansätze AQC de alta fidelidade que compactam os circuitos de Trotter profundos para execução em hardware.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "de1ad5a2",
"metadata": {},
"source": [
"\n",
"\n",
"## Requisitos\n",
"\n",
"Antes de iniciar este tutorial, certifique-se de ter os seguintes itens instalados:\n",
"\n",
"* Qiskit SDK com suporte [à visualização](/docs/api/qiskit/visualization)\n",
"* Qiskit Runtime (`pip install qiskit-ibm-runtime`)\n",
"* `qiskit-addon-aqc-tensor` com `quimb` e os recursos adicionais do JAX (`pip install 'qiskit-addon-aqc-tensor[quimb-jax]'`)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9938e4bd",
"metadata": {},
"source": [
"\n",
"\n",
"## Instalação\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",
"## Exemplo de simulador em pequena escala\n",
"\n",
"Primeiramente, demonstramos o fluxo de trabalho completo em **10 qubits**, otimizando o ansatz do estado fundamental HVA por meio da simulação MPS e utilizando o simulador de vetor de estado do Qiskit para a evolução temporal. O DMRG fornece uma energia de referência do estado fundamental e um MPS de referência. Esse exemplo em pequena escala nos permite validar cada etapa antes de ampliar a escala.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "3818ed8f",
"metadata": {},
"source": [
"\n",
"\n",
"### Etapa 1: Mapeie entradas clássicas para um problema quântico\n",
"\n",
"Começamos definindo o modelo físico e construindo os circuitos quânticos.\n",
"\n",
"**Hamiltoniano.** KCuF$_3$ é modelado pelo hamiltoniano XXZ de 1D no ponto isotrópico ( $J = 1$, $\\epsilon = 1$, definindo $J = 1$ como unidade de energia).\n",
"\n",
"**Estado fundamental.** Utilizamos um circuito de `build_ground_state_ansatz` ansatz variacional hamiltoniano (HVA) como circuito de preparação do estado fundamental. `CircuitMPS`Os parâmetros do HVA são otimizados de forma clássica, maximizando a fidelidade do estado $|\\langle\\psi_{\\mathrm{HVA}}(\\theta)|\\psi_{\\mathrm{DMRG}}\\rangle|^2$ com o MPS do estado fundamental do DMRG, em que o estado do HVA é avaliado por meio da simulação do circuito como um estado de produto matricial com o quimb. Utilizamos `scipy.optimize.minimize` para a otimização. A energia $\\langle H \\rangle$ do ansatz otimizado também é calculada para fins de referência.\n",
"\n",
"**Portões para trotadores.** Cada termo de interação de vizinho mais próximo $e^{-i\\Delta t\\, H_{\\mathrm{pair}}}$ com $H_{\\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ é construído com `PauliEvolutionGate(H_pair, time=time_step)`. O Qiskit sintetiza isso na decomposição ótima de três CNOTs durante a transpilagem.\n",
"\n",
"**Perturbação.** Um portão “ $R_z(\\pi/2)$ ” no qubit central implementa um “ $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ ”, imitando a inversão de spin local produzida por um nêutron espalhado.\n",
"\n",
"**Observáveis.** Construímos um observável de tipo “ $\\sigma^z$ ” para cada local de qubit. Esses `SparsePauliOp` objetos são passados para a primitiva Estimator na Etapa 3, a fim de extrair $\\langle\\sigma_i^z\\rangle$ em cada passo temporal.\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",
"### Etapa 2: Otimizar o problema para execução em hardware quântico\n",
"\n",
"Para um hardware real, os circuitos de Trotter acima seriam muito complexos. **A compilação quântica aproximada (AQC)** resolve essa questão substituindo as primeiras camadas de Trotter de $k$ (incluindo o circuito do estado fundamental) por um ansatz parametrizado mais curto, otimizado para maximizar a fidelidade no nível do MPS em relação ao circuito profundo original. Os passos restantes do algoritmo de Trotter são anexados exatamente, resultando em um circuito “AQC + Trotter” menos complexo.\n",
"\n",
"Para equilibrar a expressividade com a profundidade do circuito, são utilizadas duas abordagens: uma **abordagem de uma camada** (gerada a partir de um único passo de Trotter) comprime os primeiros passos temporais, e uma **abordagem** mais profunda de duas camadas (gerada a partir de dois passos de Trotter) comprime os passos seguintes, nos quais é necessária maior fidelidade.\n",
"\n",
"O fluxo de trabalho do AQC possui quatro subetapas:\n",
"\n",
"1. **Construa os circuitos-alvo** — os primeiros circuitos “ $k_1 + k_2$ ” da Etapa 1 servem diretamente como alvos para o AQC.\n",
"2. **Calcular o MPS de destino** — simular cada circuito de destino como um estado de produto matricial usando `quimb.tensor.CircuitMPS`.\n",
"3. **Gerar e otimizar ansätze** — `generate_ansatz_from_circuit` cria um ansatz parametrizado de uma camada e outro de duas camadas; os parâmetros são otimizados utilizando o algoritmo L-BFGS-B com gradientes acelerados por JAX para minimizar $1 - |\\langle\\psi_{\\mathrm{ansatz}}|\\psi_{\\mathrm{target}}\\rangle|^2$. Os primeiros passos $k_1$ utilizam o ansatz de uma camada e os passos seguintes $k_2$ utilizam o ansatz de duas camadas; dentro de cada estágio, cada passo inicia a partir dos parâmetros otimizados do passo anterior, e os parâmetros são redefinidos para os padrões do estágio no limite do estágio.\n",
"4. **Montar circuitos mistos** — para intervalos de tempo além do ponto de verificação do AQC, acrescentar camadas exatas de Trotter ao circuito AQC otimizado de duas camadas usando `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": [
""
]
},
"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",
"### Etapa 3: Executar usando Qiskit primitives\n",
"\n",
"Simulamos cada circuito compilado pelo AQC utilizando a `StatevectorEstimator` primitiva.\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",
"### Etapa 4: Realizar o pós-processamento e apresentar o resultado no formato clássico desejado\n",
"\n",
"Agora, extraímos o valor esperado $\\langle\\sigma_i^z\\rangle$ para cada qubit $i$ em cada passo temporal $t$. Esses valores formam a matriz da função de Green retardada $G^R(j, j_c, t)$. Em seguida, aplicamos a transformada de Fourier à RGF para obter o fator de estrutura dinâmica $S(q, \\omega)$ e representamos graficamente tanto a RGF quanto o DSF. `get_dsf` aplica a simetria espelhada e limita internamente os valores negativos.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "4585d139",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"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",
"## Execução em hardware em grande escala\n",
"\n",
"Agora ampliamos para **50 qubits**. Nessa escala, as otimizações atingem fidelidades menores do que no exemplo em pequena escala: a fidelidade do ansatz do estado fundamental cai para cerca de 0.65, e as fidelidades do AQC diminuem para cerca de 0.7 nos pontos de verificação posteriores. Isso é esperado, e é possível melhorar essas precisões aumentando o número de camadas do ansatz do estado fundamental (`gs_n_layers`) ou o número de iterações de otimização (`maxiter`) com um custo clássico adicional. Observe também que as fidelidades do AQC de duas camadas são menores do que as do AQC de uma camada. Isso não é uma regressão: os passos temporais posteriores geram mais entrelaçamento e são simplesmente mais difíceis de comprimir, razão pela qual se utiliza para eles o ansatz de duas camadas, que é mais expressivo.\n",
"\n",
"A otimização do AQC também pode levar várias horas de tempo de computação clássico (aproximadamente seis horas na execução mostrada aqui, sendo que a maior parte desse tempo é dedicada aos pontos de verificação de duas camadas na Etapa 2c ). Para reduzir o tempo real, considere executar este notebook em um hardware clássico mais potente, como um sistema de computação de alto desempenho (HPC). Como alternativa, você pode reduzir a escala para uma instância menor do problema (por exemplo, com menos qubits ou menos passos temporais); nesse caso, seus resultados serão diferentes dos apresentados aqui.\n",
"\n",
"O código abaixo segue a mesma estrutura de quatro etapas do exemplo em pequena escala. Os parâmetros do estado fundamental do HVA são novamente otimizados por meio da simulação MPS. Na QPU, ativamos o desacoplamento dinâmico (DD), o “Pauli twirling” e a supressão de erros por leitura com “twirling” (TREX) para suprimir e mitigar erros. A tabela a seguir resume as diferenças entre o experimento em grande escala e o de pequena escala:\n",
"\n",
"| | Pequena escala | Em grande escala |\n",
"| --------------------------------------------------- | ---------------------- | -------------------------------- |\n",
"| Qubits | 22 | 50 |\n",
"| Intervalos de tempo | 22 | 20 |\n",
"| Pontos de verificação de AQC (1 camada + 2 camadas) | 3 + 2 = 5 | 6 + 4 = 10 |\n",
"| Camadas do ansatz do estado fundamental | 3 | 5 |\n",
"| Dimensão máxima da ligação do MPS | 32 | 128 |\n",
"| Orçador | `StatevectorEstimator` | QPU com DD, giro de Pauli e TREX |\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9e85fbd1",
"metadata": {},
"source": [
"\n",
" Durante a otimização do AQC (Etapa 2c abaixo), você poderá ver `stderr` mensagens como:\n",
"\n",
" ```\n",
" [Compiling module jit_MakeArrayFn for CPU] Very slow compile? ...\n",
" The operation took 2m14s\n",
" ```\n",
"\n",
" Esses são diagnósticos benignos do XLA, o compilador que dá suporte `qiskit-addon-aqc-tensor`ao autodiff do JAX. Com 50 qubits e dimensão de ligação MPS igual a 128, o XLA leva alguns minutos para compilar a função de gradiente na primeira vez em que ela é traçada. A compilação é bem-sucedida, e os resultados da otimização não são afetados.\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": [
""
]
},
"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": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"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": [
"Os resultados do hardware reproduzem as principais características do continuum de dois spinons: a intensidade de espalhamento concentra-se próximo ao vetor de onda antiferromagnético $q = \\pi$ em baixas energias e é limitada por baixo pela dispersão sinusoidal dos spinons, com um amplo continuum de peso espectral acima dela, em vez de um único modo bem definido. Essa é a mesma estrutura medida por espalhamento inelástico de nêutrons em um KCuF$_3$, e ela valida todo o fluxo de trabalho — desde a preparação do estado fundamental, passando pela perturbação, pela evolução de Trotter comprimida por AQC, até a medição com mitigação de erros — na escala de 50 qubits.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ls-nextsteps",
"metadata": {},
"source": [
"\n",
"\n",
"## Próximas etapas\n",
"\n",
"Se você achou este trabalho interessante, talvez se interesse pelo material a seguir:\n",
"\n",
"\n",
" * [Simulação da dispersão de nêutrons com um fluxo de trabalho sem servidor (serverless) baseado em AQC + dinâmica de Trotter](/docs/tutorials/simulate-neutron-scattering-with-a-serverless-workflow) — um tutorial complementar que executa esse mesmo experimento por meio de um modelo de função do Qiskit Serverless já implantado\n",
" * [Lee et al., \"Avaliação comparativa da simulação quântica com experimentos de espalhamento de nêutrons\" ( arXiv:2603.15608 )](https://arxiv.org/abs/2603.15608) — o artigo de referência no qual este tutorial se baseia\n",
" * [Técnicas de mitigação e supressão de erros](/docs/guides/error-mitigation-and-suppression-techniques) — DD, Pauli twirling e TREX utilizadas nos experimentos com hardware\n",
" * [Compilação quântica aproximada para circuitos de evolução temporal](/docs/tutorials/approximate-quantum-compilation-for-time-evolution) — tutorial sobre o AQC-Tensor\n",
" * [Documentação do 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
}