תהודה מגנטית גרעינית עם Qiskit

במודול Qiskit in Classrooms הזה, הסטודנטים צריכים סביבת Python עובדת שבה מותקנות החבילות הבאות:

• qiskit v2.1 ומעלה
• qiskit-ibm-runtime v0.40.1 ומעלה
• qiskit-aer v0.17 ומעלה
• scipy ו-matplotlib (להתאמת עקומות ולשרטוט)
• pylatexenc (להצגת דיאגרמות של מעגלים)

כדי להגדיר ולהתקין את החבילות שלמעלה, ראו את המדריך Install Qiskit. כדי להריץ משימות על מחשבים קוונטיים אמיתיים, על הסטודנטים להקים חשבון ב-IBM Quantum® לפי השלבים במדריך Set up your IBM Cloud® account.

המודול הזה נבדק ומשתמש בכ-5 דקות של זמן QPU. זו הערכה, והשימוש בפועל שלכם עשוי להיות שונה.

מבוא

כשכימאי מריץ ניסוי NMR, הוא מניח דגימה של ספינים גרעיניים בשדה מגנטי חזק, מפעיל פולס בתדר רדיו (RF) ומודד את האות שחוזר. הטכניקה נשענת על כמה רעיונות ממכניקת הקוונטים של חלקיק בעל ספין 1/2 בשדה מגנטי. לספין יש שתי רמות אנרגיה המופרדות באנרגיה שמתאימה לתדר לרמור (Larmor) אופייני. פולס בתהודה בתדר הזה מטה את הספין למישור הרוחבי, שם הוא מבצע פרצסיה ומשרה אות שדוהה אחר כך בשני זמנים אופייניים, $T_1$ ו-$T_2$.

qubit מוליך-על הוא גרסה מהונדסת של אותו אובייקט: מערכת קוונטית דו-רמתית שניתנת לשליטה. המצבים $|0\rangle$ ו-$|1\rangle$ שלו הם האנלוג של ספין למעלה וספין למטה, השער שבו אנחנו מסובבים אותו הוא האנלוג של פולס ה-RF, והוא מאבד קוהרנטיות עם $T_1$ ו-$T_2$ משלו. זו לא רק השוואה חופשית; חלק מהמעבדים הקוונטיים הניסיוניים הראשונים, בסוף שנות ה-90 ובתחילת שנות ה-2000, השתמשו ב-NMR.

במודול הזה תשתמשו במחשב קוונטי אמיתי כדי להריץ את אותם ניסויי אפיון שמבצע ספקטרוסקופיסט NMR (או מהנדס קוונטי שמכייל חומרה):

• Rabi / נוטציה — הפעלת הספין וצפייה בו מסתובב. כך מכיילים סיבוב של $\pi$ או $\pi/2$.

• $T_1$ — הרפיית ספין-סריג, הדעיכה של ספין מעורר בחזרה לשיווי משקל.

• $T_2$ — הרפיית ספין-ספין (dephasing), נמדדת בניסוי Ramsey ומתמקדת מחדש בעזרת Hahn echo.

• ספקטרוסקופיית NMR (הסחה כימית וצימוד $J$) — סימולציה של מולקולה קטנה על כמה qubits וקריאת ספקטרום ה-multiplet שלה כדי לשחזר קישוריות באמצעות כלל $n+1$.

עבור כל מושג נציג קודם את התיאוריה המינימלית, ואז נממש את הניסוי ב-Qiskit לפי תהליך העבודה בן ארבעת השלבים של Qiskit patterns.

תיאוריה

חלקיקים בעלי ספין 1/2 בשדה מגנטי

לגרעין בעל ספין 1/2 יש רגע מגנטי $\boldsymbol{\mu} = \gamma \hbar \mathbf{I}$, כאשר $\gamma$ הוא היחס הג'ירומגנטי. בשדה סטטי $B_0$ לאורך $z$, שני מצבי הספין מתפצלים באנרגיה (אפקט זימן):

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

פער האנרגיה מגדיר את תדר לרמור (תדר התהודה) $\omega_0 = \gamma B_0$. הרגע המגנטי מבצע פרצסיה סביב $B_0$ באותו תדר.

כדי לסובב את הספין אנחנו מפעילים שדה RF מתנדנד בתדר $\omega_0$ או קרוב אליו. במערכת הייחוס המסתובבת בתדר ההפעלה (rotating frame), הפעלה בתהודה מסובבת את וקטור Bloch סביב ציר במישור $x$-$y$. כמה רחוק הוא מסתובב (זווית ההטיה) נקבע לפי עוצמת הפולס ומשכו. שתי בחירות חשובות במיוחד: פולס של $90^\circ$ מטה ספין $+1/2$ למישור הרוחבי, ופולס של $180^\circ$ מסובב אותו לגמרי אל $-1/2$.

תרגום בין הז'רגון של NMR לזה של המחשוב הקוונטי

qubit הוא מערכת דו-רמתית, זהה מתמטית לספין 1/2 בשדה מגנטי — ולכן, אף שאנחנו מריצים על מעבד קוונטי ולא על גרעין, הפיזיקה היא אותה פיזיקה שעומדת בבסיס NMR. רק המילון שונה. הטבלה הזו עוזרת לדבר עם כימאים פיזיקליים ועם מומחי מחשוב קוונטי כאחד.

NMR / ספין 1/2	מעבד קוונטי (qubit)
ספין למעלה / ספין למטה	$\ket{0}$ / $\ket{1}$
פער זימן $\Delta E = \hbar\omega_0$	אנרגיית המעבר של ה-qubit $\hbar\omega_{01}$
תדר לרמור $\omega_0$	תדר ה-qubit $\omega_{01}$
פולס RF בתהודה (זווית הטיה $\theta$)	שער סיבוב על 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$
נוטציה	תנודות Rabi
Free induction decay (FID), Hahn echo	רצפי Ramsey ו-echo

כל מה שלהלן הוא תרגיל במדידת העמודה הימנית וקריאתה כעמודה השמאלית.

הגדרת סביבת Qiskit

נעבוד לאורך כל הדרך לפי תהליך העבודה של Qiskit patterns:

1. למפות את הבעיה למעגל קוונטי.

2. לבצע אופטימיזציה של המעגל להרצה (transpile ל-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 (נוטציה)

תיאוריה

אם מפעילים ספין בתהודה, וקטור 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$ של שער $R_X(\theta)$. העקומה והזיהוי של $\pi$/$\pi/2$ זהים; ההבדל היחיד הוא שאנחנו בוחרים זווית ישירות במקום לכייל אורך פולס. בחומרה עם fractional gates, $R_X(\theta)$ עבור $\theta$ שרירותי רץ באופן native.

מבט על כדור Bloch

האנימציה שלמטה מציגה את הנוטציה: הפעלה בתהודה מסובבת בקביעות את וקטור Bloch סביב ציר $x$, ומעבירה אותו מ-$|0\rangle$ כלפי $|1\rangle$ וחזרה. זו בדיוק התנועה שעקומת Rabi מודדת.

מיפוי הבעיה למעגל

נבנה מעגל אחד עם פרמטרים, $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 למעגל ארכיטקטורת קבוצת הוראות (ISA) עבור ה-backend שנבחר.

rabi_isa = pm_q.run(rabi)

הרצה

אנחנו סורקים את $\theta$ מ-$0$ עד $4\pi$ (שני מחזורי Rabi שלמים). העברת מערך של ערכי פרמטרים כ-Primitive Unified Bloc (PUB) אחד מריצה את כולם בג'וב אחד.

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

עיבוד לאחר ההרצה וניתוח

מחלצים את $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$). עכשיו נניח שהמכשיר מכויל בצורה לא נכונה, כך שזווית נומינלית $\theta$ על ציר ה-x מייצרת בפועל סיבוב של $1.05\,\theta$.

(א) באיזו זווית נומינלית יופיע עכשיו המקסימום הראשון של $P_1$? (ב) האם הפעלת פולס $\pi$ יחיד היא דרך רגישה לזהות את השגיאה הזו של 5%?

תשובה

(א) המקסימום מופיע במקום שבו הסיבוב בפועל שווה ל-$\pi$, כלומר בזווית הנומינלית $\theta = \pi/1.05 \approx 0.95\pi$ — הזזה קטנה של השיא שמאלה.

(ב) לא. פולס $\pi$ נומינלי מייצר סיבוב בפועל של $1.05\pi$, ונותן $P_1 = \sin^2(1.05\pi/2) \approx 0.994$ במקום $1$ — שינוי קטן מדי מכדי לראות אותו על רקע רעש המדידה. לכן, כדי לכייל בדיוק, מדענים שעובדים עם ספקטרוסקופיית NMR מגבירים את השגיאה על ידי הפעלת הסיבוב פעמים רבות (למשל שיטת ה-$360^\circ$, או רצף פולסים ארוך), כך ששגיאה קטנה בכל פולס מצטברת לסטייה גדולה שאפשר למדוד. נחזור לרעיון הזה בבעיות האתגר.

ניסוי: $T_1$ (הרפיית ספין-סריג)

תיאוריה

אם משאירים ספין מעורר לנפשו, הוא מוסר אנרגיה לסביבה שלו (ה"סריג") ומתרפה בחזרה לעבר שיווי משקל תרמי. אוכלוסיית המצב המעורר דועכת אקספוננציאלית עם זמן הרפיית ספין-סריג $T_1$:

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

עבור כימאי, $T_1$ קובע את עיכוב המחזור שצריך להמתין בין סריקות ב-NMR כמותי, והוא מעיד על התהפכות מולקולרית ועל דינמיקה. בקיוביט הוא זמן הרפיית האנרגיה, ואחד משני המספרים המרכזיים שמשמשים למידול השפעות הדקוהרנטיות במחשב קוונטי.

הניסוי. מכינים את $|1\rangle$ באמצעות שער $X$ (פולס $\pi$), ממתינים השהיה משתנה $\tau$ ואז מודדים. החלק שעדיין נמצא ב-$|1\rangle$ משרטט את הדעיכה.

מבט על כדור בלוך

הרפיית $T_1$ מחזירה את המצב המעורר לשיווי משקל, ולכן וקטור בלוך גולש מ-$|1\rangle$ חזרה אל $|0\rangle$ לאורך ציר ה-$z$.

בנייה ו-Transpile של סריקת השהיות

משכי ההשהיה חייבים להיות ערכים קונקרטיים (הם מתוזמנים בזמן החומרה), ולכן אנחנו בונים מעגל אחד לכל השהיה.

# 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

הרצה

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 ל-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$ מודד אובדן קוהרנטיות פאזה במישור הרוחבי. אפשר לצפות בזה באמצעות ניסוי Ramsey: פולס $\pi/2$ מסובב את כל הספינים אל קו המשווה, וכולם מצביעים לאותו כיוון. הם עושים פרצסיית Larmor, אבל כל אחד מהם מבצע פרצסיה בתדר מעט שונה, ולכן הם מתפזרים לאורך קו המשווה והאות הרוחבי הכולל דועך, עם זמן דעיכה אופייני של $T_2^*$.

הערה לקוראים מעולם הקוונטים: תמונת ההתפזרות מראה דרך אחת לאבד קוהרנטיות: פיזור של פאזות מוגדרות שמתמצעות החוצה. כדאי לדעת שזה לא המנגנון היחיד. אפשר לאבד קוהרנטיות גם כי הפאזה הופכת למשוכלת עם הסביבה, ואז היא כבר אינה ערך מוגדר שהמערכת מחזיקה בכלל. אי אפשר להבחין בין שני המנגנונים רק על ידי התבוננות בדעיכת האות, ולכן תמונת ההתפזרות היא אינטואיציה טובה. היזהרו לא לקרוא אותה כסיפור המלא או המילולי.

שתי השפעות גורמות לספין לאבד פאזה:

• דקוהרנטיות בלתי הפיכה (או בלתי הפיכה למעשה) שנובעת מצימוד בלתי מבוקר לדרגות חופש אחרות, ומאופיינת על ידי $T_2$ האמיתי. ב-NMR זה נובע בדרך כלל מתנועה מולקולרית שמאפננת את הצימודים בין ספינים. במחשוב קוונטי זה מתאים לצימודים לא רצויים בין קיוביטים לסביבה.
• היסטי תדר הפיכים, סטטיים או משתנים לאט, כמו אי-הומוגניות של השדה, שונות בתדר בין קיוביטים ושגיאות כיול. בשילוב עם $T_2$ האמיתי הם מייצרים את זמן הדעיכה הקצר יותר, $T_2^*$, שנצפה בניסוי Ramsey.

אז איך משחזרים את $T_2$ האמיתי אם ניסוי Ramsey מספר לנו רק על $T_2^*$? הניסוי השני, שנקרא הד Hahn, מבטל את כל הרעש מהסוג השני כדי למצוא את $T_2$ האמיתי. הוא מתחיל כמו ניסוי Ramsey, עם פולס $\pi/2$, ואז הספינים מתפזרים כמקודם. אחרי זמן $\tau$ מפעילים פולס $\pi$. הפולס הופך את צבירת הפאזה שנגרמה מהיסטי תדר סטטיים, ולכן הספינים מתמקדים מחדש אחרי השהיה שנייה באותה משך, $\tau$.

$T_2^*$ ו-$T_2$ הם שניהם מספרים חשובים לדעת, הן בהקשר של הכימיה והן בהקשר של המחשוב הקוונטי. בהקשר הכימי, $T_2^*$ קובע את רוחב הקו בניסוי NMR רגיל (תראו מה זה אומר בסעיף מאוחר יותר). $T_2$ יכול לספר הרבה על המולקולות, כמו התהפכות מולקולרית, ושינוי ב-$T_2$ יכול לסמן שמולקולה נקשרה למשהו גדול יותר. בהקשר של מחשוב קוונטי, אנחנו רוצים בדרך כלל ששניהם, $T_2$ ו-$T_2^*$, יהיו ארוכים ככל האפשר, כדי שיהיה לנו יותר זמן לבצע ניסויים לפני שהדקוהרנטיות משתלטת על המערכת. משפחה שלמה של טכניקות המבוססות על הד Hahn מאריכה את הזמן הזה על ידי ביטול הרעש ההפיך, כך שהמעגל מוגבל באופן אידיאלי על ידי $T_2$ האמיתי ולא על ידי $T_2^*$. זה נקרא "דינמיקל דקאפלינג".

מימוש

כשמממשים את ניסויי Ramsey והד Hahn על קיוביט, הקריאה שונה במקצת מהמקרה של NMR. ב-NMR אפשר לצפות ישירות בשדות המתנדנדים שמושרים על ידי הספינים שמבצעים פרצסיה בתדר Larmor במישור הרוחבי. אבל בקיוביט אפשר למדוד רק את ההטלה שלו לאורך ציר $Z$, כלומר אם הוא $\vert 0 \rangle$ או $\vert 1 \rangle$. לכן, בשני הניסויים אנחנו מסיימים בפולס $\pi/2$ נוסף ואחריו מדידה. פולס $\pi/2$ הזה מתרגם את מידע הפאזה של הקיוביט לאוכלוסייה.

בהמשך אנחנו מריצים כל ניסוי כדי לקבוע את $T_2^*$ ואת $T_2$.

Ramsey ‏($T_2^{*}$): שני פולסי $\pi/2$ עם זמן השהיה $\tau$

פולס $\pi/2$ הראשון מטה את וקטור בלוך אל קו המשווה, והפולס הבא מתרגם את הפאזה של וקטור בלוך לאוכלוסייה. מכיוון שהפאזה משתנה בזמן בזמן שוקטור בלוך מבצע פרצסיה, אנחנו רואים תבנית סינוסואידלית באוכלוסייה כפונקציה של זמן הפרצסיה $\tau$. כשהקיוביט מאבד קוהרנטיות פאזה, אנחנו אמורים לראות את האמפליטודה של עקומת הסינוס הזו דועכת.

עוד כמה דקויות לקוראים מעולם הקוונטים: בניסויי מחשוב קוונטי, התנודות שאנחנו רואים אינן בעצם בתדר Larmor, שנמצא בתחום ה-GHz ומהיר מדי מכדי שנוכל לפתור אותו. ניסויים קוונטיים פועלים במערכת הייחוס המסתובבת: מערכת ייחוס שמסתובבת יחד עם הקיוביט בתדר Larmor (או קרוב אליו), כך שהפרצסיה המהירה בלתי נראית ווקטור בלוך נראה כמעט נייח. כדי לראות שוליים, אנחנו מכוונים בכוונה את מערכת הייחוס המסתובבת לתדר מעט שונה מתדר Larmor. וקטור בלוך נראה אז כמבצע פרצסיה בסטייה הקטנה הזו, ואנחנו מקבלים אות מתנודד בתדר הסטייה. בהמשך, איננו קובעים את הסטייה בערך אחד; אנחנו בוחרים אותה כך שכתריסר שוליים יתפרשו על פני חלון המדידה (שאת גודלו אנחנו קובעים לפי זמן הקוהרנטיות של הקיוביט), מה שמקל על דגימת התנודה ללא aliasing. סיבוך נוסף: ב-Qiskit אין לנו את שליטת ה-pulse-level הדרושה כדי לשנות את תדר הפולס (שהוא זה שקובע את תדר מערכת הייחוס המסתובבת). במקום זאת, אנחנו מוסיפים סטייה אפקטיבית ידנית, ומשתמשים בה כדי לקבוע את הפאזה של פולס Ramsey השני. זה קצת מלאכותי, אבל התוצאה זהה והדעיכה שאנחנו צופים בה היא עדיין $T_2^*$ אמיתי.

מבט על כדור בלוך

פולס $\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 ‏($T_2$): הוספת פולס $\pi$ מתמקד באמצע

מבט על כדור בלוך

הד Hahn: הספינים מתפזרים כמקודם, ואז פולס $\pi$ הופך אותם, כך שהפרצסיה המתמשכת שלהם ממקדת אותם מחדש יחד (הוקטור האדום גדל שוב), ומשחזרת את הקוהרנטיות שאבדה בגלל איבוד פאזה סטטי.

# 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, כמצופה. סולמות הזמן של הדעיכות האלה הם $T_2^*$ ו-$T_2$ בהתאמה.

בדקו את ההבנה שלכם

בדרך כלל תמצאו ש-$T_2 \gtrsim T_2^{*}$. הסבירו מדוע במשפט אחד, וציינו איזו השפעה פיזיקלית פולס ה-$\pi$ של הד Hahn מסיר, ומה הוא אינו יכול להסיר.

תשובה

$T_2^{*}$ כולל גם איבוד פאזה בלתי הפיך וגם פיזור סטטי הפיך, ולכן הוא דועך מהר יותר; פולס ה-$\pi$ של ההד ממקד מחדש את איבוד הפאזה הסטטי וההפיך (אי-הומוגניות של השדה / היסטי תדר קבועים), ומשאיר רק את איבוד הפאזה הבלתי הפיך של ספין-ספין, שהוא $T_2$ האמיתי. ההד אינו יכול לבטל תנודות שמשתנות בין שני חצאי הרצף.

ספקטרוסקופיית NMR של מולקולה: הזזה כימית וצימוד $J$

מספין בודד למולקולה

הניסויים עד כה אפיינו ספין בודד ומבודד. NMR כימי אמיתי קורא מולקולה (הרבה ספינים מצומדים בו זמנית) ומשתמש בספקטרום המתקבל כדי להסיק תכונות של מבנה המולקולה. זה יכול לספר לנו על שתי השפעות:

• הזזה כימית. כל גרעין שונה כימית מבצע פרצסיה בתדר מעט שונה $\nu_i$, כי האלקטרונים המקומיים שלו מסוככים על השדה המופעל בכמויות שונות. לכן מיקומי הפסגות מזהים סביבות כימיות (למשל פרוטון של CH$_3$ לעומת פרוטון ארומטי).

• צימוד סקלרי $J$. אינטראקציה דרך קשרים בין שני ספינים מפצלת כל פסגה לרב-קו. גרעין עם $n$ שכנים שקולים מתפצל ל-$n+1$ קווים בעוצמות בינומיות, כלל ה-$n+1$. תבנית הפיצול מגלה קישוריות (אילו גרעינים קשורים זה לזה), והמרווח הוא קבוע הצימוד $J$.

הניסוי

כדי למדל זאת במחשב הקוונטי שלנו, אנחנו צריכים לממש המילטוניאן ספינים על הקיוביטים. בגבול הצימוד החלש והשדה החזק (הפרשי ההזזה הכימית $\gg$ הצימודים), החלקים $XX$ ו-$YY$ של הצימוד מתמצעים החוצה ונשאר רק $ZZ$. לכן ההמילטוניאן הוא (ביחידות של תדר זוויתי, $\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}}.$

כל גרעין הוא קיוביט אחד: ההזזה הכימית היא סיבוב $Z$ של קיוביט יחיד, והצימוד הוא אינטראקציית $ZZ$ של שני קיוביטים.

הניסוי מתנהל כך. פולס $\pi/2$ מטה כל ספין אל המישור הרוחבי, והספינים מבצעים פרצסיה ומצטמדים במשך זמן $t$. אנחנו מתעדים את המגנטיזציה הרוחבית המרוכבת $s(t) = \langle X\rangle + i\langle Y\rangle$ המסוכמת על כל הספינים. כמו בניסוי Ramsey שלמעלה, המגנטיזציה הרוחבית הזו משרטטת סינוסואיד דועך כפונקציה של זמן הפרצסיה $t$. כשמבצעים התמרת פורייה של האות, מקבלים ספקטרום עם פסגות בהזזה הכימית $\nu_i$ המפוצלות על ידי הצימוד הסקלרי $J_{ij}$.

הערה

אנחנו מדמים את התפתחות המולקולה על ידי מיפוי הזמן $t$ לזוויות סיבוב: $t$ גדול יותר פירושו סיבובי $Z$ ו-$ZZ$ גדולים יותר, לא מעגל ארוך או עמוק יותר. לכן, אף שהמצב המדומה אכן מתקדם עם $t$ (וזה מה שגורם לדעיכת האינדוקציה החופשית (FID) להתנדנד), הקיוביטים הפיזיים אף פעם לא ממתינים זמן רב יותר ואין דקוהרנטיות התלויה ב-$t$. לכן ה-FID המדומה אף פעם לא דועך מעצמו. עבור FID אידיאלי ארוך לאין שיעור, קווי הספקטרום היו חדים לאין שיעור; חלון הדגימה הסופי שלנו נותן להם במקום זאת רזולוציה סופית ועלול להכניס דליפה ספקטרלית. לקווי NMR אמיתיים יש רוחב סופי שנקבע על ידי $T_2^*$. אנחנו מוסיפים אותו בחזרה על ידי הכפלת ה-FID בחלון $e^{-t/T_2^{*}}$ לפני התמרת פורייה. זה נקרא אפודיזציה להרחבת קו, שספקטרוסקופיסטים מפעילים בעיבוד כדי לקבוע את רוחב הקו (תוך החלפת רזולוציה ביחס אות לרעש). (תא אופציונלי בנספח של המודול הזה מציג חלופה מבוססת Trotter, שבה הדעיכה נובעת במקום זאת מדקוהרנטיות אמיתית שגדלה עם עומק המעגל.)

מערכת AX

נתחיל במערכת AX, כביכול, מולקולה המורכבת משני ספינים מצומדים השונים זה מזה כימית. התוויות A ו-X מציינות את שתי הסביבות הכימיות השונות של הספינים. האותיות שרירותיות, אבל העובדה שהן רחוקות זו מזו באלפבית ממחישה את השוני הכימי ביניהן.

מיפוי ואופטימיזציה

אנחנו כותבים פונקציה אחת שבהינתן היסטי ההזזה הכימית $\nu_i$ והצימודים $J_{ij}$, מחזירה מעגל הכנת מצב עם פרמטרים. פולס ה-$\pi/2$ מטה את הספינים אל קו המשווה, ואז סיבובי $Z$ מדמים את התפתחות המצב תחת איברי ההזזה הכימית והצימוד $J$. כל $\nu_i$ הוא היסט פרצסיה קטן ב-Hz ביחס להתייחסות נבחרת, או נושא של מערכת ייחוס מסתובבת. כאן אנחנו קובעים אותו ידנית; במולקולה אמיתית, סיכוך אלקטרוני משנה את השדה המגנטי המקומי ומזיז את התהודה של כל גרעין ביחס להתייחסות הזו. זה מודל על ידי שערי $R_Z$ של קיוביט יחיד, אחד לכל ספין, וכל $\nu_i$ קובע היכן נמצאת הפסגה של אותו ספין בספקטרום. שער $ZZ$ של שני קיוביטים מפעיל את הצימוד $J_{ij}$, שמפצל כל פסגה לרב-קו וקובע את המרווח בין הקווים שלו.

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")

הרצה ועיבוד לאחר ההרצה

אנחנו מעריכים את המעגל על פני רשת של זמני התפתחות, וה-Estimator מחזיר את $\langle X\rangle$ ואת $\langle Y\rangle$ בכל אחד מהם. אנחנו משלבים אותם ל-FID המרוכב $\langle X\rangle + i\langle Y\rangle$, מפעילים את חלון האפודיזציה, ומבצעים התמרת פורייה של התוצאה כדי לקבל את הספקטרום.

# --- 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 הוא רק תווית לשני הספינים השקולים; אין לה קשר לאופרטורי פאולי $X$. כל הצימודים בהמשך הם עדיין $ZZ$.) ספין A רואה אז ארבע תצורות שכנים $(\uparrow\uparrow,\ \uparrow\downarrow,\ \downarrow\uparrow,\ \downarrow\downarrow)$; שתי התצורות המעורבות מתלכדות, ולכן A מתפצל לשלשה 1:2:1 ($n=2 \Rightarrow n+1 = 3$ קווים). כל ספין X רואה רק את A ונשאר דובלט. ספינים שקולים אינם מפצלים זה את זה, ולכן אנחנו לא כוללים צימוד X–X.

מולקולה אמיתית עם התבנית הזו היא 1,1,2-טריכלורואתאן (Cl$_2$CH–CH$_2$Cl). כשמכוונים לספיני המימן (הצורה הנפוצה ביותר של NMR), פרוטון המימן הבודד ב-CH ממלא את תפקיד A, ושני פרוטוני CH$_2$ השקולים הם שני ספיני X. עכשיו פשוט ספרו שכנים: לפרוטון CH יש שני שכנים (זוג ה-CH$_2$), ולכן האות שלו מתפצל לשלשה 1:2:1 ($n=2$); לכל פרוטון CH$_2$ יש רק שכן אחד (ה-CH), ולכן אות ה-CH$_2$ הוא דובלט 1:1 ($n=1$). שימו לב שהקבוצה עם פחות פרוטונים (CH) מראה יותר קווים; מספר הקווים משקף כמה שכנים יש לאות, לא כמה פרוטונים חולקים אותו.

בהמשך אנחנו מדמים מערכת AX$_2$ עם שלושה קיוביטים.

# 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 שני הדובלטים מופרדים (מרכז למרכז) ב-12 Hz, ושני הקווים בתוך כל דובלט נמצאים במרחק 2 Hz זה מזה. איזה מספר הוא הפרש ההזזה הכימית ואיזה הוא קבוע הצימוד $J$? אם הייתם עוברים לספקטרומטר עם שדה מגנטי כפול, איזה מהם היה משתנה?

תשובה

ההפרדה של 12 Hz מרכז למרכז היא הפרש ההזזה הכימית ($\nu_1-\nu_2$); מרווח הקווים של 2 Hz הוא קבוע הצימוד $J$. הזזות כימיות נובעות מסיכוך של $B_0$, ולכן הן גדלות עם השדה. בשדה כפול ההפרדה מכפילה את עצמה ל-24 Hz. הצימוד $J$ הוא תכונה מולקולרית דרך קשרים והוא בלתי תלוי בשדה, ולכן נשאר 2 Hz. כימאים מבחינים בין הזזה לצימוד על ידי התבוננות באופן שבו פיצול משתנה עם השדה.

בעזרת כלל ה-$n+1$, נבאו את הרב-קו של A אילו הוספנו ספין X שקול שלישי (מערכת AX$_3$), והסבירו מדוע שני ספיני ה-X השקולים ב-AX$_2$ אינם מפצלים זה את זה.

תשובה

עם $n=3$ שכנים שקולים, A הופך לרבעיה 1:3:3:1. ספינים שקולים אינם מפצלים זה את זה, כי הצימוד ביניהם מזיז את רמות האנרגיה הרלוונטיות כולן באותה כמות, ולכן הפערים בין הרמות האלה — התדרים שאנחנו רואים — אינם משתנים.

סיכום

השתמשנו במחשב קוונטי אמיתי כדי לחקור כיצד ספינים גרעיניים מבצעים פרצסיה ומאבדים קוהרנטיות בניסוי NMR, וכיצד לפרש את הספקטרום שמתקבל. כדאי לעצור ולסקור מדוע זה היה אפשרי בכלל, ולמה זה יותר מסתם קוריוז.

קיוביט וגרעין ספין-1/2 הם, מבחינה מתמטית, אותו אובייקט: מערכות קוונטיות דו-רמתיות החיות על כדור בלוך, המסובבות על ידי פולסים תהודתיים, מבצעות פרצסיה סביב ציר ומאבדות קוהרנטיות לסביבתן. כל אוצר המילים של המודול הזה (תדר Larmor, מערכת הייחוס המסתובבת, $T_1$, $T_2$ ו-$T_2^*$, שולי Ramsey, הד Hahn) שייך באותה מידה ל-NMR ולמחשוב קוונטי. מה שכימאי מכנה הזזה כימית, מהנדס קוונטי מכנה סיבוב $Z$ של קיוביט יחיד; מה שכימאי מכנה צימוד $J$, מהנדס קוונטי מכנה אינטראקציית $ZZ$. שני התחומים תיארו את אותה פיזיקה במילים שונות.

הקשר הזה הוא היסטורי, לא רק מושגי: NMR היה אחת הפלטפורמות הראשונות שהריצו אלגוריתמים קוונטיים אמיתיים. עד 2001, מולקולה בתמיסה הריצה את האלגוריתם של Shor כדי לפרק את 15 לגורמים, כשכל ספין גרעיני משמש כקיוביט. NMR גם העניק למחשוב הקוונטי חלק גדול מערכת הכלים לשליטה שלו. הד הספין, רצפי פולסים של CPMG, פולסים מורכבים ו-dynamical decoupling כולם הומצאו על ידי ספקטרוסקופיסטים עשרות שנים לפני שקיוביטים היו זקוקים להם כדי להילחם בדקוהרנטיות. יש סימטריה מספקת במה שעשינו זה עתה: התחום שהוליד את המחשבים הקוונטיים הראשונים הוא עכשיו משהו שמחשב קוונטי יכול לדמות.

סקירה ושאלות

מושגים קריטיים

• קיוביט הוא מערכת דו-רמתית מהונדסת השקולה לספין-1/2; ה-$|0\rangle/|1\rangle$ שלו, התדר, השערים, $T_1$ ו-$T_2$ הם האנלוגים הישירים של כמויות ה-NMR.

• פולס תהודתי מסובב את וקטור בלוך; זווית ההטיה קובעת את הסיבוב, ועקומת Rabi/נוטציה $P_1=\sin^2(\theta/2)$ מזהה סיבובי $\pi$ ו-$\pi/2$.

• $T_1$ הוא הרפיית אנרגיה (ספין-סריג); $T_2$ הוא הרפיית פאזה (ספין-ספין); $T_2^{*}$ כולל בנוסף איבוד פאזה סטטי הפיך שהד Hahn ממקד מחדש.

• מולקולה היא קבוצה של ספינים מצומדים: ההזזה הכימית קובעת את מיקום כל פסגה, וצימוד $J$ מפצל פסגות לרב-קווים לפי כלל ה-$n+1$, ויחד הם מקודדים סביבה כימית וקישוריות.

נכון / לא נכון

1. נ/ל תדר Larmor הוא פער האנרגיה בין שני מצבי הספין חלקי $\hbar$.

2. נ/ל פולס $90^\circ$ ‏($\pi/2$) הופך לגמרי ספין בשיווי משקל.

3. נ/ל $T_2^{*}$ ארוך בדרך כלל מ-$T_2$.

4. נ/ל הד Hahn ממקד מחדש איבוד פאזה בלתי הפיך של ספין-ספין.

5. נ/ל במערכת AX$_2$ גרעין A מופיע כשלשה 1:2:1.

6. נ/ל קבוע הצימוד $J$ מוכפל כשמכפילים את השדה המגנטי של הספקטרומטר.

7. נ/ל הסרת Qiskit Pulse הופכת את מדידת $T_1$ בחומרה הנוכחית לבלתי אפשרית.

תשובה קצרה

1. הסבירו בשני משפטים מדוע פולס $\pi$ יחיד הוא דרך גרועה לזהות שגיאת זווית הטיה של $5\%$, וכיצד שרשור פולסים פותר זאת.

2. דגימה מראה קווי NMR רחבים. במונחים של $T_2$, מה זה אומר לכם, ואיזו כמות של קיוביט היא האנלוגית לכך?

3. תנו את האנלוג בקיוביט של כל אחד מאלה: (א) עיכוב המחזור, (ב) פולס ה-$90^\circ$, (ג) ההזזה הכימית.

4. בעזרת כלל ה-$n+1$, הסבירו מדוע קבוצת אתיל (CH$_3$CH$_2$–) מראה שלשה 1:2:1 ורבעיה 1:3:3:1, וזהו איזו קבוצה מייצרת איזה רב-קו.

בעיות אתגר

1. רצף Carr-Purcell. שנו את הד Hahn כך שיוכנסו $n$ פולסי $\pi$ במרווחים שווים (רצף בסגנון CPMG) ומדדו את $T_2$ כפונקציה של $n$. דונו מדוע יותר פולסי מיקוד מחדש יכולים להאריך את $T_2$ הנראה.

2. הגברת שגיאת פולס. בנו מעגל שמפעיל את $R_X(\pi)$ ‏$\,n$ פעמים על $|0\rangle$ ושרטטו את $P_1$ מול $n$ עבור זווית שהוגדרה בכוונה באופן שגוי $R_X(1.05\pi)$. הראו כיצד השגיאה הופכת לגלויה ככל ש-$n$ גדל. רמז: הוסיפו barrier אחרי כל סיבוב; אחרת ה-transpiler עלול לשלב את הסיבובים החוזרים לשער שקול יחיד.

3. קבוצת אתיל. בעזרת fid_circuit, בנו את ספקטרום הפרוטונים של שבר אתיל (CH$_3$–CH$_2$): חמישה קיוביטים, כששלושת פרוטוני CH$_3$ חולקים הזזה כימית אחת ושני פרוטוני CH$_2$ חולקים אחרת, וכל זוג CH$_3$–CH$_2$ מצומד באותו $J$ (ללא צימוד בתוך קבוצה שקולה). נבאו את התבנית בעזרת כלל ה-$n+1$, ואז דמו אותה; אתם אמורים לקבל רבעיה 1:3:3:1 (CH$_2$) ושלשה 1:2:1 (CH$_3$).

נספח (אופציונלי): התפתחות מבוססת Trotter — לתת לדקוהרנטיות לקבוע את רוחב הקו

בסעיף ספקטרוסקופיית NMR של מולקולה, ההתפתחות $e^{-iHt}$ מופעלת בשכבה אחת שעומקה אינו גדל עם $t$, ולכן ה-FID אינו דועך מעצמו ואנחנו כופים את רוחב הקו בעזרת חלון האפודיזציה $e^{-t/T_2}$. התא האופציונלי הזה נוקט בגישה ההפוכה. הוא מפרק את ההתפתחות ל-$p$ צעדי Trotter ומציב barrier אחרי כל צעד (כך שה-transpiler לא יוכל למזג אותם בחזרה לשכבה אחת). נקודת זמן $k$ משתמשת ב-$k$ צעדים, ולכן עומק המעגל גדל עם הזמן.

מכיוון שכל האיברים ב-$H$ מתחלפים, זה מדויק: בסימולטור נטול רעש מתקבלות אותן פסגות חדות. אבל ב-backend עם מודל רעש או בחומרה אמיתית, שגיאת השערים המצטברת גורמת ל-FID לדעוך מעצמו, והדעיכה הזו היא שקובעת את רוחב הקו. זהו רעיון ה-$T_2\leftrightarrow$רוחב הקו מהסעיף ניסוי $T_2$, $T_2^*$, אך עכשיו הוא נובע מדקוהרנטיות אמיתית ולא מחלון עיבוד. אנחנו מריצים כך גם את מערכת AX (שני דובלטים) וגם את מערכת AX$_2$ (שלשה 1:2:1 ועוד דובלט).

שתי הערות מעשיות לשים לב אליהן. ראשית, המחיר הוא עומק: מעגל ה-AX העמוק ביותר משתמש בכ-$\sim\!47$ שערי שני קיוביטים וגרסת AX$_2$ בערך בכפל מזה, ולכן שניהם רועשים ויקרים בהרבה מהגרסה המדויקת. שנית, כל מערכת ספינים עדיין נשלחת כג'וב אחד, וכל נקודת זמן משתמשת באותה פריסת קיוביטים פיזיים שנבחרה בסעיף הספקטרוסקופיה. אם USE_BATCH=True, הג'ובים האלה משתמשים גם ב-Batch האופציונלי של כל המחברת.

# 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()