第2章: NVセンター磁気計測

スピン依存蛍光、ODMR、フィルター関数、そして磁性体の走査プローブとしてのダイヤモンド

📖 読了時間: 45-50分 📊 難易度: 中級 💻 コード例: 7個 📝 演習問題: 5問

🌐 JP | 🇬🇧 EN | Last sync: 2026-08-13

材料科学基礎道場 > 量子センシング入門 > 第2章

炭素原子があるべき場所に窒素原子が座り、その隣の格子点が空いている。窒素-空孔中心とはこれだけのものです。ダイヤモンド中の隣接する2つの点欠陥であり、結晶育成者が一生を費やして排除しようとする類の欠陥です。第1章は量子センシングの一般的な道具立て — 位相の積分器としての2準位系、射影ノイズ、磁場毎ルートヘルツを単位とする感度 $\eta$、そして環境のどの周波数をパルス系列が見られるかを決めるフィルター関数 — を組み立てました。本章はその道具立てを欠陥の中に据え付けます。

第3章のSQUIDや第4章の原子系ではなくここから始める理由は、姉妹コース 量子ハードウェア入門 の読者ならすぐ気づくやり方で、NV中心が材料科学と量子技術の関係を反転させることにあります。あちらでは、どの方式でもコヒーレンス時間を制限していたのは欠陥でした。表面の不対スピン、酸化物中の2準位揺動子、チャネル中の $^{29}\mathrm{Si}$ 核です。ここでは欠陥そのものが装置です。欠陥を量子ビットにとって厄介者にしている物理 — 局所的な場に結合し、近傍にあるものに応答し、そのスペクトルが周囲を報告する — が、そのまま計測器としての性質になります。しかもこの計測器は結晶表面の数十ナノメートル下に埋め込まれた原子スケールの単一物体なので、他のどんな磁力計よりも試料に近づけられます。

以下は1つの問いを軸に構成されています。材料中の磁気現象 — 磁壁、単層磁性体、電流分布、ゆらぐスピンのバス — が与えられたとき、それを見るには何が必要か。答えは常に感度、帯域、そして試料までの距離であり、その3つはいずれもこの欠陥の具体的な準位構造に第1章の道具立てを適用して出てきます。7つの実例が、ハミルトニアンから画像までの計算を組み立てます。

学習目標

本章を修了すると、以下のことができるようになります：

記法

本章は第1章で固定した記法をそのまま使い、繰り返しません。以下の3つは以降のすべての数値が依拠するので、思い出しておく価値があります。

感度。 $\eta$ は帯域1 Hzにおけるノイズの磁場換算値で、単位は T/$\sqrt{\mathrm{Hz}}$、時間 $t$ 平均すると不確かさが $\eta/\sqrt{t}$ になるように定義されます。小さいほど良く、帯域を伴わない感度の数値は数値ではありません。

コヒーレンス時間。 $T_1$、$T_2$、$T_2^\ast$ は 量子ハードウェア入門第1章 と厳密に同じ意味で、限界 $T_2 \le 2T_1$、$T_2^\ast$ が遷移周波数の静的なばらつきであること、フィルター関数 $|\tilde{s}(f,t)|^2$ の定義も含みます。Ramsey、Hahnエコー、CPMG-$N$ の各系列もそこで定義したものです。再導出はせず、閉形式を引用して Code Example 4 で数値的に検証します。本コースの第1章は最適時間と最良感度を $\tau_{\mathrm{opt}} = T_2/(2p)^{1/p}$、$\eta^{\min} = \sqrt{2e}/(\gamma\sqrt{N T_2})$ の形で固定しました。第2.2節はNVの記法で両方を再導出し、項ごとに一致します。数値的には、2つの章は $1.7\times10^{-6}$ で一致し、それ以上ではありません。第1章の $T_2^\ast = 1\ \mu$s における単一スピンの値は厳密な $\gamma_e = 1.76085963\times10^{11}$ rad s$^{-1}$T$^{-1}$ から13.241510 nT/$\sqrt{\mathrm{Hz}}$、本章は本章のすべての表が引用する丸めた $\gamma_e/2\pi = 28.025$ GHz/T から13.241487 nT/$\sqrt{\mathrm{Hz}}$ を計算し、$28.025/28.024951 - 1 = 1.7\times10^{-6}$ が差の全部を説明します。Code Example 3 の中で報告される6桁の一致は別の検証 — 本章の閉形式を本章の数値と、1つの $\gamma$ の値で照合したもの — であり、章をまたぐ照合について何も述べていません。

磁気回転比。 2つの記号を混同すると $2\pi$ の因子を失います。円振動数でない比 $\gamma_e/2\pi = 28.025$ GHz/T は磁場を周波数シフトに変換し、角振動数の比 $\gamma_e = 1.7609 \times 10^{11}$ rad s$^{-1}$T$^{-1}$ は黄金律の速度に現れます。コード中では Code Example 1 から 5 で前者を GAMMA、角振動数側も必要になる Code Example 6 では GAMMA_C と呼びます。

2.1 たまたま磁力計である欠陥

その対象

負に帯電した窒素-空孔中心 NV$^-$ は、ダイヤモンド中の置換型窒素と隣接する格子空孔からなり、近傍のドナーから6番目の電子を捕獲しています。対称性は $C_{3v}$ で、対称軸は4本の $\langle 111 \rangle$ 結晶方位のいずれかを向きます。これがすぐ後で効いてきます。欠陥の6電子が残す基底状態はスピン三重項 $S = 1$ であり、$^3A_2$ と表記されます。

三重項は第1章の2準位系ではなく、その違いが有用性のすべての源です。3つの状態は1本ではなく2本の遷移を意味し、それが測定をベクトル的にします。そして三重項は2つの不対電子間の双極子相互作用によってゼロ磁場でも分裂しており、それがこの欠陥を自己参照的な周波数標準にします。

基底状態のハミルトニアン

NV軸を局所的な $z$ とする座標系で、基底状態のスピンハミルトニアンを $h$ で割ると

$$ \frac{H}{h} = D\, S_z^2 + E\left( S_x^2 - S_y^2 \right) + \frac{\gamma_e}{2\pi}\, \mathbf{B}\cdot\mathbf{S} $$

となります。$D = 2.870$ GHz がゼロ磁場分裂、$E$ が局所的な歪みまたは電場による横方向の項、最後の項がゼーマン相互作用です。それぞれの記号が物理的な物語を担っています。

$D$ は欠陥の2つの不対電子間の双極子相互作用で、その値はダイヤモンド格子の幾何で決まります。したがって、あらゆるダイヤモンド中のあらゆるNV中心で $10^4$ 分の1まで同一であり、これが測定をマイクロ波源だけを頼りに校正できる理由です。また室温付近では $\mathrm{d}D/\mathrm{d}T \approx -74$ kHz/K の弱い温度依存性をもちます。熱膨張が電子間距離を変えるからで、これはNV温度計測の基礎であり、同時に磁気計測では安定化するか測って差し引くしかない系統誤差です。$E$ は厄介者とプローブを兼ねています。ゼロ磁場で $m_s = \pm 1$ の対を $2E$ だけ分裂させるので、バイアスのないODMRスペクトルは二重線になります。良質なバルク材料では数MHz、格子が歪んだ表面近くの浅い中心では急激に大きくなります。$E$ のマップは歪みのマップです。

ゼーマン項が計測の対象であり、最低次で寄与するのはその軸方向成分だけです。$\gamma_e B \ll D$ のとき $m_s = 0$ からの2本の遷移周波数は

$$ f_{\pm} = D \pm \frac{\gamma_e}{2\pi} B_{\parallel} + \frac{3}{4}\left(\frac{\gamma_e}{2\pi}\right)^2 B_{\perp}^2 \left[ \frac{1}{D + \gamma_e B_{\parallel}/2\pi} + \frac{1}{D - \gamma_e B_{\parallel}/2\pi} \right] + \ldots $$

です。このように書いた2次項は対の中心について厳密であり、中心についてのみ厳密です。$A_\pm = [D \pm \gamma_e B_\parallel/2\pi]^{-1}$ と略記すると、個々の遷移が担うのは $(\gamma_e/2\pi)^2B_\perp^2(A_+ + A_-/2)$ と $(\gamma_e/2\pi)^2B_\perp^2(A_- + A_+/2)$ であり、その相加平均が上の $\tfrac{3}{4}(A_+ + A_-)$、その差 $\tfrac{1}{2}(\gamma_e/2\pi)^2B_\perp^2(A_+ - A_-)$ が分裂幅への小さなバイアスです。$B_\parallel = B_\perp = 1$ mT では56 MHzの分裂幅に対して $-2.7$ kHz です。したがって $f_+ - f_- = 2(\gamma_e/2\pi)B_\parallel$ が軸方向磁場に線形で他の一切に依存しないのは1次においてであり、この2次の残差がまさに Code Example 1 のベクトル復元をミリテスラ磁場で $10^3$ 分の1に制限します。（この展開の外にある項が1つあります。$B_\parallel \to 0$ では $m_s = \pm 1$ の対が縮退し、$B_\perp$ が $m_s = 0$ を介して両者を混ぜて $(\gamma_e/2\pi)^2B_\perp^2/D$ だけ分裂させます。$B_\perp = 1$ mT で274 kHzです。バイアス磁場なしのスペクトルが磁場ではなく $E$ と歪みのために読まれるのはこのためです。）いずれの場合も横磁場は $B_\perp/D$ で抑制され、ゼロバイアスから離れれば分裂幅よりも中心をはるかに大きくずらします。Code Example 1 は厳密対角化に対して両方の主張を検証し、間違えやすい中心シフトの係数 $3/4$ も確認します。

4つの配向、1つのベクトル

単一のNV中心は1つの射影 $B_\parallel$ を測ります。しかし $\langle 111 \rangle$ ファミリーには4本のメンバーがあり、巨視的なダイヤモンドは4つの部分アンサンブルを等量含みます。スペクトルには8本の共鳴が現れ、4本の分裂幅が読み取れ、既知の4本の軸への4つの射影が3成分ベクトルを過剰決定します。これがベクトル磁気計測の標準的な経路であり、Code Example 1 が実行します。

これには2つの注意が伴います。分裂幅から得られるのは $|B_\parallel|$ だけなので、相対符号は小さな組合せ探索を要し、$\mathbf{B}$ の全体符号は分裂幅だけでは決まりません。既知のバイアス磁場が縮退を破ります。また $B_\perp$ が $D/(\gamma_e/2\pi) \approx 102$ mT に近づくと線形関係は破れ、完全な対角化を数値的に反転する必要があります。

なぜ蛍光がスピン状態を知っているのか

以上はスピンを読む手段がなければ役に立ちません。NV中心の答えこそがこれが使われている理由です。光サイクルがスピン選択的なので、欠陥は自分のスピン状態を可視光で報告するのです。

慣例的に532 nm付近の緑色励起は $^3A_2$ 基底状態を励起三重項 $^3E$ へ引き上げます。この光学遷移はスピン保存であり、吸収でも続く放射崩壊でも $m_s$ は変わりません。この放射崩壊が実際に検出される広帯域の赤色蛍光を生みます。これだけなら蛍光はスピン情報を運びません。

追加の要素は2つの一重項状態 $^1A_1$ と $^1E$ で、2つの三重項の間に位置します。$^3E$ から一重項への項間交差は $m_s = 0$ よりも $m_s = \pm 1$ からのほうがはるかに速く、スピン軌道結合がそれらの状態をより強く結ぶためです。ここから2つの帰結が生じ、それがNV読み出しの両輪です。

スピン依存蛍光。 $m_s = \pm 1$ の中心は優先的に一重項多重項に格納され、検出帯域で発光せずに200 ns程度を過ごします。$m_s = 0$ の中心は放射的にサイクルします。照射の最初の数百ナノ秒で平均すると $m_s = 0$ のほうが明るいわけで、その相対差がODMRコントラスト $C$ です。良質な単一中心のパルス読み出しでは典型的に0.2から0.3、連続照射ではかなり小さくなります。

光ポンピング。 一重項は $m_s$ を保存せずに基底状態へ崩壊し、優先的に $m_s = 0$ に入ります。したがって数マイクロ秒の緑色光は、どこから始めても、そして温度にかかわらず、中心を確率0.8から0.9程度で $m_s = 0$ に置きます。これは散逸的な初期化であり、姉妹コースのDiVincenzo基準の2番であり、NV磁力計が300 Kで動くのにトランズモンが20 mKを必要とする理由です。

このコストを率直に述べておきます。本章のあらゆる感度がこれを負うからです。読み出しは射影的でも効率的でもありません。1ショットで得られる検出光子は $0.01$ から $0.1$ 程度なので、測定は $10^4$ から $10^6$ 回繰り返さねばならず、第2.2節はそれを乗数 $\sigma_R$ として定量化します。

ODMR

3つの要素を組み合わせます。連続照射で $m_s = 0$ にポンプしつつ読み出し、マイクロ波を掃引し、$f_-$ または $f_+$ に当たると占有が $m_s = \mp 1$ に移り、そちらは暗いので蛍光がへこむ。これが光検出磁気共鳴であり、その連続波形式は、誘導ではなく光学的な読み出しで単一スピンに対して行う磁気共鳴実験です。

Code Example 1: 準位構造とODMRスペクトル

上のハミルトニアンを4配向について対角化すればスペクトル全体が得られます。この例ではベクトル再構成も実行し、2次の横方向シフトも検証します。

import numpy as np

# 負に帯電したNV中心の基底状態パラメータです。D、E、gamma_e はダイヤモンド中の
# 欠陥の性質であり、装置の仕様ではありません。
D_GS = 2.870e9        # ゼロ磁場分裂、Hz
E_ST = 2.0e6          # 横方向の歪み・電場による分裂、Hz
GAMMA = 28.025e9      # gamma_e / 2 pi、Hz/T

# 基底 (|+1>, |0>, |-1>) におけるスピン1演算子です。
r2 = 1.0 / np.sqrt(2.0)
Sx = np.array([[0, r2, 0], [r2, 0, r2], [0, r2, 0]], dtype=complex)
Sy = np.array([[0, -1j * r2, 0], [1j * r2, 0, -1j * r2], [0, 1j * r2, 0]])
Sz = np.diag([1.0, 0.0, -1.0]).astype(complex)

# <111>方向を向く、結晶学的に等価な4本のNV軸です。
AXES = np.array([[1, 1, 1], [1, -1, -1], [-1, 1, -1], [-1, -1, 1]],
                dtype=float) / np.sqrt(3.0)


def nv_triad(axis):
    """NV軸を局所zとする正規直交三脚 (e1, e2, axis) を返します。"""
    seed = np.array([1.0, 0.0, 0.0])
    if abs(np.dot(seed, axis)) > 0.9:
        seed = np.array([0.0, 1.0, 0.0])
    e1 = seed - np.dot(seed, axis) * axis
    e1 /= np.linalg.norm(e1)
    e2 = np.cross(axis, e1)
    return e1, e2, axis


def transitions(B_lab, axis, D=D_GS, E=E_ST):
    """磁場 B_lab（テスラ）における1つのNV配向のODMR周波数2本（Hz）です。

    H / h = D Sz^2 + E (Sx^2 - Sy^2) + gamma (B . S)、S はNV座標系で表します。
    """
    e1, e2, e3 = nv_triad(axis)
    Bx, By, Bz = np.dot(B_lab, e1), np.dot(B_lab, e2), np.dot(B_lab, e3)
    H = (D * (Sz @ Sz) + E * (Sx @ Sx - Sy @ Sy)
         + GAMMA * (Bx * Sx + By * Sy + Bz * Sz))
    w = np.linalg.eigvalsh(H)          # 昇順
    return w[1] - w[0], w[2] - w[0]    # 最低準位からの2本の遷移


def odmr(freqs, resonances, contrast, fwhm):
    """CW-ODMRスペクトル: ベースライン1に共鳴ごとのローレンツ型のディップです。"""
    y = np.ones_like(freqs)
    a = fwhm / 2.0
    for f0 in resonances:
        y -= contrast * a**2 / ((freqs - f0) ** 2 + a**2)
    return y


# --- ゼロ磁場、続いて1本のNV軸に沿った磁場 ---------------------------------
print("NV基底三重項: D = %.3f GHz, E = %.1f MHz, "
      "gamma_e/2pi = %.3f GHz/T" % (D_GS / 1e9, E_ST / 1e6, GAMMA / 1e9))
print("\n[111] NV軸に沿った磁場での共鳴周波数:")
print(f"{'B (mT)':>8}{'f- (GHz)':>11}{'f+ (GHz)':>11}{'split (MHz)':>13}"
      f"{'2 gamma B (MHz)':>17}")
print("-" * 60)
for B_mT in [0.0, 0.5, 2.0, 5.0, 20.0]:
    B = B_mT * 1e-3 * AXES[0]
    fm, fp = transitions(B, AXES[0])
    print(f"{B_mT:>8.1f}{fm/1e9:>11.5f}{fp/1e9:>11.5f}"
          f"{(fp-fm)/1e6:>13.4f}{2*GAMMA*B_mT*1e-3/1e6:>17.4f}")

# --- 一般の磁場中における4配向すべて --------------------------------------
B_true = np.array([1.20e-3, -0.70e-3, 2.10e-3])      # テスラ、任意方向
print(f"\n印加磁場 (mT): [{B_true[0]*1e3:.3f}, {B_true[1]*1e3:.3f}, "
      f"{B_true[2]*1e3:.3f}]   |B| = {np.linalg.norm(B_true)*1e3:.4f} mT")
print(f"\n{'family':>7}{'axis':>16}{'B_par (mT)':>12}{'f- (GHz)':>11}"
      f"{'f+ (GHz)':>11}{'split (MHz)':>13}")
print("-" * 70)
res_all, splits = [], []
for k, ax in enumerate(AXES):
    fm, fp = transitions(B_true, ax)
    res_all += [fm, fp]
    splits.append(fp - fm)
    lbl = "[" + " ".join(f"{int(round(c*np.sqrt(3))):+d}" for c in ax) + "]"
    print(f"{k+1:>7}{lbl:>16}{np.dot(B_true, ax)*1e3:>12.4f}"
          f"{fm/1e9:>11.5f}{fp/1e9:>11.5f}{(fp-fm)/1e6:>13.4f}")

# --- 4本の分裂幅からのベクトル再構成 --------------------------------------
# 各分裂幅は最低次で |B . axis| = split / (2 gamma) を与えます。相対符号は
# 全組合せを試して最も整合するものを採用します。この方法ではBの全体符号は
# 観測できないので、最初の射影を正と宣言します。
proj = np.array(splits) / (2.0 * GAMMA)
best = None
for bits in range(8):
    signs = np.array([1.0] + [1.0 if (bits >> i) & 1 else -1.0
                              for i in range(3)])
    B_fit, *_ = np.linalg.lstsq(AXES, signs * proj, rcond=None)
    resid = np.linalg.norm(AXES @ B_fit - signs * proj)
    if best is None or resid < best[0]:
        best = (resid, B_fit)
resid, B_fit = best
print("\n4本の分裂幅からのベクトル再構成（最小二乗）:")
print(f"  復元された B (mT) : [{B_fit[0]*1e3:.4f}, {B_fit[1]*1e3:.4f}, "
      f"{B_fit[2]*1e3:.4f}]")
print(f"  誤差 (uT)         : [{(B_fit[0]-B_true[0])*1e6:.3f}, "
      f"{(B_fit[1]-B_true[1])*1e6:.3f}, {(B_fit[2]-B_true[2])*1e6:.3f}]")
print(f"  残差 (uT)         : {resid*1e6:.4f}")
print(f"  復元された |B|    : {np.linalg.norm(B_fit)*1e3:.4f} mT "
      f"（真値 {np.linalg.norm(B_true)*1e3:.4f} mT）")

# 各対の中心は横磁場によって押し上げられます。B_perp の2次に比例する、単一の
# 分裂幅では得られない第2の観測量です。
print(f"\n{'family':>7}{'centre - D (MHz)':>18}{'B_perp (mT)':>13}"
      f"{'2nd-order prediction':>22}")
print("-" * 60)
for k, ax in enumerate(AXES):
    fm, fp = transitions(B_true, ax)
    Bpar = np.dot(B_true, ax)
    Bperp = np.sqrt(max(np.dot(B_true, B_true) - Bpar**2, 0.0))
    shift = (0.75 * GAMMA**2 * Bperp**2
             * (1.0 / (D_GS + GAMMA * Bpar) + 1.0 / (D_GS - GAMMA * Bpar)))
    print(f"{k+1:>7}{((fp+fm)/2 - D_GS)/1e6:>18.4f}{Bperp*1e3:>13.4f}"
          f"{shift/1e6:>19.4f} MHz")

# --- CW-ODMRスペクトル ------------------------------------------------------
f = np.linspace(2.60e9, 3.14e9, 20001)
trace = odmr(f, res_all, contrast=0.10, fwhm=6.0e6)
print(f"\nODMRスペクトル: ディップ8本、各コントラスト0.10、FWHM 6 MHz")
print(f"  最深点                : ベースラインの {trace.min():.4f}")
print(f"  分離できた極小の本数  : "
      f"{int(np.sum((trace[1:-1] < trace[:-2]) & (trace[1:-1] < trace[2:])))}")

NV基底三重項: D = 2.870 GHz, E = 2.0 MHz, gamma_e/2pi = 28.025 GHz/T

[111] NV軸に沿った磁場での共鳴周波数:
  B (mT)   f- (GHz)   f+ (GHz)  split (MHz)  2 gamma B (MHz)
------------------------------------------------------------
     0.0    2.86800    2.87200       4.0000           0.0000
     0.5    2.85585    2.88415      28.3090          28.0250
     2.0    2.81391    2.92609     112.1713         112.1000
     5.0    2.72986    3.01014     280.2785         280.2500
    20.0    2.30950    3.43050    1121.0071        1121.0000

印加磁場 (mT): [1.200, -0.700, 2.100]   |B| = 2.5179 mT

 family            axis  B_par (mT)   f- (GHz)   f+ (GHz)  split (MHz)
----------------------------------------------------------------------
      1      [+1 +1 +1]      1.5011    2.82959    2.91377      84.1745
      2      [+1 -1 -1]     -0.1155    2.86880    2.87639       7.5853
      3      [-1 +1 -1]     -2.3094    2.80567    2.93516     129.4897
      4      [-1 -1 +1]      0.9238    2.84624    2.89826      52.0134

4本の分裂幅からのベクトル再構成（最小二乗）:
  復元された B (mT) : [1.1902, -0.6933, 2.1111]
  誤差 (uT)         : [-9.771, 6.691, 11.085]
  残差 (uT)         : 7.9138
  復元された |B|    : 2.5207 mT （真値 2.5179 mT）

 family  centre - D (MHz)  B_perp (mT)  2nd-order prediction
------------------------------------------------------------
      1            1.6783       2.0216             1.6779 MHz
      2            2.5959       2.5153             2.5970 MHz
      3            0.4137       1.0033             0.4134 MHz
      4            2.2502       2.3424             2.2524 MHz

ODMRスペクトル: ディップ8本、各コントラスト0.10、FWHM 6 MHz
  最深点                : ベースラインの 0.8820
  分離できた極小の本数  : 8

注目すべき点。 最初の表は計測器の校正です。NV軸に沿った20 mTでは分裂幅が1121.007 MHzで、線形の予測 $2(\gamma_e/2\pi)B = 1121.000$ MHz に対して6桁の線形性を示します。残る7 kHzは2次で寄与する $E = 2$ MHz の歪み項から来ます。ゼロ磁場では分裂幅は厳密に $2E = 4$ MHz、上で述べた二重線です。基準磁石に照らして校正することなく3桁の磁場範囲でこれだけの線形性が得られること、これが「ゼロ磁場分裂は格子の性質である」の実務的な内容です。

ベクトル再構成は4本の分裂幅から2.5 mTの磁場を1成分あたり約10 $\mu$Tで復元し、最小二乗残差は7.9 $\mu$Tです。この誤差は統計的なものではありません。この計算にはノイズがありません。無視した2次の横方向項が各分裂幅をわずかにバイアスしているのであり、したがってベクトルNV磁気計測はミリテスラ領域で $10^3$ 分の1まで線形であり、それを超えると完全な対角化が必要です。最後の表はその2次の式の検査です。測定された中心シフトはファミリーごとに $\frac{3}{4}(\gamma_e/2\pi)^2 B_\perp^2 [\,(D+\gamma_e B_\parallel/2\pi)^{-1} + (D-\gamma_e B_\parallel/2\pi)^{-1}]$ と4桁まで一致します。係数が $1/2$ でなく $3/4$ なのは、横磁場が $m_s = \pm 1$ の準位を押し上げかつ $m_s = 0$ の準位を押し下げ、両方の変位が同じ符号で遷移周波数に入るからです。ここは正しく持っておく価値があります。中心シフトは第2の独立な観測量であり、単一の配向が与える $B_\perp$ への唯一の手がかりです。

2.2 DC磁気計測

連続波ODMRと、その感度を決めるもの

CW-ODMR実験は光子計数率にローレンツ型のディップを生みます。その感度は2つの要素と1つの動作点の選び方から従います。

検出率を

$$ R(\delta) = R_0 \left[ 1 - C\, \frac{(\Gamma/2)^2}{\delta^2 + (\Gamma/2)^2} \right] $$

と書きます。$\delta$ は共鳴からのマイクロ波離調、$\Gamma$ は半値全幅、$C$ はコントラスト、$R_0$ は非共鳴での計数率です。線が最も急な場所にマイクロ波周波数を置きます。$a = \Gamma/2$、$L(\delta) = a^2/(\delta^2+a^2)$ と書けば

$$ \frac{\mathrm{d}L}{\mathrm{d}\delta} = -\frac{2a^2\delta}{(\delta^2+a^2)^2}, \qquad \frac{\mathrm{d}}{\mathrm{d}\delta}\left[\frac{\delta}{(\delta^2+a^2)^2}\right] = 0 \;\Longrightarrow\; \delta = \frac{a}{\sqrt{3}} $$

であり、その点で $|\mathrm{d}L/\mathrm{d}\delta| = 9/(8\sqrt{3}\,a) = 3\sqrt{3}/(4\Gamma)$、$L = 3/4$ です。時間 $t$ の検出光子数は $N = Rt$ でポアソン不確かさ $\sqrt{R_0 t}$、磁場変化 $\delta B$ は線を $(\gamma_e/2\pi)\delta B$ だけずらします。ノイズを傾きで割って

$$ \sigma_B = \frac{\sqrt{R_0 t}}{t\, R_0 C\, \dfrac{3\sqrt{3}}{4\Gamma}\, \dfrac{\gamma_e}{2\pi}} \qquad\Longrightarrow\qquad \boxed{\ \eta_{\mathrm{CW}} = \sigma_B \sqrt{t} = \frac{4}{3\sqrt{3}}\; \frac{2\pi}{\gamma_e}\; \frac{\Gamma}{C\sqrt{R_0}}\ } $$

3つの因子と1つの幾何学的定数です。前係数 $4/(3\sqrt{3}) = 0.7698$ はローレンツ型という線形状の代価であって、それ以外の何でもありません。物理は $\Gamma/(C\sqrt{R_0})$ にあります。これを独立な3つのつまみとして読めば、線幅は1乗で入るので最も効き、次にコントラスト、計数率は平方根でしか入らないので最も効かない、ということになります。指数のこの読み方は正しく、そして3つを独立と見なせる最後の地点でもあります。続く2つの小節でその独立性は二重に崩れます。マイクロ波パワーが $\Gamma$ と $C$ を結びつけ、$R_0$ を決めるレーザー光が $\Gamma$ をも広げるからです。その後では、この重要度の順は式についての言明であって実験の処方ではありません。

CW-ODMRの罠

この式は明らかな最適化 — 線を狭くする — を誘い、その明らかな最適化は失敗します。量子センシング全般で繰り返される理由なので述べておく価値があります。$\Gamma$ と $C$ は独立ではありません。コントラストにはスピン遷移の駆動が必要で、駆動すれば線が広がるからです。$s = \Omega^2 T_1 T_2$ をマイクロ波駆動の飽和パラメータとして、標準的な飽和共鳴の形は

$$ \Gamma(s) = \Gamma_0\sqrt{1+s}, \qquad C(s) = C_{\max}\frac{s}{1+s} $$

であり、したがって $\eta_{\mathrm{CW}} \propto (1+s)^{3/2}/s$ となって内点に最小値をもちます。微分すると $\frac{3}{2}s = 1+s$ なので

$$ s_{\mathrm{opt}} = 2, \qquad \Gamma = \sqrt{3}\,\Gamma_0, \qquad C = \tfrac{2}{3}C_{\max}, \qquad \eta_{\mathrm{CW}}^{\min} = \frac{2\,\Gamma_0}{(\gamma_e/2\pi)\, C_{\max}\sqrt{R_0}} $$

です。$\frac{4}{3\sqrt{3}}\cdot\frac{3^{3/2}}{2} = 2$ が厳密に成り立つからです。最適点は広がっていない線幅の $\sqrt{3}$ 倍、最大コントラストの3分の2にあり、結果はマイクロ波を切ったときの線幅 $\Gamma_0$ だけに依存します。

そして $\Gamma_0$ は $1/(\pi T_2^\ast)$ ではありません。光子を生み出すレーザーは同時にスピンを再分極させ、その速度は $R_0$ を決める光学励起率に比例します。この再ポンプは、光が当たっている間つねに存在する追加の均一広がりです。したがってCW-ODMRの線幅は、$T_2^\ast$ 律速の幅が数百キロヘルツであるような場合でも、マイクロ波パワーがゼロに近づいてもメガヘルツ程度あります。

このことは感度の式が隠している帰結をもち、それがこの小節の2つめの最適点です。マイクロ波を切ったときの幅を $\Gamma_0 = \Gamma_{\mathrm{int}} + \kappa R_0$ と書きます。$\Gamma_{\mathrm{int}}$ は暗状態での幅、$\kappa$ は単位計数率あたりの再ポンプ広がりです。すると

$$ \eta_{\mathrm{CW}}^{\min} \;\propto\; \frac{\Gamma_{\mathrm{int}} + \kappa R_0}{\sqrt{R_0}} $$

となり、これが $1/\sqrt{R_0}$ で下がるのは $\kappa R_0 \ll \Gamma_{\mathrm{int}}$ のあいだだけです。$\kappa R_0 = \Gamma_{\mathrm{int}}$ — レーザー起因の広がりが残留幅に等しくなる点 — で内点最小値をとり、それを越えると $\sqrt{R_0}$ で増加します。したがって最適なマイクロ波パワーと同様に最適なレーザーパワーが存在し、それを越えると光を増やすほど磁力計は悪くなります。上の重要度の順における「明るい発光体」が $\Gamma_0$ 一定のもとでしか成り立たない理由がこれです。単一中心の光学読み出しはそれを許しません。1個の欠陥では、明るさと線幅は同じ勘定から支払われるのです。

これがCW-ODMRが非効率だという正確な意味でもあり、パルス系列がより良い正確な理由です。位相が蓄積している間、レーザーもマイクロ波も切ることで、2つの勘定を完全に分離するのです。

Code Example 2: CW-ODMR感度の検証

import numpy as np

GAMMA = 28.025e9        # gamma_e / 2 pi、Hz/T


def lorentz(detuning, fwhm):
    """線中心で1に規格化したローレンツ関数です。"""
    a = fwhm / 2.0
    return a**2 / (detuning**2 + a**2)


def eta_cw(fwhm, contrast, rate, gamma=GAMMA):
    """ショットノイズ限界のCW-ODMR感度を T/sqrt(Hz) で返します。"""
    return 4.0 / (3.0 * np.sqrt(3.0)) * fwhm / (gamma * contrast * np.sqrt(rate))


# --- 手順1: ローレンツ型ディップの最大傾き ----------------------------------
fwhm = 6.0e6
d = np.linspace(-3 * fwhm, 3 * fwhm, 4000001)
slope = np.abs(np.gradient(lorentz(d, fwhm), d))
i = int(np.argmax(slope))
print("ローレンツ型ディップの最大傾き")
print(f"  FWHM                          : {fwhm/1e6:.3f} MHz")
print(f"  数値の |dL/d(delta)|_max      : {slope[i]*1e6:.6f} per MHz")
print(f"  解析式 (3 sqrt3 / 4) / FWHM   : "
      f"{3*np.sqrt(3)/4/fwhm*1e6:.6f} per MHz")
print(f"  数値の最大点における|離調|    : {abs(d[i])/1e6:.5f} MHz")
print(f"  解析式 FWHM/(2 sqrt3)         : {fwhm/(2*np.sqrt(3))/1e6:.5f} MHz")
print(f"  その離調における L            : {lorentz(d[i], fwhm):.6f} "
      f"（解析値 0.750000）")

# --- 手順2: 感度の式とモンテカルロ検証 --------------------------------------
C, R0 = 0.20, 1.0e5          # ODMRコントラスト、検出光子数（毎秒）
delta_op = fwhm / (2 * np.sqrt(3.0))
R_op = R0 * (1.0 - C * lorentz(delta_op, fwhm))
dRdB = R0 * C * (3 * np.sqrt(3) / 4) / fwhm * GAMMA     # photons/s per tesla

eta = eta_cw(fwhm, C, R0)
print(f"\nCW-ODMR感度、C = {C:.2f}、R = {R0:.0e} counts/s")
print(f"  傾き dR/dB                    : {dRdB:.6e} counts/s/T")
print(f"  式による eta                  : {eta*1e6:.4f} uT/sqrt(Hz)")

rng = np.random.default_rng(20260813)
print(f"\n{'t (ms)':>8}{'trials':>8}{'std(B_est) (uT)':>17}"
      f"{'eta/sqrt(t) (uT)':>18}{'ratio':>8}")
print("-" * 59)
for t_ms in [1.0, 10.0, 100.0]:
    t = t_ms * 1e-3
    n_trials = 20000
    counts = rng.poisson(R_op * t, n_trials)
    B_est = (counts - R_op * t) / (dRdB * t)
    pred = eta / np.sqrt(t)
    print(f"{t_ms:>8.1f}{n_trials:>8d}{B_est.std(ddof=1)*1e6:>17.4f}"
          f"{pred*1e6:>18.4f}{B_est.std(ddof=1)/pred:>8.4f}")
print(f"  モンテカルロは式より {np.sqrt(R_op/R0):.4f} 倍静かです。動作点が")
print(f"  ベースライン計数率の {R_op/R0:.4f} 倍にあるのに対し、標準的な式は")
print("  ショットノイズに sqrt(R_baseline) を使うためです。")

# --- 手順3: マイクロ波パワーに関する自己整合的な最適点 ----------------------
# 線幅とコントラストは独立ではありません。マイクロ波パワーはコントラストを買い、
# 線幅の広がりで支払います。飽和パラメータ s = Omega^2 T1 T2 として
#   FWHM = Gamma_0 sqrt(1 + s)、C = C_max s / (1 + s)
# なので eta は (1 + s)^(3/2) / s に比例し、内点に最小値をもちます。
Gamma_0, C_max = 2.0e6, 0.30
print(f"\n自己整合的なCW最適点、Gamma_0 = {Gamma_0/1e6:.1f} MHz、"
      f"C_max = {C_max:.2f}、R = 1e5/s")
print(f"{'s':>8}{'FWHM (MHz)':>12}{'contrast':>10}"
      f"{'eta (uT/sqrt(Hz))':>19}")
print("-" * 49)
for sat in [0.25, 0.5, 1.0, 2.0, 4.0, 8.0, 20.0]:
    fw = Gamma_0 * np.sqrt(1.0 + sat)
    c = C_max * sat / (1.0 + sat)
    print(f"{sat:>8.2f}{fw/1e6:>12.4f}{c:>10.4f}"
          f"{eta_cw(fw, c, 1e5)*1e6:>19.4f}")
sat = np.linspace(0.05, 50.0, 200000)
etas = eta_cw(Gamma_0 * np.sqrt(1 + sat), C_max * sat / (1 + sat), 1e5)
j = int(np.argmin(etas))
print(f"  数値の最適点 : s = {sat[j]:.4f}, "
      f"eta = {etas[j]*1e6:.4f} uT/sqrt(Hz)")
print(f"  解析の最適点 : s = 2, eta_min = 2 Gamma_0 / "
      f"(gamma C_max sqrt(R)) = "
      f"{2*Gamma_0/(GAMMA*C_max*np.sqrt(1e5))*1e6:.4f} uT/sqrt(Hz)")
print("  Gamma_0 はマイクロ波を切ったときの線幅です。1/(pi T2*) ではありません。")
print("  光子を生み出すレーザーは掃引中もスピンを再分極させ、その再ポンプ速度が")
print("  残留線幅を支配します。")

# --- 手順4: アンサンブル ----------------------------------------------------
print("\n集光効率一定でのアンサンブルのスケーリング（R は N に比例）:")
print(f"{'N spins':>10}{'R (1/s)':>12}{'eta (uT/sqrt(Hz))':>19}"
      f"{'eta * sqrt(N)':>16}")
print("-" * 57)
for N in [1, 100, 10000, 1000000]:
    r = 1e5 * N
    e = eta_cw(6.0e6, 0.20, r)
    print(f"{N:>10d}{r:>12.0e}{e*1e6:>19.6f}{e*np.sqrt(N)*1e6:>16.4f}")

ローレンツ型ディップの最大傾き
  FWHM                          : 6.000 MHz
  数値の |dL/d(delta)|_max      : 0.216506 per MHz
  解析式 (3 sqrt3 / 4) / FWHM   : 0.216506 per MHz
  数値の最大点における|離調|    : 1.73205 MHz
  解析式 FWHM/(2 sqrt3)         : 1.73205 MHz
  その離調における L            : 0.750000 （解析値 0.750000）

CW-ODMR感度、C = 0.20、R = 1e+05 counts/s
  傾き dR/dB                    : 1.213518e+08 counts/s/T
  式による eta                  : 2.6059 uT/sqrt(Hz)

  t (ms)  trials  std(B_est) (uT)  eta/sqrt(t) (uT)   ratio
-----------------------------------------------------------
     1.0   20000          76.2912           82.4050  0.9258
    10.0   20000          24.0765           26.0588  0.9239
   100.0   20000           7.6149            8.2405  0.9241
  モンテカルロは式より 0.9220 倍静かです。動作点が
  ベースライン計数率の 0.8500 倍にあるのに対し、標準的な式は
  ショットノイズに sqrt(R_baseline) を使うためです。

自己整合的なCW最適点、Gamma_0 = 2.0 MHz、C_max = 0.30、R = 1e5/s
       s  FWHM (MHz)  contrast  eta (uT/sqrt(Hz))
-------------------------------------------------
    0.25      2.2361    0.0600             3.2372
    0.50      2.4495    0.1000             2.1277
    1.00      2.8284    0.1500             1.6379
    2.00      3.4641    0.2000             1.5045
    4.00      4.4721    0.2400             1.6186
    8.00      6.0000    0.2667             1.9544
   20.00      9.1652    0.2857             2.7864
  数値の最適点 : s = 2.0001, eta = 1.5045 uT/sqrt(Hz)
  解析の最適点 : s = 2, eta_min = 2 Gamma_0 / (gamma C_max sqrt(R)) = 1.5045 uT/sqrt(Hz)
  Gamma_0 はマイクロ波を切ったときの線幅です。1/(pi T2*) ではありません。
  光子を生み出すレーザーは掃引中もスピンを再分極させ、その再ポンプ速度が
  残留線幅を支配します。

集光効率一定でのアンサンブルのスケーリング（R は N に比例）:
   N spins     R (1/s)  eta (uT/sqrt(Hz))   eta * sqrt(N)
---------------------------------------------------------
         1       1e+05           2.605876          2.6059
       100       1e+07           0.260588          2.6059
     10000       1e+09           0.026059          2.6059
   1000000       1e+11           0.002606          2.6059

注目すべき点。 最初のブロックが幾何を確認します。ローレンツ型ディップの数値的な最大傾きは0.216506毎MHzで、解析式 $3\sqrt{3}/(4\Gamma)$ と6桁まで一致し、最急点は共鳴から $\Gamma/(2\sqrt{3}) = 1.73205$ MHz にあり、そこで $L = 0.750000$ です。これらは線形状の性質であって欠陥の性質ではありません。

モンテカルロが誠実な部分です。ポアソン光子計数をシミュレートし、解析的な傾きを通じて磁場推定に変換すると、3桁の平均時間にわたって $\eta/\sqrt{t}$ を再現し、しかも一貫して7.6%下回ります。これは誤りではありません。標準的な式はショットノイズに $\sqrt{R_0}$ を使いますが、動作点はベースライン計数率の $1 - 3C/4 = 0.85$ 倍にあり、真のノイズは仮定値の $\sqrt{0.85} = 0.922$ 倍です。通常のコントラストでは式は数パーセント保守的であり、その差は $\sqrt{1-3C/4}$ の大きさです。10%の不一致を物理に帰する前に知っておくべき値です。

自己整合的な最適化は $s_{\mathrm{opt}} = 2$ を4桁まで確認し、閉形式 $2\Gamma_0/(\gamma C_{\max}\sqrt{R_0})$ を厳密に再現します。マイクロ波を切ったときの幅2 MHz、最大コントラスト0.30、毎秒 $10^5$ カウント — 装置仕様ではないので有効数字1桁で示した代表的な単一中心の値 — で答えは1.50 $\mu$T/$\sqrt{\mathrm{Hz}}$ です。最後のブロックがアンサンブルの存在理由であり、その理由は慎重に帰属させる必要があります。この表は中心あたりの照射一定のもとで $R_0$ を $N$ に比例して上げています。1個の発光体に強いビームを当てるのではなく、同じビーム強度のもとで発光体を増やすのです。だからこそ $\Gamma_0$ は変わらず、それが「$R_0$ を上げる」の唯一の版として $\eta \propto 1/\sqrt{N}$ を厳密に成立させます。ここでは実際に成立し、$\eta\sqrt{N}$ は6桁まで一定で、100万個の中心が3桁を買います。単一中心でレーザーを強くするのがもう一方の版で、上で導いた最適点を越えれば何も買えません。そこでは $R_0$ と $\Gamma_0$ が一緒に上がるからです。この表が示せないのは — 示せないのですが — その代価です。窒素濃度を上げると中心同士の双極子結合で $T_2^\ast$ が短くなるので、分子の線幅 $\Gamma_0$ が分母の $N$ と同じだけ増えます。最適濃度は材料の問題であり、アンサンブルNV磁気計測における中心的な材料の問題です。

パルスという選択肢: Ramsey

第1章のRamsey系列は、レーザー由来の広がりとパワー広がりの両方を構成上取り除きます。光学的に初期化し、レーザーを切り、$\pi/2$ を当て、マイクロ波なしの暗闇で $\tau$ 待ち、$\pi/2$ を当て、それからレーザーを入れて読み出す。暗闇で蓄積される位相は

$$ \varphi = 2\pi \frac{\gamma_e}{2\pi} B_\parallel \tau $$

であり、測定される占有数は $P = \frac{1}{2}\left[1 - C\cos\varphi\right]$、$C = C_0 \exp[-(\tau/T_2^\ast)^p]$ です。系列を最急点 $\varphi = \pi/2$ にバイアスすると $P = 1/2$ で

$$ \left|\frac{\mathrm{d}P}{\mathrm{d}B}\right| = \pi \frac{\gamma_e}{2\pi} \tau\, C $$

となります。射影的に測定された単一スピンは1ショットあたり $\sigma_P = 1/2$ の二項分布の結果を与えるので、読み出しと再分極のデッドタイムを $t_d$ として全時間 $T = N(\tau + t_d)$ の $N$ ショットは

$$ \sigma_B = \frac{1}{2\sqrt{N}}\left|\frac{\mathrm{d}P}{\mathrm{d}B}\right|^{-1} \qquad\Longrightarrow\qquad \boxed{\ \eta_{\mathrm{Ramsey}}(\tau) = \frac{\sqrt{\tau + t_d}}{2\pi (\gamma_e/2\pi)\, C_0\, \tau}\, e^{(\tau/T_2^\ast)^p}\ } $$

を与えます。この式の構造は本コースのあらゆる感度の構造です。分子は $\sqrt{\tau}$ で増えます。系列が長いと繰り返し回数が減るからです。分母は $\tau$ で増えます。長い系列はより多くの位相を蓄積するからです。そしてコヒーレンスが有限なので最終的には指数が勝ちます。$t_d = 0$ として対数を微分すると

$$ \frac{\mathrm{d}}{\mathrm{d}\tau}\left[\left(\frac{\tau}{T_2^\ast}\right)^p - \frac{1}{2}\ln\tau\right] = 0 \qquad\Longrightarrow\qquad \tau^p = \frac{(T_2^\ast)^p}{2p} $$

となり、$p = 1$ でも $p = 2$ でも $\tau_{\mathrm{opt}} = T_2^\ast/2$ です。この一致は深い理由によるものではありませんが便利であり、「コヒーレンス時間の半分で走らせる」という規則を減衰形の知識なしに使えるようにします。最小値は

$$ \eta^{\min} = \frac{e^{1/(2p)}}{2\pi(\gamma_e/2\pi) C_0}\sqrt{\frac{2}{T_2^\ast}} \;\propto\; \frac{1}{\sqrt{T_2^\ast}} $$

であり、$1/\sqrt{T_2^\ast}$ のスケーリングが最も重要な帰結です。コヒーレンス時間を倍にしても買えるのは $\sqrt{2}$ であって2ではありません。感度はコヒーレンスにおいて高価なのです。

読み出しのペナルティ

以上はすべて射影的な単一ショット測定を仮定していました。NVの光学読み出しはそのようなものではありません。$m_s = 0$ と $m_s = \pm 1$ について1ショットあたりの平均検出光子数を $\alpha_0$、$\alpha_1$ とします。標準的な結果は、感度が

$$ \sigma_R = \sqrt{1 + \frac{2(\alpha_0 + \alpha_1)}{(\alpha_0 - \alpha_1)^2}} \;\ge\; 1 $$

倍されるというものです。光子数が多くその差が1ショットで分解されるとき1に帰着します。$\alpha_0 \approx 0.03$ の単一中心では30を超え、単一NV磁力計の感度を支配する項になります。これを改善するもの — ナノ構造化した表面による集光の改善、核スピン補助による反復読み出し、スピン電荷変換 — はいずれもスピンの物理ではなく $\sigma_R$ を攻めています。

Code Example 3: 射影ノイズ込みのRamsey感度

"""第2章 Code Example 3: Ramsey系列によるDC磁気計測。
Code Example 2 の続き（同一セッション）。"""

from scipy.optimize import minimize_scalar


def eta_ramsey(tau, T2s, C0=1.0, t_dead=0.0, p=1.0, gamma=GAMMA):
    """射影ノイズ限界のRamsey感度を T/sqrt(Hz) で返します。

    eta = sqrt(tau + t_dead) * exp((tau/T2s)^p) / (2 pi gamma C0 tau)
    """
    return (np.sqrt(tau + t_dead) * np.exp((tau / T2s) ** p)
            / (2 * np.pi * gamma * C0 * tau))


def readout_factor(a0, a1):
    """雑音の多い光学読み出しのオーバーヘッドです。sigma_R >= 1 であり、
    単一ショットの測定が射影的なとき1に等しくなります。"""
    return np.sqrt(1.0 + 2.0 * (a0 + a1) / (a0 - a1) ** 2)


T2s = 1.0e-6
print("Ramsey DC磁気計測、T2* = %.2f us" % (T2s * 1e6))
print(f"\n{'decay p':>9}{'t_dead (us)':>13}{'tau_opt (us)':>14}"
      f"{'tau_opt / T2*':>15}{'eta_min (nT/sqrt(Hz))':>23}")
print("-" * 74)
for p in (1.0, 2.0):
    for td in (0.0, 0.3e-6, 3.0e-6):
        r = minimize_scalar(eta_ramsey, bracket=(1e-9, T2s, 20 * T2s),
                            args=(T2s, 1.0, td, p))
        print(f"{p:>9.1f}{td*1e6:>13.2f}{r.x*1e6:>14.5f}{r.x/T2s:>15.5f}"
              f"{r.fun*1e9:>23.4f}")
print("  デッドタイムがなければ最適値はどちらの減衰形でも厳密に T2*/2 です。")
print("  d/dtau [ (tau/T2*)^p - (1/2) ln tau ] = 0 が tau^p = T2*^p/(2p) を")
print("  与えるからです。")

# --- p = 1 と p = 2、デッドタイムなしの閉形式 -------------------------------
print("\nt_dead = 0、tau = T2*/2 における閉形式:")
for p, coef in ((1.0, np.exp(0.5)), (2.0, np.exp(0.25))):
    ana = coef / (2 * np.pi * GAMMA) * np.sqrt(2.0 / T2s)
    num = eta_ramsey(T2s / 2, T2s, 1.0, 0.0, p)
    print(f"  p = {p:.0f}: exp(1/(2p)) sqrt(2/T2*) / (2 pi gamma) = "
          f"{ana*1e9:.6f} nT/sqrt(Hz)   数値 {num*1e9:.6f}")

# --- モンテカルロ: 射影ノイズ、1スピン、最急点にバイアス --------------------
rng = np.random.default_rng(4712)
tau = T2s / 2
C = np.exp(-(tau / T2s))            # 最適点における p = 1 のコントラスト
slope = np.pi * GAMMA * tau * C     # pi/2 バイアス点での dP/dB
print(f"\nモンテカルロ、tau = T2*/2 = {tau*1e6:.3f} us、コントラスト "
      f"{C:.6f}、dP/dB = {slope:.6e} /T")
print(f"{'shots N':>10}{'T = N tau (ms)':>16}{'std(B) (nT)':>14}"
      f"{'eta/sqrt(T) (nT)':>18}{'ratio':>8}")
print("-" * 66)
for N in (1000, 10000, 100000):
    k = rng.binomial(N, 0.5, size=4000)
    B_est = (k / N - 0.5) / slope
    T_tot = N * tau
    pred = eta_ramsey(tau, T2s, 1.0, 0.0, 1.0) / np.sqrt(T_tot)
    print(f"{N:>10d}{T_tot*1e3:>16.4f}{B_est.std(ddof=1)*1e9:>14.5f}"
          f"{pred*1e9:>18.5f}{B_est.std(ddof=1)/pred:>8.4f}")

# --- 光学読み出しのペナルティ ----------------------------------------------
print("\n光学読み出しは射影的ではありません。sigma_R があらゆる eta に掛かります。")
print(f"{'alpha_0':>9}{'alpha_1':>9}{'contrast':>10}{'sigma_R':>10}"
      f"{'eta (nT/sqrt(Hz))':>19}")
print("-" * 57)
eta_pn = eta_ramsey(T2s / 2, T2s, 1.0, 0.0, 1.0)
for a0, a1 in [(0.030, 0.020), (0.30, 0.20), (3.0, 2.0), (1000.0, 700.0)]:
    sR = readout_factor(a0, a1)
    print(f"{a0:>9.3f}{a1:>9.3f}{(a0-a1)/a0:>10.4f}{sR:>10.3f}"
          f"{eta_pn*sR*1e9:>19.3f}")
print("  最後の行は反復読み出しまたはスピン電荷変換です。1ショットあたりの光子数")
print("  が多く、sigma_R は1に近づき、射影限界が射程に入ります。")

# --- 同一スピンでのパルスとCWの比較 ----------------------------------------
# CW-ODMRは連続動作なので、その eta にはすでにデューティ比100%が含まれています。
# t_dead = 0 の Ramsey の eta には含まれないので、その比較は上限にすぎません。
# 対等な比較は上の表と同じ 3 us の読み出し時間を使います。
sR = readout_factor(0.030, 0.020)
eta_R = eta_pn * sR
r3 = minimize_scalar(eta_ramsey, bracket=(1e-9, T2s, 20 * T2s),
                     args=(T2s, 1.0, 3.0e-6, 1.0))
eta_R3 = r3.fun * sR
eta_C = 2.0 * 2.0e6 / (GAMMA * 0.30 * np.sqrt(1e5))     # Code Example 2 のCW最適値
print(f"\n同一の単一NV、T2* = {T2s*1e6:.1f} us、R = 1e5 counts/s、"
      f"sigma_R = {sR:.2f}:")
print(f"  CW-ODMR（自身の最適点）    : {eta_C*1e9:9.1f} nT/sqrt(Hz)")
print(f"  Ramsey、t_dead = 0         : {eta_R*1e9:9.1f} nT/sqrt(Hz)"
      f"   比 {eta_C/eta_R:5.2f}")
print(f"  Ramsey、t_dead = 3 us      : {eta_R3*1e9:9.1f} nT/sqrt(Hz)"
      f"   比 {eta_C/eta_R3:5.2f}")
print("  対等なのは2つめのRamseyの行だけなので、1つめの比は利得の上限です。")
print("  利得の源は線幅を狭めること自体ではありません。自由発展の間はレーザーも")
print("  マイクロ波も切れているので、再ポンプもパワー広がりも位相蓄積に入らない")
print("  という点にあります。")

Ramsey DC磁気計測、T2* = 1.00 us

  decay p  t_dead (us)  tau_opt (us)  tau_opt / T2*  eta_min (nT/sqrt(Hz))
--------------------------------------------------------------------------
      1.0         0.00       0.50000        0.50000                13.2415
      1.0         0.30       0.65678        0.65678                16.3116
      1.0         3.00       0.88600        0.88600                30.6462
      2.0         0.00       0.50000        0.50000                10.3125
      2.0         0.30       0.57906        0.57906                12.8583
      2.0         3.00       0.67390        0.67390                25.4375
  デッドタイムがなければ最適値はどちらの減衰形でも厳密に T2*/2 です。
  d/dtau [ (tau/T2*)^p - (1/2) ln tau ] = 0 が tau^p = T2*^p/(2p) を
  与えるからです。

t_dead = 0、tau = T2*/2 における閉形式:
  p = 1: exp(1/(2p)) sqrt(2/T2*) / (2 pi gamma) = 13.241487 nT/sqrt(Hz)   数値 13.241487
  p = 2: exp(1/(2p)) sqrt(2/T2*) / (2 pi gamma) = 10.312480 nT/sqrt(Hz)   数値 10.312480

モンテカルロ、tau = T2*/2 = 0.500 us、コントラスト 0.606531、dP/dB = 2.670043e+04 /T
   shots N  T = N tau (ms)   std(B) (nT)  eta/sqrt(T) (nT)   ratio
------------------------------------------------------------------
      1000          0.5000     602.54211         592.17729  1.0175
     10000          5.0000     187.48253         187.26290  1.0012
    100000         50.0000      58.57488          59.21773  0.9891

光学読み出しは射影的ではありません。sigma_R があらゆる eta に掛かります。
  alpha_0  alpha_1  contrast   sigma_R  eta (nT/sqrt(Hz))
---------------------------------------------------------
    0.030    0.020    0.3333    31.639            418.942
    0.300    0.200    0.3333    10.050            133.075
    3.000    2.000    0.3333     3.317             43.917
 1000.000  700.000    0.3000     1.019             13.489
  最後の行は反復読み出しまたはスピン電荷変換です。1ショットあたりの光子数
  が多く、sigma_R は1に近づき、射影限界が射程に入ります。

同一の単一NV、T2* = 1.0 us、R = 1e5 counts/s、sigma_R = 31.64:
  CW-ODMR（自身の最適点）    :    1504.5 nT/sqrt(Hz)
  Ramsey、t_dead = 0         :     418.9 nT/sqrt(Hz)   比  3.59
  Ramsey、t_dead = 3 us      :     969.6 nT/sqrt(Hz)   比  1.55
  対等なのは2つめのRamseyの行だけなので、1つめの比は利得の上限です。
  利得の源は線幅を狭めること自体ではありません。自由発展の間はレーザーも
  マイクロ波も切れているので、再ポンプもパワー広がりも位相蓄積に入らない
  という点にあります。

注目すべき点。 デッドタイムがなければ最適点は代数の要求どおりどちらの減衰形でも厳密に $\tau = T_2^\ast/2$ にあり、閉形式 $e^{1/(2p)}\sqrt{2/T_2^\ast}/(2\pi\gamma)$ が数値の最小値を6桁まで再現します。ガウス型減衰（$p = 2$）は同じ $T_2^\ast$ で指数型より22%良く、崩れる前に長く平らでいるからです。

デッドタイムの行は誠実な補正であり、NV中心では小さくありません。読み出しと再分極は1から3 $\mu$s 程度を要し、自然同位体比のダイヤモンドでの $T_2^\ast$ は1 $\mu$s 程度です。$t_d = 3\ \mu$s では最適点が $0.886\,T_2^\ast$ まで伸び、感度は純粋にデューティ比のせいで2.3倍悪化します。$^{12}$C 同位体精製が二重に効く理由がこれです。長い $T_2^\ast$ は $1/\sqrt{T_2^\ast}$ で $\eta$ を改善し、かつデッドタイムが系列と同程度でなくなることでデューティ比も改善します。

モンテカルロは射影ノイズの計算を2桁のショット数にわたって約1%で確認します。バイアスされた縞の二項標本を解析的な傾きで換算すると、自由パラメータなしで $\eta/\sqrt{T}$ を再現します。読み出しの表がそこにペナルティを適用します。1ショットあたり0.03光子を集める単一中心は31.6倍を支払い、射影限界13.2 nT/$\sqrt{\mathrm{Hz}}$ に対して419 nT/$\sqrt{\mathrm{Hz}}$ に着地します。光子予算を2桁改善して1ショットあたり数光子にすると1桁の感度が戻り、1000光子に達すれば実質すべて戻ります。スピン電荷読み出しの論拠を数値で述べればこうなります。

最後に比較ですが、読む前に1つ補正が必要です。同じ中心、同じ光子率で、最適化したCW測定は1.50 $\mu$T/$\sqrt{\mathrm{Hz}}$、Ramsey測定は419 nT/$\sqrt{\mathrm{Hz}}$、比は3.6倍です。しかし419 nTという値は $t_d = 0$ の行のものであり、CW-ODMRは構成上連続動作なのでその数値にはすでにデューティ比100%が含まれています。上のデッドタイムの表が仮定したのと同じ3 $\mu$s の読み出し・再分極で比べると、Ramseyは970 nT/$\sqrt{\mathrm{Hz}}$ となり、比は1.55倍です。両方の数値が出力にあります。3.6倍はパルス方式の優位のデッドタイムゼロにおける上限、1.55倍が実際の単一中心実験が得る値です。したがって利得は控えめに聞こえるだけでなく実際に控えめであり、その起源は正確に述べる価値があります。Ramsey縞が共鳴線より本質的に優れているのではなく、自由発展の間はレーザーもマイクロ波も切れているので、光学的再ポンプもパワー広がりも位相蓄積に入らないということです。

2.3 AC磁気計測、動的デカップリング、ノイズ分光

なぜRamseyは振動する磁場を見られないのか

Ramsey系列は磁場を $\tau$ にわたって重み1で積分します。振動する磁場は整数周期にわたって平均するとゼロなので、$f \gg 1/\tau$ の信号のRamsey測定は何も返しません。これは、自由誘導が直流で感度最大であるという第1章の事実を逆から読んだものです。$|\tilde{s}_{\mathrm{FID}}(f,\tau)|^2 = \sin^2(\pi f\tau)/(\pi f)^2 \to \tau^2$（$f \to 0$）であり、その後は $1/f^2$ で落ちます。

対策は第1章の対策を逆の目的で使うことです。系列の中央に置いた $\pi$ パルスは蓄積位相の符号を反転し、系列をロックイン増幅器に変えます。直流で盲目になり、通過帯域で感度をもつのです。位相緩和にとってはノイズを排除するフィルターであり、単色信号にとっては復調器です。

復調器としてのエコー

第1章の変調関数 $s(t) = \pm 1$ を取り、$\pi$ パルスの前で $+1$、後で $-1$ とします。磁場 $B(t) = B_{\mathrm{AC}}\sin(2\pi f t)$ に対して蓄積位相は

$$ \varphi = 2\pi \frac{\gamma_e}{2\pi}\int_0^T s(t)\,B(t)\,\mathrm{d}t $$

です。$T/2$ に1つの $\pi$ パルスがあり $f = 1/T$ のとき — 系列あたり信号がちょうど1周期で、ゼロ交差がパルス位置にあるとき — 前半と後半はそれぞれ $+T/\pi$ を寄与するので

$$ \varphi = 2\pi \frac{\gamma_e}{2\pi} B_{\mathrm{AC}} \cdot \frac{2T}{\pi} = 4\frac{\gamma_e}{2\pi}\, B_{\mathrm{AC}}\, T $$

です。同じ振幅の直流磁場を同じ時間測った場合 $\varphi_{\mathrm{DC}} = 2\pi(\gamma_e/2\pi)BT$ と比べると、AC測定は厳密に $2/\pi = 0.6366$ 倍だけ非効率であり、これは矩形波の基本フーリエ係数です。この因子が復調の代価であり、一度だけ支払われます。

買えるものはずっと大きいのです。エコーは $T_2^\ast$ ではなく $T_2$ を測り、固体母材ではこの2つが1桁から2桁違います。$T_2^\ast$ を制限するノイズが低周波に集中しており、そこはエコーのフィルターが零点をもつ場所だからです。信号で $\pi/2$ を失う代わりにコヒーレンス時間で10倍から100倍 — したがって感度で $\sqrt{10}$ から $10$ 倍 — を得るのは価値ある取引であり、AC磁気計測がNVセンシングで最も競争力をもつ理由です。

CPMG、XY8、そしてパルスを増やすと何が買えるのか

$\pi$ パルスを増やすと通過帯域が上へ移り、コヒーレンスが伸びます。CPMG-$N$ は $N$ 個のパルスを $t_j = (j-\tfrac{1}{2})T/N$ に置くので、変調関数は周期 $2T/N$ の矩形波となり、系列は

$$ f_{\mathrm{res}} = \frac{N}{2T} $$

の信号に応答します。

XY8 は回転軸を $X, Y, X, Y, Y, X, Y, X$ と巡回させたCPMGです。この巡回は系列をパルス振幅誤差と離調誤差に対して頑健にします。数百発のパルス列は信号より速く制御誤差を蓄積するので、これは重要です。フィルター関数はまったく変わりません。$\tilde{s}(f,T)$ に入るのはパルスの時刻だけなので、XY8-$N$ のフィルターは厳密にCPMG-$8N$ のものです。

共鳴応答についての中心的な結果は、推測するより述べるほうが簡単です。$f = f_{\mathrm{res}}$ におけるCPMG-$N$ では、$N-1$ 個の内部の半周期がそれぞれ $2T/(\pi N)$ を寄与し、両端の区間がそれぞれ $T/(\pi N)$ を寄与するので

$$ \left|\tilde{s}(f_{\mathrm{res}}, T)\right| = \frac{2T}{\pi} \qquad \text{（すべての } N \text{について）} $$

です。パルスを増やしても信号は買えません。 共鳴応答は系列の全長にしか依存しません。パルスを増やして買えるのは $T$ を長くする許可です。$T_2(N)$ が $N$ とともに伸びるからで、そして $\eta \propto 1/\sqrt{T_2}$ なのでデカップリングの利得はコヒーレンスの利得の平方根にすぎません。

Code Example 4: フィルター関数を厳密に

区分定数な変調のフィルター関数には閉形式があります。各区間のフーリエ変換を足すだけです。数値積分は不要であり、結果は第1章で引用した2つの閉形式に対して検証できます。

import numpy as np
import matplotlib.pyplot as plt

GAMMA = 28.025e9        # gamma_e / 2 pi、Hz/T


def cpmg_times(n_pulses, T):
    """CPMG-Nのpiパルス時刻 t_j = (j - 1/2) T / N、j = 1 ... N です。"""
    if n_pulses == 0:
        return np.array([])
    return (np.arange(1, n_pulses + 1) - 0.5) * T / n_pulses


def s_tilde(freqs, T, pulse_times):
    """変調関数 s(t) = +-1 の厳密なフーリエ変換です。

    s(t) はパルス時刻ごとに符号を反転するので、その変換は区分定数な区間に
    ついての閉じた和になり、数値積分を必要としません。
    """
    edges = np.concatenate(([0.0], np.asarray(pulse_times, float), [T]))
    signs = (-1.0) ** np.arange(len(edges) - 1)
    f = np.atleast_1d(np.asarray(freqs, float))
    out = np.zeros(f.shape, dtype=complex)
    nz = f != 0.0
    w = -2j * np.pi * f[nz]
    for k, sgn in enumerate(signs):
        out[nz] += sgn * (np.exp(w * edges[k + 1]) - np.exp(w * edges[k])) / w
    out[~nz] = np.sum(signs * np.diff(edges))
    return out if np.ndim(freqs) else out[0]


# --- 第1章の2つの閉形式に対する検証 ----------------------------------------
T = 10.0e-6
f = np.array([0.03, 0.1, 0.37, 1.0, 2.5]) / T
fid = np.abs(s_tilde(f, T, [])) ** 2
echo = np.abs(s_tilde(f, T, cpmg_times(1, T))) ** 2
print("量子ハードウェア入門第1章の閉形式に対するフィルター関数の検証")
print(f"{'f T':>8}{'|s~_FID|^2 num':>17}{'analytic':>15}"
      f"{'|s~_echo|^2 num':>18}{'analytic':>15}")
print("-" * 73)
for fi, a, b in zip(f, fid, echo):
    ana_f = np.sin(np.pi * fi * T) ** 2 / (np.pi * fi) ** 2
    ana_e = 4 * np.sin(np.pi * fi * T / 2) ** 4 / (np.pi * fi) ** 2
    print(f"{fi*T:>8.2f}{a:>17.8e}{ana_f:>15.8e}{b:>18.8e}{ana_e:>15.8e}")

# --- CPMG / XY8 フィルターの共鳴 -------------------------------------------
# CPMG-Nの変調は周期 2T/N の矩形波なので、f_res = N/(2T) の信号は打ち消され
# ずに位相を蓄積します。XY8-N のパルス*時刻*はCPMG-8Nと厳密に同一であり、
# X軸とY軸の交替が買うのはパルス誤差への耐性であって別のフィルターではありません。
print(f"\n共鳴応答、系列全長 T = {T*1e6:.1f} us")
print(f"{'sequence':>12}{'N':>5}{'f_res (MHz)':>13}{'|s~(f_res)|/T':>15}"
       f"{'2/pi':>9}{'max|s~|/T':>12}{'f_max/f_res':>13}")
print("-" * 79)
grid = np.linspace(1e2, 40e6, 400001)
for label, N in [("FID", 0), ("Hahn echo", 1), ("CPMG-2", 2), ("CPMG-8", 8),
                 ("XY8-1", 8), ("XY8-4", 32), ("CPMG-64", 64)]:
    tp = cpmg_times(N, T)
    if N == 0:
        print(f"{label:>12}{N:>5}{'-':>13}{'-':>15}{'-':>9}"
              f"{abs(s_tilde(1e-9, T, tp))/T:>12.6f}{'DC':>13}")
        continue
    f_res = N / (2 * T)
    amp = abs(s_tilde(f_res, T, tp)) / T
    prof = np.abs(s_tilde(grid, T, tp)) / T
    j = int(np.argmax(prof))
    print(f"{label:>12}{N:>5}{f_res/1e6:>13.4f}{amp:>15.8f}{2/np.pi:>9.6f}"
          f"{prof[j]:>12.6f}{grid[j]/f_res:>13.5f}")
print("  共鳴振幅はどのNでも厳密に 2T/pi です。パルスを増やしても信号は増えず、")
print("  増えるのは T を長くする許可です。")

# --- 高調波と奇数次のみの規則 ----------------------------------------------
N = 8
tp = cpmg_times(N, T)
f_res = N / (2 * T)
print(f"\nCPMG-{N} の高調波応答（矩形波変調）:")
print(f"{'k':>4}{'k f_res (MHz)':>15}{'|s~|/T':>12}{'2/(pi k)':>11}")
print("-" * 42)
for k in range(1, 8):
    print(f"{k:>4}{k*f_res/1e6:>15.4f}{abs(s_tilde(k*f_res, T, tp))/T:>12.8f}"
          f"{2/(np.pi*k) if k % 2 else 0.0:>11.8f}")
print("  偶数次高調波には感度がなく、奇数次は 1/k で応答します。したがって")
print("  3 f_res の信号は振幅3分の1で記録され、共鳴上のより弱い信号と区別が")
print("  つきません。")

# --- Parseval検査と直流での盲点 --------------------------------------------
fine = np.linspace(0.0, 400e6, 4000001)
for label, N in [("FID", 0), ("Hahn echo", 1), ("CPMG-8", 8)]:
    prof = np.abs(s_tilde(fine, T, cpmg_times(N, T))) ** 2
    integ = np.trapezoid(prof, fine)
    print(f"  {label:<10} f>0 での |s~|^2 の積分 = {integ/T:.6f} T "
          f"(Parseval: 0.500000 T);  |s~(0)|/T = "
          f"{abs(s_tilde(0.0, T, cpmg_times(N, T)))/T:.6f}")

fig, ax = plt.subplots(figsize=(7, 4))
fg = np.linspace(1e3, 20e6, 60001)
for label, N in [("FID", 0), ("Hahn echo", 1), ("CPMG-8", 8), ("XY8-4", 32)]:
    ax.loglog(fg / 1e6, np.abs(s_tilde(fg, T, cpmg_times(N, T)))**2 / T**2,
              lw=1.0, label=label)
ax.set_xlabel("frequency (MHz)"); ax.set_ylabel(r"$|\tilde{s}(f,T)|^2/T^2$")
ax.set_ylim(1e-8, 2); ax.legend(fontsize=8)
ax.set_title(f"Filter functions, T = {T*1e6:.0f} us")
plt.tight_layout()
plt.show()

量子ハードウェア入門第1章の閉形式に対するフィルター関数の検証
     f T   |s~_FID|^2 num       analytic   |s~_echo|^2 num       analytic
-------------------------------------------------------------------------
    0.03   9.97042623e-11 9.97042623e-11    2.21737562e-13 2.21737562e-13
    0.10   9.67531209e-11 9.67531209e-11    2.42711308e-12 2.42711308e-12
    0.37   6.23375846e-11 6.23375846e-11    2.68979001e-11 2.68979001e-11
    1.00   1.51957436e-43 1.51957436e-43    4.05284735e-11 4.05284735e-11
    2.50   1.62113894e-12 1.62113894e-12    1.62113894e-12 1.62113894e-12

共鳴応答、系列全長 T = 10.0 us
    sequence    N  f_res (MHz)  |s~(f_res)|/T     2/pi   max|s~|/T  f_max/f_res
-------------------------------------------------------------------------------
         FID    0            -              -        -    1.000000           DC
   Hahn echo    1       0.0500     0.63661977 0.636620    0.724611      1.48400
      CPMG-2    2       0.1000     0.63661977 0.636620    0.663953      1.14800
      CPMG-8    8       0.4000     0.63661977 0.636620    0.638568      1.01075
       XY8-1    8       0.4000     0.63661977 0.636620    0.638568      1.01075
       XY8-4   32       1.6000     0.63661977 0.636620    0.636743      1.00068
     CPMG-64   64       3.2000     0.63661977 0.636620    0.636650      1.00018
  共鳴振幅はどのNでも厳密に 2T/pi です。パルスを増やしても信号は増えず、
  増えるのは T を長くする許可です。

CPMG-8 の高調波応答（矩形波変調）:
   k  k f_res (MHz)      |s~|/T   2/(pi k)
------------------------------------------
   1         0.4000  0.63661977 0.63661977
   2         0.8000  0.00000000 0.00000000
   3         1.2000  0.21220659 0.21220659
   4         1.6000  0.00000000 0.00000000
   5         2.0000  0.12732395 0.12732395
   6         2.4000  0.00000000 0.00000000
   7         2.8000  0.09094568 0.09094568
  偶数次高調波には感度がなく、奇数次は 1/k で応答します。したがって
  3 f_res の信号は振幅3分の1で記録され、共鳴上のより弱い信号と区別が
  つきません。
  FID        f>0 での |s~|^2 の積分 = 0.499987 T (Parseval: 0.500000 T);  |s~(0)|/T = 1.000000
  Hahn echo  f>0 での |s~|^2 の積分 = 0.499962 T (Parseval: 0.500000 T);  |s~(0)|/T = 0.000000
  CPMG-8     f>0 での |s~|^2 の積分 = 0.499785 T (Parseval: 0.500000 T);  |s~(0)|/T = 0.000000

注目すべき点。 最初の表が検証です。区間和による変換は $\sin^2(\pi fT)/(\pi f)^2$ と $4\sin^4(\pi fT/2)/(\pi f)^2$ を $fT$ の2桁にわたって8桁まで再現し、$fT = 1$ におけるFIDフィルターの厳密な零点も含みます。コードは $1.5\times10^{-43}$ を返しますが、これは $10^{-11}$ 程度の量に対する浮動小数点のゼロです。

共鳴の表が主結果であり、2度読む価値があります。共鳴振幅はCPMG-2、CPMG-8、XY8-1、XY8-4、CPMG-64 のいずれでも $0.63661977\,T$ で、$2/\pi = 0.636620$ です。XY8-1 と CPMG-8 は同一の数値を出し、軸の巡回がフィルターに見えないことを確認します。最後の2列が唯一の注意点を記録しています。$|\tilde{s}|$ の大域的最大値は共鳴値をわずかに上回り、Hahnエコーでは $f = 1.484 f_{\mathrm{res}}$ で13.8%上、$N = 64$ では0.02%まで落ちます。単一の $\pi$ パルスは貧弱な櫛であり、長いCPMG列は良い櫛なのです。

高調波の表は実務上の警告です。奇数次高調波の応答は $1/k$ で落ち — $3f_{\mathrm{res}}$ で $0.2122\,T$、$5f_{\mathrm{res}}$ で $0.1273\,T$ — 偶数次では恒等的に消えるので、共鳴の3倍にある実信号は振幅3分の1で記録され、$T$ を変えない限り共鳴上のより弱い信号と区別がつきません。XY8で行われるナノスケールNMR測定はすべてこれと戦わなければなりません。Parseval検査が輪を閉じます。正の周波数における $|\tilde{s}|^2$ の積分はどの系列でも5桁まで $T/2$ であり、$|\tilde{s}(0)|$ は自由誘導で厳密に $T$、$\pi$ パルスを含むものでは厳密にゼロです。総感度は保存され、系列はそれをどこに置くかだけを決めます。

AC感度

共鳴時の位相は $\varphi = 2\pi(\gamma_e/2\pi)B_{\mathrm{AC}}\cdot 2T/\pi = 4(\gamma_e/2\pi)B_{\mathrm{AC}}T$ なので、この傾きで第2.2節の導出を繰り返すと

$$ \eta_{\mathrm{AC}}(T) = \frac{\sqrt{T + t_d}}{4(\gamma_e/2\pi)\, C_0\, T}\, e^{(T/T_2)^p} $$

を得ます。Ramseyの式と違うのは分母の $2\pi \to 4$ — $\pi/2$ の復調ペナルティ — と指数の $T_2^\ast \to T_2$ だけであり、利得のすべては後者にあります。

指数 $p$ は自由ではありません。量子ハードウェア入門第1章は、スペクトル密度 $S(f) = A/f^\alpha$ のノイズがガウス型のエコー減衰を生むことを確立しました。Code Example 5 はそれを定量化します。自己相似なフィルターの下での位相の平均二乗は

$$ \left\langle \varphi^2 \right\rangle = \left(2\pi\frac{\gamma_e}{2\pi}\right)^2 A\, T^{1+\alpha}\, I_N(\alpha), \qquad I_N(\alpha) = \int_0^\infty u^{-\alpha}\,\frac{\left|\tilde{s}(u/T, T)\right|^2}{T^2}\,\mathrm{d}u $$

であり、$u = fT$、$I_N(\alpha)$ は純粋な数です。$|\tilde{s}(f,T)|^2/T^2$ が $fT$ だけの関数だからです。したがってコヒーレンスは $\exp[-(T/T_2)^{1+\alpha}]$ で減衰し、$1/f$ ノイズでは $p = 2$ です。そして $T_2^{1+\alpha} \propto 1/I_N$ なので、コヒーレンス時間のパルス数依存性は1つの積分のパルス数依存性になります。

ノイズ分光: 同じフィルターを逆に読む

フィルターの2つの用途を並べて述べる価値があります。環境の記述を変えただけの同じ計算だからです。

周波数 $f$ の単色信号 ノイズスペクトル $S(f)$
磁場が何か $B_{\mathrm{AC}}\sin(2\pi f t + \phi)$ 定常、$S(f)$ で記述される
系列が返すもの 位相 $\varphi = 2\pi\gamma B_{\mathrm{AC}}|\tilde{s}(f,T)|$ 減衰 $\exp(-\langle\varphi^2\rangle/2)$
関わる積分 1つの周波数における $\tilde{s}$ $\int S(f)|\tilde{s}(f,T)|^2\mathrm{d}f$
掃引するもの $f$ を探すための $T$ $f_{\mathrm{res}}$ を動かすための $N$
分かること 単色成分の振幅 母材のスペクトル密度

第2列が姉妹コースのノイズ分光の手法であり、磁力計が同時に特性評価の道具である理由です。フィルターを $f_{\mathrm{res}}$ のデルタ関数で近似すると実用的な反転式

$$ \left\langle \varphi^2 \right\rangle \approx \left(2\pi\frac{\gamma_e}{2\pi}\right)^2 S_B(f_{\mathrm{res}})\, \kappa\, T $$

が得られます。$\kappa$ は Code Example 5 が評価する純粋な数です。測定された減衰の族を数桁にわたって $S_B(f)$ に反転するのは、そのあとは算術です。

Code Example 5: AC感度とノイズ分光の双対性

"""第2章 Code Example 5: AC感度と、同じフィルターをノイズ分光計として読む
双対性。Code Example 4 の続き（同一セッション）。"""

from scipy.optimize import minimize_scalar


def eta_ac(T_seq, T2, C0=1.0, t_dead=0.0, p=2.0, gamma=GAMMA):
    """継続時間 T_seq の系列に対するAC感度を T/sqrt(Hz) で返します。

    共鳴時の位相は phi = 2 pi gamma B |s~| = 4 gamma B T_seq なので、分母は
    Ramseyの 2 pi gamma tau ではなく 4 gamma T_seq になります。
    """
    return (np.sqrt(T_seq + t_dead) * np.exp((T_seq / T2) ** p)
            / (4.0 * gamma * C0 * T_seq))


# --- 無次元のフィルター積分 -------------------------------------------------
# S_B(f) = A / f^alpha に対して位相の平均二乗は
#   <phi^2> = (2 pi gamma)^2 A T^(1+alpha) I_N(alpha)、
#   I_N(alpha) = u^(-alpha) |s~(u/T, T)|^2 / T^2 の u についての積分、u = f T
# であり、これは純粋な数です。フィルター関数が u について自己相似だからです。
# したがって減衰は exp(-(T/T2)^(1+alpha)) となり、1/fノイズではガウス型で、
# T2 は I_N から決まります。
def filter_integrals(N, alphas=(1.0, 2.0), T=1.0):
    u_lo = np.logspace(-6, -2, 4001)
    u_hi = np.arange(1e-2, 100.0 * max(N, 1) + 1e-2, 5e-3)
    u = np.concatenate((u_lo, u_hi))
    prof = np.abs(s_tilde(u / T, T, cpmg_times(N, T))) ** 2 / T**2
    return {a: np.trapezoid(u**(-a) * prof, u) for a in alphas}


print("無次元のフィルター積分 I_N(alpha) とそれが導く T2(N)")
print(f"{'N':>5}{'I_N(1)':>12}{'N I_N(1)':>11}{'I_N(2)':>12}{'N^2 I_N(2)':>13}"
      f"{'T2(N)/T2(1) a=1':>17}{'sqrt(N)':>10}")
print("-" * 80)
Ns = [1, 2, 4, 8, 16, 32]
I1, I2 = {}, {}
for N in Ns:
    d = filter_integrals(N)
    I1[N], I2[N] = d[1.0], d[2.0]
    # T2 ~ I_N^(-1/(1+alpha))
    print(f"{N:>5}{I1[N]:>12.6f}{N*I1[N]:>11.6f}{I2[N]:>12.6f}"
          f"{N**2*I2[N]:>13.6f}{(I1[1]/I1[N])**0.5:>17.5f}"
          f"{np.sqrt(N):>10.5f}")
c1 = np.polyfit(np.log(Ns), np.log([I1[N] for N in Ns]), 1)[0]
c2 = np.polyfit(np.log(Ns), np.log([I2[N] for N in Ns]), 1)[0]
print(f"  I_N の N に対するフィット指数 : alpha=1 -> {c1:.4f}（漸近値 -1）、"
      f"alpha=2 -> {c2:.4f}（漸近値 -2）")
print(f"  したがって T2(N) ~ N^(alpha/(alpha+1)): "
      f"{-c1/2:.4f} 対 0.5000、{-c2/3:.4f} 対 0.6667")
print(f"  N^2 I_N(2) はどのNでも厳密に pi^2/6 = {np.pi**2/6:.6f} です。奇数次")
print("  高調波の櫛から I_N(2) = (16/(pi^2 N^2)) sum_odd 1/k^4")
print("  = (16/(pi^2 N^2))(pi^4/96) = pi^2/(6 N^2) となり、T2 は N^(2/3) で伸びます。")

# --- ノイズ分光: 同じ積分を逆に使う ----------------------------------------
# フィルターを f_res = N/(2T) のデルタ関数で近似すると
#   <phi^2> = (2 pi gamma)^2 S_B(f_res) kappa T、kappa = N I_N(1) / 2
# となります。
kappa_inf = 4 / np.pi**2 * (7.0 / 8.0) * 1.2020569031595943      # (4/pi^2) zeta(3) 7/8
print("\n減衰を S_B(f) に反転するために使うデルタフィルター近似:")
print(f"{'N':>5}{'kappa = N I_N(1)/2':>21}{'large-N limit':>15}")
print("-" * 41)
for N in Ns:
    print(f"{N:>5}{N*I1[N]/2:>21.6f}{kappa_inf:>15.6f}")
print("  フィルターが真の櫛になると kappa は N に依存しなくなります。極限は")
print("  (4/pi^2) x (7/8) zeta(3) です。Parseval重み T の半分、矩形波パワーの")
print("  基本波成分 8/pi^2、そして 1/f ノイズが上乗せする奇数次高調波の 1/k^3 の")
print("  和の積です。Hahnエコー（N=1）は単一パルスが櫛として貧弱なため19%低い値です。")

# --- AC感度 ----------------------------------------------------------------
T2_1 = 100e-6                      # Hahnエコーの T2、1/f律速
sigma_R = np.sqrt(1.0 + 2.0 * (0.030 + 0.020) / (0.030 - 0.020) ** 2)
print(f"\nAC磁気計測、エコー T2 = {T2_1*1e6:.0f} us、sigma_R = "
      f"{sigma_R:.2f}、1/fノイズなので p = 2")
print(f"{'sequence':>10}{'N':>5}{'T2(N) (us)':>12}{'T_opt (us)':>12}"
      f"{'f_res (kHz)':>13}{'eta_PN (nT)':>13}{'eta (nT)':>11}")
print("-" * 76)
etas = []
for label, N in [("echo", 1), ("CPMG-2", 2), ("XY8-1", 8), ("XY8-2", 16),
                 ("XY8-4", 32)]:
    T2N = T2_1 * (I1[1] / I1[N]) ** 0.5
    r = minimize_scalar(eta_ac, bracket=(1e-9, T2N, 20 * T2N),
                        args=(T2N, 1.0, 0.0, 2.0))
    etas.append((N, r.fun * sigma_R))
    print(f"{label:>10}{N:>5}{T2N*1e6:>12.3f}{r.x*1e6:>12.3f}"
          f"{N/(2*r.x)/1e3:>13.3f}{r.fun*1e9:>13.4f}"
          f"{r.fun*sigma_R*1e9:>11.3f}")
sl = np.polyfit(np.log([e[0] for e in etas]), np.log([e[1] for e in etas]), 1)[0]
print(f"  どの行でも T_opt = T2(N)/2 であり、eta_min は "
      f"N^({sl:.4f}) でスケールします")
print(f"  予測値 -alpha/(2(alpha+1)) = {-1/4:.4f}。パルス数を倍にしても感度は")
print(f"  {2**0.25:.3f} 倍しか買えず、f_res は2倍上へ移ります。")

# --- 1つのフィルター、2つの使い方 ------------------------------------------
N, T_seq = 8, 50e-6
f_res = N / (2 * T_seq)
A_1f = 2.0 / ((2 * np.pi * GAMMA) ** 2 * T2_1**2 * I1[1])
phi_coh = 2 * np.pi * GAMMA * 100e-9 * abs(s_tilde(f_res, T_seq,
                                                   cpmg_times(N, T_seq)))
phi2_inc = (2 * np.pi * GAMMA) ** 2 * A_1f * T_seq**2 * I1[N]
print(f"\n1つのフィルター、2つの使い方。N = {N}、T = {T_seq*1e6:.0f} us "
      f"(f_res = {f_res/1e3:.1f} kHz):")
print(f"  共鳴上の 100 nT の単色信号     -> 位相 {phi_coh:.6f} rad")
print(f"  1/fバス、S_B = {A_1f/f_res:.4e} T^2/Hz -> rms位相 "
      f"{np.sqrt(phi2_inc):.6f} rad")
print(f"  位相における信号対雑音比       : "
      f"{phi_coh/np.sqrt(phi2_inc):.4f}")
print(f"  バスを磁場換算すると sqrt(<phi^2>)/(4 gamma T) = "
      f"{np.sqrt(phi2_inc)/(4*GAMMA*T_seq)*1e9:.2f} nT")
print("  単色信号は位相を与え、スペクトルは減衰を与えます。N を掃引すれば f_res")
print("  が走査され、後者の読み方が母材のスペクトルになります。")

無次元のフィルター積分 I_N(alpha) とそれが導く T2(N)
    N      I_N(1)   N I_N(1)      I_N(2)   N^2 I_N(2)  T2(N)/T2(1) a=1   sqrt(N)
--------------------------------------------------------------------------------
    1    0.693134   0.693134    1.644932     1.644932          1.00000   1.00000
    2    0.392433   0.784866    0.411233     1.644934          1.32900   1.41421
    4    0.204399   0.817596    0.102808     1.644934          1.84149   2.00000
    8    0.104363   0.834906    0.025702     1.644934          2.57712   2.82843
   16    0.052732   0.843707    0.006426     1.644934          3.62554   4.00000
   32    0.026504   0.848126    0.001606     1.644934          5.11391   5.65685
  I_N の N に対するフィット指数 : alpha=1 -> -0.9486（漸近値 -1）、alpha=2 -> -2.0000（漸近値 -2）
  したがって T2(N) ~ N^(alpha/(alpha+1)): 0.4743 対 0.5000、0.6667 対 0.6667
  N^2 I_N(2) はどのNでも厳密に pi^2/6 = 1.644934 です。奇数次
  高調波の櫛から I_N(2) = (16/(pi^2 N^2)) sum_odd 1/k^4
  = (16/(pi^2 N^2))(pi^4/96) = pi^2/(6 N^2) となり、T2 は N^(2/3) で伸びます。

減衰を S_B(f) に反転するために使うデルタフィルター近似:
    N   kappa = N I_N(1)/2  large-N limit
-----------------------------------------
    1             0.346567       0.426278
    2             0.392433       0.426278
    4             0.408798       0.426278
    8             0.417453       0.426278
   16             0.421853       0.426278
   32             0.424063       0.426278
  フィルターが真の櫛になると kappa は N に依存しなくなります。極限は
  (4/pi^2) x (7/8) zeta(3) です。Parseval重み T の半分、矩形波パワーの
  基本波成分 8/pi^2、そして 1/f ノイズが上乗せする奇数次高調波の 1/k^3 の
  和の積です。Hahnエコー（N=1）は単一パルスが櫛として貧弱なため19%低い値です。

AC磁気計測、エコー T2 = 100 us、sigma_R = 31.64、1/fノイズなので p = 2
  sequence    N  T2(N) (us)  T_opt (us)  f_res (kHz)  eta_PN (nT)   eta (nT)
----------------------------------------------------------------------------
      echo    1     100.000      50.000       10.000       1.6199     51.251
    CPMG-2    2     132.900      66.450       15.049       1.4051     44.457
     XY8-1    8     257.712     128.856       31.042       1.0091     31.925
     XY8-2   16     362.554     181.277       44.131       0.8507     26.916
     XY8-4   32     511.391     255.696       62.574       0.7163     22.663
  どの行でも T_opt = T2(N)/2 であり、eta_min は N^(-0.2367) でスケールします
  予測値 -alpha/(2(alpha+1)) = -0.2500。パルス数を倍にしても感度は
  1.189 倍しか買えず、f_res は2倍上へ移ります。

1つのフィルター、2つの使い方。N = 8、T = 50 us (f_res = 80.0 kHz):
  共鳴上の 100 nT の単色信号     -> 位相 0.560500 rad
  1/fバス、S_B = 1.1632e-19 T^2/Hz -> rms位相 0.274378 rad
  位相における信号対雑音比       : 2.0428
  バスを磁場換算すると sqrt(<phi^2>)/(4 gamma T) = 48.95 nT
  単色信号は位相を与え、スペクトルは減衰を与えます。N を掃引すれば f_res
  が走査され、後者の読み方が母材のスペクトルになります。

注目すべき点。 フィルター積分がこの例の数値的な核心です。$\alpha = 1$ について $N\,I_N(1)$ は $N = 1$ の0.693から $N = 32$ の0.848へ這い上がるので、漸近的に $I_N \propto N^{-1}$ であり $T_2(N) \propto N^{1/2}$ です。この範囲でのフィット指数は $0.4743$ で、量子ハードウェア入門の Code Example 7 がまったく別の方法 — フィルター積分ではなく明示的なノイズ軌跡 — で同じ予測値 $0.5$ に対して得た $0.476$ と比べるべき値です。2つの独立な経路が同じ有限 $N$ の不足を与えることは、どちらも正しいという良い兆候です。

$\alpha = 2$ では厳密なことが起きます。$N^2 I_N(2) = 1.644934 = \pi^2/6$ が $N = 1$ を含めてどの $N$ でも6桁まで成り立ちます。奇数次高調波の櫛がこれを見通しよくします — $I_N(2) = (16/\pi^2 N^2)\sum_{k\ \mathrm{odd}}k^{-4} = (16/\pi^2N^2)(\pi^4/96) = \pi^2/(6N^2)$ — そしてフィットなしに $T_2(N) \propto N^{2/3}$ を与えます。演習3で $N = 1$ の場合を解析的に検証してもらいます。

$\kappa$ の表はデルタフィルター近似を定量化します。$\kappa = N I_N(1)/2$ は $(4/\pi^2)\cdot\frac{7}{8}\zeta(3) = 0.426278$ に収束します。Parseval重み $T$ の半分、矩形波パワーの基本波成分 $8/\pi^2$、そして $1/f$ ノイズがなお届く高調波についての $\sum_{k\ \mathrm{odd}}k^{-3}$ の積です。Hahnエコーは極限より19%低く、XY8-4 では誤差が0.6%を切ります。単一の $\pi$ パルスによるノイズ分光には厳密な積分が必要で、長い列ならデルタ近似で十分です。感度の表が成果であり、意図的に地味です。

Hahnエコーから XY8-4 へ移ると $T_2$ は5.1倍伸び、$\eta$ は約束どおりその平方根の2.26倍改善して51.3から22.7 nT/$\sqrt{\mathrm{Hz}}$ になり、通過帯域は10 kHzから63 kHzへ移ります。フィットしたスケーリングは $N^{-0.237}$、予測は $-\alpha/(2(\alpha+1)) = -0.25$ です。32発のパルスが買うのは2倍です。動的デカップリングだけからそれ以上を約束する者は、フィルター関数が与えないものを約束しています。

最後のブロックは同じ系列に2つの読み方を並べます。80 kHzの100 nTの単色信号は0.5605 radの位相を与え、同じ母材の $1/f$ バスは0.2744 radのrms位相を寄与し、信号対雑音比は2.04です。磁場に換算すればバスはその系列にわたって49 nTに相当します。計測器とそれ自身のノイズ床が同じ積分で計算されるのです。

2.4 $T_1$ リラクソメトリ: GHzの窓

パルス系列が届かない周波数

第2.3節のすべては数十メガヘルツ以下に住んでいます。理由は初等的です。CPMG列が作れる最も高い通過帯域は $\pi$ パルスをどれだけ密に置けるかで決まり、$f_{\max} \approx 1/(2t_\pi)$ です。20 nsの $\pi$ パルスは系列を25 MHzで打ち止めにします。マグノンが伝播し、強磁性膜のスピン波が共鳴し、金属中のジョンソンノイズがパワーの大半を置き、ナノ秒の相関時間をもつ常磁性スピンがノイズを置くギガヘルツ帯には、パルス系列の道具箱では届きません。

NV中心はパルスをまったく使わずにその帯域に届きます。自身の遷移周波数が2.87 GHzだからです。$f_0$ の磁気ノイズは $m_s = 0$ からの遷移を駆動し、生じる緩和は $f_0$ におけるノイズパワーの直接測定です。これがリラクソメトリであり、本章で最も地味かつ最も応用範囲の広いNV手法です。

緩和速度からスペクトル密度へ

横磁場の1成分が片側パワースペクトル密度 $S_\perp(f)$（T$^2$/Hz、$\langle B_\perp^2\rangle = \int_0^\infty S_\perp\,\mathrm{d}f$ で規格化）をもち、$x$ と $y$ が無相関で等パワーだとします。相互作用は $H_{\mathrm{int}} = \hbar\gamma_e(B_x S_x + B_y S_y)$ で、$S_x$、$S_y$ はスピン1演算子、$m_s = 0$ と $m_s = \mp 1$ の間の行列要素の絶対値は $1/\sqrt{2}$ です。フェルミの黄金律から

$$ \Gamma_{0\to\mp1} = \int_{-\infty}^{\infty}\mathrm{d}t\; e^{i\omega_0 t}\left\langle \left(H_{\mathrm{int}}/\hbar\right)_{\mp1,0}(t)\, \left(H_{\mathrm{int}}/\hbar\right)^\ast_{\mp1,0}(0)\right\rangle = \gamma_e^2\left[\tfrac{1}{2}\tfrac{S_\perp(f_0)}{2} + \tfrac{1}{2}\tfrac{S_\perp(f_0)}{2}\right] = \frac{\gamma_e^2 S_\perp(f_0)}{2} $$

を得ます。$S_\perp(f_0)/2$ という因子は、実の定常過程の自己相関の $\omega_0 = 2\pi f_0$ におけるフーリエ変換を片側密度で表したものです。上下の速度 $\Gamma$ が対称なとき — 2.87 GHzと300 Kでは高温極限が優れた近似です — $m_s = 0$ の占有数は $\dot{P_0} = -2\Gamma P_0 + \Gamma(1 - P_0) = -3\Gamma(P_0 - \tfrac{1}{3})$ に従うので

$$ \boxed{\ \frac{1}{T_1} = 3\Gamma = \frac{3\gamma_e^2}{2}\, S_\perp(f_0) \qquad\Longleftrightarrow\qquad S_\perp(f_0) = \frac{2}{3\gamma_e^2 T_1}\ } $$

です。測定された $T_1$ は測定されたスペクトル密度です。ここでの $\gamma_e$ は角振動数の磁気回転比 $1.7609\times10^{11}$ rad s$^{-1}$T$^{-1}$ であり、円振動数側を使うと39倍間違えます。

固有の $T_1$ を決めるもの、そこに加わるもの

室温でのバルクNV中心の固有 $T_1$ はミリ秒程度で、2フォノンのラマン過程で決まります。これは格子の性質で強い温度依存性をもち、磁気的ではまったくありません。約100 K以下では桁で伸びます。リラクソメトリ実験が検出するものはすべてその固有速度への付加として現れるので、固有 $T_1$ が手法のノイズ床であり、冷却が改善する理由です。

プローブ周波数は同調可能です。$f_\pm = D \pm (\gamma_e/2\pi)B_\parallel$ なので、バイアス磁場は下側の枝を2.87 GHzから $B = D/(\gamma_e/2\pi) = 102.4$ mT のゼロまで、上側の枝を5.7 GHz以上まで掃引します。数ギガヘルツの同調範囲をもつスペクトラムアナライザであり、下側の枝がゼロを横切る点が基底状態レベル反交差です。ここでNVは環境中の裸の電子スピンと縮退し、交差緩和が最も強くなります。

Code Example 6: 分光としてのリラクソメトリ

import numpy as np

GAMMA_C = 28.025e9              # gamma_e / 2 pi、Hz/T
GAMMA = 2 * np.pi * GAMMA_C     # gamma_e、rad/s/T
D_GS = 2.870e9                  # ゼロ磁場分裂、Hz
MU_B = 9.2740100783e-24         # ボーア磁子、J/T
MU0_4PI = 1.0e-7                # mu_0 / 4 pi、T m / A


def s_perp_from_T1(T1):
    """横磁場1成分の片側PSDを T^2/Hz で返します。

    上下の遷移速度が対称なスピン1の m_s = 0 の占有数について
    1/T1 = 3 gamma^2 S_perp(f_0) / 2 が成り立ちます。
    """
    return 2.0 / (3.0 * GAMMA**2 * T1)


def T1_from_s_perp(S_perp):
    return 2.0 / (3.0 * GAMMA**2 * S_perp)


print("T1リラクソメトリ: NVをGHz帯のスペクトラムアナライザとして使う")
print(f"  gamma_e = {GAMMA:.6e} rad/s/T、3 gamma^2 / 2 = "
      f"{1.5*GAMMA**2:.6e} (rad/s)^2/T^2")
print(f"\n{'T1':>10}{'1/T1 (1/s)':>13}{'S_perp (T^2/Hz)':>18}"
      f"{'sqrt(S_perp)':>16}")
print("-" * 57)
for label, T1 in [("6 ms", 6.0e-3), ("1 ms", 1.0e-3), ("100 us", 1.0e-4),
                  ("10 us", 1.0e-5)]:
    S = s_perp_from_T1(T1)
    print(f"{label:>10}{1/T1:>13.4f}{S:>18.4e}"
          f"{np.sqrt(S)*1e12:>13.3f} pT/rtHz")

# --- 通過帯域の位置と、その同調のしかた ------------------------------------
print("\nプローブ周波数はNV自身の遷移であり、バイアス磁場で同調します:")
print(f"{'B (mT)':>9}{'f- = D - gamma B (GHz)':>25}"
      f"{'f+ = D + gamma B (GHz)':>25}")
print("-" * 59)
for B_mT in [0.0, 25.0, 50.0, 90.0, 102.4]:
    print(f"{B_mT:>9.1f}{(D_GS - GAMMA_C*B_mT*1e-3)/1e9:>25.4f}"
          f"{(D_GS + GAMMA_C*B_mT*1e-3)/1e9:>25.4f}")
print("  102.4 mT で下側の枝がゼロに達します。これが基底状態レベル反交差であり、")
print("  NVが環境中の裸の電子スピンと共鳴して交差緩和が最も強くなる点です。")

# --- 表面のゆらぐスピンバス ------------------------------------------------
# 桁数レベルの順問題モデルです。面密度 sigma_s の常磁性モーメントのシートが
# 距離 d にあり、各モーメントは相関時間 tau_c をもちます。磁場の平均二乗は
# (mu_0 mu / 4 pi)^2 sigma_s をシート上で積分したものであり、積分
# 2 pi rho drho / (rho^2 + d^2)^3 = pi / (2 d^4) を使います。
def bath_noise(sigma_s, d, tau_c, moment=MU_B, f=D_GS):
    b2 = (MU0_4PI * moment) ** 2 * sigma_s * np.pi / (2.0 * d**4)
    lorentz = 4.0 * tau_c / (1.0 + (2 * np.pi * f * tau_c) ** 2)
    return b2, b2 * lorentz


print("\nゆらぐ表面スピンのシートを f_0 = 2.870 GHz で見る")
print(f"{'sigma_s (1/nm^2)':>18}{'d (nm)':>8}{'tau_c (ns)':>12}"
      f"{'B_rms (uT)':>12}{'S_perp (T^2/Hz)':>18}{'added 1/T1 (1/s)':>18}")
print("-" * 86)
for sig_nm, d_nm, tc_ns in [(0.1, 5.0, 1.0), (0.1, 10.0, 1.0),
                            (0.1, 5.0, 0.1), (1.0, 5.0, 1.0),
                            (0.1, 20.0, 1.0)]:
    b2, S = bath_noise(sig_nm * 1e18, d_nm * 1e-9, tc_ns * 1e-9)
    print(f"{sig_nm:>18.2f}{d_nm:>8.1f}{tc_ns:>12.2f}"
          f"{np.sqrt(b2)*1e6:>12.4f}{S:>18.4e}{1.5*GAMMA**2*S:>18.4f}")
print("  1/d^4 依存性がリラクソメトリの設計則そのものです。距離を半分にすると")
print("  検出されるノイズパワーは16倍になります。")

# --- 緩和測定が分解できるもの ----------------------------------------------
T1_0 = 6.0e-3
G0 = 1.0 / T1_0
print(f"\n固有の T1 = {T1_0*1e3:.0f} ms（Gamma_0 = {G0:.2f} /s）に対して"
      f"付加された速度を分解する")
print(f"{'rate resolution':>17}{'delta Gamma (1/s)':>19}"
      f"{'delta S_perp (T^2/Hz)':>23}{'field equivalent':>19}")
print("-" * 78)
for frac in [0.10, 0.01, 0.001]:
    dG = frac * G0
    dS = 2.0 * dG / (3.0 * GAMMA**2)
    print(f"{frac*100:>15.1f} %{dG:>19.4f}{dS:>23.4e}"
          f"{np.sqrt(dS)*1e12:>16.3f} pT/rtHz")

T1リラクソメトリ: NVをGHz帯のスペクトラムアナライザとして使う
  gamma_e = 1.760863e+11 rad/s/T、3 gamma^2 / 2 = 4.650956e+22 (rad/s)^2/T^2

        T1   1/T1 (1/s)   S_perp (T^2/Hz)    sqrt(S_perp)
---------------------------------------------------------
      6 ms     166.6667        3.5835e-21       59.862 pT/rtHz
      1 ms    1000.0000        2.1501e-20      146.632 pT/rtHz
    100 us   10000.0000        2.1501e-19      463.691 pT/rtHz
     10 us  100000.0000        2.1501e-18     1466.320 pT/rtHz

プローブ周波数はNV自身の遷移であり、バイアス磁場で同調します:
   B (mT)   f- = D - gamma B (GHz)   f+ = D + gamma B (GHz)
-----------------------------------------------------------
      0.0                   2.8700                   2.8700
     25.0                   2.1694                   3.5706
     50.0                   1.4688                   4.2713
     90.0                   0.3478                   5.3922
    102.4                   0.0002                   5.7398
  102.4 mT で下側の枝がゼロに達します。これが基底状態レベル反交差であり、
  NVが環境中の裸の電子スピンと共鳴して交差緩和が最も強くなる点です。

ゆらぐ表面スピンのシートを f_0 = 2.870 GHz で見る
  sigma_s (1/nm^2)  d (nm)  tau_c (ns)  B_rms (uT)   S_perp (T^2/Hz)  added 1/T1 (1/s)
--------------------------------------------------------------------------------------
              0.10     5.0        1.00     14.7024        2.6508e-21          123.2878
              0.10    10.0        1.00      3.6756        1.6568e-22            7.7055
              0.10     5.0        0.10     14.7024        2.0336e-20          945.8116
              1.00     5.0        1.00     46.4930        2.6508e-20         1232.8782
              0.10    20.0        1.00      0.9189        1.0355e-23            0.4816
  1/d^4 依存性がリラクソメトリの設計則そのものです。距離を半分にすると
  検出されるノイズパワーは16倍になります。

固有の T1 = 6 ms（Gamma_0 = 166.67 /s）に対して付加された速度を分解する
  rate resolution  delta Gamma (1/s)  delta S_perp (T^2/Hz)   field equivalent
------------------------------------------------------------------------------
           10.0 %            16.6667             3.5835e-22          18.930 pT/rtHz
            1.0 %             1.6667             3.5835e-23           5.986 pT/rtHz
            0.1 %             0.1667             3.5835e-24           1.893 pT/rtHz

注目すべき点。 最初の表がこの手法の換算表です。6 msの $T_1$ は2.87 GHzにおける横磁場ノイズ $3.58\times10^{-21}$ T$^2$/Hz、振幅60 pT/$\sqrt{\mathrm{Hz}}$ に対応します。この数値は第2.3節で最良のパルス系列が与えた22 nT/$\sqrt{\mathrm{Hz}}$ と比べるべきです。リラクソメトリは、動的デカップリングが磁場単色成分に対して持つ感度より、磁場ノイズに対してほぼ3桁高い感度をもちます。数百マイクロ秒ではなくミリ秒積分するからであり、緩和速度が2次ではなく1次の効果だからです。その代価は位相情報の完全な喪失です。リラクソメトリはパワーを測り、振幅も方向も決して測りません。

表面スピンの表が設計則を示します。距離 $d$ にあるモーメントのシートによる磁場の平均二乗は $1/d^4$ でスケールし — 双極子磁場が $1/d^3$ で落ち、届く範囲のスピン数が $d^2$ で増えるのでパワーは $d^{-6}\cdot d^2$ となります — 距離を半分にすると検出ノイズパワーは16倍になります。1平方ナノメートルあたり0.1個の不対スピンのシートが5 nmにあり相関時間1 nsでゆらぐと、固有値が167 s$^{-1}$ の緩和速度に123 s$^{-1}$ を加えます。74%の変化で、検出は容易です。同じシートを20 nmに動かすと0.48 s$^{-1}$、0.3%の変化で、注意深い測定が分解できる限界付近です。意図するかどうかにかかわらず、リラクソメトリは表面の手法なのです。

モデル中のローレンツ関数が実際に働いています。0.1 nsの行は1 nsの行と同じrms磁場をもちながら8倍近い緩和を生みます。相関時間を短くするとスペクトル重みが2.87 GHzまで上がるからです。リラクソメトリはバスの大きさよりその動力学に感度をもつので、静的あるいは非常に遅いバスはどれだけ強くても見えません。静的な磁場を見て速いものに盲目な第2.2節のRamsey測定の、ちょうど相補です。分解能の表は $\Gamma$ の相対精度を検出可能なスペクトル密度に換算します。6 msの $T_1$ に対する1%の分解能は2.87 GHzで6 pT/$\sqrt{\mathrm{Hz}}$ の付加磁場ノイズに相当します。これが本節冒頭で述べた分業を締めくくります。パルス系列がキロヘルツから数十メガヘルツを、リラクソメトリがギガヘルツを担い、材料の磁気ノイズの完全な特性評価は同じ $S(f)$ への両方の窓を使います。

2.5 磁性体をイメージングする

2つの配置

ここまではセンサ1つと磁場1つの話でした。イメージングには空間的な配置が必要で、それは2種類あります。

走査NV顕微鏡は単一の中心をダイヤモンド探針の先端に置き、試料上をラスタ走査します。分解能は光の波長ではなく距離 $d$ で決まります。磁場を局所的に測るからです。浅いイオン注入と探針加工が与えるのは20から100 nmです。代価はスループットで、単一中心の感度で一度に1ピクセルです。

広視野アンサンブルイメージングは表面近傍のNV層をもつダイヤモンドの上に試料を置き、蛍光をカメラに結像します。どのピクセルも $\eta$ が $\sqrt{N}$ 改善した独立の磁力計なので視野全体が一度に取得できますが、分解能は回折限界の数百ナノメートルで、距離はNV層の厚みです。トレードオフは第1章が地図にしたものです。

何を見るべきか

以下は、NV磁力計が扱う材料科学の磁気現象を、生じる磁場の大きさのおおよその順に並べたものです。

面直磁化膜における磁区と磁壁。 一様に磁化した無限シートは外部磁場を作らず、漏れ磁場はすべて端と壁から来ます。面モーメント $m_s = M_s t$ の膜中の壁は磁荷の線として働き、距離 $d$ で $\mu_0 m_s/(2\pi d)$ を生みます。典型的な遷移金属強磁性体の1ナノメートルなら数十ナノメートルでミリテスラであり、NVの尺度では巨大で、第2.1節の線形ゼーマン近似が破れるほど大きい値です。したがって磁壁イメージングは信号としては容易で解釈としては難しく、有用な測定は壁の構造 — 内部の磁化プロファイル、Dzyaloshinskii-Moriya壁のカイラリティ、スキルミオンのプロファイル — についてのものです。

2次元ファンデルワールス磁性体。 単層磁性体は1平方ナノメートルあたり10から20ボーア磁子程度をもち、金属膜より1桁から2桁小さく、数十ナノメートルで数百マイクロテスラの端磁場を生みます。代替手段が貧弱なので科学的価値は高いのです。試料は数マイクロメートルのフレークで、その総モーメントはバルク磁力計の分解能をはるかに下回りますが、NV磁気計測は単層の絶対磁化を基準なしに測ります。

反強磁性体。 補償された反強磁性体は正味の漏れ磁場を作らず、まさにそれが難しい理由であり興味深い理由です。NV磁気計測が見るのは不完全さ、すなわち表面や磁壁における非補償モーメントと、マルチフェロイックの磁気電気応答です。信号は強磁性体より1桁から3桁小さい。反強磁性共鳴が第2.4節のギガヘルツの窓に入りうるので、リラクソメトリが第2の経路を提供します。

電流分布。 磁気画像は反転を経れば電流の画像であり、動作中の素子の内部でどこを電流が流れているかを見る他の方法はありません。グラフェン中の流体力学的な電子流、コンタクトでの電流集中、抵抗変化素子中のフィラメント伝導、電池電極中の不均一な電流です。難しいのは反転で、Code Example 7 がそれを定量化します。

ゆらぐスピンとGHzノイズ。 第2.4節のリラクソメトリです。常磁性ラベル、表面スピンバス、隣接膜中のマグノン、金属中のジョンソンノイズ。

その一覧の試料の多くはスピントロニクスが供給します。本道場の スピントロニクス入門 は、ここで測る磁場を生む磁性 — 磁壁、Dzyaloshinskii-Moriya相互作用、スピン軌道トルク、反強磁性秩序、マグノン輸送 — を扱います。両者を併せて読めば、センサと試料が同じ絵の中に入ります。

距離が分解能を指数関数的に決める

すべてを支配する結果が1つあります。平面状のソースの上には電流も電荷もないので磁場はラプラス方程式を満たし、面内波数 $k$ のフーリエ成分は $e^{-kd}$ で減衰します。周期 $L$ のソース構造は距離 $d$ で $e^{-2\pi d/L}$ だけ減衰します。$d$ より細かい構造は $e^{-2\pi} \approx 2\times10^{-3}$ 以上の割合で抑制されます。

その帰結は両方向に働きます。順方向では、走査磁力計の分解能は走査格子をどれだけ細かくしても距離であることを意味します。逆方向では、磁場マップからソースを復元するには測定されたフーリエ成分に $e^{+kd}$ を掛ける必要があり、これはノイズを同じ倍率で掛けることを意味します。逆問題は $kd$ について指数関数的に悪条件であり、平均化では指数は変わりません。$k \approx 1/d$ でのカットオフやTikhonov罰則という形の正則化は数値上の便宜ではなく、$k \approx 1/d$ を超える情報は測定されていないという主張なのです。

Code Example 7: 漏れ磁場と、それを反転するコスト

import numpy as np

MU0_4PI = 1.0e-7                # mu_0 / 4 pi、T m / A
MU0 = 4 * np.pi * 1e-7          # T m / A
MU_B = 9.2740100783e-24         # J/T
MU_P = 1.41060679736e-26        # J/T

ETA_DC = 419e-9                 # 単一NVのRamsey感度、Code Example 3
ETA_AC = 51.3e-9                # 単一NVのHahnエコーAC感度、Code Example 5


def b_dipole(moment, d):
    """点モーメントの軸上、距離 d における磁場 mu_0 m / (2 pi d^3) です。"""
    return 2.0 * MU0_4PI * moment / d**3


def b_edge(sheet_moment, d):
    """面直磁化したシートの端から距離 d の位置における磁場です。

    端は線密度 lambda = m_s（面モーメント×単位長さ）の磁荷の線とみなせ、
    mu_0 lambda / (2 pi d) を与えます。
    """
    return MU0 * sheet_moment / (2.0 * np.pi * d)


def b_wire(current, d):
    """直線電流の磁場 mu_0 I / (2 pi d) です。"""
    return MU0 * current / (2.0 * np.pi * d)


d = 20e-9
print(f"距離 {d*1e9:.0f} nm における漏れ磁場と、その検出に要するコスト")
print(f"  単一NVのDC感度 {ETA_DC*1e9:.0f} nT/sqrt(Hz)、"
      f"AC {ETA_AC*1e9:.1f} nT/sqrt(Hz)")

# 面モーメントです。Co相当: M_s = 1.4e6 A/m を厚さ1 nmに。単層vdW磁性体:
# 1 nm^2あたり約15 mu_B。反強磁性体の非補償表面: (0.4 nm)^2 の表面セルあたり
# 0.01 mu_B。
m_co = 1.4e6 * 1e-9
m_vdw = 15.0 * MU_B / 1e-18
m_afm = 0.01 * MU_B / (0.4e-9) ** 2
sources = [
    ("single electron spin",       b_dipole(MU_B, d)),
    ("single proton",              b_dipole(MU_P, d)),
    ("1 nm Co-like film, edge",    b_edge(m_co, d)),
    ("monolayer vdW magnet, edge", b_edge(m_vdw, d)),
    ("uncompensated AFM surface",  b_edge(m_afm, d)),
    ("1 uA current line",          b_wire(1e-6, d)),
    ("1 nA current line",          b_wire(1e-9, d)),
]
print(f"\n{'source':>28}{'B (T)':>13}{'B':>14}{'t for SNR 1, DC':>18}"
      f"{'AC':>14}")
print("-" * 87)
for name, B in sources:
    unit = ("%9.3f mT" % (B * 1e3)) if B > 1e-4 else (
        ("%9.3f uT" % (B * 1e6)) if B > 1e-7 else ("%9.3f nT" % (B * 1e9)))
    t_dc, t_ac = (ETA_DC / B) ** 2, (ETA_AC / B) ** 2
    print(f"{name:>28}{B:>13.4e}{unit:>14}{t_dc:>15.3e} s{t_ac:>11.3e} s")

# --- 距離が分解能を指数関数的に決める ---------------------------------------
# 平面状のソースの上ではラプラス方程式が成り立つので、ソース面内の波数 k の
# フーリエ成分は exp(-k d) で減衰して届きます。1/d より細かい空間周波数は、
# 走査のピクセルサイズをどれだけ細かくしても失われています。
print(f"\nラプラスフィルター exp(-k d): 周期 L のソース構造を距離 d で見る")
print(f"{'L / d':>8}{'k d = 2 pi d / L':>19}{'exp(-k d)':>13}"
      f"{'exp(+k d) needed to invert':>29}")
print("-" * 69)
for ratio in [10.0, 4.0, 2.0, 1.0, 0.5, 0.25]:
    kd = 2 * np.pi / ratio
    print(f"{ratio:>8.2f}{kd:>19.5f}{np.exp(-kd):>13.4e}"
          f"{np.exp(kd):>29.4e}")
print("  距離程度より細かい構造は e^{-2 pi} = 1.9e-3 以上の割合で減衰します。")
print("  それを復元するには測定マップに e^{+k d} を掛けることになり、ノイズも")
print("  同じ倍率で増えます。逆問題は指数関数的に悪条件であり、正則化は選択肢")
print("  ではなく必須です。")

# --- 正則化のある場合とない場合の電流再構成 --------------------------------
# 幅 200 nm のストリップが 1 uA を運ぶとシート電流密度は K = I/w = 5 A/m です。
rng = np.random.default_rng(881)
n, dx = 512, 5e-9
x = (np.arange(n) - n // 2) * dx
w, K0 = 200e-9, 1.0e-6 / 200e-9
Ky = np.where(np.abs(x) < w / 2, K0, 0.0)
k = 2 * np.pi * np.fft.fftfreq(n, dx)
for standoff in (20e-9, 50e-9):
    # シート電流の順問題: Bz(k) = (mu_0/2) i sign(k) K(k) e^{-|k| d}
    kernel = 0.5 * MU0 * 1j * np.sign(k) * np.exp(-np.abs(k) * standoff)
    Bz = np.real(np.fft.ifft(kernel * np.fft.fft(Ky)))
    Bm = Bz + ETA_AC * rng.standard_normal(n)        # 1ピクセルあたり1秒の平均
    good = np.abs(kernel) > 0.0
    inv = np.zeros(n, dtype=complex)
    inv[good] = np.fft.fft(Bm)[good] / kernel[good]
    K_naive = np.real(np.fft.ifft(inv))
    K_reg = np.real(np.fft.ifft(np.where(np.abs(k) < 1.0 / standoff, inv, 0.0)))
    inside = np.abs(x) < w / 2
    print(f"\n  距離 {standoff*1e9:.0f} nm: |Bz| の最大 = "
          f"{np.abs(Bz).max()*1e6:.3f} uT、ノイズ {ETA_AC*1e9:.1f} nT/pixel、"
          f"SNR {np.abs(Bz).max()/ETA_AC:.0f}")
    print(f"    真のシート電流密度            : {K0:.4f} A/m")
    print(f"    素朴な反転、ストリップ上の平均: "
          f"{np.mean(K_naive[inside]):.4e} A/m")
    print(f"    k = 1/d で打ち切り、平均      : "
          f"{np.mean(K_reg[inside]):.4f} A/m")
    print(f"    k = 1/d で打ち切り、rms誤差   : "
          f"{np.sqrt(np.mean((K_reg - Ky)**2)):.4f} A/m")
    print(f"    k = 0 の残差（カーネルの零空間）: "
          f"{np.abs(kernel[0]):.1e} - 一様なシート電流は見えません")

距離 20 nm における漏れ磁場と、その検出に要するコスト
  単一NVのDC感度 419 nT/sqrt(Hz)、AC 51.3 nT/sqrt(Hz)

                      source        B (T)             B   t for SNR 1, DC            AC
---------------------------------------------------------------------------------------
        single electron spin   2.3185e-07      0.232 uT      3.266e+00 s  4.896e-02 s
               single proton   3.5265e-10      0.353 nT      1.412e+06 s  2.116e+04 s
     1 nm Co-like film, edge   1.4000e-02     14.000 mT      8.957e-10 s  1.343e-11 s
  monolayer vdW magnet, edge   1.3911e-03      1.391 mT      9.072e-08 s  1.360e-09 s
   uncompensated AFM surface   5.7963e-06      5.796 uT      5.226e-03 s  7.833e-05 s
           1 uA current line   1.0000e-05     10.000 uT      1.756e-03 s  2.632e-05 s
           1 nA current line   1.0000e-08     10.000 nT      1.756e+03 s  2.632e+01 s

ラプラスフィルター exp(-k d): 周期 L のソース構造を距離 d で見る
   L / d   k d = 2 pi d / L    exp(-k d)   exp(+k d) needed to invert
---------------------------------------------------------------------
   10.00            0.62832   5.3349e-01                   1.8745e+00
    4.00            1.57080   2.0788e-01                   4.8105e+00
    2.00            3.14159   4.3214e-02                   2.3141e+01
    1.00            6.28319   1.8674e-03                   5.3549e+02
    0.50           12.56637   3.4873e-06                   2.8675e+05
    0.25           25.13274   1.2162e-11                   8.2226e+10
  距離程度より細かい構造は e^{-2 pi} = 1.9e-3 以上の割合で減衰します。
  それを復元するには測定マップに e^{+k d} を掛けることになり、ノイズも
  同じ倍率で増えます。逆問題は指数関数的に悪条件であり、正則化は選択肢
  ではなく必須です。

  距離 20 nm: |Bz| の最大 = 2.280 uT、ノイズ 51.3 nT/pixel、SNR 44
    真のシート電流密度            : 5.0000 A/m
    素朴な反転、ストリップ上の平均: 4.3526e+01 A/m
    k = 1/d で打ち切り、平均      : 4.2810 A/m
    k = 1/d で打ち切り、rms誤差   : 0.5168 A/m
    k = 0 の残差（カーネルの零空間）: 0.0e+00 - 一様なシート電流は見えません

  距離 50 nm: |Bz| の最大 = 1.411 uT、ノイズ 51.3 nT/pixel、SNR 27
    真のシート電流密度            : 5.0000 A/m
    素朴な反転、ストリップ上の平均: -1.8706e+09 A/m
    k = 1/d で打ち切り、平均      : 3.9271 A/m
    k = 1/d で打ち切り、rms誤差   : 0.6404 A/m
    k = 0 の残差（カーネルの零空間）: 0.0e+00 - 一様なシート電流は見えません

注目すべき点。 最初の表が保存すべきものです。距離20 nmで単一電子スピンは232 nTを与え、信号対雑音比1にはDCで3 s、ACで49 msの平均を要します。単一陽子は0.35 nTで $1.4\times10^6$ s あるいは $2.1\times10^4$ s です。これが、単一電子スピン検出が日常的で単一核スピン検出が可能性の縁にあるというよく知られた事実の定量的な表現です。2つのモーメントは657倍違うので時間は $4\times10^5$ 倍違います。金属膜中の磁壁は14 mTで、ピコ秒の平均で検出されます。この測定はまったく感度律速ではなく、実験上の困難は安定な20 nmの距離を保つことに完全に尽きています。

非補償の反強磁性表面は5.8 $\mu$T でその間に位置します。余裕をもって検出でき、強磁性体より3桁小さく、この差こそが反強磁性体のイメージングがナノメートルの距離でナノテスラの感度をもつ手法を待たねばならなかった理由です。ラプラスの表はなぜ距離が最も効くパラメータなのかを説明します。$L = 10d$ から $L = d/4$ へ移ると信号は10桁失われ、本章でここまで論じたどの感度因子よりはるかに急なペナルティであり、走査NV顕微鏡の設計を浅いイオン注入、探針の鋭さ、表面の清浄さが支配する理由です。それらが $d$ を決め、$d$ は指数で入るのです。

再構成のブロックが悪条件性を具体的にします。1 $\mu$A を運ぶ幅200 nmのストリップは20 nmで2.3 $\mu$T を生み、1ピクセルあたり1秒のAC平均に対して信号対雑音比44です。素朴に反転する — すべてのフーリエ成分をカーネルで割る — と、真値5 A/mに対して20 nmで平均 $4\times10^{1}$ A/m、50 nmで $-2\times10^{9}$ A/m を返します。答えは不正確なだけでなく無意味であり、距離が増えるほど悪化します。$k = 1/d$ で打ち切ると4.28 A/mと3.93 A/mを返し、14%と21%低い。ストリップの鋭い端が、実際に測定されていない高い $k$ に住んでいるからです。このバイアスが正則化の誠実な代価であり、カットオフを述べずに再構成された電流密度を報告することは誤差範囲不明の数値を報告することです。

最後の行は知っておく価値のある細部を記録しています。順方向のカーネルは $k = 0$ で消えるので、空間的に一様なシート電流は面直磁場をまったく作らず $B_z$ マップでは見えません。総電流は面内磁場の線積分か2次元マップから得るしかなく、$B_z$ 画像だけからは得られません。

演習

演習1: ODMRスペクトルから磁場を読む

単一NV中心のODMRスペクトルに2本のディップが2.6980 GHzと3.0460 GHzに見えました。試料を交換したあとの2枚目のスペクトルではディップが2.7150 GHzと3.0350 GHzにあります。

  1. 各スペクトルについて軸方向磁場 $B_\parallel$ をミリテスラで求めてください。
  2. 各スペクトルについて対の中心、したがって見かけのゼロ磁場分裂を求めてください。
  3. 2枚目の中心が動いています。物理的に異なる2つの説明を挙げ、それらを区別する追加測定を1つ述べてください。
  4. このシフトがすべて横磁場によるものだとすると、$B_\perp$ はどれだけの大きさですか。
解答

1. 分裂幅は \(2(\gamma_e/2\pi)B_\parallel\) です。1枚目: \(3.0460 - 2.6980 = 0.3480\) GHz なので \(B_\parallel = 0.3480\times10^9/(2\times28.025\times10^9) = 6.208\) mT。2枚目: \(3.0350 - 2.7150 = 0.3200\) GHz なので \(B_\parallel = 5.7092\) mT。

2. 中心は \((2.6980+3.0460)/2 = 2.8720\) GHz と \((2.7150+3.0350)/2 = 2.8750\) GHz。見かけの \(D\) は3.0 MHz上がっています。

3. 温度が下がったか（\(\mathrm{d}D/\mathrm{d}T \approx -74\) kHz/K なので \(-3.0\ \mathrm{MHz}/(-74\ \mathrm{kHz/K}) = 40.5\) K の冷却で説明できます）、横磁場が現れて中心を2次で押し上げたかです。第3の可能性は局所歪みの変化ですが、これは \(E\) に入り中心をずらすのは \(E/D\) の2次なのではるかに弱い効果です。区別する測定: 意図的に別のバイアス磁場で繰り返すこと（横磁場の寄与は \(B_\perp^2/D\) でスケールし、熱的な寄与はしません）、あるいは第2のNV配向を使うこと（ある配向にとって横向きの磁場は別の配向にとって横向きではない一方、温度は4配向に共通です）。

4. 小さい \(B_\parallel\) での \(\Delta f_{\mathrm{centre}} \approx \tfrac{3}{2}(\gamma_e/2\pi)^2B_\perp^2/D\) から、2つの分母をもう少し丁寧に扱うと \(B_\perp^2 = \Delta f \cdot D / (\tfrac{3}{2}(\gamma_e/2\pi)^2)\)。\(\Delta f = 3.0\) MHz、\(D = 2.872\) GHz として \(B_\perp^2 = 3.0\times10^6\times2.872\times10^9/(1.5\times(2.8025\times10^{10})^2) = 7.31\times10^{-6}\) T\(^2\)、すなわち \(B_\perp = 2.7043\) mT。これは \(B_\parallel\) と同程度なので少なくとも自己整合的です。ただしその大きさの横磁場もどこかから来なければならず、通常は40 Kの冷却のほうがありそうな話です。

import numpy as np
G = 28.025e9
for f1, f2 in [(2.6980e9, 3.0460e9), (2.7150e9, 3.0350e9)]:
    print(f"B_par = {(f2-f1)/(2*G)*1e3:.4f} mT   centre = {(f1+f2)/2/1e9:.4f} GHz")
# B_par = 6.2087 mT   centre = 2.8720 GHz
# B_par = 5.7092 mT   centre = 2.8750 GHz
D = 2.872e9
print(round(float(np.sqrt(3.0e6 * D / (1.5 * G**2)) * 1e3), 4))   # 2.7043  mT
print(round(-3.0e6 / -74e3, 4))                                   # 40.5405 K

演習2: 感度はどこへ消えるのか

単一NV中心の $T_2^\ast = 1.5\ \mu$s、読み出しデッドタイム $t_d = 2\ \mu$s、光学読み出しは1ショットあたり $\alpha_0 = 0.045$、$\alpha_1 = 0.030$ 光子です。Ramsey減衰はガウス型（$p = 2$）とします。

  1. デッドタイムを無視して $\tau = T_2^\ast/2$ における射影ノイズ限界の $\eta$ を計算してください。
  2. 読み出し因子 $\sigma_R$ と、その結果の $\eta$ を計算してください。
  3. デッドタイムを入れてください。$\eta$ は何倍悪化し、最適な $\tau$ はどこへ動きますか。
  4. 同僚が同位体精製で $T_2^\ast$ を $15\ \mu$s にすることを提案しました。$\eta$ の改善を見積もり、2つの効果 — コヒーレンスの伸びとデューティ比の改善 — のどちらが大きいか判定してください。
解答

1. \(\eta = e^{1/4}\sqrt{2/T_2^\ast}/(2\pi\gamma_c)\)、\(\gamma_c = 28.025\times10^9\) Hz/T。\(\sqrt{2/1.5\times10^{-6}} = 1154.7\)、\(2\pi\gamma_c = 1.7609\times10^{11}\) なので \(\eta = 1.2840\times1154.7/1.7609\times10^{11} = 8.421\times10^{-9}\) T/\(\sqrt{\mathrm{Hz}}\)、すなわち 8.42 nT/\(\sqrt{\mathrm{Hz}}\)。

2. \((\alpha_0-\alpha_1)^2 = 2.25\times10^{-4}\)、\(2(\alpha_0+\alpha_1) = 0.15\) なので \(\sigma_R = \sqrt{1+0.15/2.25\times10^{-4}} = \sqrt{667.7} = 25.84\)。したがって \(\eta = 217.6\) nT/\(\sqrt{\mathrm{Hz}}\)。

3. \(\sqrt{\tau+t_d}\,e^{(\tau/T_2^\ast)^2}/\tau\) を数値的に最小化すると \(\tau_{\mathrm{opt}} = 0.970\ \mu\mathrm{s} = 0.647\,T_2^\ast\)、射影限界の \(\eta\) は 15.33 nT/\(\sqrt{\mathrm{Hz}}\) で1.82倍の悪化です。\(\sigma_R\) を含めると \(\eta = 396\) nT/\(\sqrt{\mathrm{Hz}}\)。コヒーレンス時間より長いデッドタイムは単一NVにとって大きな損失であり、物理の問題ではなくデューティ比の問題です。

4. 2つの効果が重なります。\(1/\sqrt{T_2^\ast}\) のスケーリングだけで \(\sqrt{10} = 3.16\)。デューティ比も改善します。\(t_d = 2\ \mu\)s が \(\tau_{\mathrm{opt}}\) に比べて小さくなるからです。\(T_2^\ast = 15\ \mu\)s、\(t_d = 2\ \mu\)s での数値最小値は射影限界で 2.98 nT/\(\sqrt{\mathrm{Hz}}\)、以前の 15.33 に対して合計5.14倍で、うち3.16がコヒーレンス、残る1.63がデューティ比です。コヒーレンスの項が大きいのですが、デューティ比の項もさらに60%に相当し、しばしば忘れられます。

import numpy as np
from scipy.optimize import minimize_scalar
gc = 28.025e9
sR = np.sqrt(1 + 2 * (0.045 + 0.030) / (0.045 - 0.030) ** 2)
def eta(tau, T2s, td):
    return np.sqrt(tau + td) * np.exp((tau / T2s) ** 2) / (2 * np.pi * gc * tau)
print(round(sR, 2))                                            # 25.84
for T2s, td in [(1.5e-6, 0.0), (1.5e-6, 2e-6), (15e-6, 2e-6)]:
    r = minimize_scalar(eta, bracket=(1e-9, T2s, 20 * T2s), args=(T2s, td))
    print(f"T2*={T2s*1e6:5.1f} us td={td*1e6:3.1f} us  tau_opt={r.x*1e6:.3f} us "
          f"eta={r.fun*1e9:7.3f} nT  sigma_R込み {r.fun*sR*1e9:7.1f} nT")
# T2*=  1.5 us td=0.0 us  tau_opt=0.750 us eta=  8.420 nT  sigma_R込み   217.6 nT
# T2*=  1.5 us td=2.0 us  tau_opt=0.970 us eta= 15.328 nT  sigma_R込み   396.1 nT
# T2*= 15.0 us td=2.0 us  tau_opt=8.202 us eta=  2.982 nT  sigma_R込み    77.1 nT

演習3: エコーのフィルター積分を解析的に

第1章は $|\tilde{s}_{\mathrm{echo}}(f,T)|^2 = 4\sin^4(\pi f T/2)/(\pi f)^2$ を与えています。

  1. $\int_0^\infty \sin^4(au)/u^4\,\mathrm{d}u = \pi a^3/3$ を使って $I_1(\alpha=2) \equiv \int_0^\infty u^{-2}\,|\tilde{s}(u/T,T)|^2/T^2\,\mathrm{d}u = \pi^2/6$ を示してください。
  2. Code Example 5 はどの $N$ でも $N^2 I_N(2) = \pi^2/6$ と報告します。$\sum_{k\ \mathrm{odd}}k^{-4} = \pi^4/96$ を使って奇数次高調波の櫛からこれを導いてください。
  3. これは $S(f) = A/f^2$ のノイズについてコヒーレンス時間の $N$ 依存性をどう予言しますか。
  4. 同じ議論が $\alpha = 1$ では有限 $N$ の補正を与えるのに $\alpha = 2$ では与えないのはなぜですか。
解答

1. \(u = fT\) として \(|\tilde{s}|^2/T^2 = 4\sin^4(\pi u/2)/(\pi u)^2\)。したがって \(I_1(2) = \int_0^\infty u^{-2}\cdot 4\sin^4(\pi u/2)/(\pi u)^2\,\mathrm{d}u = (4/\pi^2)\int_0^\infty \sin^4(\pi u/2)/u^4\,\mathrm{d}u\)。\(a = \pi/2\) として引用した積分は \(\pi(\pi/2)^3/3 = \pi^4/24\) なので \(I_1(2) = (4/\pi^2)(\pi^4/24) = \pi^2/6 = 1.644934\)、数値と6桁一致します。

2. 大きな \(N\) ではフィルターは \(f_{\mathrm{res}} = N/2T\) の奇数次高調波、すなわち \(u_k = kN/2\) における櫛です。Parsevalは正の周波数に総重み \(T/2\) を置き、矩形波はそのうち \((8/\pi^2)k^{-2}\) を第 \(k\) 高調波に置きます。したがって \(I_N(2) = \sum_{k\ \mathrm{odd}} \tfrac{1}{2}(8/\pi^2)k^{-2}u_k^{-2} = (4/\pi^2)(4/N^2)\sum_{k\ \mathrm{odd}}k^{-4} = (16/\pi^2N^2)(\pi^4/96) = \pi^2/(6N^2)\)。

3. \(\langle\varphi^2\rangle \propto T^{1+\alpha}I_N = T^3\pi^2/(6N^2)\) なので \(T_2^3 \propto N^2\)、すなわち \(T_2 \propto N^{2/3}\) です。これは \(\alpha = 2\) における \(N^{\alpha/(\alpha+1)}\) であり、第1章の予測どおりです。\(\eta \propto 1/\sqrt{T_2}\) なので感度の利得は \(N^{1/3}\) にすぎないことに注意してください。

4. 櫛近似は各通過帯域を同じ重みのデルタ関数に置き換えます。その誤差は通過帯域の有限の幅と、帯域内での \(S(f)\) の変化から来ます。\(\alpha = 2\) では厳密な結果が \(N = 1\) でも櫛の結果と一致します（1で明示的な積分により示したとおり）。\(\alpha = 1\) では重み \(u^{-1}\) がより平坦なので、櫛が無視する第1通過帯域の低周波側の裾が相対的に大きな割合を占め、その割合は \(N\) が増えて帯域が狭くなるにつれ縮みます。だから \(N I_N(1)\) は0.693から0.848へ漂う一方で \(N^2 I_N(2)\) はまったく漂わないのです。

演習4: リラクソメトリ対デカップリング

表面から10 nm下にNV中心をもつダイヤモンドの上に磁性膜を置きます。膜のスピン波はNVの位置に横磁場ノイズを生み、そのスペクトル密度は1から4 GHzでほぼ平坦な $S_\perp = 2\times10^{-20}$ T$^2$/Hz、100 MHz以下では無視できるとします。

  1. これはどれだけの緩和速度を加え、固有値が5 msなら $T_1$ はいくらになりますか。
  2. XY8系列は同じノイズを検出できますか。10 nsの $\pi$ パルスが許す最高の通過帯域を根拠に答えてください。
  3. 膜を2.87 GHzで駆動し、5 nTの単色横磁場がNVに届くようにします。どの手法が検出し、どれだけの時間を要しますか。
  4. 同僚がスピン波を $T_2$ への影響で検出することを提案しました。フィルター関数の議論を使って、これが失敗する理由を説明してください。
解答

1. \(\Delta(1/T_1) = \tfrac{3}{2}\gamma_e^2 S_\perp = 1.5\times(1.7609\times10^{11})^2\times2\times10^{-20} = 930.2\ \mathrm{s}^{-1}\)。固有の \(200\ \mathrm{s}^{-1}\) を加えて \(1130.2\ \mathrm{s}^{-1}\)、すなわち \(T_1 = 0.885\) ms で固有値の5.65分の1です。容易に測定できます。

2. できません。10 nsの \(\pi\) パルスは通過帯域を \(1/(2\times10\ \mathrm{ns}) = 50\) MHz に制限し、ノイズは100 MHz以下では無視できると指定されています。パルス系列は何もない周波数を見ているのです。2.87 GHzに届くには0.2 ns未満の \(\pi\) パルス、すなわち2.9 GHzのRabi周波数が必要で、これは遷移周波数自身と同程度なので回転座標系の記述の外に出ます。

3. 本章の2手法はどちらもそのままでは検出しません。リラクソメトリは応答しますが、厳密に共鳴上の単色成分は緩和ではなくRabi振動を駆動するので、正しい取り扱いは黄金律の速度ではなく駆動2準位系です。自然な測定は膜の磁場を駆動として使うNVのRabi実験で、\(\Omega/2\pi = (\gamma_e/2\pi)B/\sqrt{2} = 28.025\times10^9\times5\times10^{-9}/1.414 = 99.1\) Hz なので1Rabi周期は10 msかかり、\(T_1\) を超えます。この単色成分はこの方法で見るには弱すぎ、誠実な答えは2.87 GHzの5 nTの単色磁場は単一NVにとって検出可能性の縁にあるということです。（\(1/\sqrt{2}\) はスピン1の行列要素です。）

4. 位相の平均二乗は \(\int S(f)|\tilde{s}(f,T)|^2\mathrm{d}f\) であり、\(|\tilde{s}|^2\) は第1通過帯域を越えると \(1/f^2\) で落ちます。したがって2.87 GHzのノイズは帯域内のノイズに対して \((f_{\mathrm{res}}/f)^2\) だけ抑制され、\(f_{\mathrm{res}} = 3\) MHz なら6桁のパワーであり、\(T_2\) に測定可能な寄与をしません。唯一の注意点は \(\omega_0\) のノイズが \(T_1\) を短くし、\(T_2 \le 2T_1\) なので十分大きな緩和速度は最終的に \(T_2\) を制限することです。それは緩和の帰結であって位相緩和の測定ではなく、そう読むことは \(T_1\) 過程を誤ったスペクトル領域に帰することになります。

import numpy as np
g = 2 * np.pi * 28.025e9
dG = 1.5 * g**2 * 2e-20
print(round(dG, 1), round(1 / (dG + 200) * 1e3, 4))     # 930.2 0.8848  (1/s, ms)
print(round(28.025e9 * 5e-9 / np.sqrt(2), 1))           # 99.1  Hz Rabi

演習5: イメージング実験を設計する

幅2 $\mu$m の金属細線が100 $\mu$A を運んでおり、その電流密度の100 nmの構造を分解することを目標に電流分布をイメージングするよう求められました。

  1. どれだけの距離が必要ですか。その2倍の距離でラプラスフィルターは100 nmの構造に何をしますか。
  2. 細線を一様な電流密度のシートとみなして、必要な距離での磁場を見積もってください。
  3. 第2.3節の単一NVのAC感度51 nT/$\sqrt{\mathrm{Hz}}$ と、200×200ピクセルを前提に、1ピクセルあたり信号対雑音比100のための総取得時間を見積もってください。次に、1ピクセルあたり $\eta = 1$ nT/$\sqrt{\mathrm{Hz}}$ の広視野アンサンブル測定について繰り返してください。
  4. どちらを選びますか。決め手となる1つの物理的な仮定は何ですか。
解答

1. 分解能は距離なので、よくても \(d \approx 100\) nm、余裕をとれば \(d = 50\) nm です。\(d = 200\) nm では100 nmの構造は \(kd = 2\pi\times200/100 = 12.57\) で \(e^{-12.57} = 3.5\times10^{-6}\) だけ減衰します。消えており、平均化では回復しません。指数は物理の中にあってノイズの中にはないからです。

2. シート電流密度 \(K = 100\ \mu\mathrm{A}/2\ \mu\mathrm{m} = 50\) A/m。細線の内側で近距離では磁場は無限シートの値 \(B = \mu_0 K/2 = 1.2566\times10^{-6}\times50/2 = 31.4\ \mu\)T に近づきます。これは面内成分で、空間情報を運ぶ面直成分は端の近くで同程度、中央ではより小さいので、実用的な目安としては 10 \(\mu\)T です。

3. 1ピクセルあたり、10 \(\mu\)T に対する信号対雑音比100には \(\sigma_B = 100\) nT が必要なので \(t = (\eta/\sigma_B)^2\)。単一NV: \((51/100)^2 = 0.26\) s/ピクセル、40 000ピクセルで \(1.04\times10^4\) s、約2.9時間です。走査のオーバーヘッドを加えると実際にはそれが支配します。広視野: \((1/100)^2 = 10^{-4}\) s で、全ピクセルが同時に取得されるので合計 \(10^{-4}\) s。差は9桁で、その全部が並列性です。

4. 広視野測定はどんな尺度でも速いのですが、ここでは誤った選択です。分解能が回折限界の数百ナノメートルであり、距離はNV層の厚みに試料までの隙間を加えたものだからです。要求は100 nmでした。この問いを決めるのは感度ではまったくなく、要求される空間分解能が光学回折限界の上か下かです。上なら広視野を使って距離を受け入れ、下なら走査探針を使って何時間かを受け入れます。第1章1.5節のトレードオフ地図を具体的な要求に適用したものです。

まとめ

要点

1. 点欠陥が良い磁力計になる理由は3つ同時にある

2. DC感度は3つの因子をもち、それらは独立ではない

3. 1つのフィルター関数、2つの用途

4. $T_1$ はGHzの窓であり、他に届くものはない

5. イメージングは距離によって指数関数的に制限される

6. 本章のあらゆる限界は材料の限界である

実務上の含意

第3章はトレードオフ地図の反対の端に移ります。SQUIDは空間分解能4桁と10 K以上で動作する能力を手放し、その代わりに磁場感度5桁から8桁を受け取ります。8桁になるのは大きなピックアップループを単一NVと比べたときだけです。磁束を測り、磁束が $\Phi_0 = h/2e$ の単位で来るからです。そして最良のSQUIDを制限するノイズが超伝導量子ビットを制限するのと同じ欠陥集団から来ることも分かります。それが本コースを貫く材料の糸です。

← シリーズトップ 第3章: SQUID →

免責事項

← 前の章 シリーズ目次に戻る 次の章 →