第3章: イオントラップ

⚛️ Paulトラップ、共有フォノン、Mølmer-Sørensenゲート、そして異常加熱

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

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

基礎数理道場 > 量子ハードウェア入門 > 第3章

第2章は量子ビットを作りました。イオントラップは逆のことをします。量子ビットを見つけるのです。$^{40}$Ca$^+$ イオンは寿命の長い2つの内部状態をもち、その分裂は元素の性質であって、宇宙のどのイオンでも同一で、作製公差もなく素子ごとのばらつきもありません。量子ビット自身については設計するものが何もなく、間違えるものも何もありません。エンジニアリングはすべて、イオンを静止させ、アドレスし、2つを結合させることにあります。

これは困難の在り処を移動させます。その移動をはっきり述べておく価値があります。両方式がいまなお追求されている理由がそこにあるからです。超伝導量子ビットは保持が容易で、同一に作るのが難しい。イオンは同一であることが自明で、保持が難しい。超高真空、安定性が自明でない力学的問題であるラジオ波トラップ、運動の基底状態までのレーザー冷却、そして量子ゲートの位相を定義できるだけ安定なレーザー光学系が必要です。見返りはマイクロ秒ではなく秒のコヒーレンス時間、そして — イオンが振動モードを共有しているために — 隣同士だけでなく任意のペア間の2量子ビットゲートです。

本章は実験と同じ順序で物理を追います。イオンを捕まえ(3.1〜3.2)、冷やし(3.3)、ゲートをかけ(3.4〜3.5)、そして実際に何が律速しているのかに向き合います(3.6)。数値計算は第2章より重くなります。対象がなじみの薄いものだからです。Mathieu安定図を表から引かずにFloquet理論で計算し、イオン結晶の基準モードを導いてジグザグ不安定性を数値的に見つけ、Mølmer-Sørensenゲートを厳密なMagnus展開に照らして検証し — そしてLamb-Dicke近似を破って、それが何を隠していたのかを見ます。

単位と規約。 運動の周波数はMHz、ゲートの離調はkHz、レーザー線幅はMHzで表記します。ハミルトニアン中のエネルギーは周波数($E/h$)であり、シミュレーション内部では $\hbar = 1$、実験室の量が必要なときにSI単位に戻します。量子ビットの順序とゲート記号は量子コンピューティング入門に従います(big-endian、$X$、$Y$、$Z$、$H$、CNOT)。$T_1$、$T_2$、$T_2^\ast$ は第1章の定義をそのまま使います。3.3節以降で使う調和振動子と量子化された運動の道具は量子力学入門で展開されています。

学習目標

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


3.1 Paulトラップの物理

Earnshawの障害と、その回り道

静電場だけで荷電粒子を捕捉することはできません。Laplace方程式が禁じています。$\nabla^2\phi = 0$ はポテンシャルが自由空間で極小をもたないことを意味するので、2方向で閉じ込める静的な配置はどれも第3の方向で追い出します。これがEarnshawの定理であり、素直なやり方を完全に排除します。

Wolfgang Paulの解決策は、鞍点を回転させることです。四重極電極形状にラジオ波電圧をかけ、ポテンシャルを

$$ \phi(x, y, t) = \frac{U + V\cos\Omega t}{2 r_0^2}\left(x^2 - y^2\right) $$

とします。任意の瞬間には一方の軸で閉じ込め、他方で追い出しています。しかし追い出しの向きはRF半周期ごとに反転し、その反転がイオンの応答より速ければ、正味の効果は両方向での閉じ込めになります。イオンは $\Omega$ での小さく速い揺れ — マイクロモーション — を、実効的な時間平均ポテンシャル中の遅く大きな secular 振動に重ねて実行します。それが実際に働くのか、どのパラメータで働くのかは、本物の力学的問題であって本物の答えがあります。

Mathieu方程式

径方向運動についてのNewton方程式は、$\xi = \Omega t/2$ とおくと標準的なMathieu方程式になります:

$$ \frac{d^2u}{d\xi^2} + \left(a - 2q\cos 2\xi\right)u = 0 $$

無次元のトラップパラメータは

$$ q = \frac{2eV}{m r_0^2 \Omega^2}, \qquad a = -\frac{4eU}{m r_0^2 \Omega^2} $$

です(符号の規約は文献によって異なります。本章では上記を通して使います)。上に書いたのは、RFポテンシャルが $x^2 - y^2$ に比例する線形トラップの $x$ 方向のパラメータです。$y$ 方向の式は $q \to -q$、$a \to -a$ とした同じ式になります。安定性は $q$ には $q^2$ を通してしか依存しないので以下の図は $q$ について対称ですが、$a$ については対称ではありません。一方の径方向に有利なDC電圧は他方に不利であり、だから $|a|$ は常に小さく保たれます。

これは周期的な係数をもつ線形方程式なので、Floquet理論が使えます。係数の1周期($\cos 2\xi$ の周期は $\pi$ なので $\xi$ が0から $\pi$)にわたって独立な2つの解を積分し、$2\times2$ のモノドロミー行列 $M$ を作ります。減衰項がないので $\det M = 1$ が厳密に成り立ち、固有値は $\lambda^{\pm1}$ で $\lambda\lambda^{-1} = 1$ です。2つの場合があります:

捕捉される場合には $\cos\pi\beta = \mathrm{tr}\,M/2$ で特性指数 $\beta$ が定義され、secular周波数は

$$ \omega_\mathrm{sec} = \frac{\beta\Omega}{2} $$

です。$q$ が小さいとき擬ポテンシャル近似は $\beta \approx \sqrt{a + q^2/2}$、$a = 0$ では単に $\beta \approx q/\sqrt{2}$ を与えます。

Code Example 1: 安定図を計算する

Paulトラップの安定図は普通、引いてくる図として提示されます。実際には20行の計算であり、自分でやってみることで物理が — そして有名な数値0.908の意味が — 具体的になります。

"""第3章 Code Example 1: Mathieu安定図を、表を引かずにFloquet理論で計算する。

標準形(xi = Omega t / 2 とする):
    d^2 u / d xi^2 + (a - 2 q cos 2 xi) u = 0
係数の xi についての周期は pi なので、1周期分の写像で足りる。"""
import numpy as np
from scipy.integrate import solve_ivp


def monodromy(a, q):
    """Mathieu方程式の1周期 xi: 0 -> pi にわたる転送行列。"""
    def rhs(xi, y):
        u1, v1, u2, v2 = y
        k = a - 2.0 * q * np.cos(2.0 * xi)
        return [v1, -k * u1, v2, -k * u2]
    sol = solve_ivp(rhs, [0.0, np.pi], [1.0, 0.0, 0.0, 1.0],
                    rtol=1e-11, atol=1e-13, dense_output=False)
    u1, v1, u2, v2 = sol.y[:, -1]
    return np.array([[u1, u2], [v1, v2]])


def stable(a, q):
    """解が有界であることと |モノドロミー行列のトレース| < 2 は同値。"""
    return abs(np.trace(monodromy(a, q))) < 2.0


def beta_exact(a, q):
    """cos(pi beta) = tr(M)/2 から特性指数を求める。"""
    t = np.trace(monodromy(a, q)) / 2.0
    return np.arccos(np.clip(t, -1.0, 1.0)) / np.pi


# --- 検査:行列式は1でなければならない(Liouville) --------------------
M = monodromy(0.0, 0.5)
print(f"(a, q) = (0, 0.5) でのモノドロミー行列の行列式: "
      f"{np.linalg.det(M):.12f}   (Liouvilleより厳密に1)")
print()

# --- a = 0 における最初の安定境界 --------------------------------------
lo, hi = 0.5, 1.5
for _ in range(60):
    mid = 0.5 * (lo + hi)
    lo, hi = (mid, hi) if stable(0.0, mid) else (lo, mid)
print(f"a = 0: q < {0.5 * (lo + hi):.6f} でイオンは捕捉される")
print()

# --- 安定領域 ----------------------------------------------------------
print(f"{'q':>7}{'a_min':>11}{'a_max':>11}{'beta at a=0':>14}"
      f"{'sqrt(q^2/2)':>14}{'error':>10}")
for q in [0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 0.9]:
    def scan(sign):
        """a = 0 から安定性が失われるまで進み、そこから二分法にかける。

        [0, +-1.2] で素朴に二分法を使うと失敗する。最初の不安定帯の上には
        第2の安定帯があり、二分法がそこへ入り込んで、いま居る領域の縁では
        ない境界を報告してしまう。
        """
        step = 1.0e-2
        lo2 = 0.0
        while stable(lo2 + sign * step, q) and abs(lo2) < 1.5:
            lo2 += sign * step
        hi2 = lo2 + sign * step
        for _ in range(40):
            mid = 0.5 * (lo2 + hi2)
            lo2, hi2 = (mid, hi2) if stable(mid, q) else (lo2, mid)
        return 0.5 * (lo2 + hi2)
    amin, amax = scan(-1.0), scan(+1.0)
    b = beta_exact(0.0, q)
    approx = np.sqrt(q ** 2 / 2.0)
    print(f"{q:>7.2f}{amin:>11.5f}{amax:>11.5f}{b:>14.6f}"
          f"{approx:>14.6f}{abs(b - approx) / b * 100:>9.3f}%")
print()

# --- 同じ数値を実験室の単位で ------------------------------------------
u = 1.66053906660e-27
e = 1.602176634e-19
mass = 40.078 * u                # Ca-40
r0 = 0.5e-3                      # m、イオン-電極間距離
frf = 20.0e6                     # Hz、RF駆動周波数
Om = 2.0 * np.pi * frf
print(f"Ca+(m = {mass / u:.1f} u)のリニアPaulトラップ、r0 = "
      f"{r0 * 1e3:.1f} mm、Omega/2pi = {frf / 1e6:.0f} MHz。")
print(f"{'V_RF (V)':>10}{'q':>8}{'beta (exact)':>14}"
      f"{'f_sec (MHz)':>13}{'q/sqrt(2) est':>15}{'depth (eV)':>12}")
for VRF in [100.0, 245.0, 400.0, 600.0, 742.0]:
    q = 2.0 * e * VRF / (mass * r0 ** 2 * Om ** 2)
    if q >= 0.908:
        print(f"{VRF:>10.0f}{q:>8.4f}{'    不安定':>14}")
        continue
    b = beta_exact(0.0, q)
    wsec = b * Om / 2.0
    depth = 0.5 * mass * wsec ** 2 * r0 ** 2 / e
    print(f"{VRF:>10.0f}{q:>8.4f}{b:>14.6f}{wsec / 2 / np.pi / 1e6:>13.4f}"
          f"{q / np.sqrt(2.0):>15.6f}{depth:>12.2f}")
print()
q = 0.3
b = beta_exact(0.0, q)
wsec = b * Om / 2.0
depth = 0.5 * mass * wsec ** 2 * r0 ** 2
kB = 1.380649e-23
print(f"q = {q:.1f} では f_sec = {wsec / 2 / np.pi / 1e6:.4f} MHz、トラップ深さ "
      f"{depth / e:.2f} eV = {depth / kB:.2e} K。")
(a, q) = (0, 0.5) でのモノドロミー行列の行列式: 1.000000000000   (Liouvilleより厳密に1)

a = 0: q < 0.908046 でイオンは捕捉される

      q      a_min      a_max   beta at a=0   sqrt(q^2/2)     error
   0.05   -0.00125    0.94969      0.035373      0.035355    0.049%
   0.10   -0.00499    0.89877      0.070850      0.070711    0.196%
   0.20   -0.01991    0.79512      0.142551      0.141421    0.793%
   0.30   -0.04457    0.68917      0.216059      0.212132    1.818%
   0.50   -0.12177    0.47065      0.373744      0.353553    5.402%
   0.70   -0.23317    0.24391      0.563066      0.494975   12.093%
   0.90   -0.37456    0.00958      0.915911      0.636396   30.518%

Ca+(m = 40.1 u)のリニアPaulトラップ、r0 = 0.5 mm、Omega/2pi = 20 MHz。
  V_RF (V)       q  beta (exact)  f_sec (MHz)  q/sqrt(2) est  depth (eV)
       100  0.1220      0.086493       0.8649       0.086240        1.53
       245  0.2988      0.215168       2.1517       0.211289        9.49
       400  0.4878      0.363556       3.6356       0.344961       27.09
       600  0.7318      0.599230       5.9923       0.517442       73.60
       742  0.9050      0.947949       9.4795       0.639903      184.20

q = 0.3 では f_sec = 2.1606 MHz、トラップ深さ 9.57 eV = 1.11e+05 K。

着目点。 モノドロミー行列の行列式は12桁で1です。これは数値的な偶然ではありません。ハミルトニアン系についてのLiouvilleの定理であり、積分が正しいことを確かめる最も安価な検査です。ここにドリフトがあれば許容誤差が緩すぎます。

$a = 0$ での境界は $q = 0.908046$ と出て、教科書の0.90800に対応します。この数値はもう借り物ではなく自分のものです。これはPaulトラップをどこまで強く駆動できるかの実用的な上限です。それを越えるとイオンは弱く閉じ込められるのではなく、指数関数的に追い出されます

安定領域は $q$ が大きくなると狭まります。$q = 0.05$ では直流パラメータ $a$ は $[-0.001, 0.950]$ の範囲を取れますが、$q = 0.9$ では $[-0.375, 0.010]$ だけです。したがってトラップは控えめな $q$、典型的には0.1から0.3で運転され、軸方向の閉じ込めとマイクロモーション補償を担う直流電極のための余地を残します。

$\beta$ と $q/\sqrt{2}$ の比較は、擬ポテンシャル描像をどこまで信用できるかを述べています。$q = 0.05$ で0.05%、$q = 0.3$ で1.8%、$q = 0.9$ で30%の誤差です。動作点 $q = 0.3$ では近似は有効数字2桁まで良く、だからこそ本章の残りではRF駆動を忘れてイオンが静的な調和井戸にいるものとして扱えるのです。

実験室の表は電圧に換算します。20 MHzで駆動する0.5 mmのトラップは $q = 0.30$ に245 VのRF振幅を必要とし、2.16 MHzのsecular周波数と9.6 eVのトラップ深さ — すなわち $1.1 \times 10^5$ K — を与えます。$10^5$ Kの井戸の中のミリケルビンのイオンはどこにも行きません。保持時間を制限するのは深さではなく化学です。背景ガスとの衝突、そしてCa$^+$をCaH$^+$に変えてしまう反応です。だから超高真空であり、だから1個のイオンを数日保持できるのです。

マイクロモーション

擬ポテンシャルは近似であり、そこで捨てられるものが物理的に重要です。厳密なFloquet解は $u(\xi) = e^{i\beta\xi}P(\xi)$($P$ は周期関数)なので、secular周波数 $\beta\Omega/2$ だけでなく $\beta\Omega/2 \pm n\Omega$ のサイドバンドも含みます。$n = \pm1$ の成分がマイクロモーションで、secular振幅に対する相対振幅は $q/2$ です。

イオンがRFのゼロ点にちょうど座っている場合、これは小さくてほぼ無害な効果です。問題は過剰マイクロモーションです。迷い直流電場 — 電極上の帯電した絶縁体パッチから、真空容器内のどこかの迷い電子から — がイオンをゼロ点から押し出し、RF電場が消えない領域に入れてしまいます。すると振幅は温度ではなく変位が決めるかたちで、$\Omega$ で駆動されます。冷却は助けになりません。

帰結は分光的なものです。レーザー方向に沿って $\Omega$ で振動するイオンは位相変調された場を見て、変調指数は $k \cdot x_\mathrm{micro}$ です。スペクトルは $\pm\Omega$ にマイクロモーションサイドバンドを獲得し、キャリアとサイドバンドの比はBessel関数で与えられ、キャリア自身は弱まります。それらのサイドバンドをゼロにすることが迷い電場の測定と補償の方法であり — そして繰り返さなければなりません。迷い電場はドリフトするからです。

Code Example 2: マイクロモーションを厳密に

"""第3章 Code Example 2: マイクロモーション。Mathieu方程式の厳密なFloquet解を、
本章の残りで用いる擬ポテンシャル(secular)近似と比較する。"""
import numpy as np
from scipy.integrate import solve_ivp
from scipy.special import jv

TWOPI = 2.0 * np.pi


def monodromy(q, a=0.0):
    def rhs(xi, y):
        k = a - 2.0 * q * np.cos(2.0 * xi)
        return [y[1], -k * y[0], y[3], -k * y[2]]
    s = solve_ivp(rhs, [0.0, np.pi], [1.0, 0.0, 0.0, 1.0],
                  rtol=1e-12, atol=1e-14)
    u1, v1, u2, v2 = s.y[:, -1]
    return np.array([[u1, u2], [v1, v2]])


def floquet_harmonics(q, a=0.0, nmax=4, npts=4001):
    """Floquet解 u = exp(i beta xi) P(xi) のFourier成分を求める。

    P は周期 pi をもつので u = sum_n C_n exp(i(beta + 2n) xi) と書ける。
    実時間では成分 n の周波数は beta*Omega/2 + n*Omega である。n = 0 が
    secular運動、n = +-1 が駆動周波数のマイクロモーションにあたる。
    """
    M = monodromy(q, a)
    w, V = np.linalg.eig(M)
    k = int(np.argmax(np.imag(w)))
    beta = np.angle(w[k]) / np.pi
    y0 = V[:, k]

    def rhs(xi, y):
        kk = a - 2.0 * q * np.cos(2.0 * xi)
        return [y[1], -kk * y[0]]
    xi = np.linspace(0.0, np.pi, npts)
    s = solve_ivp(rhs, [0.0, np.pi], y0, t_eval=xi, rtol=1e-12, atol=1e-14)
    P = s.y[0] * np.exp(-1j * beta * xi)
    C = {n: np.trapezoid(P * np.exp(-2j * n * xi), xi) / np.pi
         for n in range(-nmax, nmax + 1)}
    return beta, C


print(f"{'q':>7}{'beta':>10}{'|C1|/|C0|':>12}{'|C-1|/|C0|':>12}"
      f"{'micro/secular':>15}{'q/2':>8}{'|C2|/|C0|':>12}")
for q in [0.05, 0.1, 0.2, 0.3, 0.5]:
    beta, C = floquet_harmonics(q)
    c0 = abs(C[0])
    ratio = (abs(C[1]) + abs(C[-1])) / c0
    print(f"{q:>7.2f}{beta:>10.6f}{abs(C[1]) / c0:>12.6f}"
          f"{abs(C[-1]) / c0:>12.6f}{ratio:>15.6f}{q / 2.0:>8.4f}"
          f"{abs(C[2]) / c0:>12.3e}")
print()

# --- 実験室の数値 ------------------------------------------------------
u_amu = 1.66053906660e-27
e = 1.602176634e-19
kB = 1.380649e-23
hbar = 1.054571817e-34
mass = 40.078 * u_amu
r0 = 0.5e-3
frf = 20.0e6
Om = TWOPI * frf
q = 0.3
beta, _ = floquet_harmonics(q)
wsec = beta * Om / 2.0
lam = 729e-9
kvec = TWOPI / lam
print(f"Ca+ at q = {q}: f_sec = {wsec / TWOPI / 1e6:.4f} MHz, "
      f"f_RF = {frf / 1e6:.0f} MHz")
x0 = np.sqrt(hbar / (2.0 * mass * wsec))
print(f"  基底状態の広がり x0 = sqrt(hbar / 2 m omega) = {x0 * 1e9:.3f} nm")
print(f"  基底状態のイオンが本来もつマイクロモーション: "
      f"{q / 2.0 * x0 * 1e9:.3f} nm、"
      f"変調指数 k x = {kvec * q / 2 * x0:.4f}")
print()
print(f"{'E_stray (V/m)':>15}{'offset (nm)':>13}{'micromotion (nm)':>18}"
      f"{'v (m/s)':>10}{'mod. index k x':>16}{'sideband/carrier':>18}")
for Estray in [0.1, 1.0, 10.0, 100.0]:
    d = e * Estray / (mass * wsec ** 2)
    amp = q / 2.0 * d
    v = amp * Om
    mi = kvec * amp
    sb = (jv(1, mi) / jv(0, mi)) ** 2
    print(f"{Estray:>15.1f}{d * 1e9:>13.2f}{amp * 1e9:>18.3f}{v:>10.3f}"
          f"{mi:>16.4f}{sb:>18.3e}")
print()
      q      beta   |C1|/|C0|  |C-1|/|C0|  micro/secular     q/2   |C2|/|C0|
   0.05  0.035373    0.012070    0.012955       0.025024  0.0250   3.706e-05
   0.10  0.070850    0.023322    0.026875       0.050197  0.0500   1.407e-04
   0.20  0.142551    0.043590    0.058014       0.101604  0.1000   5.081e-04
   0.30  0.216059    0.061151    0.094454       0.155605  0.1500   1.032e-03
   0.50  0.373744    0.088943    0.190427       0.279369  0.2500   2.325e-03

Ca+ at q = 0.3: f_sec = 2.1606 MHz, f_RF = 20 MHz
  基底状態の広がり x0 = sqrt(hbar / 2 m omega) = 7.640 nm
  基底状態のイオンが本来もつマイクロモーション: 1.146 nm、変調指数 k x = 0.0099

  E_stray (V/m)  offset (nm)  micromotion (nm)   v (m/s)  mod. index k x  sideband/carrier
            0.1         1.31             0.196     0.025          0.0017         7.131e-07
            1.0        13.06             1.959     0.246          0.0169         7.131e-05
           10.0       130.63            19.595     2.462          0.1689         7.182e-03
          100.0      1306.32           195.948    24.624          1.6889         2.036e+00

着目点。 厳密なFloquet解を高調波に分解すると、フィッティングなしで標準的な描像が確認されます。$n = \pm1$ 成分の $n = 0$ に対する合計重みは $q = 0.05$ で0.02502(予測 $q/2 = 0.025$)、$q = 0.3$ で0.1556(予測0.15)— 3次の補正が4%です。第2高調波は $q^2$ の次数でしか現れず、そうあるべきです。

$|C_{+1}|$ と $|C_{-1}|$ の非対称性($q = 0.3$ で0.061対0.094)は実在するもので、1次の解では見えません。上側と下側のマイクロモーションサイドバンドは同じ強さではないのです。

実験室の数値がスケールを与えます。2.16 MHzの井戸中の基底状態のCa$^+$イオンは7.6 nmに広がり、本来のマイクロモーションは1.1 nm — 729 nm遷移に対して変調指数0.01で、無視できます。ここに迷い電場を加えます。1 V/mではイオンは13 nm変位し、振幅2.0 nmで駆動され、変調指数0.017、サイドバンド/キャリア比 $7\times10^{-5}$ — 検出できて、小さい値です。10 V/mでは比は0.7%、100 V/mでは変調指数が1を超え、マイクロモーションサイドバンドがキャリアより強くなります。

数十V/mの迷い電場は、誘電体表面が露出したトラップではまったくありふれています。したがってマイクロモーション補償は任意の洗練ではありません。物理的な起源に注意してください。電極近傍の絶縁パッチに蓄積する電荷であり、それらのパッチが帯電・放電するにつれてドリフトします。これは第2章の2準位欠陥と同じ種類の問題 — 遅いダイナミクスをもつ制御されていない表面状態 — であり、3.6節で別の名前をまとって再登場します。


3.2 量子ビットを担うもの

トラップはイオンを保持し、量子ビットはその内部状態に住みます。使われている系統は2つあり、その選択は実験全体に波及します。

超微細量子ビットは電子基底状態の2つの副準位を使い、電子と核スピンの相互作用による分裂 — $^{171}$Yb$^+$ で 12.6 GHz、$^9$Be$^+$ で 1.25 GHz — を利用します。両状態が基底多重項にあるので自発放出はまったくなく、自然寿命は本質的に無限です。コヒーレンスは代わりに磁場ノイズが制限します。2つの準位の磁気モーメントが異なるからです。対処法はクロック遷移です。ゼロ磁場または「マジック」磁場でエネルギー差が磁場について停留する、$\partial\nu/\partial B = 0$ となる状態の対です。そこでは磁場に対する1次感度が消え、数秒 — 専用の実験ではさらに桁違いに長い — のコヒーレンス時間が従います。代価は、分裂がマイクロ波周波数であることです。マイクロ波を直接かけても空間分解能はなく運動との結合もないので、2量子ビットゲートには誘導Raman遷移を駆動する2本のレーザーが必要で、実効的な波数ベクトルは2本のになります。

光学量子ビットは基底状態と禁制遷移でつながる準安定励起状態を使います。$^{40}$Ca$^+$ の729 nmにある $S_{1/2} \to D_{5/2}$ 四重極子線で、自然寿命は1秒程度です。狭線幅レーザー1本で直接駆動でき、単一イオンのアドレッシングが素直になり、1本のビームで運動サイドバンドに直接アクセスできます。代価は、励起状態が実際に崩壊するので $T_1$ が準安定寿命で打ち止めになること、そしてレーザー自身の周波数安定性が量子ビットのコヒーレンスに入ってくることです。1 Hzの線幅のレーザーは装置であって、買ってくるものではありません。

どちらが支配的ということはありません。本章の残りにとって重要なのは、どちらもレーザーの波数ベクトルを通じてイオンの運動に結合するということ、そしてその結合が2量子ビットゲートを組み立てる資源であるということです。

イオン結晶

同じトラップ中の複数のイオンは互いに反発し、弱く閉じ込められた軸に沿って直線結晶に落ち着きます。平衡位置まわりの微小振動は全イオンが共有する基準モードであり — その共有こそが任意ペア間の2量子ビットゲートを可能にするものです。

軸方向の調和井戸では、平衡位置は

$$ V = \sum_i \frac{u_i^2}{2} + \sum_{i<j} \frac{1}{|u_i - u_j|} $$

を $\ell = \left(e^2/4\pi\epsilon_0 m\omega_z^2\right)^{1/3}$ の単位で最小化することで得られ、モード周波数はHessianの固有値の平方根です。概念的に重要なモードが2つあります。全イオンが一緒に動くちょうど $\omega_z$ の重心モードと、2イオンで $\sqrt{3}\omega_z$ の伸縮モードです。

径方向モードは違うふるまいをし、実際にゲートで使われるのはそちらです。Coulomb反発は径方向の運動を軟化させます — 横に押されたイオンは隣のイオンによってさらに横に押されます — ので、径方向スペクトルは単一イオンの径方向周波数 $\omega_r$ よりにあり、狭い帯に押し込められます。軟化が閉じ込めを超えると最低の径方向モードが虚数になり、鎖はジグザグに座屈します。これは純粋にCoulomb行列についての幾何学的な言明であり、Code Example 3 が境界を数値的に見つけます。

Code Example 3: イオン結晶とそのモード

"""第3章 Code Example 3: イオン結晶。平衡位置、軸方向と径方向の基準モード、
ジグザグ不安定性、そして長い鎖のゲート速度を縛るモードの混雑。"""
import numpy as np

u_amu = 1.66053906660e-27
e = 1.602176634e-19
eps0 = 8.8541878128e-12
TWOPI = 2.0 * np.pi
KE = e ** 2 / (4.0 * np.pi * eps0)


def equilibrium(N, iters=300):
    """調和的な軸方向ポテンシャル中のN個のイオンの無次元平衡位置。

    長さの単位は l = (e^2 / 4 pi eps0 m omega_z^2)^(1/3)、ポテンシャルは
    V = sum_i u_i^2/2 + sum_{i<j} 1/|u_i - u_j|。Newton法で解く。
    """
    if N == 1:
        return np.zeros(1)
    u = np.linspace(-1.0, 1.0, N) * (0.5 * N ** 0.6)
    for _ in range(iters):
        d = u[:, None] - u[None, :]
        np.fill_diagonal(d, np.inf)
        inv3 = 1.0 / np.abs(d) ** 3
        grad = u - np.sum(np.sign(d) / d ** 2, axis=1)
        H = -2.0 * inv3
        np.fill_diagonal(H, 1.0 + 2.0 * np.sum(inv3, axis=1))
        u = u - np.linalg.solve(H, grad)
    return np.sort(u)


def modes(N, ratio=None):
    """軸方向モード周波数(omega_z単位)と、径方向/軸方向の周波数比が
    与えられた場合は径方向モードも返す。戻り値は (u, w_ax, V_ax, w_tr)。"""
    u = equilibrium(N)
    d = u[:, None] - u[None, :]
    np.fill_diagonal(d, np.inf)
    inv3 = 1.0 / np.abs(d) ** 3
    A = -2.0 * inv3
    np.fill_diagonal(A, 1.0 + 2.0 * np.sum(inv3, axis=1))
    w2, V = np.linalg.eigh(A)
    w_ax = np.sqrt(np.maximum(w2, 0.0))
    w_tr = None
    if ratio is not None:
        B = inv3.copy()
        np.fill_diagonal(B, ratio ** 2 - np.sum(inv3, axis=1))
        t2 = np.linalg.eigvalsh(B)
        w_tr = np.sign(t2) * np.sqrt(np.abs(t2))    # 負なら不安定
    return u, w_ax, V, w_tr


print(f"{'N':>4}{'length (l)':>12}{'min spacing (l)':>17}"
      f"   mode frequencies (omega_z)")
for N in [1, 2, 3, 5, 8, 10]:
    u, w, V, _ = modes(N)
    sp = np.diff(u).min() if N > 1 else np.nan
    print(f"{N:>4}{(u[-1] - u[0]):>12.4f}{sp:>17.4f}   "
          + "  ".join(f"{x:.4f}" for x in w[:6]))
print()
u2, w2, V2, _ = modes(2)
print(f"  イオン間隔      {u2[1] - u2[0]:.9f} l   "
      f"(2 x 4^(-1/3) = {2.0 * 4.0 ** (-1.0 / 3.0):.9f})")
print(f"  モード0 (COM)   {w2[0]:.9f} omega_z   (厳密に1)")
print(f"  モード1 (伸縮)  {w2[1]:.9f} omega_z   "
      f"(sqrt(3) = {np.sqrt(3.0):.9f})")
print(f"  COMモードベクトル   {np.round(np.abs(V2[:, 0]), 6)}  "
      f"(一様、1/sqrt(2) = {1 / np.sqrt(2):.6f})")
print(f"  伸縮モードベクトル  {np.round(V2[:, 1], 6)}")
print()

# --- 実験室の単位 ------------------------------------------------------
mass = 40.078 * u_amu
fz = 1.0e6
wz = TWOPI * fz
ell = (KE / (mass * wz ** 2)) ** (1.0 / 3.0)
print(f"f_z = {fz / 1e6:.1f} MHz の Ca+: 長さの単位 l = "
      f"{ell * 1e6:.3f} um。イオンは数ミクロン間隔で並ぶ。")
print(f"{'N':>4}{'chain length (um)':>19}{'min spacing (um)':>18}"
      f"{'f_COM (MHz)':>13}{'f_ax,max (MHz)':>16}")
for N in [2, 6, 10, 20, 30]:
    u, w, V, _ = modes(N)
    print(f"{N:>4}{(u[-1] - u[0]) * ell * 1e6:>19.2f}"
          f"{np.diff(u).min() * ell * 1e6:>18.3f}"
          f"{w[0] * fz / 1e6:>13.4f}{w[-1] * fz / 1e6:>16.4f}")
print()

# --- 径方向モードとジグザグ不安定性 ------------------------------------
print(f"{'omega_r/omega_z':>16}{'largest linear N':>18}"
      f"{'0.73 N^0.86 at that N':>24}")
for ratio in [3.0, 5.0, 8.0, 12.0, 20.0]:
    Nmax = 1
    for N in range(2, 200):
        _, _, _, wt = modes(N, ratio)
        if wt.min() <= 0.0:
            break
        Nmax = N
    print(f"{ratio:>16.1f}{Nmax:>18}{0.73 * Nmax ** 0.86:>24.2f}")
print()

print(f"{'N':>4}{'bandwidth (kHz)':>17}{'min spacing (kHz)':>19}"
      f"{'implied gate time (us)':>24}")
fr = 5.0e6
rows = []
for N in [2, 6, 10, 14, 20, 30]:
    _, _, _, wt = modes(N, fr / fz)
    if wt.min() <= 0:
        # 径方向の閉じ込めを強めて直線配置の安定性を回復させる
        ratio = 0.9 * N ** 0.86 + 1.0
        _, _, _, wt = modes(N, ratio)
        tag = f"  (omega_r/omega_z = {ratio:.1f} が必要)"
    else:
        tag = ""
    f = np.sort(wt) * fz
    gap = np.diff(f).min()
    rows.append((N, gap / fz))
    print(f"{N:>4}{(f[-1] - f[0]) / 1e3:>17.2f}{gap / 1e3:>19.3f}"
          f"{1e6 / gap:>24.2f}{tag}")
print()
sp = []
for N in [8, 16, 32, 64]:
    _, _, _, wt = modes(N, 0.9 * N ** 0.86 + 1.0)
    sp.append((N, np.diff(np.sort(wt)).min()))
c = np.polyfit(np.log([s[0] for s in sp]), np.log([s[1] for s in sp]), 1)
print(f"安定限界でのフィッティング: 最小モード間隔 / omega_z ~ "
      f"N^{c[0]:.2f}")
print()
u, w, V, _ = modes(5)
print("N = 5 の軸方向モード参加率(行=イオン、列=モード)")
print("  周波数 (omega_z): " + "  ".join(f"{x:.4f}" for x in w))
for i in range(5):
    print(f"   イオン{i}: " + "  ".join(f"{V[i, m]:+.4f}" for m in range(5)))
print()
   N  length (l)  min spacing (l)   mode frequencies (omega_z)
   1      0.0000              nan   1.0000
   2      1.2599           1.2599   1.0000  1.7321
   3      2.1544           1.0772   1.0000  1.7321  2.4083
   5      3.4858           0.8221   1.0000  1.7321  2.4120  3.0549  3.6708
   8      4.9516           0.6360   1.0000  1.7321  2.4153  3.0632  3.6847  4.2859
  10      5.7417           0.5642   1.0000  1.7321  2.4168  3.0672  3.6914  4.2955

  イオン間隔      1.259921050 l   (2 x 4^(-1/3) = 1.259921050)
  モード0 (COM)   1.000000000 omega_z   (厳密に1)
  モード1 (伸縮)  1.732050808 omega_z   (sqrt(3) = 1.732050808)
  COMモードベクトル   [0.707107 0.707107]  (一様、1/sqrt(2) = 0.707107)
  伸縮モードベクトル  [-0.707107  0.707107]

f_z = 1.0 MHz の Ca+: 長さの単位 l = 4.445 um。イオンは数ミクロン間隔で並ぶ。
   N  chain length (um)  min spacing (um)  f_COM (MHz)  f_ax,max (MHz)
   2               5.60             5.600       1.0000          1.7321
   6              17.89             3.288       1.0000          4.2738
  10              25.52             2.508       1.0000          6.5758
  20              38.47             1.708       1.0000         11.9280
  30              47.75             1.355       1.0000         16.9902

 omega_r/omega_z  largest linear N   0.73 N^0.86 at that N
             3.0                 6                    3.41
             5.0                11                    5.74
             8.0                18                    8.77
            12.0                30                   13.60
            20.0                53                   22.19

   N  bandwidth (kHz)  min spacing (kHz)  implied gate time (us)
   2           101.02            101.021                    9.90
   6           954.37            101.021                    9.90
  10          3030.29            101.021                    9.90
  14          2203.44             51.642                   19.36  (omega_r/omega_z = 9.7 が必要)
  20          3134.88             39.019                   25.63  (omega_r/omega_z = 12.8 が必要)
  30          4656.95             28.158                   35.51  (omega_r/omega_z = 17.8 が必要)

安定限界でのフィッティング: 最小モード間隔 / omega_z ~ N^-0.80

N = 5 の軸方向モード参加率(行=イオン、列=モード)
  周波数 (omega_z): 1.0000  1.7321  2.4120  3.0549  3.6708
   イオン0: +0.4472  +0.6395  +0.5377  -0.3017  +0.1045
   イオン1: +0.4472  +0.3017  -0.2805  +0.6395  -0.4704
   イオン2: +0.4472  -0.0000  -0.5143  +0.0000  +0.7318
   イオン3: +0.4472  -0.3017  -0.2805  -0.6395  -0.4704
   イオン4: +0.4472  -0.6395  +0.5377  +0.3017  +0.1045

着目点。 2イオンの場合は厳密で、すべてが合います。間隔 $2\times4^{-1/3} = 1.259921$、モードはちょうど1と $\sqrt{3} = 1.7320508$、重心モードベクトルは $1/\sqrt{2}$ で一様、伸縮モードは反対称です。閉じた形の結果と9桁一致することが、閉じた形が存在しない $N = 30$ でコードを信用する免許になります。

実験室の単位では、$f_z = 1$ MHzのCa$^+$ で $\ell = 4.4\ \mu$m なのでイオンは数ミクロン間隔で並びます — 光学的に十分分解でき、それが単一イオンのアドレッシングと単一イオンのイメージングを可能にしています。間隔は鎖が伸びるほど縮みます。2イオンで5.6 $\mu$m、30イオンで1.4 $\mu$m、そしてどこかでアドレッシングビームを隣同士のあいだに絞り込めなくなります。

ジグザグの表が印象的な結果です。最低の径方向固有値が正であり続けることを要求すると、経験的な基準 $\omega_r/\omega_z \approx 0.73\,N^{0.86}$ が $N$ の10倍の範囲にわたって再現されます(比3で6イオン、比20で53イオン、基準の値は3.4と22.2)。径方向の閉じ込めは、鎖をまっすぐに保つだけのためにイオン数に対してほぼ線形に強くしなければなりません。この導出に量子力学は一切出てきません。静電気学とHessianだけです。

モード混雑の表は、全結合性の代価が現れる場所です。$\omega_r/\omega_z = 5$ 固定では最小の径方向モード間隔は101 kHz でNに依存しません — 最小のギャップが最上位2モードのあいだにあり、そのギャップは $\approx \omega_z^2/2\omega_r$ で $N$ を知らないからです。しかしその比では鎖が直線であるのは $N = 11$ までです。より大きな $N$ で直線性を安定性基準が要求するように $\omega_r$ を上げて強制すると、同じギャップは $1/\omega_r$ で縮みます。$N = 14$ で51 kHz、$N = 30$ で28 kHzです。安定限界でのフィッティングされたスケーリングは $N^{-0.80}$ です。

1つのモードをアドレスするゲートは、隣接モードからスペクトル的に分解されていなければならないので、離調 — したがってゲート時間 $\sim 2\pi/\delta$ — はその間隔で縛られます。含意される最小ゲート時間は $N = 2$ の10 $\mu$sから $N = 30$ の36 $\mu$sへ、$N^{0.8}$ で伸びます。これが長い鎖の本当の代価であり、電極の数でもレーザー出力でもありません。だからこの分野は1つの鎖を長くするのではなく、短い鎖を多数つなぐアーキテクチャ — 分割トラップ上でイオンをゾーン間に輸送する、あるいは別々のトラップを光子でつなぐ — に向かったのです。

$N = 5$ の参加率行列が結合性を具体的にします。どのイオンもほぼすべてのモードで非ゼロの振幅をもつので、イオン $i$ へのビームとイオン $j$ へのビームは、どれだけ離れていても同じ振動子に結合します。そこが全結合性の由来です。そして同時にクロストークの由来でもあります。同じことが他のすべてのペアについても同時に成り立つからです。


3.3 レーザー冷却とサイドバンド分光

2段階、2つの限界

冷却は2段階で行われ、その2つの限界はまったく別の性質のものです。

Doppler冷却は、共鳴から赤方離調した広い双極子許容遷移を使います。レーザーに向かって動くイオンはDopplerシフトで共鳴に入り優先的に吸収するので、光子散乱が運動に逆らいます。限界は自発放出光子のランダムな反跳と冷却力の釣り合いが決め、

$$ T_\mathrm{Doppler} = \frac{\hbar\Gamma}{2k_B} $$

となります。これが何を述べているかに注意してください。広い遷移は速く冷やし、早く止まるのです。$\Gamma/2\pi = 21.6$ MHzのCa$^+$の397 nm線なら限界は0.5 mKで、2 MHzのトラップでは $\bar{n} \approx 5$ を意味します。捕捉と撮像には十分冷たく、ゲートにははるかに暖かすぎます。Code Example 4 が示すように、サイドバンドのRabi周波数が $n$ に依存するからです。

分解サイドバンド冷却は逆を要求します。トラップ周波数に比べて狭い遷移であり、キャリアも青サイドバンドも触らずに赤サイドバンド($|g, n\rangle \to |e, n-1\rangle$)を駆動できることです。1サイクルごとに量子が1つ取り除かれます。残留する加熱は青サイドバンドの非共鳴励起から来て、線幅とサイドバンド間隔の比で抑えられ、定常状態は

$$ \bar{n}_\mathrm{ss} \approx \left(\frac{\Gamma}{4\omega}\right)^2 $$

の熱的状態になります。

Lamb-Dickeパラメータ

内部状態と運動の結合は、1つの無次元数

$$ \eta = k x_0 = \frac{2\pi}{\lambda}\sqrt{\frac{\hbar}{2m\omega}} $$

が支配します。$k = 2\pi/\lambda$ は(実効的な)レーザー波数ベクトルなので、$\eta$ はイオンの基底状態の広がりとレーザー波長の比そのものではなく、その $2\pi$ 倍です。相互作用演算子は $e^{i\eta(a + a^\dagger)}$ で、その行列要素は

$$ \langle n + s|e^{i\eta(a+a^\dagger)}|n\rangle = e^{-\eta^2/2}\,(i\eta)^s\sqrt{\frac{n!}{(n+s)!}}\,L_n^{s}(\eta^2) $$

で、$L_n^s$ は一般化Laguerre多項式です。Lamb-Dicke領域とは条件 $\eta^2(2n+1) \ll 1$ のことです。$\eta$ が小さいことイオンが冷たいことの両方が要ります。波長に比べて小さくなければならないのは基底状態の広がりではなく波束の広がりだからです。そこではキャリアが強く、第1サイドバンドは $\eta\sqrt{n+1}$(青)と $\eta\sqrt{n}$(赤)で効き、標準的なゲート理論はすべてこの極限で導かれています。条件がぎりぎりしか満たされないときの代価は Code Example 6 が測ります。

Code Example 4: サイドバンドを厳密に、そして近似的に

"""第3章 Code Example 4: Lamb-Dickeパラメータとサイドバンドスペクトル。
exp(i eta (a + a^dag)) の厳密な行列要素、熱的状態に対するサイドバンド
スペクトルのシミュレーション、そしてサイドバンド比温度測定。"""
import numpy as np
from scipy.linalg import expm
from scipy.special import eval_genlaguerre, factorial

hbar = 1.054571817e-34
u_amu = 1.66053906660e-27
kB = 1.380649e-23
TWOPI = 2.0 * np.pi


def lamb_dicke(mass_u, f_trap, wavelength, n_beams=1):
    """eta = k x0、x0 = sqrt(hbar / 2 m omega)。対向伝播するRaman2本の場合は
    実効的な k が2倍になるので n_beams = 2 とする。"""
    m = mass_u * u_amu
    w = TWOPI * f_trap
    x0 = np.sqrt(hbar / (2.0 * m * w))
    return n_beams * TWOPI / wavelength * x0, x0


print(f"{'ion':>6}{'m (u)':>9}{'lambda (nm)':>13}{'beams':>7}"
      f"{'f_trap (MHz)':>14}{'x0 (nm)':>10}{'eta':>9}   transition")
cases = [
    ("Ca+", 40.078, 729e-9, 1, 2.0e6, "四重極子光学量子ビット"),
    ("Ca+", 40.078, 729e-9, 1, 1.0e6, "同じ、より弱いトラップ"),
    ("Sr+", 87.62, 674e-9, 1, 1.0e6, "四重極子光学量子ビット"),
    ("Yb+", 171.0, 355e-9, 2, 3.0e6, "Raman対、超微細量子ビット"),
    ("Be+", 9.012, 313e-9, 2, 3.0e6, "Raman対、超微細量子ビット"),
]
for name, m_u, lam, nb, ft, trans in cases:
    eta, x0 = lamb_dicke(m_u, ft, lam, nb)
    print(f"{name:>6}{m_u:>9.3f}{lam * 1e9:>13.0f}{nb:>7}{ft / 1e6:>14.1f}"
          f"{x0 * 1e9:>10.3f}{eta:>9.4f}   {trans}")
print()

# --- 厳密な行列要素 ----------------------------------------------------
def rabi_matrix(eta, nmax=90):
    """行列指数関数による |<m| exp(i eta (a + a^dag)) |n>|。"""
    a = np.diag(np.sqrt(np.arange(1, nmax + 1)), 1)
    return np.abs(expm(1j * eta * (a + a.conj().T)))


def rabi_analytic(eta, n, s):
    """s >= 0 に対する解析形 |<n+s| exp(i eta (a+a^dag)) |n>|。"""
    return (np.exp(-eta ** 2 / 2.0) * eta ** s
            * np.sqrt(factorial(n) / factorial(n + s))
            * abs(eval_genlaguerre(n, s, eta ** 2)))


eta = 0.0684
M = rabi_matrix(eta)
print(f"eta = {eta:.4f} でのサイドバンド行列要素 "
      "(Ca+、729 nm、2 MHzトラップ)")
print(f"{'n':>4}{'carrier':>12}{'red (n-1)':>12}{'blue (n+1)':>12}"
      f"{'analytic blue':>15}{'LD: eta sqrt(n+1)':>19}{'LD error':>10}")
for n in [0, 1, 5, 10, 20, 40, 80]:
    car, blue = M[n, n], M[n + 1, n]
    red = M[n - 1, n] if n > 0 else 0.0
    ld = eta * np.sqrt(n + 1)
    print(f"{n:>4}{car:>12.6f}{red:>12.6f}{blue:>12.6f}"
          f"{rabi_analytic(eta, n, 1):>15.6f}{ld:>19.6f}"
          f"{(ld / blue - 1) * 100:>9.2f}%")
print()

# --- サイドバンドスペクトルのシミュレーション --------------------------
def thermal(nbar, nmax):
    """熱的なFock分布。オーバーフローを避けるため対数空間で計算する。"""
    if nbar == 0.0:
        p = np.zeros(nmax)
        p[0] = 1.0
        return p
    n = np.arange(nmax)
    p = np.exp(n * np.log(nbar) - (n + 1) * np.log1p(nbar))
    return p / p.sum()


def excitation(eta, nbar, om_car, t_pulse, f_trap, detuning, nmax=80):
    """1つの離調における励起確率。

    各Fock状態は、それぞれのサイドバンド上で離調した2準位Rabi振動を行う。
    遷移は独立に扱うが、これはRabi周波数がトラップ周波数より小さいときに
    正当化される。非共鳴のキャリアは残してある。この手法の限界を決めるのが
    まさにそれだからである。
    """
    Mm = rabi_matrix(eta, nmax + 6)
    P = thermal(nbar, nmax)
    pe = 0.0
    for s in range(-3, 4):
        off = detuning - s * f_trap
        for n in range(nmax):
            if n + s < 0:
                continue
            om = om_car * Mm[n + s, n]
            geff = np.hypot(om, off)
            if geff == 0.0:
                continue
            pe += P[n] * (om / geff) ** 2 * np.sin(np.pi * geff * t_pulse) ** 2
    return pe


f_trap, om_car = 2.0e6, 50.0e3
t_pi = 1.0 / (2.0 * om_car)
print(f"スペクトルのシミュレーション: f_trap = {f_trap / 1e6:.1f} MHz、"
      f"キャリアRabi {om_car / 1e3:.0f} kHz、パルス長 = キャリアのπ時間 = "
      f"{t_pi * 1e6:.1f} us")
print(f"{'line':>10}{'detuning':>12}{'nbar=0':>10}{'nbar=1':>10}"
      f"{'nbar=5':>10}{'nbar=20':>10}")
for s, lab in [(-2, "第2赤"), (-1, "第1赤"), (0, "キャリア"),
               (1, "第1青"), (2, "第2青")]:
    vals = [excitation(eta, nb, om_car, t_pi, f_trap, s * f_trap)
            for nb in [0.0, 1.0, 5.0, 20.0]]
    print(f"{lab:>10}{s:>+9d} x f" + "".join(f"{v:>10.5f}" for v in vals))
print()
off_bg = excitation(eta, 0.0, om_car, t_pi, f_trap, -f_trap)
print(f"  nbar = 0 の赤サイドバンドに残る {off_bg:.2e} は非共鳴キャリアである:")
print(f"  上限は (Omega/2 f_trap)^2 = "
      f"{(om_car / (2 * f_trap)) ** 2:.2e} で、パルス長とともに振動する。")
print()

# --- サイドバンド比温度測定 --------------------------------------------
print(f"{'nbar (true)':>13}{'red strength':>15}{'blue strength':>15}"
      f"{'ratio r':>11}{'r/(1-r)':>11}{'error':>9}")
Mt = rabi_matrix(eta, 400)
for nb in [0.02, 0.2, 1.0, 5.0, 20.0]:
    nmx = min(380, int(40 + 12 * nb))
    P = thermal(nb, nmx)
    red = sum(P[n] * Mt[n - 1, n] ** 2 for n in range(1, nmx))
    blue = sum(P[n] * Mt[n + 1, n] ** 2 for n in range(nmx))
    r = red / blue
    est = r / (1.0 - r)
    print(f"{nb:>13.2f}{red:>15.4e}{blue:>15.4e}{r:>11.6f}{est:>11.4f}"
          f"{(est / nb - 1) * 100:>8.3f}%")
print()
for nb in [0.02, 5.0]:
    T = TWOPI * f_trap * hbar / (kB * np.log(1.0 + 1.0 / nb))
    print(f"  f_trap = {f_trap / 1e6:.1f} MHz で nbar = {nb:>5.2f} は "
          f"T = {T * 1e6:>6.1f} uK に相当")
   ion    m (u)  lambda (nm)  beams  f_trap (MHz)   x0 (nm)      eta   transition
   Ca+   40.078          729      1           2.0     7.940   0.0684   四重極子光学量子ビット
   Ca+   40.078          729      1           1.0    11.229   0.0968   同じ、より弱いトラップ
   Sr+   87.620          674      1           1.0     7.595   0.0708   四重極子光学量子ビット
   Yb+  171.000          355      2           3.0     3.139   0.1111   Raman対、超微細量子ビット
   Be+    9.012          313      2           3.0    13.672   0.5489   Raman対、超微細量子ビット

eta = 0.0684 でのサイドバンド行列要素 (Ca+、729 nm、2 MHzトラップ)
   n     carrier   red (n-1)  blue (n+1)  analytic blue  LD: eta sqrt(n+1)  LD error
   0    0.997663    0.000000    0.068240       0.068240           0.068400     0.23%
   1    0.992996    0.068240    0.096280       0.096280           0.096732     0.47%
   5    0.974434    0.151165    0.165205       0.165205           0.167545     1.42%
  10    0.951476    0.211279    0.221070       0.221070           0.226857     2.62%
  20    0.906366    0.291804    0.298300       0.298300           0.313448     5.08%
  40    0.819309    0.393364    0.397289       0.397289           0.437974    10.24%
  80    0.657389    0.504221    0.506095       0.506095           0.615600    21.64%

スペクトルのシミュレーション: f_trap = 2.0 MHz、キャリアRabi 50 kHz、パルス長 = キャリアのπ時間 = 10.0 us
      line    detuning    nbar=0    nbar=1    nbar=5   nbar=20
       第2赤       -2 x f   0.00000   0.00003   0.00064   0.00704
       第1赤       -1 x f   0.00000   0.01126   0.05279   0.16363
      キャリア       +0 x f   0.99999   0.99977   0.99689   0.96993
       第1青       +1 x f   0.01145   0.02251   0.06334   0.17233
       第2青       +2 x f   0.00003   0.00011   0.00092   0.00790

  nbar = 0 の赤サイドバンドに残る 2.37e-07 は非共鳴キャリアである:
  上限は (Omega/2 f_trap)^2 = 1.56e-04 で、パルス長とともに振動する。

  nbar (true)   red strength  blue strength    ratio r    r/(1-r)    error
         0.02     9.3117e-05     4.7490e-03   0.019608     0.0200  -0.000%
         0.20     9.2961e-04     5.5776e-03   0.166667     0.2000  -0.000%
         1.00     4.6135e-03     9.2269e-03   0.500000     1.0000  -0.000%
         5.00     2.2227e-02     2.6672e-02   0.833333     5.0000  -0.000%
        20.00     7.7594e-02     8.1474e-02   0.952381    19.9999  -0.000%

  f_trap = 2.0 MHz で nbar =  0.02 は T =   24.4 uK に相当
  f_trap = 2.0 MHz で nbar =  5.00 は T =  526.5 uK に相当

着目点。 $\eta$ の表が範囲を示します。Ca$^+$の光学量子ビットで0.068、Yb$^+$のRaman駆動超微細量子ビットで0.11(Raman対は実効 $k$ を2倍にします)、そしてBe$^+$では0.55 — 軽いので非局在しているのです。小さな $\eta$ は弱いサイドバンド、したがって遅いゲートを意味します。大きな $\eta$ は速いゲートと、Lamb-Dicke極限で導かれたすべてへの大きな補正を意味します。Be$^+$は本当に別の領域です。

行列要素の表はLaguerreの式を行列指数関数と印字された全桁で照合し、続いてLamb-Dicke近似が破れていくのを示します。$n = 0$ で0.2%、$n = 5$ で1.4%、$n = 40$ で10%、$n = 80$ で22%の誤差で、$\eta^2(n+1)$ が支配します。これはノイズの問題ではありません。熱いイオンは間違ったRabi周波数をもつので、校正された $\pi$ パルスが $\pi$ パルスではなく、しかもその誤差は $n$ が揺らぐのでショットごとに異なります。それはコヒーレントな誤差であり、平均化では消えません。

シミュレートしたスペクトルは実験が見る構造を示します。$\bar{n} = 0$ では赤サイドバンドが厳密に消えます — 移るべき下のフォノン状態がないのです — そしてこの消失が運動の基底状態についての最も清潔に得られる署名です。$\bar{n}$ が上がると赤サイドバンドは青に向かって育ち、2次のサイドバンドが現れます。$\bar{n} = 0$ の赤サイドバンドに残る $2\times10^{-7}$ は非共鳴キャリアで、$(\Omega/2f_\mathrm{trap})^2 = 1.6\times10^{-4}$ で上から抑えられ、パルス長とともに振動します。これは実在する背景であり、基底状態冷却が短く強いサイドバンドパルスではなく長く弱いパルスで検証される理由です。

温度測定の表が嬉しい驚きです。逆変換 $\bar{n} = r/(1-r)$ は3桁にわたって印字された全桁で厳密です — $\eta$ の1次までではなく、厳密に。理由はLamb-Dicke近似ではなく詳細釣り合いです。熱的状態は $P(n) = P(n+1)e^{\hbar\omega/k_BT}$ を満たし、行列要素は $|\langle n|\ldots|n+1\rangle| = |\langle n+1|\ldots|n\rangle|$ を満たすので、2つの重み付き和の比は行列要素が何であれBoltzmann因子になります。サイドバンド温度測定は何の絶対校正も必要としないのであり、だからこの分野の冷却の成果はすべて $\bar{n}$ で報告されるのです。

Code Example 5: 基底状態まで冷やす

""" 第3章 Code Example 5: Doppler限界から運動の基底状態まで。
Doppler冷却限界を数値で押さえ、続いて分解サイドバンド冷却の速度方程式を
積分する。"""
import numpy as np
from scipy.linalg import expm

hbar = 1.054571817e-34
h = 2.0 * np.pi * hbar
u_amu = 1.66053906660e-27
kB = 1.380649e-23
TWOPI = 2.0 * np.pi

# --- Doppler限界 -------------------------------------------------------
print(f"{'ion':>6}{'transition':>14}{'Gamma/2pi (MHz)':>17}{'T_D (uK)':>11}"
      f"{'v_rms (m/s)':>13}{'nbar at 2 MHz':>15}")
f_trap = 2.0e6
for name, m_u, lab, gam in [("Ca+", 40.078, "397 nm S-P", 21.6e6),
                            ("Sr+", 87.62, "422 nm S-P", 21.5e6),
                            ("Yb+", 171.0, "369 nm S-P", 19.6e6),
                            ("Be+", 9.012, "313 nm S-P", 19.4e6)]:
    T = hbar * TWOPI * gam / (2.0 * kB)
    m = m_u * u_amu
    vrms = np.sqrt(kB * T / m)
    nbar = 1.0 / (np.exp(h * f_trap / (kB * T)) - 1.0)
    print(f"{name:>6}{lab:>14}{gam / 1e6:>17.1f}{T * 1e6:>11.1f}"
          f"{vrms:>13.4f}{nbar:>15.2f}")
print()

# --- 狭い線は「量子ビット」であって第2の冷却段ではない -------------------
# Ca+ 729 nm S1/2-D5/2 は電気四重極遷移で、D5/2 の寿命は 1.168 s なので
# Gamma/2pi = 0.136 Hz。ここから2つのことが従い、どちらも冷却ではない。
# (i) Doppler の式 T_D = hbar Gamma / 2 kB はこれよりずっと手前で下限では
# なくなる。単一光子の反跳エネルギーが hbar Gamma より大きいからである。
# (ii) 散乱率が双極子線より9桁小さく、実験室の時間スケールでは何も冷えない。
tau_D5 = 1.168                              # s, Ca+ D5/2 の寿命
gam_q = 1.0 / tau_D5                        # s^-1
lam_q = 729.147e-9
k_q = TWOPI / lam_q
m_ca = 40.078 * u_amu
E_r = hbar ** 2 * k_q ** 2 / (2.0 * m_ca)   # 単一光子の反跳エネルギー
print(f"Ca+ 729 nm S1/2-D5/2 四重極線(光学量子ビット):")
print(f"  寿命 {tau_D5:.3f} s -> Gamma/2pi = {gam_q / TWOPI:.3f} Hz, "
      f"Q = f/(Gamma/2pi) = {(2.998e8 / lam_q) / (gam_q / TWOPI):.2e}")
print(f"  Doppler の式なら T_D = {hbar * gam_q / (2 * kB):.2e} K だが、"
      f"反跳限界は")
print(f"  T_r = 2 E_r/kB = {2 * E_r / kB * 1e9:.0f} nK(E_r/h = "
      f"{E_r / h / 1e3:.2f} kHz)。hbar*Gamma は E_r の "
      f"{(E_r / h) / (gam_q / TWOPI):.1e} 分の1なので、")
print(f"  Doppler の式はここでは無意味であり、散乱率は 397 nm 線の "
      f"{gam_q / (TWOPI * 21.6e6):.1e} 倍しかない。")
print()

# --- 分解サイドバンド冷却を速度方程式として扱う ------------------------
def rates(om_car, gamma_eff, f_trap, eta):
    """量子1個あたりの冷却率と加熱率(1/s)。

    キャリアRabi周波数 Omega の赤サイドバンド駆動と、実効線幅 Gamma の
    励起状態のもとで、n -> n-1 の速度は量子あたり
        R(delta) = eta^2 Omega^2 Gamma / [(Gamma/2)^2 + delta^2]
    となる。冷却は共鳴(delta = 0)、加熱は非共鳴(delta = 2 omega、
    すなわち青サイドバンド)で評価する。Omega < Gamma で有効。
    """
    Om, G, w = TWOPI * om_car, TWOPI * gamma_eff, TWOPI * f_trap
    Rc = eta ** 2 * Om ** 2 * G / (G / 2.0) ** 2
    Rh = eta ** 2 * Om ** 2 * G / ((G / 2.0) ** 2 + (2.0 * w) ** 2)
    return Rc, Rh


def cool(nbar0, Rc, Rh, nmax=60, tmax=None, nt=400):
    """冷却率 Rc*n、加熱率 Rh*(n+1) として dP_n/dt を積分する。"""
    n = np.arange(nmax)
    P = np.exp(n * np.log(nbar0) - (n + 1) * np.log1p(nbar0))
    P /= P.sum()
    A = np.zeros((nmax, nmax))
    for k in range(nmax):
        if k > 0:
            A[k - 1, k] += Rc * k
            A[k, k] -= Rc * k
        if k < nmax - 1:
            A[k + 1, k] += Rh * (k + 1)
            A[k, k] -= Rh * (k + 1)
    tmax = tmax if tmax else 25.0 / Rc
    out = []
    for t in np.linspace(0.0, tmax, nt):
        Pt = expm(A * t) @ P
        out.append((t, float(np.sum(n * Pt)), float(Pt[0])))
    return out


eta = 0.0684
om_car = 50.0e3
T_D = hbar * TWOPI * 21.6e6 / (2.0 * kB)
nbar_D = 1.0 / (np.exp(h * f_trap / (kB * T_D)) - 1.0)
print(f"分解サイドバンド冷却、f_trap = {f_trap / 1e6:.1f} MHz の Ca+、"
      f"eta = {eta:.4f}、")
print(f"キャリアRabi周波数 {om_car / 1e3:.0f} kHz、Doppler限界の "
      f"nbar = {nbar_D:.2f} から出発する。")
print(f"{'Gamma_eff/2pi':>16}{'R_c (1/s)':>12}{'r = R_h/R_c':>14}"
      f"{'nbar steady':>14}{'(Gamma/4 omega)^2':>20}{'P(n=0)':>10}")
for gam in [10.0e6, 2.0e6, 200.0e3, 100.0e3]:
    Rc, Rh = rates(om_car, gam, f_trap, eta)
    r = Rh / Rc
    nss = r / (1.0 - r)
    approx = (gam / (4.0 * f_trap)) ** 2
    print(f"{gam / 1e6:>14.3f} M{Rc:>12.3e}{r:>14.3e}{nss:>14.3e}"
          f"{approx:>20.3e}{1.0 / (1.0 + nss):>10.6f}")
print()

Rc, Rh = rates(om_car, 200.0e3, f_trap, eta)
traj = cool(nbar_D, Rc, Rh)
print(f"Gamma_eff/2pi = 200 kHz での冷却の時間発展 "
      f"(量子あたり R_c = {Rc:.1f} 1/s):")
print(f"{'t (ms)':>10}{'nbar':>12}{'P(n=0)':>12}")
for k in [0, 12, 25, 50, 100, 200, 399]:
    t, nb, p0 = traj[k]
    print(f"{t * 1e3:>10.3f}{nb:>12.6f}{p0:>12.6f}")
print()

# --- 残留 nbar がゲートに与える代価 ------------------------------------
eta = 0.0684
for nb in [0.01, 0.05, 0.5, nbar_D]:
    nmx = min(300, int(40 + 12 * nb))
    n = np.arange(nmx)
    P = np.exp(n * np.log(nb) - (n + 1) * np.log1p(nb))
    P /= P.sum()
    om = eta * np.sqrt(n + 1)                 # Lamb-Dicke近似のサイドバンド強度
    mean = float(np.sum(P * om))
    rms = float(np.sqrt(np.sum(P * om ** 2) - mean ** 2))
    print(f"  nbar = {nb:>5.2f}: サイドバンド強度の平均 {mean:.5f}、"
          f"広がり {rms:.5f}(平均の {rms / mean * 100:>5.2f}%)")
print()
   ion    transition  Gamma/2pi (MHz)   T_D (uK)  v_rms (m/s)  nbar at 2 MHz
   Ca+    397 nm S-P             21.6      518.3       0.3279           4.92
   Sr+    422 nm S-P             21.5      515.9       0.2213           4.89
   Yb+    369 nm S-P             19.6      470.3       0.1512           4.42
   Be+    313 nm S-P             19.4      465.5       0.6554           4.37

Ca+ 729 nm S1/2-D5/2 四重極線(光学量子ビット):
  寿命 1.168 s -> Gamma/2pi = 0.136 Hz, Q = f/(Gamma/2pi) = 3.02e+15
  Doppler の式なら T_D = 3.27e-12 K だが、反跳限界は
  T_r = 2 E_r/kB = 899 nK(E_r/h = 9.36 kHz)。hbar*Gamma は E_r の 6.9e+04 分の1なので、
  Doppler の式はここでは無意味であり、散乱率は 397 nm 線の 6.3e-09 倍しかない。

分解サイドバンド冷却、f_trap = 2.0 MHz の Ca+、eta = 0.0684、
キャリアRabi周波数 50 kHz、Doppler限界の nbar = 4.92 から出発する。
   Gamma_eff/2pi   R_c (1/s)   r = R_h/R_c   nbar steady   (Gamma/4 omega)^2    P(n=0)
        10.000 M   2.940e+01     6.098e-01     1.563e+00           1.562e+00  0.390244
         2.000 M   1.470e+02     5.882e-02     6.250e-02           6.250e-02  0.941176
         0.200 M   1.470e+03     6.246e-04     6.250e-04           6.250e-04  0.999375
         0.100 M   2.940e+03     1.562e-04     1.562e-04           1.563e-04  0.999844

Gamma_eff/2pi = 200 kHz での冷却の時間発展 (量子あたり R_c = 1469.8 1/s):
    t (ms)        nbar      P(n=0)
     0.000    4.914527    0.169052
     0.512    2.318517    0.301306
     1.066    1.027613    0.493152
     2.131    0.215262    0.822853
     4.263    0.010000    0.990099
     8.526    0.000643    0.999358
    17.009    0.000625    0.999375

  nbar =  0.01: サイドバンド強度の平均 0.06868、広がり 0.00283(平均の  4.13%)
  nbar =  0.05: サイドバンド強度の平均 0.06980、広がり 0.00635(平均の  9.10%)
  nbar =  0.50: サイドバンド強度の平均 0.08122、広がり 0.02050(平均の 25.24%)
  nbar =  4.92: サイドバンド強度の平均 0.15192、広がり 0.06778(平均の 44.62%)

着目点。 Doppler の表は、この段階についてはイオンの選択がほとんど関係ないことを示します。使われている双極子許容の冷却線はどれも $\Gamma/2\pi \approx 20$ MHz なので、どの種でも0.5 mKで止まり、2 MHzのトラップで $\bar{n} \approx 4.5$ です。Doppler冷却は既製品であり、そしてそれでは足りません。$\bar{n} = 4.9$ は基底状態ではないのです。

その下のブロックが、より狭い線を選べば解決するという道がない理由を述べます。Ca$^+$ での明らかな候補は 729 nm の $S_{1/2}$-$D_{5/2}$ 四重極遷移で、上位状態の寿命は 1.168 s、したがって $\Gamma/2\pi = 0.136$ Hz、線のQ値は $3\times10^{15}$ です。これを $T_D = \hbar\Gamma/2k_B$ に代入すると $3\times10^{-12}$ K が返りますが、これはどのイオンも到達しない温度です。この式は $\hbar\Gamma$ が単一光子の反跳エネルギーより小さくなった時点で下限ではなくなり、ここでは $E_r/h = 9.4$ kHz が $\Gamma/2\pi$ の $7\times10^4$ 倍もあります。狭い線に対する本当の下限は反跳限界 $T_r = 2E_r/k_B = 0.9\ \mu$K であり、そもそもこの線での散乱率は 397 nm 線の $6\times10^{-9}$ 倍 — 毎秒およそ1光子 — なので、これでは何も冷えません。0.136 Hz の線幅は冷却の資源ではなく、その遷移を量子ビットにしているものです。したがって深い冷却には別の線ではなく別の発想が必要で、その発想が分解された運動サイドバンドで冷やすこと、すなわちこの例の後半です。

(Ca$^+$ で狭い冷却線と間違われやすい遷移が 866 nm です。これは $D_{3/2}$-$P_{1/2}$ の再励起線で、397 nm 冷却線と $P_{1/2}$ 準位を共有しているので、その自然幅はその準位の $\approx 22$ MHz です。ときに引用される $\sim$1.7 MHz は $P_{1/2}$ からの分岐の部分率であって線幅ではなく、この遷移を狭くはしません。)

サイドバンド冷却の表は $\bar{n}_\mathrm{ss} = (\Gamma/4\omega)^2$ を4桁にわたって3桁の精度で確認します。設計則として読んでください。冷却限界はサイドバンドがどれだけよく分解されているかが決めるのであって、いかなる温度でもありません。$\Gamma/2\pi = 10$ MHz — 双極子線 — では限界は $\bar{n} = 1.6$ で、Doppler冷却と大差ありません。100 kHzでは $1.6\times10^{-4}$、すなわち確率0.9998で基底状態です。$\omega$ が分母に入っているので、トラップ周波数自体が冷却パラメータです。硬いトラップはより深く冷える。これがゲートが高周波の径方向モードを使うもう1つの理由です。

時間発展は接近の形を示します。冷却率は $n$ に比例するので、高いFock状態が先に空になり、最後の2分の1がそれ以前のすべてと同じだけ時間を食います。$\bar{n}$ は1 msで4.9から1.0に落ち、$10^{-2}$ に達するのにさらに3 msかかります。実際のシーケンスは数百回のサイドバンドパルスをリポンプと交互に挟んだもので、ミリ秒かかります。この準備時間がイオントラップ型プロセッサのデューティサイクルの大部分を占めます。この機械のクロック速度は、ゲートと同じくらい冷却と読み出しが決めているのです。

最後の表が残留分の値段を付けます。$\bar{n} = 0.5$ でサイドバンドRabi周波数の広がりは平均の25%、Doppler限界では45%です。Rabi周波数の4分の1は回転角の4分の1です。だから第1世代のCirac-Zollerゲート — サイドバンド $\pi$ パルスを直接使います — は $\bar{n} \ll 1$ を要求し、そしてそれに代わったゲートは要求しないのです。


3.4 ゲートの機構

Cirac-Zoller:直接的な道筋

もとの提案(CiracとZoller、1995)は3段階のシーケンスです。赤サイドバンドの $\pi$ パルスでイオン1の内部状態を共有フォノンモードに写し、モードにフォノンがあるかどうかに依存する操作をイオン2にかけ、それからフォノンをイオン1に写し戻します。バスは運動であり、ゲートは原理的に厳密です。

実際には2つの厳しい要求があります。モードは $|n = 0\rangle$ から始まらなければなりません。サイドバンドRabi周波数が $n$ に依存するので、そうでなければ $\pi$ パルスがCode Example 5 の表にある量だけ狂うからです。そしてシーケンス中にフォノンが失われたり加熱されたりしてはなりません。中間状態では量子情報が本当に運動に格納されているからです。両方とも英雄的な努力で達成されており、どちらもスケールしません。

Mølmer-Sørensen:幾何学的な道筋

Mølmer-Sørensenゲートは、運動に実在の励起をまったく入れないことで両方の要求を取り除きます。両方のイオンを二色の場で駆動し、量子ビット周波数の両側に $\omega_0 \pm (\omega - \delta)$ で対称に配置して、どちらのトーンもどちらのサイドバンドとも共鳴しないようにします — それぞれが $\delta$ だけ離調しています。どちらのイオンも単独では光子を吸収できません。しかしなら吸収できます。2光子過程 $|00\rangle \to |11\rangle$ が各トーンから1光子ずつを使い、フォノン数の正味の変化なしにエネルギーを保存するのです。運動は仮想的に使われます。

Lamb-Dicke領域で回転波近似を行うと、相互作用は

$$ H(t) = g\left(a e^{-i\delta t} + a^\dagger e^{i\delta t}\right)\left(\sigma_y^{(1)} + \sigma_y^{(2)}\right), \qquad g = \frac{\eta\Omega}{2} $$

となります($\Omega$ はキャリアRabi周波数)。このハミルトニアンは注目すべき性質をもちます。Magnus展開が有限項で終わるのです。$H$ が $a$ と $a^\dagger$ について線形なので、交換子 $[H(t_1), H(t_2)]$ はc数と $S_y^2$ の積になり、これは $H$ と可換です。したがって高次のMagnus項はすべて消え、時間発展は厳密に

$$ U(\tau) = e^{-i\Phi S_y^2}\,D!\left(\alpha(\tau)S_y\right), \qquad S_y = \sigma_y^{(1)} + \sigma_y^{(2)} $$

$$ \alpha(\tau) = -\frac{g\left(e^{i\delta\tau} - 1\right)}{\delta}, \qquad \Phi(\tau) = -\frac{g^2}{\delta}\left(\tau - \frac{\sin\delta\tau}{\delta}\right) $$

となります。これをゲートにする条件が2つあります。ループ閉条件:$\alpha(\tau) = 0$ には整数 $K$ について $\delta\tau = 2\pi K$ が必要です — 運動の位相空間中の軌跡が出発点に戻らなければならず、さもなければ運動がスピンとエンタングルしたまま残り、スピンだけのユニタリでは結果を記述できません。ゲート角:ループが閉じていれば $\Phi = -2\pi K g^2/\delta^2$ であり、$S_y^2 = 2 + 2\sigma_y^{(1)}\sigma_y^{(2)}$ なので操作は大域位相を除いて $e^{-2i\Phi\sigma_y\sigma_y}$ です。$2|\Phi| = \pi/4$ を要求すると

$$ g = \frac{\delta}{4\sqrt{K}} \quad\Longleftrightarrow\quad \eta\Omega = \frac{\delta}{2\sqrt{K}} $$

が得られます。

そしてここが見返りです。$\Phi$ に $n$ が入っていません。ハミルトニアンはフォノン数に依存しないので、ゲートは冷たいイオンと暖かいイオンで同一に働きます。原理的にはDoppler冷却で足りるのです。この1つの事実が、以後建設されたすべてのイオントラップ型プロセッサがこのゲートを使う理由です。

Code Example 6: Mølmer-Sørensenゲートを検証する

"""第3章 Code Example 6: Mølmer-Sørensenゲートを数値的に検証する。

Lamb-Dicke形(g = eta Omega / 2、Omega はキャリアRabi周波数):
    H(t) = g (a exp(-i delta t) + a^dag exp(+i delta t)) (sy_1 + sy_2)
このMagnus展開は有限項で終わり、厳密に
    U(tau) = exp(-i Phi Sy^2) D(alpha Sy)
    alpha(tau) = -g (exp(i delta tau) - 1) / delta
    Phi(tau)   = -(g^2/delta) (tau - sin(delta tau)/delta)
となる。Schrödinger方程式を積分して各主張を確認し、続いてLamb-Dicke近似を
外して同じ計算を繰り返す。"""
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import expm

TWOPI = 2.0 * np.pi

# --- 演算子 ------------------------------------------------------------
NMAX = 30                                    # Fock空間のカットオフ
sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
sp = np.array([[0, 1], [0, 0]], dtype=complex)     # 本書の順序では |0><1|
I2 = np.eye(2, dtype=complex)
a = np.diag(np.sqrt(np.arange(1, NMAX + 1)), 1).astype(complex)
Im = np.eye(NMAX + 1, dtype=complex)


def spin(op, j):
    return np.kron(op, I2) if j == 0 else np.kron(I2, op)


SY = np.kron(spin(sy, 0) + spin(sy, 1), Im)
A = np.kron(np.eye(4, dtype=complex), a)
AD = A.conj().T
NOP = AD @ A
DIM = 4 * (NMAX + 1)

# --- パラメータ --------------------------------------------------------
eta = 0.0684
f_mode = 2.0e6            # Hz、共有モード
K = 1                     # 位相空間ループの回数
tau = 100e-6              # s、ゲート時間
delta = K / tau           # Hz。delta*tau = K サイクルとなるように
g = TWOPI * delta / 4.0 / np.sqrt(K)      # rad/s。下記のゲート条件
om_car = 2.0 * g / (TWOPI * eta)          # Hz、キャリアRabi周波数
wdelta = TWOPI * delta
print(f"eta = {eta:.4f}、モード f = {f_mode / 1e6:.1f} MHz、"
      f"ゲート時間 tau = {tau * 1e6:.0f} us、ループ回数 K = {K}")
print(f"サイドバンドからの離調 delta/2pi = {delta / 1e3:.3f} kHz")
print(f"結合 g/2pi = {g / TWOPI / 1e3:.4f} kHz  "
      f"-> キャリアRabi Omega/2pi = {om_car / 1e3:.3f} kHz")
Phi = -(g ** 2 / wdelta) * (tau - np.sin(wdelta * tau) / wdelta)
print(f"解析解 Phi(tau) = {Phi:.9f} rad、2 Phi = {2 * Phi:.9f}、"
      f"-pi/4 = {-np.pi / 4:.9f}")
print(f"解析解 |alpha(tau)| = "
      f"{abs(-g * (np.exp(1j * wdelta * tau) - 1) / wdelta):.3e} "
      f"(delta tau が整数サイクルのとき0)")
print()


def evolve_ld(psi0, tau, delta_hz, gcoup):
    """Lamb-Dicke近似のMølmer-Sørensenハミルトニアンを積分する。"""
    wd = TWOPI * delta_hz

    def rhs(t, y):
        H = gcoup * (A * np.exp(-1j * wd * t) + AD * np.exp(1j * wd * t)) @ SY
        return -1j * (H @ y)
    s = solve_ivp(rhs, [0.0, tau], psi0, rtol=1e-10, atol=1e-12,
                  method="DOP853")
    return s.y[:, -1]


def spin_block(nfock, tau, delta_hz, gcoup):
    """4x4のスピン時間発展演算子。Fock状態 nfock から出発し、同じ状態へ
    射影して取り出す。そのユニタリ性が、運動がどれだけ完全に解けたかを
    測る指標になる。"""
    ev = evolve_ld
    U = np.zeros((4, 4), dtype=complex)
    for s0 in range(4):
        psi0 = np.zeros(DIM, dtype=complex)
        psi0[s0 * (NMAX + 1) + nfock] = 1.0
        out = ev(psi0, tau, delta_hz, gcoup)
        for s1 in range(4):
            U[s1, s0] = out[s1 * (NMAX + 1) + nfock]
    return U


# --- 解析的な目標 ------------------------------------------------------
SY4 = (np.kron(sy, I2) + np.kron(I2, sy))
U_target = expm(-1j * Phi * (SY4 @ SY4))
U0 = spin_block(0, tau, delta, g)
print("スピン時間発展演算子、行・列の順序は 00, 01, 10, 11")
print("  数値積分                        |  解析解 exp(-i Phi Sy^2)")
for r0, r1 in zip(np.round(U0, 4), np.round(U_target, 4)):
    fmt = lambda row: " ".join(f"{z.real:+.3f}{z.imag:+.3f}j" for z in row)
    print(f"  {fmt(r0)}  |  {fmt(r1)}")
print(f"  ユニタリ性の検査 ||U^dag U - I|| = "
      f"{np.linalg.norm(U0.conj().T @ U0 - np.eye(4)):.3e}")
print(f"  解析解との一致 ||U - U_target|| = "
      f"{np.linalg.norm(U0 - U_target):.3e}")
print()

# --- Bell状態 ----------------------------------------------------------
psi0 = np.zeros(DIM, dtype=complex)
psi0[0] = 1.0                                  # |00> |n=0>
out = evolve_ld(psi0, tau, delta, g)
bell = np.zeros(DIM, dtype=complex)
bell[0 * (NMAX + 1)] = 1.0 / np.sqrt(2.0)
bell[3 * (NMAX + 1)] = -1j / np.sqrt(2.0)      # (|00> - i|11>)/sqrt(2)
print("from |00>|n=0>:")
pops = [float(np.sum(np.abs(out[s * (NMAX + 1):(s + 1) * (NMAX + 1)]) ** 2))
        for s in range(4)]
print(f"  スピンの占有数 00, 01, 10, 11 = "
      + ", ".join(f"{p:.6f}" for p in pops))
print(f"  終了時の平均フォノン数          = "
      f"{float(np.real(out.conj() @ (NOP @ out))):.3e}")
print(f"  目標Bell状態 x |n=0> との忠実度 = "
      f"{abs(bell.conj() @ out) ** 2:.9f}")
print()

# --- 初期フォノン数に対する鈍感さ --------------------------------------
print(f"{'n initial':>11}{'||U(n) - U(0)||':>18}{'Bell fidelity':>16}")
for n0 in [0, 1, 3, 8]:
    Un = spin_block(n0, tau, delta, g)
    psi0 = np.zeros(DIM, dtype=complex)
    psi0[n0] = 1.0
    o = evolve_ld(psi0, tau, delta, g)
    b = np.zeros(DIM, dtype=complex)
    b[n0] = 1.0 / np.sqrt(2.0)
    b[3 * (NMAX + 1) + n0] = -1j / np.sqrt(2.0)
    print(f"{n0:>11}{np.linalg.norm(Un - U0):>18.3e}"
          f"{abs(b.conj() @ o) ** 2:>16.9f}")
print()

# --- 離調がずれたゲートに何が起きるか ----------------------------------
print(f"{'delta tau (cycles)':>19}{'|alpha| analytic':>18}"
      f"{'spin-block norm defect':>24}{'Bell fidelity':>16}")
for cycles in [1.0, 1.02, 1.05, 1.1, 1.5]:
    d = cycles / tau
    wd = TWOPI * d
    al = abs(-g * (np.exp(1j * wd * tau) - 1) / wd)
    Um = spin_block(0, tau, d, g)
    psi0 = np.zeros(DIM, dtype=complex)
    psi0[0] = 1.0
    o = evolve_ld(psi0, tau, d, g)
    print(f"{cycles:>19.2f}{al:>18.4f}"
          f"{np.linalg.norm(Um.conj().T @ Um - np.eye(4)):>24.3e}"
          f"{abs(bell.conj() @ o) ** 2:>16.6f}")
print()

# --- Lamb-Dicke近似の外へ ----------------------------------------------
D0 = expm(1j * eta * (a + a.conj().T))
nvec = np.arange(NMAX + 1)
SPS = np.array(spin(sp, 0) + spin(sp, 1))          # 4x4のスピン上昇演算子の和
SPSd = SPS.conj().T


def evolve_exact(psi0, tau, delta_hz, gcoup):
    """Lamb-Dicke展開をせずに、二色駆動の完全なハミルトニアンを積分する。

    D を1次まで展開するとLamb-Dicke形(結合 Omega eta / 2)が再現される。
    したがって与えられた g に対応するキャリアRabi周波数は
    Omega = 2 g / eta である。
    """
    wm = TWOPI * f_mode
    wbi = TWOPI * (f_mode - delta_hz)
    Om = 2.0 * gcoup / eta

    def rhs(t, y):
        ph = np.exp(1j * wm * t * nvec)
        Dt = (ph[:, None] * D0) * np.conj(ph)[None, :]
        Y = y.reshape(4, NMAX + 1)
        out = SPS @ Y @ Dt.T + SPSd @ Y @ Dt.conj()
        return (-1j * Om * np.cos(wbi * t) * out).ravel()
    s = solve_ivp(rhs, [0.0, tau], psi0, rtol=1e-9, atol=1e-11,
                  method="DOP853")
    return s.y[:, -1]


print()
print(f"{'n initial':>11}{'Bell fidelity (LD)':>21}{'Bell fidelity (exact)':>23}"
      f"{'infidelity (exact)':>20}")
for n0 in [0, 1, 3, 8]:
    b = np.zeros(DIM, dtype=complex)
    b[n0] = 1.0 / np.sqrt(2.0)
    b[3 * (NMAX + 1) + n0] = -1j / np.sqrt(2.0)
    p0 = np.zeros(DIM, dtype=complex)
    p0[n0] = 1.0
    f_ld = abs(b.conj() @ evolve_ld(p0, tau, delta, g)) ** 2
    oe = evolve_exact(p0, tau, delta, g)
    # 幾何位相の符号は駆動位相に依存するので、-i と +i の与える2つのBell
    # 状態のどちらでも受け入れる
    bp = np.zeros(DIM, dtype=complex)
    bp[n0] = 1.0 / np.sqrt(2.0)
    bp[3 * (NMAX + 1) + n0] = +1j / np.sqrt(2.0)
    f_ex = max(abs(b.conj() @ oe) ** 2, abs(bp.conj() @ oe) ** 2)
    print(f"{n0:>11}{f_ld:>21.9f}{f_ex:>23.9f}{1 - f_ex:>20.3e}")
print()
print(f"    n = 0 でのDebye-Waller減衰: 結合の {eta ** 2 / 2 * 100:.3f}%")
eta = 0.0684、モード f = 2.0 MHz、ゲート時間 tau = 100 us、ループ回数 K = 1
サイドバンドからの離調 delta/2pi = 10.000 kHz
結合 g/2pi = 2.5000 kHz  -> キャリアRabi Omega/2pi = 73.099 kHz
解析解 Phi(tau) = -0.392699082 rad、2 Phi = -0.785398163、-pi/4 = -0.785398163
解析解 |alpha(tau)| = 1.608e-16 (delta tau が整数サイクルのとき0)

スピン時間発展演算子、行・列の順序は 00, 01, 10, 11
  数値積分                        |  解析解 exp(-i Phi Sy^2)
  +0.500+0.500j +0.000+0.000j +0.000+0.000j +0.500-0.500j  |  +0.500+0.500j +0.000+0.000j +0.000+0.000j +0.500-0.500j
  +0.000+0.000j +0.500+0.500j -0.500+0.500j +0.000+0.000j  |  +0.000+0.000j +0.500+0.500j -0.500+0.500j +0.000+0.000j
  +0.000+0.000j -0.500+0.500j +0.500+0.500j +0.000+0.000j  |  +0.000+0.000j -0.500+0.500j +0.500+0.500j +0.000+0.000j
  +0.500-0.500j +0.000+0.000j +0.000+0.000j +0.500+0.500j  |  +0.500-0.500j +0.000+0.000j +0.000+0.000j +0.500+0.500j
  ユニタリ性の検査 ||U^dag U - I|| = 3.493e-12
  解析解との一致 ||U - U_target|| = 5.771e-12

from |00>|n=0>:
  スピンの占有数 00, 01, 10, 11 = 0.500000, 0.000000, 0.000000, 0.500000
  終了時の平均フォノン数          = 8.475e-25
  目標Bell状態 x |n=0> との忠実度 = 1.000000000

  n initial   ||U(n) - U(0)||   Bell fidelity
          0         0.000e+00     1.000000000
          1         3.104e-12     1.000000000
          3         8.131e-12     1.000000000
          8         1.832e-11     1.000000000

 delta tau (cycles)  |alpha| analytic  spin-block norm defect   Bell fidelity
               1.00            0.0000               3.493e-12        1.000000
               1.02            0.0308               5.349e-03        0.997182
               1.05            0.0745               3.104e-02        0.983814
               1.10            0.1405               1.073e-01        0.945000
               1.50            0.3333               5.074e-01        0.757025


  n initial   Bell fidelity (LD)  Bell fidelity (exact)  infidelity (exact)
          0          1.000000000            0.999943995           5.600e-05
          1          1.000000000            0.999760759           2.392e-04
          3          1.000000000            0.999090978           9.090e-04
          8          1.000000000            0.995730633           4.269e-03

    n = 0 でのDebye-Waller減衰: 結合の 0.234%

着目点。 解析解の $2\Phi$ は $-0.785398163$ で、9桁で $-\pi/4$ です。そして $|\alpha(\tau)| = 1.6\times10^{-16}$ — ループは浮動小数点精度で閉じています。数値積分したスピン時間発展演算子は $e^{-i\Phi S_y^2}$ と $6\times10^{-12}$ で一致し、$3.5\times10^{-12}$ でユニタリです。これは運動が本当に因子化されたことを意味します。$|00\rangle|n=0\rangle$ からゲートは $|00\rangle$ と $|11\rangle$ にちょうど半分ずつのpopulationを作り、$8\times10^{-25}$ 個のフォノンを残し、忠実度1.000000000でBell状態に到達します。

熱的鈍感さの表が重要なもので、数値はほとんど清潔すぎるほどです。$n = 0$、1、3、8 から出発しても同じスピン時間発展演算子が $2\times10^{-11}$ で得られます。これは近似的な相殺ではなく、$a$ について線形で $\hat{n}$ に依存しないハミルトニアンの厳密な性質であって、数値はそれを発見したのではなく確認しているのです。

離調ずれの表は「ループ閉」が実務上何を意味するかを示します。$\delta\tau = 1.02$ サイクルではループが $|\alpha| = 0.031$ だけ外れ、スピンブロックは $5\times10^{-3}$ だけユニタリ性を失い、Bell忠実度は0.997に落ちます。1.10サイクルでは0.945です。この破綻は校正で消せる位相誤差ではありません。残留変位はスピン状態についての経路情報を運動が持ち去ることを意味するので、2量子ビット操作はそもそもユニタリではないのです。2%のタイミング誤差が $3\times10^{-3}$ の忠実度を食い、これがパルスとモード周波数への安定性要求を定めます。

最後の表はLamb-Dicke近似を破っており、それに要した追加の仕掛けに見合う価値があります。変位演算子を $D(t) = R(\omega t)D(i\eta)R^\dagger(\omega t)$($R$ は対角)と書けば厳密な二色ハミルトニアンの積分が安価になり、結果として理想的なゲートは $n = 0$ で $5.6\times10^{-5}$ — $\eta^4$ の程度 — を失い、$n = 8$ では $4.3\times10^{-3}$ になります。

この損失を2つの機構が分け合っており、壊れ方が違うので分けておく価値があります。第一はDebye-Waller因子 $e^{-\eta^2(2n+1)/2}$ で、サイドバンド結合を $n = 0$ で0.23%、$n = 8$ で3.9%減らします。幾何位相は結合の2乗で効くので $\Phi$ が $\pi/4$ に届かなくなり、$n = 8$ では7.9%足りません。Debye-Waller因子の予測7.8%に対してこれは $\eta^4$ の見積もりとしては十分な一致です。その不足がどう現れるかに注意してください。$|00\rangle$ と $|11\rangle$ の相対位相が狂うのではありません(厳密積分はそれを $2.6\times10^{-5}$ rad の精度で当てています)。現れるのは移行の不完全さです。占有は半分ずつではなく $\cos^2\Phi$ と $\sin^2\Phi$ になり — $n = 0$ で1.1%、$n = 8$ で12.4%の不均衡 — 結果として生じる非忠実度 $\cos^2 2\Phi/4$ は $n = 0$ で $3.0\times10^{-5}$、$n = 8$ で $3.9\times10^{-3}$ です。つまりこのチャネルは $n = 0$ で損失の53%、$n = 8$ で90%を供給し、回転の行き過ぎ/不足なので、$\bar{n}$ が既知なら原理的には校正で除けます。

第二は残留するスピン-運動エンタングルメントです。ループがもはや厳密に閉じないので、$n = 0$ で占有の $2.6\times10^{-5}$($n = 8$ では $4.0\times10^{-4}$)が $|01\rangle$ と $|10\rangle$ に、しかも運動がフォノン1個ぶん変位した状態で残ります。この部分はパルス面積の調整では校正できません。運動がどちらの経路を通ったかの情報を持ち去るからです — 離調の表と同じ破綻モードの縮小版です。$n = 0$ では2つのチャネルは同程度で、イオンが暖まると校正可能なほうが速く増えます。

つまり熱的鈍感さは主要項の言明であって厳密なものではなく、補正は $\bar{n}$ と $\eta$ の両方とともに増えます。だからMølmer-Sørensenゲートは厳密には必要としないにもかかわらず基底状態冷却が先立つのであり、強い結合を与える軽いイオンが最大の非Lamb-Dicke代価を払うのです。厳密な計算はどちらかを教えてくれますが、近似的な計算は教えてくれません。


3.5 全結合性とその代償

共有バスが買うもの

どのイオンもどのモードにも参加しているので、任意のペアを直接エンタングルできます。超伝導チップ上では離れた量子ビット間のゲートにSWAPの連鎖が必要で、一般の回路に対する深さのオーバーヘッドは大きくなりえます。イオン鎖ではゲート1個です。自然な相互作用グラフが平面的でないアルゴリズム — 全結合のCoulomb項をもつ量子化学が標準的な例です — に対してこれは本物のアーキテクチャ上の利点であり、単一ゲートの数値ではなく実効的な回路の深さに現れます。

同じ物理からさらに2つの利点が従います。ある種のイオンはすべて同一なので周波数衝突の問題がなく作製歩留まりもありません。イオンの鎖には、多量子ビットチップの混雑して個別に校正された周波数地形に相当するものが存在しないのです。そして量子ビットのコヒーレンス時間は秒であり、10〜100 $\mu$s のゲート時間の4〜5桁上です。

代償

ゲート速度。 バスは数MHzの力学的振動子であり、ゲートの離調 $\delta$ はモード間隔より小さくなければなりません。ゲート時間は数十から数百マイクロ秒で、超伝導回路の数十ナノ秒に対するものです。秒のコヒーレンスに対してこの比はなお有利ですが、深い回路の実時間はそうではありません。

モードの混雑。 Code Example 3 が定量化しました。直線安定性の限界では最小の径方向モード間隔が $N^{-0.8}$ で落ちるので、最小ゲート時間は $N^{0.8}$ で伸びます。30イオンの鎖は、他の一切を考える前に2イオンの鎖より4倍遅いゲートを必要とします。

傍観者モード。 ゲートは1つのモードをアドレスしますが、駆動はスペクトル的に清潔ではなく、隣接モードへの残留結合が残留スピン-運動エンタングルメントを残します — Code Example 6 の離調ずれの破綻が、望まないモード1つごとに1回起きるのです。すべての位相空間ループを同時に閉じる多トーン・振幅整形パルスが標準的な対処法で、パルス整形の帯域と校正の手間を代価に要します。

冷却と読み出しのオーバーヘッド。 実験サイクルごとに鎖を再冷却し(ミリ秒、Code Example 5)、状態依存蛍光で読み出します(数百マイクロ秒)。どちらも $N$ に対して厄介にスケールします。

アーキテクチャ側の応答

モード混雑への応答は、それを直すのではなく回避することです。鎖を短く保って、それらをつなぐのです。

QCCD(quantum charge-coupled device)は多数の電極をもつ分割トラップを使い、いくつかの短い鎖を別々のゾーンに保持し、直流電圧をランプさせてイオンをゾーン間で物理的に輸送します。ゲートは常に短い鎖で起こるのでモードスペクトルは清潔に保たれ、結合性は輸送が与えます。原理的にはイオンは断熱的に動かされて運動状態は保たれますが、実際には輸送はイオンを加熱し、その加熱率は3.6節のすべてと同じ電場ノイズが決めます。輸送はゲートに比べて遅くもあるので、実効的な結合性は時間で買われます。

光子インターコネクトは、物理的に別のトラップにあるイオンが放出した光子を干渉させて同時計数を検出することでエンタングルします。試行あたりの成功確率は小さい — 立体角、ファイバ結合、検出器効率 — のでリンク率は低いのですが、heraldedです。検出パターンが正しければ、エンタングルメントはそこにあります。これは難しいスケーリングの問題を率の問題に変換し、率の問題はより良い種類の問題です。

どちらも解決された工学ではなく原理であり、そしてどちらも同じ根底の材料の問題に律速されています。


3.6 異常加熱:表面材料の問題

観測されること

$|n = 0\rangle$ に冷やして暗闇で放置したイオンは暖まります。その率は、可変の時間だけ待ってサイドバンド温度測定を行うことで測られ、トラップ回路のJohnsonノイズで説明できる値より1桁から数桁大きく出ます。これが異常加熱であり、最初に測定されて以来この方式の中心的な技術的障害です。

物理は述べるのは簡単です。イオン位置での揺らぐ電場が運動を駆動し、$\mathrm{V}^2\mathrm{m}^{-2}\mathrm{Hz}^{-1}$ 単位の電場ノイズスペクトル密度 $S_E(\omega)$ について

$$ \frac{d\bar{n}}{dt} = \frac{e^2 S_E(\omega)}{4m\hbar\omega} $$

(1価イオンの電荷はここでは $q$ ではなく $e$ と書きます。$q$ は3.1節のMathieuパラメータに使われているためです)

です。困難のすべてが $S_E$ にあり、そして $S_E$ は材料量です。

スケーリングが語ること

経験的事実と、それぞれが含意するもの:

まとめると、ノイズは電極表面の最後の数ナノメートルにある吸着物とパッチ電位の熱活性化的な運動から来ます。それはPaulトラップの性質ではありません。

Code Example 7: 加熱のバジェット

"""第3章 Code Example 7: 異常加熱を表面材料の問題として扱う。

イオン位置での電場ノイズが運動を駆動する:
    dn/dt = e^2 S_E(omega) / (4 m hbar omega),  e はイオンの電荷
したがって加熱率は、電極表面のノイズスペクトル密度 S_E という1つの材料量と、
それを通じてトラップ形状によって決まる。"""
import numpy as np

hbar = 1.054571817e-34
e = 1.602176634e-19
u_amu = 1.66053906660e-27
kB = 1.380649e-23
TWOPI = 2.0 * np.pi


def heating_rate(S_E, mass_u, f_trap):
    """V^2 m^-2 Hz^-1 単位の電場ノイズ密度 S_E から毎秒の量子数を返す。"""
    m = mass_u * u_amu
    w = TWOPI * f_trap
    return e ** 2 * S_E / (4.0 * m * hbar * w)


# 基準点:室温の表面電極トラップ
D_REF, F_REF, S_REF = 50e-6, 1.0e6, 1.0e-11
print(f"基準となる表面トラップ: d = {D_REF * 1e6:.0f} um、"
      f"f = {F_REF / 1e6:.1f} MHz、S_E = {S_REF:.0e} V^2 m^-2 Hz^-1")
print(f"Ca+ ではこれは dn/dt = {heating_rate(S_REF, 40.078, F_REF):.1f} "
      f"量子/s に相当する。")
print()

print(f"{'d (um)':>9}{'S_E (V^2/m^2/Hz)':>19}{'dn/dt (quanta/s)':>19}"
      f"{'quanta in 100 us':>19}")
for d in [20e-6, 50e-6, 100e-6, 500e-6]:
    S = S_REF * (D_REF / d) ** 4
    r = heating_rate(S, 40.078, F_REF)
    print(f"{d * 1e6:>9.0f}{S:>19.3e}{r:>19.3f}{r * 100e-6:>19.3e}")
print()

print(f"{'f (MHz)':>9}{'S_E (1/f model)':>18}{'dn/dt (quanta/s)':>19}"
      f"{'relative':>11}")
base = None
for f in [0.5e6, 1.0e6, 2.0e6, 5.0e6]:
    S = S_REF * (F_REF / f)
    r = heating_rate(S, 40.078, f)
    base = base or r
    print(f"{f / 1e6:>9.1f}{S:>18.3e}{r:>19.3f}{r / base:>11.3f}")
print()

print(f"{'electrode T':>13}{'suppression':>13}{'S_E':>13}"
      f"{'dn/dt (quanta/s)':>19}{'quanta in 100 us':>19}")
for T, supp in [(300.0, 1.0), (77.0, 10.0), (10.0, 100.0), (4.0, 200.0)]:
    S = S_REF / supp
    r = heating_rate(S, 40.078, F_REF)
    print(f"{T:>11.0f} K{supp:>13.0f}x{S:>13.3e}{r:>19.4f}"
          f"{r * 100e-6:>19.3e}")
print()

print(f"{'ion':>7}{'m (u)':>9}{'dn/dt (quanta/s)':>19}{'relative':>11}")
ref = heating_rate(S_REF, 40.078, F_REF)
for name, m_u in [("Be+", 9.012), ("Ca+", 40.078), ("Yb+", 171.0)]:
    r = heating_rate(S_REF, m_u, F_REF)
    print(f"{name:>7}{m_u:>9.3f}{r:>19.3f}{r / ref:>11.3f}")
print()

# --- ゲート誤差のバジェット --------------------------------------------
print(f"{'trap':>34}{'dn/dt':>12}{'tau = 30 us':>14}{'100 us':>12}{'1 ms':>12}")
configs = [
    ("マクロ型, 300 K, d = 500 um", 500e-6, 1.0, 1.0e6),
    ("表面型, 300 K, d = 50 um", 50e-6, 1.0, 1.0e6),
    ("表面型, 10 K, d = 50 um", 50e-6, 100.0, 1.0e6),
    ("表面型, 10 K, d = 50 um, 3 MHz", 50e-6, 100.0, 3.0e6),
    ("表面型, 10 K, d = 30 um, 3 MHz", 30e-6, 100.0, 3.0e6),
]
for lab, d, supp, f in configs:
    S = S_REF * (D_REF / d) ** 4 * (F_REF / f) / supp
    r = heating_rate(S, 40.078, f)
    print(f"{lab:>34}{r:>12.3f}"
          + "".join(f"{r * t:>12.2e}" for t in [30e-6, 100e-6, 1e-3]))
print()

# --- 表面で何が起きているのか ------------------------------------------
print()
S_needed = S_REF * 1e-4 / (heating_rate(S_REF, 40.078, 1e6) * 100e-6)
print(f"1 MHz の Ca+ で100 usのゲート中の吸収量子数を1e-4にするには")
print(f"S_E = {S_needed:.3e} V^2 m^-2 Hz^-1、すなわち上の基準より "
      f"{S_REF / S_needed:.0f} 倍低い値が必要")
基準となる表面トラップ: d = 50 um、f = 1.0 MHz、S_E = 1e-11 V^2 m^-2 Hz^-1
Ca+ ではこれは dn/dt = 1455.3 量子/s に相当する。

   d (um)   S_E (V^2/m^2/Hz)   dn/dt (quanta/s)   quanta in 100 us
       20          3.906e-10          56847.277          5.685e+00
       50          1.000e-11           1455.290          1.455e-01
      100          6.250e-13             90.956          9.096e-03
      500          1.000e-15              0.146          1.455e-05

  f (MHz)   S_E (1/f model)   dn/dt (quanta/s)   relative
      0.5         2.000e-11           5821.161      1.000
      1.0         1.000e-11           1455.290      0.250
      2.0         5.000e-12            363.823      0.062
      5.0         2.000e-12             58.212      0.010

  electrode T  suppression          S_E   dn/dt (quanta/s)   quanta in 100 us
        300 K            1x    1.000e-11          1455.2903          1.455e-01
         77 K           10x    1.000e-12           145.5290          1.455e-02
         10 K          100x    1.000e-13            14.5529          1.455e-03
          4 K          200x    5.000e-14             7.2765          7.276e-04

    ion    m (u)   dn/dt (quanta/s)   relative
    Be+    9.012           6471.940      4.447
    Ca+   40.078           1455.290      1.000
    Yb+  171.000            341.083      0.234

                              trap       dn/dt   tau = 30 us      100 us        1 ms
           マクロ型, 300 K, d = 500 um       0.146    4.37e-06    1.46e-05    1.46e-04
             表面型, 300 K, d = 50 um    1455.290    4.37e-02    1.46e-01    1.46e+00
              表面型, 10 K, d = 50 um      14.553    4.37e-04    1.46e-03    1.46e-02
       表面型, 10 K, d = 50 um, 3 MHz       1.617    4.85e-05    1.62e-04    1.62e-03
       表面型, 10 K, d = 30 um, 3 MHz      12.477    3.74e-04    1.25e-03    1.25e-02


1 MHz の Ca+ で100 usのゲート中の吸収量子数を1e-4にするには
S_E = 6.871e-15 V^2 m^-2 Hz^-1、すなわち上の基準より 1455 倍低い値が必要

着目点。 基準点 — $d = 50\ \mu$m の室温表面トラップ — はCa$^+$を1 MHzで毎秒1455量子です。これは敵対的な数値です。イオンは700 $\mu$sごとに量子1個を得ており、ゲート時間と同程度です。

形状の走査が、微細化が無償でない理由を示します。$d = 500\ \mu$m から $d = 20\ \mu$m へ行くと加熱率は毎秒0.15から57000量子へ — 距離25分の1から $4\times10^5$ 倍、すなわち $25^4$ です。トラップを小さくするあらゆる論拠(電極を増やす、閉じ込めを強める、光学系と電子回路を集積する、輸送を速くする)が真っ向から $d^{-4}$ にぶつかります。これがこの方式の中心的な工学的緊張であり、しかも定量的です。

周波数の走査は無償の得点です。$S_E \sim 1/f$ で加熱率が明示的な $1/\omega$ をもつので、率は $1/f^2$ で落ちます。5 MHzのモードは0.5 MHzのそれより100倍ゆっくり暖まります。Code Example 5 のより深いサイドバンド冷却限界と、より小さい $\eta$ と合わせて、これがゲートが高周波の径方向モードを使う理由です。

温度の走査が、極低温イオントラップが存在する理由です。10 Kでの2桁が基準を毎秒1455から14.6量子にし、その代価は研究プログラムではなくクライオスタットです。

誤差バジェットの表がすべてをまとめます。ゲート誤差の目標 $10^{-4}$ には、ゲート中の吸収量子数がその程度である必要があり、表はどの構成がそれを供給できるかを述べています。50 $\mu$mの室温表面トラップは無理(100 $\mu$sで0.15量子を吸収、3桁多すぎ)、50 $\mu$mで3 MHzの極低温トラップならぎりぎり可能($1.6\times10^{-4}$)、その同じ極低温トラップを30 $\mu$mに縮めると8倍損してまた手が届かなくなります。これはスローガンではなく条件つきの言明として読んでください。目標が数十マイクロメートルの電極距離で $10^{-4}$ のゲート誤差であるなら、極低温運転は便宜ではなく要件です — 一方で、より控えめな誤差目標であれば、より大きな室温トラップ中の小規模レジスタは十分に動作し、この分野の仕事の多くはそこで行われてきました。無条件に成り立つのは緊張のほうです。トラップの微細化とゲート忠実度は逆方向に引き合います。

最後の行は材料の目標を数値として述べます。1 MHzのCa$^+$で100 $\mu$sのゲート中の吸収量子数を $10^{-4}$ にするには $S_E = 7\times10^{-15}\ \mathrm{V^2m^{-2}Hz^{-1}}$、室温表面トラップの基準の1455分の1が必要です。その一部は極低温から、一部は表面処理から、一部は形状とモードの選択から来ます。そのすべてが原理的に入手可能であり、そのどれも回路設計の改良では入手できません。

また同じパターン

第2章と第3章を並べてみます。

超伝導トランズモン イオントラップ
ハミルトニアン 厳密に既知、1%まで設計 周期表の性質
何が律速するか 誘電損失と準粒子による $T_1$ 電場ノイズによる $\bar{n}$ の増大
損失の住処 表面の3 nmのアモルファス酸化膜 電極の最後の数ナノメートル
スペクトルの署名 $1/f$、活性化エネルギーの広い分布 $1/f$、活性化エネルギーの広い分布
洗浄への応答 in-situ洗浄、封止、別の酸化膜 アルゴンイオンミリング、数日で劣化
冷却への応答 20 mKが必要、準粒子は非熱的 4〜10 Kで $10^2$ の抑制
形状のてこ $p \propto t/w$:電場を薄める $S_E \propto d^{-4}$:表面から離れる
名目上同一の素子の差 $T_1$ で1桁 $S_E$ で数桁

物理的に何ひとつ共通点のない2つの方式 — 一方は20ミリケルビンの1センチメートルのアルミニウム、他方は真空中の1個の原子 — が同じものに律速されています。数ナノメートルの表面における制御されていないダイナミクスであり、$1/f$ スペクトルをもち、洗浄と冷却に応答し、名目上同一の素子間で再現しない。これは偶然ではなく、本コースの主張です。量子ハードウェアのボトルネックは材料の問題であり、しかも同じ材料の問題が2度現れているのです。


演習

本章のコードを手元に置いて取り組んでください。各問のあとに解答があります。

演習1: トラップを設計する

$^{171}$Yb$^+$($m = 171$ u)用のトラップで、$r_0 = 200\ \mu$m を使い $q = 0.25$ で3 MHzの径方向secular周波数を与えたいとします。(a) 必要なRF周波数と振幅はいくらですか。(b) トラップ深さは何eVですか。(c) 同じ電圧のまま $^{40}$Ca$^+$ を入れると $q$ はどうなりますか。

解答

(a) \(q = 0.25\) では \(\beta \approx q/\sqrt{2} = 0.177\)(Code Example 1 の厳密値は0.178)なので \(\Omega = 2\omega_\mathrm{sec}/\beta = 2\times2\pi\times3\times10^6/0.178\)、すなわち \(\Omega/2\pi = 33.7\) MHz です。次に \(V = q m r_0^2\Omega^2/2e\):\(m = 171\times1.66\times10^{-27} = 2.84\times10^{-25}\) kg、\(r_0^2 = 4\times10^{-8}\) m\(^2\)、\(\Omega = 2.12\times10^8\) s\(^{-1}\) から \(V = 0.25\times2.84\times10^{-25}\times4\times10^{-8}\times4.49\times10^{16}/(2\times1.602\times10^{-19}) = 398\) V。

(b) \(U = \frac{1}{2}m\omega_\mathrm{sec}^2 r_0^2 = 0.5\times2.84\times10^{-25}\times(1.885\times10^7)^2\times4\times10^{-8} = 2.02\times10^{-18}\) J = 12.6 eV。

(c) \(q \propto 1/m\) なので同じ電圧では Ca\(^+\) の \(q = 0.25\times171/40 = 1.07\) となり、安定境界0.908を超えます — 軽いイオンはまったく捕捉されません。これは利用されています。質量選択的な捕捉は望まない種を追い出す方法であり、2種を同時に入れるトラップではパラメータを両方について検査しなければならない理由でもあります。

演習2: マイクロモーションを補償する

Code Example 2 の $q = 0.3$ のCa$^+$のモデルを使って、(a) サイドバンド/キャリア比をちょうど1%にする迷い電場はいくらですか。(b) 補償電極が電場を0.2 V/mまで打ち消せるとき、残る比はいくらですか。(c) 答えがレーザー波長に依存するのはなぜで、どちらの量子ビット型がより寛容ですか。

解答

(a) 表が挟んでいます。10 V/m で \(7.2\times10^{-3}\) です。比は変調指数が小さいとき \((J_1/J_0)^2 \approx (\mathrm{mi}/2)^2\) で、mi は \(E_\mathrm{stray}\) に比例するので、10 V/m からのスケーリングで \(\sqrt{0.01/0.00718} = 1.18\) 倍、すなわち約11.8 V/m です。

(b) 比は \(E^2\) でスケールします:\((0.2/11.8)^2\times0.01 = 2.9\times10^{-6}\)。V/mの端数まで補償すればマイクロモーションサイドバンドは無視できるようになり、まさにそれがこの手続きに手間をかける価値がある理由です。

(c) 変調指数は \(k x_\mathrm{micro}\) なので \(1/\lambda\) に比例します — Raman駆動の超微細量子ビットならの波数ベクトルに比例します。マイクロ波駆動の超微細遷移は \(k\) が4桁小さく、マイクロモーションに本質的に免疫です。729 nmの光学量子ビットはそうではありません。その免疫は超微細量子ビットの過小評価されている利点の1つであり、ゲートにRaman対を使った瞬間に失われます。

演習3: 鎖はいつ座屈するか

(a) Code Example 3 を使って、$\omega_r/\omega_z = 10$ での最大の直線鎖はいくらですか。(b) $f_z = 200$ kHz で20イオンの鎖が欲しいとします。直線性が要求する径方向周波数はいくらで、その結果の径方向モード間隔はいくらですか。(c) それはどんなゲート時間を含意し、1本の鎖のスケーリングについて何を語りますか。

解答

(a) 印字された表の比8(\(N = 18\))と12(\(N = 30\))のあいだを補間すると、比10では約 \(N = 24\) です。\(\mathrm{ratio} = 10\) でループを走らせると確認できます。

(b) 基準 \(\omega_r/\omega_z \gtrsim 0.73 N^{0.86}\) は \(N = 20\) で9.6 を与えるので \(f_r \gtrsim 1.9\) MHz、2 MHzを取ります。最小の径方向間隔は \(\approx \omega_z^2/2\omega_r = (0.2)^2/(2\times2) = 0.01\) MHz = 10 kHz です。

(c) ゲートの離調は10 kHzの十分内側でなければならないので \(\tau \gtrsim 1/10\) kHz \(= 100\ \mu\)s、実際にはループ閉条件 \(\delta\tau = 2\pi K\) と傍観者モードを勘定に入れるとその数倍になります。一方200 kHzでの加熱率は1 MHzでの25倍悪く(Code Example 7)、ゲートは10倍長いので、吸収量子数は2桁以上増えます。弱い軸方向閉じ込めは鎖を長くし、同時に3つの意味で物理を悪くします。これが短い鎖+輸送という選択の定量的な論拠です。

演習4: 冷却のバジェット

(a) Code Example 5 を使って、1 MHzのトラップで $\bar{n} = 0.01$ に到達するには実効線幅はいくら必要ですか。(b) そのときの $R_c$ でDoppler限界からの冷却にどれだけかかりますか。(c) 自然線幅を使うのではなく狭い遷移をクエンチする(意図的に広げる)のが標準的な技法である理由は何ですか。

解答

(a) \(\bar{n}_\mathrm{ss} = (\Gamma/4\omega)^2 = 0.01\) には \(\Gamma = 0.4\omega\)、すなわち \(f = 1\) MHz で \(\Gamma/2\pi = 400\) kHz が必要です。

(b) 印字された式 \(R_c = 4\eta^2\Omega^2/\Gamma\) に \(\Omega/2\pi = 50\) kHz、\(\eta = 0.097\)(1 MHzでの値、Code Example 4)、\(\Gamma/2\pi = 400\) kHz を入れると \(R_c = 4\times0.0094\times(3.14\times10^5)^2/(2.51\times10^6) = 1.48\times10^3\) s\(^{-1}\)。Code Example 5 の時間発展から接近には \(10/R_c\) 程度、すなわち約7 msかかります。

(c) 2つの要求が逆を向いているからです。狭い線は深い限界 \((\Gamma/4\omega)^2\) を与えますが、\(\Omega\) 固定では \(R_c \propto \Gamma^{-1}\) で小さくなり、さらに重要なことにサイクル率が励起状態をどれだけ速く再利用できるかで縛られ、1秒の準安定状態では絶望的です。補助レーザーによるクエンチは \(\Gamma_\mathrm{eff}\) を深さと速さの積が最適になる値、典型的には数百kHzに設定できます。これは可変パラメータであり、自然線幅はそうではありません。

演習5: 2ループのMSゲート

同じゲート時間 $\tau = 100\ \mu$s で $K = 2$ としてCode Example 6 を繰り返します。(a) $\delta$ と $g$ はいくらですか。(b) Lamb-Dicke極限でBell忠実度がなお1であることを数値的に確認してください。(c) $K > 1$ の実務上の論拠と、それに反対する論拠は何ですか。

解答

(a) \(\delta = K/\tau = 20\) kHz、\(g = 2\pi\delta/(4\sqrt{K}) = 2\pi\times20\,\mathrm{kHz}/5.657 = 2\pi\times3.54\) kHz、すなわち \(\eta\Omega/2\pi = 7.07\) kHz、\(\Omega/2\pi = 103\) kHz です。スクリプトで K = 2 とすればすべて自動で行われます。

(b) 9桁で1です。解析解の \(2\Phi\) はやはり \(-\pi/4\)、\(|\alpha(\tau)|\) はやはり計算機精度でゼロ、スピンブロックはやはりユニタリです。導出のどこにも \(K = 1\) を特別扱いした箇所はありませんでした。

(c) 賛成:\(\delta\) が大きいほど望まないモードから遠いので、傍観者モードへの残留結合が抑えられ、それらのモードの位相空間ループも閉に近づきます。またモード周波数の安定性への要求も緩みます。反対:\(\tau\) 固定では \(g\) が \(\sqrt{K}\) で増えるのでレーザー強度が \(K\) に線形に増え、非共鳴キャリア励起、光子散乱、校正すべきACStarkシフトが増えます。\(K = 1\) か2が通常の妥協点です。

演習6: 加熱とイオンの選択

Code Example 7 を使って、(a) 同じ $S_E$ と周波数で Be$^+$ が Ca$^+$ の4.4倍速く暖まるのはなぜですか。(b) Be$^+$ は $\eta = 0.55$(Ca$^+$ は0.068)でもあるので、ゲートをはるかに速くできます。吸収量子数のバジェットではどちらの効果が勝ちますか。(c) これは種の選択について何を語りますか。

解答

(a) \(d\bar{n}/dt \propto 1/m\) であり、\(40.078/9.012 = 4.45\) です。軽いイオンは同じ電場ノイズで揺らしやすいのです。

(b) ゲート条件は \(\eta\Omega = \delta/2\sqrt{K}\)、\(\tau = 2\pi K/\delta\) なので、レーザー強度固定では \(\eta\) が大きいほど \(\delta\) を大きく取れ、したがって \(\tau\) が短くなり、\(\Omega\) 固定でおおよそ \(\tau \propto 1/\eta\) です。Be\(^+\) はゲート時間で8倍得をし、加熱率で4.4倍損をするので、吸収量子数 \(\dot{\bar{n}}\tau\) は約2倍改善します。軽いイオンの勝ちですが — ぎりぎりで、その余裕は利用できるレーザー出力に完全に依存します。

(c) 種の選択が支配項をもたない多方面の妥協であるということです。質量は加熱率と \(\eta\) を決め、準位構造は超微細か光学かを、したがってコヒーレンスの機構を決め、波長はレーザーが市販品か英雄的努力かを決め(Be\(^+\) は313 nmで周波数混合が必要、Ca\(^+\) と Sr\(^+\) は便利なダイオード波長の近くにあります)、そして2種を同時捕捉して同情冷却と回路中途の読み出しを行うなら質量比が再び効きます。研究グループごとに違う項を最適化してきたので、この分野は1つのイオンに収束していないのです。


まとめ

要点

1. Paulトラップは鞍点が回転するから働く

2. マイクロモーションは形を変えた表面の問題である

3. イオン結晶はバスを供給し、速度を縛る

4. 冷却は無関係な限界をもつ2段階である

5. Mølmer-Sørensenゲートは幾何位相であり、厳密である

6. 異常加熱がボトルネックであり、それは表面である

実務上の含意

次章へ

第4章は原子を残して電荷を捨てます。光ツイーザーに保持された中性原子は任意の2次元・3次元ジオメトリに配置でき、ショット間に再配置できるので、3.5節の直線鎖の制約は完全に消えます。相互作用は共有フォノンではなくRydberg封鎖 — 1個の励起原子が隣接原子の励起そのものを禁じるほど強い双極子-双極子結合 — であり、Rydberg状態の物性の $n^{11}$ スケーリングは計算するまで極端に見えるでしょう。

材料の物語も性格を変えます。中性原子の近くに電極はないので異常加熱はありません。限界は原子損失、Rydberg状態の有限寿命、そして光学系から来ます。それを表面の問題からの脱出とみなすのか、単なる移転とみなすのかは、第4章を通して抱えておき、第5章で決着させるべき問いです。

← 第2章: 超伝導量子ビット 第4章: 中性原子 →

免責事項