นิวเคลียร์แมกเนติกเรโซแนนซ์ด้วย Qiskit

สำหรับโมดูล Qiskit in Classrooms นี้ นักเรียนต้องมีสภาพแวดล้อม Python ที่ใช้งานได้ พร้อมติดตั้งแพ็กเกจต่อไปนี้

• qiskit v2.1 หรือใหม่กว่า
• qiskit-ibm-runtime v0.40.1 หรือใหม่กว่า
• qiskit-aer v0.17 หรือใหม่กว่า
• scipy และ matplotlib (สำหรับการฟิตและการพลอต)
• pylatexenc (สำหรับเรนเดอร์ไดอะแกรม Circuit)

หากต้องการตั้งค่าและติดตั้งแพ็กเกจข้างต้น ดูคู่มือ Install Qiskit หากต้องการรันงานบนคอมพิวเตอร์ควอนตัมจริง นักเรียนต้องตั้งค่าบัญชีกับ IBM Quantum® ตามขั้นตอนในคู่มือ Set up your IBM Cloud® account

โมดูลนี้ผ่านการทดสอบแล้วและใช้เวลา QPU ประมาณ 5 นาที นี่เป็นเพียงค่าประมาณ การใช้งานจริงของคุณอาจต่างออกไป

บทนำ

เมื่อนักเคมีทำการทดลอง NMR พวกเขาวางตัวอย่างสปินของนิวเคลียสในสนามแม่เหล็กแรงสูง ส่งพัลส์ความถี่วิทยุ (RF) แล้ววัดสัญญาณที่ตอบกลับมา เทคนิคนี้ตั้งอยู่บนแนวคิดไม่กี่ข้อจากกลศาสตร์ควอนตัมของอนุภาคสปิน 1/2 ในสนามแม่เหล็ก สปินมีระดับพลังงานสองระดับที่ห่างกันด้วยพลังงานที่สอดคล้องกับความถี่ Larmor เฉพาะตัว พัลส์เรโซแนนซ์ที่ความถี่นั้นพลิกสปินไปยังระนาบขวาง ซึ่งสปินจะหมุนควง (precess) และเหนี่ยวนำสัญญาณที่จากนั้นค่อยๆ จางลงด้วยเวลาลักษณะเฉพาะสองค่า คือ $T_1$ และ $T_2$

qubit แบบตัวนำยิ่งยวดคือเวอร์ชันที่ออกแบบขึ้นของวัตถุเดียวกัน คือระบบควอนตัมสองระดับที่ควบคุมได้ สถานะ $|0\rangle$ และ $|1\rangle$ เปรียบได้กับสปินขึ้นและสปินลง Gate ที่เราใช้หมุนมันเปรียบได้กับพัลส์ RF และมันสูญเสียความต่อเนื่องเชิงควอนตัมด้วย $T_1$ และ $T_2$ ของตัวเอง นี่ไม่ใช่แค่การเปรียบเทียบหลวมๆ โปรเซสเซอร์ควอนตัมเชิงทดลองยุคแรกๆ บางตัวในช่วงปลายทศวรรษ 1990 และต้นทศวรรษ 2000 ใช้ NMR

ในโมดูลนี้ คุณจะใช้คอมพิวเตอร์ควอนตัมจริงรันการทดลองระบุลักษณะแบบเดียวกับที่นักสเปกโทรสโกปี NMR (หรือวิศวกรควอนตัมที่สอบเทียบฮาร์ดแวร์) ทำ

• Rabi / nutation — ขับสปินแล้วดูมันหมุน นี่คือวิธีสอบเทียบการหมุน $\pi$ หรือ $\pi/2$

• $T_1$ — การคลายตัวแบบสปิน-แลตทิซ คือการที่สปินกระตุ้นสลายกลับสู่สมดุล

• $T_2$ — การคลายตัวแบบสปิน-สปิน (dephasing) วัดด้วยการทดลอง Ramsey และรีโฟกัสด้วย Hahn echo

• สเปกโทรสโกปี NMR (chemical shift และ $J$-coupling) — จำลองโมเลกุลขนาดเล็กบน qubit หลายตัวและอ่านสเปกตรัม multiplet เพื่อกู้คืนการเชื่อมต่อผ่านกฎ $n+1$

สำหรับแต่ละแนวคิด เราให้ทฤษฎีขั้นต่ำก่อน แล้วจึงสร้างการทดลองใน Qiskit ตามเวิร์กโฟลว์ Qiskit patterns สี่ขั้นตอน

ทฤษฎี

อนุภาคสปิน 1/2 ในสนามแม่เหล็ก

นิวเคลียสสปิน 1/2 มีโมเมนต์แม่เหล็ก $\boldsymbol{\mu} = \gamma \hbar \mathbf{I}$ โดยที่ $\gamma$ คืออัตราส่วนไจโรแมกเนติก ในสนามสถิต $B_0$ ตามแกน $z$ สถานะสปินทั้งสองแยกออกจากกันในพลังงาน (ปรากฏการณ์ Zeeman)

$E_{\pm} = \mp \tfrac{1}{2}\gamma \hbar B_0, \qquad \Delta E = \gamma \hbar B_0 = \hbar \omega_0 .$

ช่องว่างพลังงานนิยามความถี่ Larmor (เรโซแนนซ์) $\omega_0 = \gamma B_0$ โมเมนต์แม่เหล็กหมุนควงรอบ $B_0$ ที่ความถี่เดียวกันนี้

เพื่อหมุนสปิน เราใช้สนาม RF ที่แกว่งที่ (หรือใกล้) $\omega_0$ ในเฟรมที่หมุนตามความถี่การขับ (เฟรมหมุน) การขับที่เรโซแนนซ์จะหมุนเวกเตอร์ Bloch รอบแกนในระนาบ $x$-$y$ มุมที่หมุนไป (มุมพลิก) ถูกกำหนดโดยความแรงและระยะเวลาของพัลส์ สองกรณีที่สำคัญที่สุดคือ พัลส์ $90^\circ$ พลิกสปิน $+1/2$ ไปยังระนาบขวาง และพัลส์ $180^\circ$ หมุนมันไปที่ $-1/2$ เต็มที่

การแปลระหว่างศัพท์ NMR กับศัพท์คอมพิวเตอร์ควอนตัม

qubit คือระบบสองระดับ ซึ่งในทางคณิตศาสตร์เหมือนกับสปิน 1/2 ในสนามแม่เหล็กทุกประการ — ดังนั้นแม้เราจะรันบนโปรเซสเซอร์ควอนตัมแทนที่จะเป็นนิวเคลียส ฟิสิกส์ก็เป็นฟิสิกส์เดียวกับที่เป็นพื้นฐานของ NMR ต่างกันแค่คำศัพท์ ตารางนี้ช่วยให้คุณคุยกับทั้งนักเคมีกายภาพและผู้เชี่ยวชาญคอมพิวเตอร์ควอนตัมได้

NMR / สปิน 1/2	โปรเซสเซอร์ควอนตัม (qubit)
สปินขึ้น / สปินลง	$\ket{0}$ / $\ket{1}$
ช่องว่าง Zeeman $\Delta E = \hbar\omega_0$	พลังงานการเปลี่ยนสถานะของ qubit $\hbar\omega_{01}$
ความถี่ Larmor $\omega_0$	ความถี่ของ qubit $\omega_{01}$
พัลส์ RF เรโซแนนซ์ (มุมพลิก $\theta$)	Gate หมุน qubit เดี่ยว $R_X(\theta)$
พัลส์ $90^\circ$ / $180^\circ$	$R_X(\pi/2)$ / $R_X(\pi)$
การคลายตัวแบบสปิน-แลตทิซ $T_1$	เวลาคลายพลังงานของ qubit $T_1$
การคลายตัวแบบสปิน-สปิน $T_2$	เวลา dephasing ของ qubit $T_2$
Nutation	การแกว่งแบบ Rabi
Free induction decay (FID), Hahn echo	ลำดับ Ramsey, echo

ทุกอย่างต่อจากนี้คือแบบฝึกหัดในการวัดคอลัมน์ขวามือ แล้วอ่านมันเป็นคอลัมน์ซ้ายมือ

ตั้งค่าสภาพแวดล้อม Qiskit

เราใช้เวิร์กโฟลว์ Qiskit patterns ตลอด

1. แมปปัญหาเป็น Circuit ควอนตัม

2. ปรับให้เหมาะสม Circuit สำหรับการรัน (ทรานสไพล์ไปยัง Backend)

3. รันด้วย primitive ของ IBM Quantum

4. ประมวลผลภายหลังและวิเคราะห์

ก่อนอื่นเราโหลดบริการ runtime และเลือก Backend จริงที่คิวว่างที่สุด หากคุณยังไม่ได้บันทึกข้อมูลรับรอง ดูหน้าเอกสารเรื่องวิธีบันทึกข้อมูลเข้าสู่ระบบของคุณ

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

# Load the IBM Quantum Compute Service and select the least-busy real backend
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler, Batch

# After setup, this picks up your saved credentials automatically:
service = QiskitRuntimeService()

backend = service.least_busy(operational=True, simulator=False, min_num_qubits=156)
print(f"Selected backend: {backend.name}")

Selected backend: ibm_kingston

ใช้เซลล์ด้านล่างเพื่อสลับระหว่างฮาร์ดแวร์จริงกับตัวจำลองในเครื่องสำหรับทั้งโน้ตบุ๊ก การรันบนฮาร์ดแวร์คือจุดมุ่งหมายของโมดูลนี้ แต่หากมีปัญหาด้านการเชื่อมต่อหรือการใช้งาน ตัวจำลองใช้เป็นตัวสำรองได้ ตั้งค่า USE_HARDWARE = True เพื่อรันบน QPU ที่เลือกไว้ข้างต้น

import numpy as np
from itertools import combinations
from qiskit_aer import AerSimulator
from qiskit.transpiler import generate_preset_pass_manager
from qiskit_ibm_runtime import EstimatorV2 as Estimator

USE_HARDWARE = True  # <-- set True to run on the QPU selected above
USE_BATCH = False  # Safer for a class with pauses; set True for an uninterrupted run

if USE_HARDWARE:
    run_backend = backend  # the real QPU selected above
else:
    # Local testing mode: an Aer simulator that mimics the selected device -- its noise model,
    # coupling map, and qubit count -- while layout choices still use real calibration data.
    run_backend = AerSimulator.from_backend(backend)

# --- Choose connected qubits with low readout error. --------------------------------------
def readout_error(q):
    props = backend.target["measure"].get((q,))
    return props.error if props is not None and props.error is not None else 1.0

coupling_map = backend.coupling_map or backend.target.build_coupling_map()
edges = {tuple(sorted(edge)) for edge in coupling_map.get_edges() if edge[0] != edge[1]}

# AX needs one connected pair.
AX_LAYOUT = list(min(edges, key=lambda edge: sum(readout_error(q) for q in edge)))

# AX2 needs one qubit connected to two neighbours; the first entry is the central A qubit.
neighbors = {q: set() for q in range(backend.num_qubits)}
for q0, q1 in edges:
    neighbors[q0].add(q1)
    neighbors[q1].add(q0)
ax2_layouts = (
    [a, x1, x2]
    for a, adjacent in neighbors.items()
    for x1, x2 in combinations(adjacent, 2)
)
AX2_LAYOUT = min(ax2_layouts, key=lambda layout: sum(readout_error(q) for q in layout))

# Use a long-lived qubit for the single-qubit experiments.
def best_qubit(backend):
    return max(
        range(backend.num_qubits), key=lambda q: backend.qubit_properties(q).t2 or 0.0
    )

QUBIT = best_qubit(backend)
qp = backend.qubit_properties(QUBIT)
print(
    f"Using qubit {QUBIT} for single-qubit experiments: "
    f"T1 = {qp.t1 * 1e6:.0f} us, T2 = {qp.t2 * 1e6:.0f} us"
)
print(f"Using physical qubits {AX_LAYOUT} for AX")
print(f"Using physical qubits {AX2_LAYOUT} for AX2 (A first, then the two X spins)")

pm_q = generate_preset_pass_manager(
    optimization_level=1, backend=run_backend, initial_layout=[QUBIT]
)
pm_ax = generate_preset_pass_manager(
    optimization_level=1, backend=run_backend, initial_layout=AX_LAYOUT
)
pm_ax2 = generate_preset_pass_manager(
    optimization_level=1, backend=run_backend, initial_layout=AX2_LAYOUT
)

# Job mode is the default because a classroom pause can deactivate a Batch. On the Open plan,
# a Batch has a 10-minute maximum lifetime and deactivates after 1 minute with no queued job.
batch = Batch(backend=run_backend) if USE_BATCH else None
execution_mode = batch if batch is not None else run_backend
sampler = Sampler(mode=execution_mode)
estimator = Estimator(mode=execution_mode)
estimator.options.resilience_level = 0  # raw expectation values, no error mitigation
estimator.options.default_shots = 1000

mode_name = "Batch" if USE_BATCH else "individual jobs"
device_name = "hardware" if USE_HARDWARE else "local simulator (runtime local mode)"
print(f"Ready. Running on {device_name} using {mode_name}.")

การทดลอง: การแกว่งแบบ Rabi (nutation)

ทฤษฎี

ขับสปินที่เรโซแนนซ์ แล้วเวกเตอร์ Bloch ของมันจะหมุนสม่ำเสมอจากขั้วหนึ่งไปอีกขั้ว — มันนูเทต ถ้าเราใช้การหมุนด้วยมุมพลิก $\theta$ รอบแกน $x$ กับ qubit ที่เริ่มที่ $|0\rangle$ ความน่าจะเป็นที่จะพบมันใน $|1\rangle$ คือ

$P_1(\theta) = \sin^2\!\left(\tfrac{\theta}{2}\right) = \tfrac{1}{2}\bigl(1 - \cos\theta\bigr).$

นี่คือเส้นโค้ง Rabi ค่าสูงสุดแรกคือพัลส์ $\pi$ (การกลับสถานะเต็มที่ คือพัลส์ $180^\circ$) ส่วนพัลส์ $\pi/2$ คือการหมุนหนึ่งในสี่รอบ (พัลส์ $90^\circ$ ที่พลิกสปินไปยังระนาบขวาง) ใน NMR เส้นโค้งนี้คือวิธีที่คุณสอบเทียบพัลส์ของคุณ

ข้อควรระวังของการทดลอง

การทดลอง Rabi ตามตัวอักษร กวาดแอมพลิจูดหรือระยะเวลาของการขับตามเวลา ซึ่งต้องใช้การควบคุมระดับพัลส์ ที่ถูกถอดออกใน Qiskit v2.0 ที่นี่เราจึงกวาดมุมหมุน $\theta$ ของ Gate $R_X(\theta)$ แทน เส้นโค้งและการระบุ $\pi$/$\pi/2$ เหมือนกันทุกประการ ต่างกันแค่เราเลือกมุมโดยตรงแทนการสอบเทียบความยาวพัลส์ บนฮาร์ดแวร์ที่มี fractional gate $R_X(\theta)$ สำหรับ $\theta$ ใดๆ จะรันแบบเนทีฟ

มุมมองทรงกลม Bloch

แอนิเมชันด้านล่างแสดงนูเทชัน: การขับที่เรโซแนนซ์หมุนเวกเตอร์ Bloch อย่างสม่ำเสมอรอบแกน $x$ พามันจาก $|0\rangle$ ลงไปทาง $|1\rangle$ แล้วกลับมา นี่คือการเคลื่อนที่ที่เส้นโค้ง Rabi วัดอยู่

แมปปัญหาเป็น Circuit

เราสร้าง Circuit ที่มีพารามิเตอร์หนึ่งตัว คือ $R_X(\theta)$ ตามด้วยการวัด แล้วประเมินที่หลายมุมพร้อมกันในภายหลัง

from qiskit import QuantumCircuit
from qiskit.circuit import Parameter

theta = Parameter("theta")
rabi = QuantumCircuit(1, 1)
rabi.rx(theta, 0)
rabi.measure(0, 0)
rabi.draw("mpl")

ปรับให้เหมาะสมสำหรับการรัน

Transpile เป็นวงจร Instruction Set Architecture (ISA) สำหรับ backend ที่เลือก

rabi_isa = pm_q.run(rabi)

Execute

เรา sweep ค่า $\theta$ จาก $0$ ถึง $4\pi$ (สองรอบ Rabi เต็ม) การส่งอาร์เรย์ของค่าพารามิเตอร์เป็น Primitive Unified Bloc (PUB) เดียว จะรันทั้งหมดในหนึ่ง job

angles = np.linspace(0, 4 * np.pi, 41)
job = sampler.run([(rabi_isa, angles.reshape(-1, 1))], shots=1000)
rabi_result = job.result()

Post-process และวิเคราะห์

ดึงค่า $P_1$ ที่แต่ละมุมออกมา แล้วเปรียบเทียบกับ $\sin^2(\theta/2)$

import matplotlib.pyplot as plt

ba = rabi_result[0].data.c  # BitArray, one entry per angle
shots = ba.num_shots
p1 = np.array([ba.get_counts(i).get("1", 0) / shots for i in range(len(angles))])

plt.figure(figsize=(7, 4))
plt.plot(angles, p1, "o", label="measured")
plt.plot(angles, np.sin(angles / 2) ** 2, "-", label=r"$\sin^2(\theta/2)$")
for k, name in [(np.pi / 2, r"$\pi/2$ (90°)"), (np.pi, r"$\pi$ (180°)")]:
    plt.axvline(k, color="gray", ls="--", lw=1)
    plt.text(k, 1.02, name, ha="center", fontsize=9)
plt.xlabel(r"flip angle $\theta$ (rad)")
plt.ylabel(r"$P_1$")
plt.title("Rabi / nutation curve")
plt.legend()
plt.tight_layout()
plt.show()

ตรวจสอบความเข้าใจของคุณ

บนเส้นโค้งที่คุณวัดได้ $P_1$ ถึงค่าสูงสุดครั้งแรกที่พัลส์ $\pi$ ($\theta = \pi$) ทีนี้สมมติว่าอุปกรณ์ถูกปรับเทียบผิด ทำให้มุม nominal $\theta$ บนแกน x สร้างการหมุนจริงเท่ากับ $1.05\,\theta$

(a) ค่าสูงสุดแรกของ $P_1$ จะปรากฏที่มุม nominal เท่าไร (b) การใช้พัลส์ $\pi$ เพียงครั้งเดียวเป็นวิธีที่ไวพอในการตรวจจับความคลาดเคลื่อน 5% นี้หรือไม่

คำตอบ

(a) ค่าสูงสุดปรากฏตรงที่การหมุน จริง เท่ากับ $\pi$ นั่นคือที่มุม nominal $\theta = \pi/1.05 \approx 0.95\pi$ — ยอดกราฟเลื่อนไปทางซ้ายเล็กน้อย

(b) ไม่ พัลส์ $\pi$ แบบ nominal ให้การหมุนจริง $1.05\pi$ ซึ่งได้ $P_1 = \sin^2(1.05\pi/2) \approx 0.994$ แทนที่จะเป็น $1$ — ความเปลี่ยนแปลงที่เล็กเกินกว่าจะมองเห็นท่ามกลางสัญญาณรบกวนจากการวัด นี่คือเหตุผลที่ในการ calibrate ให้แม่นยำ นักสเปกโทรสโกปี NMR จะขยายความคลาดเคลื่อนโดยทำการหมุนซ้ำหลายครั้ง (เช่น วิธี $360^\circ$ หรือขบวนพัลส์ยาว) เพื่อให้ความคลาดเคลื่อนต่อพัลส์เล็กๆ สะสมเป็นความเบี่ยงเบนใหญ่ที่วัดได้ เราจะกลับมาที่แนวคิดนี้ในโจทย์ท้าทาย

การทดลอง: $T_1$ (การคลายตัวแบบ spin-lattice)

ทฤษฎี

เมื่อปล่อยไว้เฉยๆ สปินที่ถูกกระตุ้นจะคายพลังงานให้สิ่งแวดล้อมรอบตัว ("lattice") และคลายตัวกลับสู่สมดุลทางความร้อน ประชากรของสถานะกระตุ้นลดลงแบบเอกซ์โพเนนเชียลด้วย เวลาคลายตัวแบบ spin-lattice $T_1$:

$P_1(\tau) = e^{-\tau/T_1}.$

สำหรับนักเคมี $T_1$ กำหนด recycle delay ที่คุณต้องรอระหว่างการสแกนใน NMR เชิงปริมาณ และยังบอกถึง การหมุนควงของโมเลกุล และพลวัตด้วย สำหรับ qubit มันคือเวลาคลายตัวของพลังงาน และเป็นหนึ่งในสองค่าหลักที่ใช้สร้างแบบจำลองผลของ decoherence ในคอมพิวเตอร์ควอนตัม

การทดลอง เตรียม $|1\rangle$ ด้วยเกต $X$ (พัลส์ $\pi$) รอเวลาหน่วงที่เปลี่ยนค่าได้ $\tau$ แล้วจึงวัด สัดส่วนที่ยังอยู่ใน $|1\rangle$ จะวาดเส้นการสลายตัวออกมา

มุมมองบน Bloch sphere

การคลายตัวแบบ $T_1$ พาสถานะกระตุ้นกลับสู่สมดุล ดังนั้นเวกเตอร์ Bloch จึงเลื่อนจาก $|1\rangle$ กลับขึ้นไปที่ $|0\rangle$ ตามแกน $z$

สร้างและ transpile การสแกนเวลาหน่วง

ระยะเวลาหน่วงต้องเป็นค่าที่แน่นอน (เพราะถูกจัดตารางตามเวลาของฮาร์ดแวร์) ดังนั้นเราจึงสร้างหนึ่ง circuit ต่อหนึ่งค่าหน่วง

# Scale the T1 delay sweep to the qubit's reported coherence, and snap each delay to an integer
# number of samples (dt) so there are no "duration rounded" warnings (as in the Ramsey/echo cells).
T1_us = (backend.qubit_properties(QUBIT).t1 or 100e-6) * 1e6  # reported T1 (us)
delays_us = np.linspace(0, 5 * T1_us, 20)  # 0 .. ~5 T1 captures the full decay
print(f"T1 sweep: 0..{5*T1_us:.0f} us in {len(delays_us)} points")

dt = backend.dt
t1_circuits = []
for tau in delays_us:
    n_dt = int(round(tau * 1e-6 / dt))  # integer samples -> no rounding warning
    qc = QuantumCircuit(1, 1)
    qc.x(0)  # pi pulse: prepare |1>
    qc.delay(n_dt, 0, unit="dt")  # wait tau
    qc.measure(0, 0)
    t1_circuits.append(qc)

t1_isa = pm_q.run(t1_circuits)

T1 sweep: 0..1393 us in 20 points

Execute

t1_job = sampler.run([(c,) for c in t1_isa], shots=1000)
t1_res = t1_job.result()

วิเคราะห์ — ฟิตเอกซ์โพเนนเชียล

from scipy.optimize import curve_fit

p1_t1 = np.array(
    [r.data.c.get_counts().get("1", 0) / r.data.c.num_shots for r in t1_res]
)

def decay(t, A, T1, B):
    return A * np.exp(-t / T1) + B

p0 = [1.0, 100.0, 0.0]
popt, _ = curve_fit(decay, delays_us, p1_t1, p0=p0, maxfev=10000)
T1_fit = popt[1]

plt.figure(figsize=(7, 4))
plt.plot(delays_us, p1_t1, "o", label="measured")
tt = np.linspace(0, delays_us.max(), 200)
plt.plot(tt, decay(tt, *popt), "-", label=rf"fit: $T_1$ = {T1_fit:.1f} µs")
plt.xlabel("delay τ (µs)")
plt.ylabel(r"$P_1$")
plt.title("$T_1$ relaxation")
plt.legend()
plt.tight_layout()
plt.show()
print(f"Measured T1 = {T1_fit:.1f} microseconds")

Measured T1 = 225.3 microseconds

ค่า $T_1$ ที่แน่นอนจะต่างกันไปในแต่ละการทดลองและแต่ละ backend แต่คุณควรเห็นผลลัพธ์ราวๆ ไม่กี่ร้อยไมโครวินาที

ตรวจสอบความเข้าใจของคุณ

ใน NMR เชิงปริมาณ เรามักถูกบอกให้รอประมาณ $5\,T_1$ ระหว่างการสแกน จากกฎการสลายตัว แมกเนไทเซชันสมดุลฟื้นคืนกลับมากี่เปอร์เซ็นต์หลังผ่านไป $5\,T_1$

คำตอบ

หลัง $5\,T_1$ ส่วนที่ยังเบี่ยงเบนจากสมดุลเหลือ $e^{-5}\approx 0.0067$ นั่นคือฟื้นคืนแล้วประมาณ 99.3%

การทดลอง: $T_2$, $T_2^*$

ทฤษฎี

$T_1$ วัดการสูญเสีย พลังงาน ส่วน $T_2$ วัดการสูญเสีย phase coherence ในระนาบตามขวาง เราสังเกตสิ่งนี้ได้ด้วย การทดลอง Ramsey: พัลส์ $\pi/2$ หมุนสปินทั้งหมดลงมาที่เส้นศูนย์สูตร โดยทุกตัวชี้ไปทิศเดียวกัน พวกมันหมุนควงแบบ Larmor แต่แต่ละตัวหมุนควงด้วยความถี่ต่างกันเล็กน้อย จึงกระจายออกเป็นพัดตามเส้นศูนย์สูตร และสัญญาณตามขวางรวมจะสลายตัวด้วยเวลาสลายตัวลักษณะเฉพาะ $T_2^*$

ข้อสังเกตสำหรับผู้อ่านสายควอนตัม: ภาพการกระจายเป็นพัดแสดงวิธีหนึ่งของการสูญเสีย coherence คือการที่เฟสที่แน่นอนหลายค่ากระจายตัวแล้วหักล้างกันไป ควรรู้ว่านี่ไม่ใช่กลไกเดียว coherence อาจสูญเสียได้เพราะเฟสพันกัน (entangle) กับสิ่งแวดล้อม จนมันไม่ใช่ค่าที่แน่นอนซึ่งระบบถืออยู่อีกต่อไปเลย เราแยกสองกลไกนี้ออกจากกันไม่ได้ด้วยการดูแค่การสลายตัวของสัญญาณ ดังนั้นภาพการกระจายเป็นพัดจึงเป็นสัญชาตญาณที่ดี แต่ระวังอย่าตีความว่าเป็นเรื่องราวที่สมบูรณ์หรือตามตัวอักษร

มีสองผลที่ทำให้สปินเสียเฟส:

• decoherence ที่ย้อนกลับไม่ได้ (หรือแทบย้อนกลับไม่ได้) ซึ่งเกิดจากการเชื่อมโยงที่ควบคุมไม่ได้กับองศาอิสระอื่น ลักษณะเฉพาะด้วย $T_2$ ที่แท้จริง ใน NMR ส่วนใหญ่มาจากการเคลื่อนที่ของโมเลกุลที่ปรับการเชื่อมโยงระหว่างสปิน ในคอมพิวเตอร์ควอนตัม มันสอดคล้องกับการเชื่อมโยงที่ไม่ต้องการระหว่าง qubit กับสิ่งแวดล้อม
• ออฟเซตความถี่ที่คงที่หรือเปลี่ยนช้า ซึ่งย้อนกลับได้ เช่น ความไม่สม่ำเสมอของสนาม ความต่างของความถี่ระหว่าง qubit และข้อผิดพลาดจากการ calibrate เมื่อรวมกับ $T_2$ ที่แท้จริง จะให้เวลาสลายตัวปรากฏที่สั้นลงคือ $T_2^*$ ที่สังเกตได้ในการทดลอง Ramsey

แล้วเราจะหา $T_2$ ที่แท้จริงได้อย่างไร ถ้าการทดลอง Ramsey บอกได้แค่ $T_2^*$ การทดลองที่สองเรียกว่า Hahn echo ซึ่งหักล้างสัญญาณรบกวนประเภทที่สองทั้งหมดเพื่อหา $T_2$ ที่แท้จริง มันเริ่มเหมือนการทดลอง Ramsey ด้วยพัลส์ $\pi/2$ แล้วสปินกระจายออกเหมือนเดิม หลังเวลาผ่านไป $\tau$ จะใส่พัลส์ $\pi$ พัลส์นี้ย้อนการสะสมเฟสที่เกิดจากออฟเซตความถี่คงที่ ทำให้สปินกลับมารวมกันอีกครั้งหลังการหน่วงครั้งที่สองที่ยาวเท่ากันคือ $\tau$

ทั้ง $T_2^*$ และ $T_2$ เป็นตัวเลขสำคัญที่ควรรู้ ทั้งในบริบทเคมีและบริบทคอมพิวเตอร์ควอนตัม ในบริบทเคมี $T_2^*$ กำหนดความกว้างเส้น (linewidth) ในการทดลอง NMR มาตรฐาน (คุณจะเห็นความหมายในส่วนถัดไป) $T_2$ บอกเรื่องโมเลกุลได้มาก เช่น การหมุนควงของโมเลกุล และการเปลี่ยนแปลงของ $T_2$ อาจบ่งว่าโมเลกุลไปจับกับสิ่งที่ใหญ่กว่า ในบริบทคอมพิวเตอร์ควอนตัม โดยทั่วไปเราต้องการให้ทั้ง $T_2$ และ $T_2^*$ ยาวที่สุดเท่าที่จะทำได้ เพื่อให้มีเวลาทำการทดลองมากขึ้นก่อนที่ decoherence จะครอบงำระบบ เทคนิคทั้งตระกูลที่อิงกับ Hahn echo ช่วยยืดเวลานี้ด้วยการหักล้างสัญญาณรบกวนที่ย้อนกลับได้ ทำให้ circuit ถูกจำกัดด้วย $T_2$ ที่แท้จริงในอุดมคติ ไม่ใช่ $T_2^*$ เทคนิคนี้เรียกว่า "dynamical decoupling"

การทำให้เป็นจริง

เมื่อเราทำการทดลอง Ramsey และ Hahn echo บน qubit การอ่านค่าจะต่างจากกรณี NMR เล็กน้อย ใน NMR คุณสังเกตสนามที่แกว่งซึ่งเหนี่ยวนำโดยสปินที่หมุนควงที่ความถี่ Larmor ในระนาบตามขวางได้โดยตรง แต่กับ qubit เราวัดได้เพียงการฉายลงบนแกน $Z$ คือเป็น $\vert 0 \rangle$ หรือ $\vert 1 \rangle$ ดังนั้นในทั้งสองการทดลอง เราจึงปิดท้ายด้วยพัลส์ $\pi/2$ เพิ่มอีกหนึ่งครั้งแล้วจึงวัด พัลส์ $\pi/2$ นี้แปลงข้อมูลเฟสของ qubit ให้เป็นประชากร

ด้านล่างเรารันแต่ละการทดลองเพื่อหา $T_2^*$ และ $T_2$

Ramsey ($T_2^{*}$): พัลส์ $\pi/2$ สองครั้งโดยมีเวลาหน่วง $\tau$

พัลส์ $\pi/2$ แรกพาเวกเตอร์ Bloch ไปที่เส้นศูนย์สูตร และพัลส์ถัดไปแปลงเฟสของเวกเตอร์ Bloch ให้เป็นประชากร เนื่องจากเฟสเปลี่ยนตามเวลาขณะที่เวกเตอร์ Bloch หมุนควง เราจึงเห็นรูปแบบไซน์ในประชากรเป็นฟังก์ชันของเวลาหมุนควง $\tau$ เมื่อ qubit สูญเสีย phase coherence เราควรเห็นแอมพลิจูดของเส้นไซน์นี้สลายตัวลง

ความละเอียดอ่อนเพิ่มเติมสำหรับผู้อ่านสายควอนตัม: ในการทดลองคอมพิวเตอร์ควอนตัม การแกว่งที่เราเห็นไม่ใช่ที่ความถี่ Larmor จริงๆ ซึ่งอยู่ในช่วง GHz และเร็วเกินกว่าที่เราจะแยกแยะได้ การทดลองควอนตัมทำงานใน rotating frame: กรอบอ้างอิงที่หมุนไปพร้อมกับ qubit ที่ความถี่ (ใกล้) Larmor ทำให้การหมุนควงที่เร็วมองไม่เห็น และเวกเตอร์ Bloch ดูเหมือนแทบนิ่ง เพื่อให้เห็นริ้ว (fringe) เราจงใจตั้ง rotating frame ให้เพี้ยนจากความถี่ Larmor เล็กน้อย เวกเตอร์ Bloch จึงดูเหมือนหมุนควงที่ความเพี้ยนเล็กๆ นั้น และเราได้สัญญาณแกว่งที่ความถี่ของความเพี้ยน ด้านล่างนี้ เราไม่ได้ตรึงความเพี้ยนไว้ที่ค่าเดียว แต่เลือกให้มีริ้วประมาณหนึ่งโหลครอบคลุมช่วงการวัด (ซึ่งเรากำหนดขนาดตามเวลา coherence ของ qubit) ซึ่งทำให้เก็บตัวอย่างการแกว่งได้ง่ายโดยไม่เกิด aliasing อีกเรื่องที่ยุ่งยาก: ใน Qiskit เราไม่มีการควบคุมระดับพัลส์ที่จำเป็นต่อการเปลี่ยนความถี่ของพัลส์ (ซึ่งเป็นตัวกำหนดความถี่ของ rotating frame) แทนที่จะทำเช่นนั้น เราเพิ่มความเพี้ยน ประสิทธิผล ด้วยมือ โดยใช้มันกำหนดเฟสของพัลส์ Ramsey ครั้งที่สอง มันดูประดิษฐ์ไปหน่อย แต่ผลลัพธ์เหมือนกัน และการสลายตัวที่เราสังเกตก็ยังคงเป็น $T_2^*$ ของจริง

มุมมองบน Bloch sphere

พัลส์ $\pi/2$ พาสปินไปที่เส้นศูนย์สูตร จากนั้นแต่ละตัวหมุนควงด้วยอัตราต่างกันเล็กน้อยและ กระจายเป็นพัด ทำให้เวกเตอร์สุทธิ (สีแดง) หดสั้นลงเมื่อพวกมันเสียเฟส — การสลายตัวนี้คือ $T_2^*$

# Adaptive Ramsey: scale the window to the reported coherence and choose the detuning so a fixed
# number of fringes span it. This decouples the point count from the window -- a long window stays
# cheap and cannot alias (we keep 8 samples per fringe).
T2_us = (backend.qubit_properties(QUBIT).t2 or 100e-6) * 1e6  # reported (echo) T2 (us)
ramsey_span = 0.5 * T2_us  # window = half the reported T2
n_fringes = 12  # fringes across the window
f_detune = n_fringes / (
    ramsey_span * 1e-6
)  # adaptive detuning (Hz) -> low, so a long window is cheap
samp_per_fringe = 8
n_ramsey = n_fringes * samp_per_fringe + 1  # ~97 points, fixed regardless of T2
ramsey_delays_us = np.linspace(0, ramsey_span, n_ramsey)
print(
    f"Ramsey sweep: 0..{ramsey_span:.0f} us in {n_ramsey} points "
    f"(detuning {f_detune/1e3:.1f} kHz, {samp_per_fringe} samples/fringe)"
)

dt = backend.dt  # seconds per sample
delays_dt = [int(round(t * 1e-6 / dt)) for t in ramsey_delays_us]  # integer samples
actual_delays_us = np.array(delays_dt) * dt * 1e6  # true (snapped) delays for fitting

ramsey_circ = []
for n_dt, tau_us in zip(delays_dt, actual_delays_us):
    qc = QuantumCircuit(1, 1)
    qc.rx(np.pi / 2, 0)
    qc.delay(n_dt, 0, unit="dt")  # integer duration -> no rounding warning
    qc.rz(
        2 * np.pi * f_detune * tau_us * 1e-6, 0
    )  # detuning phase (uses the true delay)
    qc.rx(np.pi / 2, 0)
    qc.measure(0, 0)
    ramsey_circ.append(qc)

ramsey_isa = pm_q.run(ramsey_circ)
ramsey_res = sampler.run([(c,) for c in ramsey_isa], shots=1000).result()
p1_ramsey = np.array(
    [r.data.c.get_counts().get("1", 0) / r.data.c.num_shots for r in ramsey_res]
)

ramsey_circ[len(ramsey_circ) // 2].draw("mpl")  # draw a sample circuit (mid-sweep)

Ramsey sweep: 0..248 us in 97 points (detuning 48.4 kHz, 8 samples/fringe)

Hahn echo ($T_2$): เพิ่มพัลส์ $\pi$ ที่ทำให้กลับมารวมกันตรงกลาง

มุมมองบน Bloch sphere

Hahn echo: สปินกระจายเป็นพัดเหมือนเดิม จากนั้นพัลส์ $\pi$ พลิกพวกมัน ทำให้การหมุนควงต่อเนื่องของพวกมัน รวมกลับมา อีกครั้ง (เวกเตอร์สีแดงโตขึ้นอีก) ฟื้นคืน coherence ที่เสียไปจากการเสียเฟสแบบคงที่

# Scale the echo sweep to the reported (echo) T2. The echo signal is monotonic (no fringes),
# so only the range matters here; 0 .. ~2.5 T2 shows the full decay for a clean fit.
T2_us = (backend.qubit_properties(QUBIT).t2 or 100e-6) * 1e6
echo_delays_us = np.linspace(0, 2.5 * T2_us, 24)  # total free-evolution time
print(f"Echo sweep: 0..{2.5*T2_us:.0f} us in {len(echo_delays_us)} points")

dt = backend.dt
half_dt = [
    int(round(t / 2 * 1e-6 / dt)) for t in echo_delays_us
]  # samples per half-delay
actual_echo_us = np.array(half_dt) * 2 * dt * 1e6  # true total delay for fitting

echo_circ = []
for n_half in half_dt:
    qc = QuantumCircuit(1, 1)
    qc.rx(np.pi / 2, 0)
    qc.delay(n_half, 0, unit="dt")  # half the free-evolution, integer samples
    qc.rx(np.pi, 0)  # refocusing pi pulse
    qc.delay(n_half, 0, unit="dt")
    qc.rx(-np.pi / 2, 0)
    qc.measure(0, 0)
    echo_circ.append(qc)

echo_isa = pm_q.run(echo_circ)
echo_res = sampler.run([(c,) for c in echo_isa], shots=1000).result()
p1_echo = np.array(
    [r.data.c.get_counts().get("1", 0) / r.data.c.num_shots for r in echo_res]
)

echo_circ[len(echo_circ) // 2].draw("mpl")  # draw a sample circuit

Echo sweep: 0..1240 us in 24 points

วิเคราะห์: ฟิตและเปรียบเทียบการสลายตัวทั้งสอง

# Initial guesses / bounds derived from the reported coherence so the fits track the scaled axes.
T2_us = (backend.qubit_properties(QUBIT).t2 or 100e-6) * 1e6

# Ramsey: decaying oscillation -> extract T2*
def ramsey_fn(t, A, T2s, f, phi, B):
    return A * np.exp(-t / T2s) * np.cos(2 * np.pi * f * t + phi) + B

rp0 = [0.4, T2_us / 2, f_detune * 1e-6, 0.0, 0.5]  # T2* guess ~ half the echo T2
rbounds = (
    [0, 0.1, 0, -np.pi, 0],
    [1, 2 * T2_us, 5, np.pi, 1],
)  # keep the fit from collapsing to ~0
rpopt, _ = curve_fit(
    ramsey_fn, ramsey_delays_us, p1_ramsey, p0=rp0, bounds=rbounds, maxfev=20000
)
T2_star = abs(rpopt[1])

# Echo: simple exponential envelope -> extract T2
def env(t, A, T2, B):
    return A * np.exp(-t / T2) + B

epopt, _ = curve_fit(env, echo_delays_us, p1_echo, p0=[0.5, T2_us, 0.0], maxfev=10000)
T2 = abs(epopt[1])

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(ramsey_delays_us, p1_ramsey, "o")
tt = np.linspace(0, ramsey_delays_us.max(), 400)
ax[0].plot(tt, ramsey_fn(tt, *rpopt), "-")
ax[0].set_title(rf"Ramsey: $T_2^*$ = {T2_star:.1f} µs")
ax[0].set_xlabel("τ (µs)")
ax[0].set_ylabel(r"$P_1$")
ax[1].plot(echo_delays_us, p1_echo, "o")
tt2 = np.linspace(0, echo_delays_us.max(), 400)
ax[1].plot(tt2, env(tt2, *epopt), "-")
ax[1].set_title(rf"Hahn echo: $T_2$ = {T2:.1f} µs")
ax[1].set_xlabel("τ (µs)")
plt.tight_layout()
plt.show()
print(f"T2* (Ramsey) = {T2_star:.1f} µs ;  T2 (echo) = {T2:.1f} µs")

T2* (Ramsey) = 171.2 µs ;  T2 (echo) = 326.0 µs

ตรงนี้เราเห็นทั้งการแกว่งแบบ Ramsey และสัญญาณ Hahn echo สลายตัวตามที่คาด สเกลเวลาของการสลายตัวทั้งสองคือ $T_2^*$ และ $T_2$ ตามลำดับ

ตรวจสอบความเข้าใจของคุณ

โดยทั่วไปคุณควรพบว่า $T_2 \gtrsim T_2^{*}$ จงอธิบายเหตุผลในหนึ่งประโยค และบอกว่าพัลส์ $\pi$ ของ Hahn echo ขจัดผลทางกายภาพอะไรได้ และอะไรที่ขจัดไม่ได้

คำตอบ

$T_2^{*}$ รวมทั้งการเสียเฟสที่ย้อนกลับไม่ได้ และ การกระจายแบบคงที่ที่ย้อนกลับได้ จึงสลายตัวเร็วกว่า พัลส์ $\pi$ ของ echo ทำให้การเสียเฟสแบบ คงที่และย้อนกลับได้ (ความไม่สม่ำเสมอของสนาม / ออฟเซตความถี่คงที่) กลับมารวมกัน เหลือเพียงการเสียเฟสแบบ spin-spin ที่ ย้อนกลับไม่ได้ ซึ่งคือ $T_2$ ที่แท้จริง echo ไม่สามารถยกเลิกความผันผวนที่เปลี่ยนไประหว่างสองครึ่งของลำดับได้

NMR spectroscopy ของโมเลกุล: Chemical shift และ $J$-coupling

จากสปินเดี่ยวสู่โมเลกุล

การทดลองที่ผ่านมาศึกษาสปินเดี่ยวที่แยกอิสระ NMR ทางเคมีจริงๆ อ่าน โมเลกุล (สปินจำนวนมากที่เชื่อมโยงกันพร้อมกัน) และใช้สเปกตรัมที่ได้เพื่ออนุมานลักษณะของโครงสร้างโมเลกุล สิ่งนี้บอกเราได้สองผล:

• Chemical shift นิวเคลียสที่ต่างกันทางเคมีแต่ละตัวหมุนควงที่ความถี่ $\nu_i$ ต่างกันเล็กน้อย เพราะอิเล็กตรอนในบริเวณใกล้เคียงกำบังสนามที่ใช้ในปริมาณต่างกัน ดังนั้น ตำแหน่ง ของพีคจึงระบุสภาพแวดล้อมทางเคมีได้ (เช่น โปรตอนของ CH$_3$ เทียบกับโปรตอนอะโรมาติก)

• Scalar $J$-coupling ปฏิสัมพันธ์ผ่านพันธะระหว่างสองสปินแยกแต่ละพีคเป็น multiplet นิวเคลียสที่มีเพื่อนบ้าน เทียบเท่ากัน $n$ ตัวจะแยกเป็น $n+1$ เส้นที่มีความเข้มตามทวินาม เรียกว่า กฎ $n+1$ รูปแบบ การแยกเผยการเชื่อมต่อ (นิวเคลียสใดเชื่อมพันธะกัน) และระยะห่างคือค่าคงที่การเชื่อมโยง $J$

การทดลอง

เพื่อจำลองสิ่งนี้บนคอมพิวเตอร์ควอนตัมของเรา เราต้องสร้าง spin Hamiltonian บน qubit ในลิมิตการเชื่อมโยงอ่อนและสนามสูง (ความต่างของ chemical shift $\gg$ การเชื่อมโยง) ส่วน $XX$ และ $YY$ ของการเชื่อมโยงเฉลี่ยหักล้างกันไป เหลือเพียง $ZZ$ ดังนั้น Hamiltonian คือ (หน่วยความถี่เชิงมุม, $\hbar=1$):

$H = \underbrace{\sum_i \pi\,\nu_i\, Z_i}_{\text{chemical shift}} \;+\; \underbrace{\sum_{i<j} \tfrac{\pi}{2}\,J_{ij}\, Z_iZ_j}_{J\text{-coupling}}.$

นิวเคลียสแต่ละตัวคือ หนึ่ง qubit: chemical shift คือการหมุน $Z$ บน qubit เดี่ยว และการเชื่อมโยงคือปฏิสัมพันธ์ $ZZ$ ระหว่างสอง qubit

การทดลองดำเนินดังนี้ พัลส์ $\pi/2$ พาทุกสปินเข้าสู่ระนาบตามขวาง และสปินหมุนควงและเชื่อมโยงกันเป็นเวลา $t$ เราบันทึกแมกเนไทเซชันตามขวางเชิงซ้อน $s(t) = \langle X\rangle + i\langle Y\rangle$ ที่รวมจากทุกสปิน เช่นเดียวกับการทดลอง Ramsey ข้างต้น แมกเนไทเซชันตามขวางนี้วาดเป็นไซน์ที่สลายตัวเป็นฟังก์ชันของเวลาหมุนควง $t$ เมื่อเราทำการแปลงฟูเรียร์ของสัญญาณ เราจะได้สเปกตรัมที่มีพีคที่ chemical shift $\nu_i$ ซึ่งถูกแยกโดยการเชื่อมโยงสเกลาร์ $J_{ij}$

หมายเหตุ

เราจำลองวิวัฒนาการของโมเลกุลโดยแมปเวลา $t$ ไปเป็น มุมการหมุน: $t$ ที่มากขึ้นหมายถึงการหมุน $Z$ และ $ZZ$ ที่ใหญ่ขึ้น ไม่ใช่ circuit ที่ยาวหรือลึกขึ้น ดังนั้นแม้สถานะที่จำลองจะก้าวหน้าไปตาม $t$ (ซึ่งทำให้ free induction decay (FID) แกว่ง) แต่ qubit จริงไม่เคยอยู่เฉยนานขึ้น และไม่มี decoherence ที่ขึ้นกับ $t$ FID ที่จำลองจึงไม่สลายตัวเองเลย สำหรับ FID ในอุดมคติที่ยาวไม่สิ้นสุด เส้นสเปกตรัมจะคมไม่สิ้นสุด แต่หน้าต่างการสุ่มตัวอย่างที่จำกัดของเรากลับให้ความละเอียดจำกัดและอาจทำให้เกิดการรั่วของสเปกตรัม เส้น NMR จริงมีความกว้างจำกัดที่กำหนดโดย $T_2^*$ เราเติมมันกลับโดยคูณ FID ด้วยหน้าต่าง $e^{-t/T_2^{*}}$ ก่อนการแปลงฟูเรียร์ ซึ่งเรียกว่า line-broadening apodization ที่นักสเปกโทรสโกปีใช้ในการประมวลผลเพื่อกำหนดความกว้างเส้น (แลกความละเอียดกับอัตราส่วนสัญญาณต่อสัญญาณรบกวน) (เซลล์ทางเลือกในภาคผนวกของโมดูลนี้แสดงทางเลือกแบบ Trotterized ที่การสลายตัวมาจาก decoherence จริงซึ่งเพิ่มตามความลึกของ circuit)

ระบบ AX

เริ่มจากระบบ AX ซึ่งเป็นโมเลกุลที่ประกอบด้วยสองสปินที่เชื่อมโยงกันและต่างกันทางเคมี ป้ายกำกับ A และ X แทนสภาพแวดล้อมสปินที่ต่างกันทางเคมีสองแบบ ตัวอักษรเลือกได้ตามใจ แต่การที่มันอยู่ห่างกันในตัวอักษรสื่อถึงความไม่เหมือนกันทางเคมี

แมปและปรับให้เหมาะสม

เราเขียนฟังก์ชันหนึ่งที่รับออฟเซต chemical shift $\nu_i$ และการเชื่อมโยง $J_{ij}$ แล้วคืนค่า circuit เตรียมสถานะแบบ มีพารามิเตอร์ พัลส์ $\pi/2$ พาสปินไปที่เส้นศูนย์สูตร จากนั้นการหมุน $Z$ จำลองวิวัฒนาการของสถานะภายใต้พจน์ chemical shift และ $J$-coupling แต่ละ $\nu_i$ คือ ออฟเซต การหมุนควงเล็กๆ เป็น Hz เทียบกับค่าอ้างอิงที่เลือก หรือคลื่นพาห์ของ rotating frame ตรงนี้เรากำหนดเอง ในโมเลกุลจริง การกำบังของอิเล็กตรอนเปลี่ยนสนามแม่เหล็กเฉพาะที่และเลื่อนเรโซแนนซ์ของนิวเคลียสแต่ละตัวเทียบกับค่าอ้างอิงนั้น สิ่งนี้จำลองโดยเกต $R_Z$ บน qubit เดี่ยว หนึ่งตัวต่อหนึ่งสปิน และแต่ละ $\nu_i$ กำหนดตำแหน่งพีคของสปินนั้นในสเปกตรัม เกต $ZZ$ สอง qubit ใช้การเชื่อมโยง $J_{ij}$ ซึ่งแยกแต่ละพีคเป็น multiplet และกำหนดระยะห่างระหว่างเส้น

from qiskit import QuantumCircuit
from qiskit.circuit import Parameter
from qiskit.quantum_info import SparsePauliOp

t = Parameter(
    "t"
)  # evolution time (arbitrary units; the frequencies below are in matching units)

def fid_circuit(offsets, couplings):
    """State-preparation circuit for a set of spins: a pi/2 pulse, then evolution under the
    chemical-shift + J-coupling Hamiltonian for time t. No measurement -- the Estimator reads
    <X> and <Y> directly.

    offsets   : list of chemical-shift frequencies nu_i, one per qubit
    couplings : dict {(i, j): J_ij}
    """
    n = len(offsets)
    qc = QuantumCircuit(n)
    qc.ry(np.pi / 2, range(n))  # pi/2 pulse: tip every spin into the transverse plane
    for i, nu in enumerate(offsets):
        qc.rz(2 * np.pi * nu * t, i)  # chemical-shift precession -> line at nu_i
    for (i, j), J in couplings.items():
        qc.rzz(np.pi * J * t, i, j)  # J-coupling -> splitting of size J
    return qc

def transverse_obs(n):
    """Total transverse observables (sum over spins): (sum_i X_i, sum_i Y_i)."""
    ox = SparsePauliOp.from_sparse_list(
        [("X", [i], 1.0) for i in range(n)], num_qubits=n
    )
    oy = SparsePauliOp.from_sparse_list(
        [("Y", [i], 1.0) for i in range(n)], num_qubits=n
    )
    return ox, oy

# AX system: two distinct chemical shifts, one coupling
ax_offsets = [6.0, -6.0]  # chemical shifts nu_1, nu_2 (Hz)
ax_couplings = {(0, 1): 2.0}  # J_12 = 2 Hz

ax_isa = pm_ax.run(fid_circuit(ax_offsets, ax_couplings))
fid_circuit(ax_offsets, ax_couplings).draw("mpl")

Execute และ post-process

เราประเมิน circuit บนกริดของเวลาวิวัฒนาการ และ Estimator คืนค่า $\langle X\rangle$ และ $\langle Y\rangle$ ที่แต่ละจุด เรารวมเป็น FID เชิงซ้อน $\langle X\rangle + i\langle Y\rangle$ ใช้หน้าต่าง apodization และแปลงฟูเรียร์ผลลัพธ์เพื่อให้ได้สเปกตรัม

# --- time grid ---
n_t = 128  # number of FID samples
dt = 1 / 32  # time step -> Nyquist frequency = 1/(2 dt) = 16 Hz
times = np.arange(n_t) * dt

# --- Step 3-4: Estimator returns <X> and <Y> directly (both quadratures in one job) ---
def run_fid(isa_circ, n, est):
    """Complex FID s(t) = sum_i <X_i> + i <Y_i> over the time grid, via the Estimator."""
    ox, oy = transverse_obs(n)
    ox = ox.apply_layout(isa_circ.layout)  # map observables onto the transpiled qubits
    oy = oy.apply_layout(isa_circ.layout)
    params = times.reshape(-1, 1)
    res = est.run(
        [(isa_circ, ox, params), (isa_circ, oy, params)]
    ).result()  # both quadratures, one job
    return res[0].data.evs + 1j * res[1].data.evs

def spectrum(fid, T2=1.0):
    sig = fid * np.exp(-times / T2)  # apodization -> Lorentzian linewidth ~ 1/(pi T2)
    spec = np.fft.fftshift(np.fft.fft(sig))
    freq = np.fft.fftshift(np.fft.fftfreq(len(times), dt))
    return freq, np.abs(spec)

fid_ax = run_fid(ax_isa, n=2, est=estimator)
freq, spec = spectrum(fid_ax, T2=1.0)

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(times, fid_ax.real, label=r"Re $s(t)=\langle X\rangle$")
ax[0].plot(times, fid_ax.imag, label=r"Im $s(t)=\langle Y\rangle$")
ax[0].set_xlabel("evolution time $t$")
ax[0].set_ylabel("transverse signal")
ax[0].set_title("AX free induction decay")
ax[0].legend()
ax[1].plot(freq, spec)
for nu in ax_offsets:
    for s in (-1, 1):
        ax[1].axvline(nu + s * ax_couplings[(0, 1)] / 2, color="gray", ls="--", lw=1)
ax[1].set_xlim(-12, 12)
ax[1].set_xlabel(r"frequency $\nu$ (Hz)")
ax[1].set_ylabel("intensity")
ax[1].set_title("AX spectrum: two doublets, split by $J$")
plt.tight_layout()
plt.show()
print(
    "Expected AX lines (Hz):",
    sorted({nu + s * ax_couplings[(0, 1)] / 2 for nu in ax_offsets for s in (-1, 1)}),
)

Expected AX lines (Hz): [-7.0, -5.0, 5.0, 7.0]

ส่วนขยาย: ระบบ AX$_2$ และกฎ $n+1$

ทีนี้เพิ่มเพื่อนบ้านที่เทียบเท่ากันตัวที่สอง ในระบบ AX$_2$ สปิน A หนึ่งตัวเชื่อมโยงด้วย $J$ เท่ากันกับสปิน X ที่เทียบเท่ากัน สองตัว (จำไว้ว่า X เป็นเพียงป้ายกำกับของสองสปินที่เทียบเท่ากัน ไม่เกี่ยวกับตัวดำเนินการ Pauli $X$ การเชื่อมโยงทุกตัวด้านล่างยังคงเป็น $ZZ$) สปิน A จึงเห็นรูปแบบเพื่อนบ้านสี่แบบ $(\uparrow\uparrow,\ \uparrow\downarrow,\ \downarrow\uparrow,\ \downarrow\downarrow)$ สองแบบที่ผสมกันทับกัน A จึงแยกเป็น triplet 1:2:1 ($n=2 \Rightarrow n+1 = 3$ เส้น) สปิน X แต่ละตัวเห็นแค่ A และยังเป็น doublet สปินที่เทียบเท่ากันไม่แยกกันเอง เราจึงไม่ใส่การเชื่อมโยง X–X

โมเลกุลจริงที่มีรูปแบบนี้คือ 1,1,2-trichloroethane (Cl$_2$CH–CH$_2$Cl) เมื่อมุ่งที่สปินของไฮโดรเจน (ซึ่งเป็นรูปแบบ NMR ที่พบบ่อยที่สุด) โปรตอนไฮโดรเจนเดี่ยวใน CH ทำหน้าที่เป็น A และโปรตอน CH$_2$ ที่เทียบเท่ากันสองตัวคือสปิน X สองตัว ทีนี้ลองนับเพื่อนบ้าน: โปรตอน CH มีเพื่อนบ้าน สอง ตัว (คู่ CH$_2$) สัญญาณจึงแยกเป็น triplet 1:2:1 ($n=2$) โปรตอน CH$_2$ แต่ละตัวมีเพื่อนบ้านเพียง หนึ่ง ตัว (CH) สัญญาณ CH$_2$ จึงเป็น doublet 1:1 ($n=1$) สังเกตว่ากลุ่มที่มีโปรตอน น้อยกว่า (CH) กลับแสดงเส้น มากกว่า จำนวนเส้นสะท้อนว่าสัญญาณมีเพื่อนบ้านกี่ตัว ไม่ใช่มีโปรตอนกี่ตัวใช้สัญญาณนั้นร่วมกัน

ด้านล่างเราจำลองระบบ AX$_2$ ด้วยสาม qubit

# AX2: spin A (qubit 0) coupled equally to two equivalent X spins (qubits 1, 2)
ax2_offsets = [6.0, -6.0, -6.0]  # A at +6 Hz; the two X both at -6 Hz
ax2_couplings = {
    (0, 1): 2.0,
    (0, 2): 2.0,
}  # equal A-X couplings; no coupling between the X spins

ax2_isa = pm_ax2.run(fid_circuit(ax2_offsets, ax2_couplings))
fid_ax2 = run_fid(ax2_isa, n=3, est=estimator)
freq2, spec2 = spectrum(fid_ax2, T2=1.0)

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(times, fid_ax2.real, label=r"Re $s(t)=\langle X\rangle$")
ax[0].plot(times, fid_ax2.imag, label=r"Im $s(t)=\langle Y\rangle$")
ax[0].set_xlabel("evolution time $t$")
ax[0].set_ylabel("transverse signal")
ax[0].set_title("AX$_2$ free induction decay")
ax[0].legend()
ax[1].plot(freq2, spec2)
ax[1].set_xlim(-12, 12)
ax[1].set_xlabel(r"frequency $\nu$ (Hz)")
ax[1].set_ylabel("intensity")
ax[1].set_title(r"AX$_2$ spectrum: A is a 1:2:1 triplet, X is a doublet")
plt.tight_layout()
plt.show()
print(
    "A triplet at (Hz):",
    [ax2_offsets[0] + k * ax2_couplings[(0, 1)] for k in (-1, 0, 1)],
)
print(
    "X doublet at (Hz):",
    [ax2_offsets[1] + s * ax2_couplings[(0, 1)] / 2 for s in (-1, 1)],
)

ตรวจสอบความเข้าใจของคุณ

อ่านแต่ละคำถาม คิดคำตอบของคุณ แล้วคลิกสามเหลี่ยมเพื่อดูเฉลย

ในสเปกตรัม AX doublet ทั้งสองห่างกัน (จากกึ่งกลางถึงกึ่งกลาง) 12 Hz และสองเส้นภายในแต่ละ doublet ห่างกัน 2 Hz ตัวเลขใดคือความต่างของ chemical shift และตัวใดคือค่าคงที่การเชื่อมโยง $J$ ถ้าเปลี่ยนไปใช้สเปกโตรมิเตอร์ที่มีสนามแม่เหล็กเป็นสองเท่า ตัวใดจะเปลี่ยน

คำตอบ

ระยะห่างกึ่งกลางถึงกึ่งกลาง 12 Hz คือ ความต่างของ chemical shift ($\nu_1-\nu_2$) ระยะห่างเส้น 2 Hz คือ ค่าคงที่การเชื่อมโยง $J$ Chemical shift มาจากการกำบัง $B_0$ จึงแปรผันตามสนาม ที่สนามสองเท่า ระยะห่างเพิ่มเป็นสองเท่าคือ 24 Hz ส่วนการเชื่อมโยง $J$ เป็นสมบัติของโมเลกุลผ่านพันธะและไม่ขึ้นกับสนาม จึงยังเป็น 2 Hz นักเคมีแยก shift ออกจาก coupling โดยดูว่าการแยกเส้นเปลี่ยนตามสนามอย่างไร

จงใช้กฎ $n+1$ ทำนาย multiplet ของ A ถ้าเพิ่มสปิน X ที่เทียบเท่ากันตัวที่ สาม (ระบบ AX$_3$) และอธิบายว่าทำไมสปิน X ที่เทียบเท่ากันสองตัวใน AX$_2$ จึงไม่แยกกันเอง

คำตอบ

เมื่อมีเพื่อนบ้านที่เทียบเท่ากัน $n=3$ ตัว A จะกลายเป็น quartet 1:3:3:1 สปินที่เทียบเท่ากันไม่แยกกันเองเพราะการเชื่อมโยงของพวกมันเลื่อนระดับพลังงานที่เกี่ยวข้องทั้งหมดไปเท่ากัน ช่องว่างระหว่างระดับเหล่านั้น — ความถี่ที่เราเห็น — จึงไม่เปลี่ยน

บทสรุป

เราใช้คอมพิวเตอร์ควอนตัมจริงสำรวจว่าสปินนิวเคลียสหมุนควงและเสีย coherence อย่างไรในการทดลอง NMR และตีความสเปกตรัมที่ได้อย่างไร ควรถอยออกมาทบทวนว่าทำไมสิ่งนี้เป็นไปได้ตั้งแต่แรก และทำไมมันจึงมากกว่าเรื่องน่าสนใจ

qubit และนิวเคลียสสปิน-1/2 ในทางคณิตศาสตร์คือสิ่งเดียวกัน: ระบบควอนตัมสองระดับบน Bloch sphere ที่ถูกหมุนด้วยพัลส์เรโซแนนซ์ หมุนควงรอบแกนหนึ่ง และสูญเสีย coherence ให้สิ่งแวดล้อม คำศัพท์ทั้งหมดของโมดูลนี้ (ความถี่ Larmor, rotating frame, $T_1$, $T_2$ และ $T_2^*$, ริ้ว Ramsey, Hahn echo) เป็นของทั้ง NMR และคอมพิวเตอร์ควอนตัมเท่าๆ กัน สิ่งที่นักเคมีเรียกว่า chemical shift วิศวกรควอนตัมเรียกว่า การหมุน $Z$ บน qubit เดี่ยว สิ่งที่นักเคมีเรียกว่า $J$-coupling วิศวกรควอนตัมเรียกว่า ปฏิสัมพันธ์ $ZZ$ ทั้งสองสาขาบรรยายฟิสิกส์เดียวกันด้วยถ้อยคำต่างกัน

ความเชื่อมโยงนี้เป็นเรื่องในประวัติศาสตร์ ไม่ใช่แค่แนวคิด: NMR เป็นหนึ่งในแพลตฟอร์มแรกที่รันอัลกอริทึมควอนตัมจริง ภายในปี 2001 โมเลกุลในสารละลายรันอัลกอริทึมของ Shor เพื่อแยกตัวประกอบ 15 โดยสปินนิวเคลียสแต่ละตัวทำหน้าที่เป็น qubit NMR ยังมอบเครื่องมือควบคุมจำนวนมากให้คอมพิวเตอร์ควอนตัม spin echo, ขบวนพัลส์ CPMG, composite pulse และ dynamical decoupling ล้วนถูกคิดค้นโดยนักสเปกโทรสโกปีหลายสิบปีก่อนที่ qubit จะต้องใช้มันต่อสู้กับ decoherence มีความสมมาตรที่น่าพอใจในสิ่งที่เราเพิ่งทำ: สาขาที่ให้กำเนิดคอมพิวเตอร์ควอนตัมยุคแรก ตอนนี้กลายเป็นสิ่งที่คอมพิวเตอร์ควอนตัมจำลองได้

ทบทวนและคำถาม

แนวคิดสำคัญ

• qubit คือระบบสองระดับที่ออกแบบขึ้นซึ่งเทียบเท่าสปิน-1/2 ค่า $|0\rangle/|1\rangle$, ความถี่, เกต, $T_1$ และ $T_2$ ของมันคือปริมาณที่เทียบได้โดยตรงกับปริมาณใน NMR

• พัลส์เรโซแนนซ์หมุนเวกเตอร์ Bloch มุมพลิกกำหนดการหมุน และเส้นโค้ง Rabi/nutation $P_1=\sin^2(\theta/2)$ ระบุการหมุน $\pi$ และ $\pi/2$

• $T_1$ คือการคลายตัวของพลังงาน (spin-lattice) $T_2$ คือการคลายตัวของเฟส (spin-spin) $T_2^{*}$ รวมการเสียเฟสแบบคงที่ที่ย้อนกลับได้เพิ่มเติม ซึ่ง Hahn echo ทำให้กลับมารวมกัน

• โมเลกุลคือชุดของสปินที่เชื่อมโยงกัน: chemical shift กำหนดตำแหน่งของแต่ละพีค และ $J$-coupling แยกพีคเป็น multiplet ตามกฎ $n+1$ ซึ่งรวมกันเข้ารหัสสภาพแวดล้อมทางเคมีและการเชื่อมต่อ

จริง / เท็จ

1. T/F ความถี่ Larmor คือช่องว่างพลังงานระหว่างสองสถานะสปินหารด้วย $\hbar$

2. T/F พัลส์ $90^\circ$ ($\pi/2$) กลับสปินสมดุลได้เต็มที่

3. T/F โดยทั่วไป $T_2^{*}$ ยาวกว่า $T_2$

4. T/F Hahn echo ทำให้การเสียเฟสแบบ spin-spin ที่ย้อนกลับไม่ได้กลับมารวมกัน

5. T/F ในระบบ AX$_2$ นิวเคลียส A ปรากฏเป็น triplet 1:2:1

6. T/F ค่าคงที่การเชื่อมโยง $J$ เพิ่มเป็นสองเท่าเมื่อคุณเพิ่มสนามแม่เหล็กของสเปกโตรมิเตอร์เป็นสองเท่า

7. T/F การเอา Qiskit Pulse ออกทำให้วัด $T_1$ บนฮาร์ดแวร์ปัจจุบันไม่ได้

คำตอบสั้น

1. จงอธิบายในสองประโยคว่าทำไมพัลส์ $\pi$ เพียงครั้งเดียวจึงเป็นวิธีที่แย่ในการตรวจจับความคลาดเคลื่อนของมุมพลิก $5\%$ และการต่อพัลส์หลายครั้งแก้ปัญหานี้อย่างไร

2. ตัวอย่างหนึ่งแสดงเส้น NMR กว้าง ในแง่ของ $T_2$ สิ่งนี้บอกอะไรคุณ และปริมาณของ qubit ที่เทียบเท่าคืออะไร

3. จงให้ปริมาณของ qubit ที่เทียบเท่ากับแต่ละข้อ: (a) recycle delay, (b) พัลส์ $90^\circ$, (c) chemical shift

4. จงใช้กฎ $n+1$ อธิบายว่าทำไมหมู่เอทิล (CH$_3$CH$_2$–) จึงแสดง triplet 1:2:1 และ quartet 1:3:3:1 และระบุว่าหมู่ใดให้ multiplet ใด

โจทย์ท้าทาย

1. ขบวน Carr-Purcell ปรับ Hahn echo ให้แทรกพัลส์ $\pi$ ที่เว้นระยะเท่ากัน $n$ ครั้ง (ขบวนแบบ CPMG) และวัด $T_2$ เป็นฟังก์ชันของ $n$ อภิปรายว่าทำไมพัลส์ที่ทำให้กลับมารวมกันมากขึ้นจึงทำให้ $T_2$ ปรากฏยาวขึ้นได้

2. การขยายความคลาดเคลื่อนของพัลส์ สร้าง circuit ที่ใช้ $R_X(\pi)$ ซ้ำ $\,n$ ครั้งกับ $|0\rangle$ และพล็อต $P_1$ เทียบกับ $n$ สำหรับมุมที่ตั้งผิดโดยตั้งใจ $R_X(1.05\pi)$ แสดงให้เห็นว่าความคลาดเคลื่อนมองเห็นได้ชัดขึ้นเมื่อ $n$ โตขึ้น คำใบ้: ใส่ barrier หลังการหมุนทุกครั้ง มิฉะนั้น transpiler อาจรวมการหมุนซ้ำๆ เป็นเกตเทียบเท่าเกตเดียว

3. หมู่เอทิล ใช้ fid_circuit สร้างสเปกตรัมโปรตอนของชิ้นส่วนเอทิล (CH$_3$–CH$_2$): ห้า qubit โดยโปรตอน CH$_3$ สามตัวใช้ chemical shift ร่วมกัน และโปรตอน CH$_2$ สองตัวอีกค่าหนึ่ง และทุกคู่ CH$_3$–CH$_2$ เชื่อมโยงด้วย $J$ เดียวกัน (ไม่มีการเชื่อมโยงภายในกลุ่มที่เทียบเท่ากัน) ทำนายรูปแบบด้วยกฎ $n+1$ แล้วจำลอง คุณควรได้ quartet 1:3:3:1 (CH$_2$) และ triplet 1:2:1 (CH$_3$)

ภาคผนวก (ทางเลือก): วิวัฒนาการแบบ Trotterized — ให้ decoherence กำหนดความกว้างเส้น

ในส่วนNMR spectroscopy ของโมเลกุล วิวัฒนาการ $e^{-iHt}$ ถูกใช้ในชั้นเดียวที่ความลึก ไม่ เพิ่มตาม $t$ ดังนั้น FID จึงไม่สลายตัวเอง และเรากำหนดความกว้างเส้นด้วยหน้าต่าง apodization $e^{-t/T_2}$ เซลล์ทางเลือกนี้ใช้แนวทางตรงข้าม มันแบ่งวิวัฒนาการเป็น $p$ ขั้น Trotter และใส่ barrier หลังแต่ละขั้น (เพื่อไม่ให้ transpiler รวมกลับเป็นชั้นเดียว) จุดเวลา $k$ ใช้ $k$ ขั้น ความลึกของ circuit จึงเพิ่มตามเวลา

เนื่องจากทุกพจน์ใน $H$ สลับที่กันได้ วิธีนี้จึง แม่นยำ: บนซิมูเลเตอร์ไร้สัญญาณรบกวนจะได้พีคคมเหมือนกัน แต่บน backend แบบมีโมเดลสัญญาณรบกวนหรือฮาร์ดแวร์จริง ข้อผิดพลาดของเกตที่สะสมทำให้ FID สลายตัวเอง และการสลายตัวนั้นเองกำหนดความกว้างเส้น นี่คือแนวคิด $T_2\leftrightarrow$ความกว้างเส้น ของส่วนการทดลอง $T_2$, $T_2^*$ ที่ตอนนี้เกิดจาก decoherence จริงแทนหน้าต่างการประมวลผล เรารันทั้งระบบ AX (สอง doublet) และระบบ AX$_2$ (triplet 1:2:1 กับ doublet) ด้วยวิธีนี้

มีข้อสังเกตเชิงปฏิบัติสองข้อ ข้อแรก ต้นทุนคือความลึก: circuit AX ที่ลึกที่สุดใช้เกตสอง qubit ประมาณ $\sim\!47$ ตัว และเวอร์ชัน AX$_2$ ราวสองเท่า ทั้งสองจึงมีสัญญาณรบกวนและแพงกว่าเวอร์ชันแม่นยำมาก ข้อสอง ระบบสปินแต่ละระบบยังถูกส่งเป็นหนึ่ง job และทุกจุดเวลาใช้เลย์เอาต์ qubit กายภาพเดิมที่เลือกไว้ในส่วนสเปกโทรสโกปี ถ้า USE_BATCH=True job เหล่านี้จะใช้ Batch ทางเลือกทั่วทั้ง notebook ด้วย

# Trotterized FID: depth grows with time, so real decoherence (not apodization) sets the linewidth.
estimator.options.default_shots = 2000  # deep circuits need more shots in this appendix
n_tr, dt_tr = (
    48,
    1 / 20,
)  # Nyquist = 1/(2 dt) = 10 Hz (covers the AX2 triplet's 8 Hz line)
times_tr = np.arange(n_tr) * dt_tr
APODIZE_TROTTER = False  # keep False so the linewidth comes from decoherence

def trotter_fid_circuit(offsets, couplings, p, step):
    n = len(offsets)
    qc = QuantumCircuit(n)
    qc.ry(np.pi / 2, range(n))  # 90-degree pulse
    for _ in range(p):  # p Trotter steps, each of duration `step`
        for i, nu in enumerate(offsets):
            qc.rz(2 * np.pi * nu * step, i)  # chemical-shift evolution for one step
        for (i, j), J in couplings.items():
            qc.rzz(np.pi * J * step, i, j)  # J-coupling evolution for one step
        qc.barrier()  # <-- stops the transpiler merging the steps
    return qc

def run_trotter(offsets, couplings, n, est, pass_manager):
    # Time point k uses k Trotter steps (its own circuit). Submit every time point and both
    # quadratures together as ONE Estimator job, using one fixed physical-qubit layout.
    ox, oy = transverse_obs(n)
    pubs = []
    for k in range(n_tr):
        isa = pass_manager.run(trotter_fid_circuit(offsets, couplings, k, dt_tr))
        pubs.append((isa, ox.apply_layout(isa.layout)))  # <X> at time k
        pubs.append((isa, oy.apply_layout(isa.layout)))  # <Y> at time k
    res = est.run(pubs).result()
    return np.array(
        [res[2 * k].data.evs + 1j * res[2 * k + 1].data.evs for k in range(n_tr)]
    )

def trotter_spectrum(fid):
    sig = fid * np.exp(-times_tr / 1.0) if APODIZE_TROTTER else fid
    spec = np.abs(np.fft.fftshift(np.fft.fft(sig)))
    freq = np.fft.fftshift(np.fft.fftfreq(n_tr, dt_tr))
    return freq, spec

fid_tr = run_trotter(
    ax_offsets, ax_couplings, 2, estimator, pm_ax
)  # AX  -> two doublets
fid2_tr = run_trotter(
    ax2_offsets, ax2_couplings, 3, estimator, pm_ax2
)  # AX2 -> triplet + doublet

f_ax, s_ax = trotter_spectrum(fid_tr)
f_ax2, s_ax2 = trotter_spectrum(fid2_tr)

fig, ax = plt.subplots(1, 3, figsize=(15, 4))
ax[0].plot(times_tr, np.abs(fid_tr), "o-", label="AX")
ax[0].plot(times_tr, np.abs(fid2_tr), "s-", label="AX$_2$")
ax[0].set_xlabel("evolution time $t$")
ax[0].set_ylabel("|FID|")
ax[0].set_title("Trotter FID envelope (decays only on a noisy backend)")
ax[0].legend()
ax[1].plot(f_ax, s_ax)
ax[1].set_xlim(-12, 12)
ax[1].set_xlabel(r"frequency $\nu$ (Hz)")
ax[1].set_ylabel("intensity")
ax[1].set_title("AX: two doublets")
ax[2].plot(f_ax2, s_ax2)
ax[2].set_xlim(-12, 12)
ax[2].set_xlabel(r"frequency $\nu$ (Hz)")
ax[2].set_title("AX$_2$: 1:2:1 triplet + doublet")
plt.tight_layout()
plt.show()
print(
    f"Deepest AX circuit: {n_tr-1} two-qubit gates;  deepest AX2 circuit: {2*(n_tr-1)} two-qubit gates."
)
print(
    "Noiseless backend: |FID| stays flat, peaks razor-sharp (set APODIZE_TROTTER=True to tame ripples)."
)
print(
    "Noise-model backend / hardware: |FID| decays on its own and the peaks gain a real linewidth."
)

# Close the optional runtime Batch, if one was opened.
if batch is not None:
    batch.close()