第1章: 何が良い量子ビットを作るのか

DiVincenzo基準、比較の軸、そしてデコヒーレンスの共通言語

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

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

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

量子コンピュータは一つの装置であり、本章はそれがどのような装置でなければならないかを扱います。姉妹コース量子コンピューティング入門ではハードウェアを一組のパラメータ — コヒーレンス時間、エラー率 — として扱い、その抽象化のままかなり遠くまで進みました。いまその抽象化を開く必要があります。それらのパラメータの背後には、ハミルトニアンと環境と作製履歴をもつ物理系があり、この3つ目こそ、計算機科学者にはできない貢献を材料研究者ができる場所です。

本章の役割は2つあります。第一は「良い」とは何かを確立することで、これは見かけよりも難しい仕事です。量子ビット方式を評価する軸は少なくとも6本あり、それらは物理的に互いに独立ではなく、そしてすべての軸で最良となる方式は存在しません。第二は、$T_1$、$T_2$、$T_2^\ast$、ノイズパワースペクトル密度、Ramsey測定とエコー、という共通言語を固定することです。第2章から第5章はこれを再定義せずに用います。本章の7つのコード例のうち5つ — 1.6節にまとめたもの — で小さなBloch方程式の実験室を構築します。標準的なコヒーレンス測定を第一原理から再現するものであり、これらはどの方式でも必ず報告される測定であると同時に、正しく読めば材料の特性評価データでもあるからです。

学習目標

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

記法と単位

本コースは以下の3つの約束を例外なく守ります。後で戸惑わないよう、いま一度目を通しておく価値があります。

換算単位。 ハミルトニアンは $\hbar = 1$ で書き、エネルギーと角周波数を同一の対象として扱います。物理的なエネルギーを意味する箇所では係数を明示的に復元します。たとえば $\omega_q/2\pi = 5$ GHz のトランズモン遷移のエネルギーは $\hbar\omega_q = h \times 5\ \mathrm{GHz} = 20.7\ \mu\mathrm{eV}$ です。

角周波数と周期周波数。 $\Omega$(Rabi周波数)や $\Delta$(離調)などの記号は rad/s 単位の周波数です。数値として引用する際は $\Omega/2\pi$ や $\Delta/2\pi$ のように Hz 単位の周期周波数を用います。コード中の $2\pi$ をすべて明示しているのはこのためです。この2つを混同することは本分野で最も多い数値的誤りであり、6.28倍の誤差を生みながら一見もっともらしいグラフを与えます。

分野ごとの慣用エネルギー単位。 超伝導回路は GHz、イオントラップと中性原子は MHz、スピン量子ビットは MHz または印加磁場のテスラで表します。Code Example 1 で換算表を一度示し、以降は当該分野の文献が用いる単位に従います。

量子ビットの順序とゲート記号は姉妹コースと完全に同一です。量子ビット0が左端かつ最上位ビットであり、$X, Y, Z, H$、CNOT はそちらと同じ意味です。第1章では1量子ビットの言語しか必要としませんが、第2章から第5章はこの規約を前提とします。


1.1 なぜこれほど多種類の量子ビットがあるのか

それ自身と矛盾する要求

量子ビットが満たすべき条件を、物理学者が書く順に書き出してみます。

  1. 2つの準位が操作可能であること。すなわち系の他のすべての準位から分離でき、制御場で到達できること。
  2. その2準位間の相対位相が、計算を終えるまで生き残ること。
  3. 制御場がその生存時間よりずっと短い時間で状態を $\pi$ だけ回転できること。

条件2と3は逆方向に引き合います。印加場に強く応答する系は、頼んでもいない揺動場にも強く応答します。同一の結合定数が両者を支配するという事実の形式的な表現が、揺動散逸定理です。完全に孤立した2準位系は無限のコヒーレンス時間をもち、そして操作できません。これは設計で回避すべき技術的な不便ではなく、本分野の中心的な設計上の緊張であり、あらゆる方式はその緊張に対する特定の妥協です。

部分的な逃げ道はあり、それがこの緊張が量子計算を単に禁止してしまわない理由です。結合を周波数選択的にすればよいのです。$\omega_q$ の場には強く結合し他の周波数の場には弱く結合する量子ビットは、環境が $\omega_q$ にスペクトル重みをほとんど持たないかぎり、制御可能かつ保護されています。本コースのすべての方式がこれを利用しており、第1.4節でそれを定量化します。同時にこれは、環境ノイズスペクトルの大きさだけでなくそのこそが問題になる理由でもあります。

2つの系統の答え

提案は2つの系統に明快に分かれ、その区分は本質的に「誰がハミルトニアンを選んだか」の違いです。

自然の量子ビットは、自然がすでに提供している準位を使います。原子・イオンの超微細遷移や光学遷移、核スピンや電子スピン、光子の偏光です。その長所は出自から直接に従います。宇宙のあらゆる $^{171}\mathrm{Yb}^{+}$ イオンは15桁まで同一の超微細分裂をもちますから、作製由来の乱れはなく、原理的に素子ごとの校正問題も存在しません。これらの準位はほとんど何にも弱く結合しており、だからこそ秒を超えるコヒーレンス時間が達成できます。短所も同様に従います。ハミルトニアンは調整できず、量子ビットは超高真空中に捕獲・保持しなければならず、制御装置 — レーザー、光学系、磁気シールド — は量子ビット数が増えても小さくなりません。

人工の量子ビットは、固体の中に人工的な2準位系を作ります。超伝導回路、ゲート形成量子ドット、結晶中の欠陥中心です。ここではハミルトニアンが設計パラメータであり、自然ではなくリソグラフィと膜厚が決めます。結合を大きく、ゲートを速くでき、素子全体をすでに存在するウェハスケールのプロセスで作製できます。その代償として、量子ビットは母材を引き継ぎます。2準位欠陥に満ちたアモルファス酸化膜、ダングリングボンドをもつ界面、核スピンをもつ基板から数十ナノメートルの距離に置かれ、そのコヒーレンス時間は設計中の何かではなくそれらによって決まります。作製はパラメータのばらつきも導入します。名目上同一の素子は同一ではなく、その分布は裾を引きます。

自然の量子ビット 人工の量子ビット
ハミルトニアン 自然が決める 設計が決める
再現性 定義上、厳密に同一 プロセス制御に律速される
コヒーレンス 長い。装置側が限界を決める 短い。母材が限界を決める
ゲート速度 遅い(kHz-MHz) 速い(MHz-GHz)
集積 光学系と真空は小型化しない ウェハスケールのプロセスが既にある
主要な研究課題 制御の複雑さ 材料

この表が、本コースがこの形をとる理由です。右列では律速要因が材料問題であり、しかも難しい材料問題です。ミリケルビン・ギガヘルツにおけるアモルファス誘電体の損失、界面欠陥密度、同位体純度、超伝導体-半導体エピタキシー。左列では律速要因は別のところにありますが、そこでもトラップ表面、真空、光学コーティングは材料です。

簡単な系譜

方式の乱立は物理的でもあり歴史的でもあります。本コースの各方式は、それ以前の方式がもつ特定の困難への解答として提案されました。

提案 物理的なアイデア
1995 Cirac と Zoller: イオントラップゲート イオン結晶の共有振動モードを内部状態間のバスとして使う
1997 Gershenfeld, Chuang, Cory: NMRアンサンブル 分子中の核スピンはすでに長い $T_2$ をもつ。アンサンブルを扱う代償を払ってそれを使う
1997-2003 Kitaev: トポロジカル量子計算 情報を非局所的な自由度に蓄え、局所ノイズが読めず壊せないようにする
1998 Loss と DiVincenzo: 量子ドットスピン 半導体中の電子スピンと交換相互作用によるゲート。既存のリソグラフィと整合する
1999 中村・Pashkin・蔡: 超伝導電荷量子ビット ジョセフソン接合により巨視的な回路が人工原子として振る舞う
2000 Jaksch ら: Rydberg封鎖 1個の原子を巨大双極子のRydberg状態に励起すると、近傍の原子の励起が禁止される
2001 Knill, Laflamme, Milburn: 線形光学 光子は相互作用しないが、測定と事後選択が相互作用の代わりになる
2007 Koch ら: トランズモン 電荷量子ビットを大きな容量で分流し、非調和性の弱さを受け入れて電荷ノイズに不感にする

このうち2つの項目は立ち止まる価値があります。それが本分野全体のパターンだからです。トランズモンが存在するのは、電荷ノイズ — 酸化膜中や表面上の可動電荷という材料問題 — がその前身を殺したからで、対処法は量子ビット周波数が電荷に依存しないようハミルトニアンを再設計することでした。トポロジカル量子ビットが存在するのは、あらゆる局所ノイズが材料問題だからで、提案されている対処法は、どの局所演算子も届かない場所に情報を蓄えることです。どちらもデコヒーレンスへの応答であり、方式とは特定のノイズチャネルへの解答として理解するのが最もよいことを示しています。

姉妹コースが空けておいた場所

量子コンピューティング入門は、自らがハードウェアの講座ではないと明言しています。超伝導量子ビット、イオントラップ、中性原子はそちらでは誤差特性が関係する箇所にのみ登場します。本コースはまさにその空白を埋めるものであり、両者の対応は直接的です。

アルゴリズム編では 本コースでは
量子ビットは規格化された複素ベクトル 実在するスペクトルの2準位部分空間。そこからの漏れ(リーケージ)がある
ゲートはユニタリ行列 整形された制御パルス。校正誤差と有限のパルス長をもつ
$T_1$ と $T_2$ はノイズモデルの入力パラメータ 欠陥・界面・同位体が決める従属変数
接続性はグラフ コンデンサのネットワーク、共有フォノンモード、あるいは封鎖半径
誤り訂正にはオーバーヘッドがある そのオーバーヘッドこそが材料問題を切迫させる

2つのコースを併せて読むと、有用な問いが「量子コンピュータは何ができるか」から「それを実現するには物理的に何が改善されなければならないか」へと変わります。


1.2 DiVincenzo基準

2000年、David DiVincenzo は物理系が量子計算を実装するために備えるべき能力の最小限のリストを書き下しました。いまなお標準的なチェックリストであり続けているのは、それが技術ではなく物理の水準で述べられているからです。5つが計算に関わり、さらに2つが通信に関わります。

5つの基準

1. スケーラブルで、よく特性づけられた量子ビットからなる物理系。 一文の中に2つの要求が隠れています。よく特性づけられているとは、ハミルトニアンが既知であることです。準位間隔、制御場との結合、他の量子ビットとの結合、そして決定的に重要な、計算部分空間の外側の準位との結合。トランズモンは2準位系ではありません。第3準位が数百MHzしか離れていない弱い非調和振動子であり、速すぎるパルスはその準位を占有してしまいます。このリーケージはそもそも量子ビットの誤りとして記述できず、だからこそ非調和性が第2章で第一級の設計パラメータになります。スケーラブルとは、量子ビットを増やすたびに新しい発明を必要としないことであり、実際の制約は量子ビットそのものよりも制御チャネル数、配線の熱流入、レーザー出力にあります。

2. 単純な基準状態へ初期化する能力。 通常は $|00\ldots0\rangle$ です。物理的な経路は2つあります。熱的初期化は $k_BT \ll \hbar\omega_q$ であれば働き、5 GHz のマイクロ波量子ビットでは希釈冷凍機を意味します。要求は単に「冷たい」ことではなく $T \lesssim 30$ mK であり、Code Example 1 がそれを示します。散逸的初期化 — 光ポンピングや能動リセット — は温度によらず働きます。平衡を待つのではなく系を暗状態へ駆動するからです。だからこそ 12.6 GHz の超微細量子ビットをもつイオンが室温の真空槽で動作し、より低い周波数のトランズモンがミリケルビンを必要とするのです。誤り訂正はこの要求をかなり厳しくします。新鮮で高忠実度の補助量子ビットが連続的に消費されるため、初期化は「可能」であるだけでなく速く、かつ再現的でなければなりません。

3. 関連するデコヒーレンス時間が、ゲート操作時間よりずっと長いこと。 効いている語は関連するです。重要なのは無次元比 $T_2/t_\mathrm{gate}$ であって、どちらか一方の数値ではありません。$T_2 = 1$ s で 100 $\mu$s のゲートをもつ方式と、$T_2 = 100\ \mu$s で 10 ns のゲートをもつ方式は、数倍の違いで比較可能な予算をもちます。第1.3節でこれを計算します。ハードウェアに関する主張から取り出すべき最も有用な単一の数値です。

4. 普遍的な量子ゲート集合。 $SU(2^n)$ を生成する任意の集合、実際には任意の1量子ビット回転と1つのエンタングリング2量子ビットゲートです。難しいのは2量子ビットゲートの側で、その物理機構は方式間で完全に異なります。回路では直接の容量結合、イオントラップでは共有運動モード、中性原子ではRydberg封鎖、スピンでは交換相互作用です。第2章から第5章はそれぞれ1節をこれに充てます。ゲート機構が接続グラフとゲート時間を同時に決めるからです。

5. 量子ビットごとの測定能力。 指定した1つの量子ビットを、他を乱すことなく高い忠実度で読み出すこと。実務では忠実度以外に2つの性質が効きます。速さ: 誤り訂正はコヒーレンス時間内に多数の測定ラウンドを要求しますから、$T_2$ の大きな割合を消費する読み出しは忠実度が完璧でも役に立ちません。非破壊性: 測定でトラップから失われる原子や、再冷却が必要なイオンは、稼働率を支配する再装填サイクルを課します。

さらに2つ

6. 静止量子ビットと飛行量子ビットを相互変換する能力。 ネットワークはモジュール間で量子情報を運ぶ必要があり、その担い手はほぼ常に光子です。物理的な課題は、GHz帯の固体量子ビットと数百THzの光子との間のインピーダンス整合です。

7. 飛行量子ビットを地点間で忠実に伝送する能力。 ファイバ中の損失、自由空間の損失、そして未知の量子状態を増幅できないこと。

この2つがモジュラーアーキテクチャを可能にするものであり、計算機として劣る場合でも光子方式が依然として興味深い理由です。光子を相互作用しにくくしているのと同じ物理が、光子を移動に卓越させています。

この基準が語っていないこと

このリストは必要条件であって十分条件ではなく、これを採点表として読むのは3つの意味で誤りです。

それでもなお、どの基準も材料科学者が影響を与えうる物理量に対応しており、それを明示的に表にする価値があります。

基準 物理的内容 材料科学が入る場所
1. 特性づけとスケーラビリティ 既知のハミルトニアン、小さなパラメータばらつき、制御されたリーケージ 膜厚と接合面積の均一性、同位体・化学的純度
2. 初期化 $k_BT \ll \hbar\omega_q$、あるいは速い散逸チャネル 熱アンカー、残留熱流入、トラップ表面の清浄度
3. $T_2 \gg t_\mathrm{gate}$ 直流付近および $\omega_q$ におけるノイズパワースペクトル密度 欠陥密度、誘電損失、核スピン浴、準粒子
4. 普遍ゲート 制御可能な相互作用と非調和性 接合材料、基板誘電率、表面電場ノイズ
5. 測定 測定器への強く速く量子ビット選択的な結合 共振器のQ値、検出器効率、増幅器ノイズ
6-7. 飛行量子ビット コヒーレントな変換と低損失伝送 電気光学・圧電材料、ファイバとコーティングの損失

1.3 6本の軸、そしてそれを1本に潰せない理由

DiVincenzo基準は合否判定です。すべて合格する方式同士を比較するには連続的な軸が必要で、本コースは6本を用います。

測るもの 限界の物理的起源
コヒーレンス時間 $T_1$、$T_2$、$T_2^\ast$ 関連周波数における環境のノイズパワースペクトル密度
ゲート忠実度と速度 ゲートあたり誤り、$t_\mathrm{gate}$ 非調和性やモード構造が速度の上限を決め、校正とデコヒーレンスが誤りを決める
接続性 どのペアが直接相互作用できるか 2量子ビットゲートの物理機構
再現性と歩留まり パラメータのばらつき、動作素子の割合 作製由来の乱れ、あるいは自然の量子ビットではその不在
動作温度 必要な $T$ と冷却の予算 $k_BT$ に対する量子ビットのエネルギースケール
スケーラビリティ $n$ を増やしたときに壊れるもの 制御チャネル、配線の熱流入、レーザー出力、占有面積

軸が互いに争う理由

これらは独立ではなく、軸の間の結びつきは偶然ではなく物理的です。

Code Example 1: エネルギースケール、単位、そして誰が希釈冷凍機を必要とするか

比較の前に換算表を1つ。量子ビットの遷移周波数がそのエネルギーと等価温度を決め、したがって冷却によって初期化できるかどうかを決めます。

import numpy as np

h = 6.62607015e-34      # プランク定数, J s
kB = 1.380649e-23       # ボルツマン定数, J/K
eV = 1.602176634e-19    # 電気素量, C(1 eV のジュール換算値)


def thermal_population(f_hz: float, T_kelvin: float) -> float:
    """熱平衡状態にある2準位系の励起状態占有数。"""
    x = h * f_hz / (kB * T_kelvin)
    return np.exp(-x) / (1.0 + np.exp(-x))


# 単位換算表は一度だけ示し、本コース全体でこれを用います。
print("遷移周波数 h f の単位換算:")
print(f"{'f':>10}{'energy':>14}{'h f / kB':>14}{'wavenumber':>14}")
print("-" * 52)
for label, f in [("1 MHz", 1e6), ("1 GHz", 1e9), ("1 THz", 1e12),
                 ("500 THz", 5e14)]:
    E_ueV = h * f / eV * 1e6
    T_K = h * f / kB
    nu_cm = f / 2.99792458e10
    print(f"{label:>10}{E_ueV:>11.4g} ueV{T_K:>11.4g} K{nu_cm:>11.4g} 1/cm")

# 代表的な量子ビット遷移。ここに挙げた周波数は各準位構造の*物理*
# (超微細分裂、ジョセフソン回路のプラズマ周波数、Zeeman分裂)であり、
# 装置のスペックではありません。
qubits = [
    ("Transmon, microwave",        5.0e9,   "thermal"),
    ("Electron spin, B = 1 T",     28.0e9,  "thermal"),
    ("NV centre, zero field",      2.87e9,  "optical pumping"),
    ("Rb-87 hyperfine",            6.835e9, "optical pumping"),
    ("Yb-171+ hyperfine",          12.64e9, "optical pumping"),
    ("Ca-40+ optical, 729 nm",     4.11e14, "optical pumping"),
]

print("\n2準位系の熱励起状態占有数:")
header = (f"{'qubit transition':<26}{'f (GHz)':>12}{'hf/kB (K)':>11}"
          f"{'300 K':>10}{'4 K':>10}{'100 mK':>10}{'10 mK':>10}")
print(header)
print("-" * len(header))
for name, f, _ in qubits:
    pops = [thermal_population(f, T) for T in (300.0, 4.0, 0.1, 0.010)]
    print(f"{name:<26}{f/1e9:>12.4g}{h*f/kB:>11.4g}"
          + "".join(f"{p:>10.2e}" for p in pops))

print("\n各方式が実際に初期状態へ到達する経路:")
for name, f, route in qubits:
    if route == "thermal":
        T_needed = h * f / (kB * np.log(1.0 / 1e-3 - 1.0))   # P_exc = 1e-3
        print(f"  {name:<26} 熱平衡        P_exc < 1e-3 には "
              f"T < {T_needed*1e3:6.1f} mK が必要")
    else:
        print(f"  {name:<26} 光ポンピング  装置の温度は無関係")

# 同じエネルギースケールを駆動側から見る: 1 Rabi周期に要する時間。
print("\n駆動強度とパルス長(Omega = 2 pi f_Rabi):")
print(f"{'f_Rabi':>10}{'Omega (rad/s)':>16}{'pi pulse':>14}")
print("-" * 40)
for label, f_rabi in [("1 kHz", 1e3), ("100 kHz", 1e5),
                      ("1 MHz", 1e6), ("50 MHz", 5e7)]:
    Omega = 2 * np.pi * f_rabi
    print(f"{label:>10}{Omega:>16.4g}{np.pi/Omega*1e6:>11.4g} us")
遷移周波数 h f の単位換算:
         f        energy      h f / kB    wavenumber
----------------------------------------------------
     1 MHz   0.004136 ueV  4.799e-05 K  3.336e-05 1/cm
     1 GHz      4.136 ueV    0.04799 K    0.03336 1/cm
     1 THz       4136 ueV      47.99 K      33.36 1/cm
   500 THz  2.068e+06 ueV    2.4e+04 K  1.668e+04 1/cm

2準位系の熱励起状態占有数:
qubit transition               f (GHz)  hf/kB (K)     300 K       4 K    100 mK     10 mK
-----------------------------------------------------------------------------------------
Transmon, microwave                  5       0.24  5.00e-01  4.85e-01  8.32e-02  3.79e-11
Electron spin, B = 1 T              28      1.344  4.99e-01  4.17e-01  1.46e-06  4.37e-59
NV centre, zero field             2.87     0.1377  5.00e-01  4.91e-01  2.01e-01  1.04e-06
Rb-87 hyperfine                  6.835      0.328  5.00e-01  4.80e-01  3.63e-02  5.67e-15
Yb-171+ hyperfine                12.64     0.6066  4.99e-01  4.62e-01  2.31e-03  4.51e-27
Ca-40+ optical, 729 nm        4.11e+05  1.972e+04  2.79e-29  0.00e+00  0.00e+00  0.00e+00

各方式が実際に初期状態へ到達する経路:
  Transmon, microwave        熱平衡        P_exc < 1e-3 には T <   34.7 mK が必要
  Electron spin, B = 1 T     熱平衡        P_exc < 1e-3 には T <  194.6 mK が必要
  NV centre, zero field      光ポンピング  装置の温度は無関係
  Rb-87 hyperfine            光ポンピング  装置の温度は無関係
  Yb-171+ hyperfine          光ポンピング  装置の温度は無関係
  Ca-40+ optical, 729 nm     光ポンピング  装置の温度は無関係

駆動強度とパルス長(Omega = 2 pi f_Rabi):
    f_Rabi   Omega (rad/s)      pi pulse
----------------------------------------
     1 kHz            6283        500 us
   100 kHz       6.283e+05          5 us
     1 MHz       6.283e+06        0.5 us
    50 MHz       3.142e+08       0.01 us

注目すべき点。 最初の表が覚えるべき換算です。1 GHz は 4.14 $\mu$eV に、そして 48 mK に対応します。この最後の数値こそ、超伝導量子計算が希釈冷凍機の中で行われる理由そのものです。5 GHz の量子ビットは $\hbar\omega_q/k_B = 0.24$ K ですから、4 K のヘリウム浴では48%が励起しており、残留励起が $10^{-3}$ を下回るのは約35 mK 以下でのみです。光学遷移の行は逆の極端です。$\hbar\omega/k_B \approx 2\times10^4$ K で、室温での熱占有数は $10^{-29}$、4 K 以下では $e^{-5000}$ が倍精度でアンダーフローするため厳密な0として表示されます。これは直すべき丸めの副作用ではありません。光学量子ビットは実験室が経験するどんな温度でも絶対的に基底状態にある、という主張です。

2番目のブロックが重要な留保です。12.6 GHz の $^{171}\mathrm{Yb}^{+}$ 超微細量子ビットは 4 K で46%が励起しているはずですが、イオントラップは室温で動作します。熱平衡はそもそも初期化の手段ではないのです。光ポンピングがマイクロ秒でイオンを特定の超微細状態へ駆動し、装置の温度はどこにも入りません。DiVincenzoの基準2はある機構についての主張であって冷却についての主張ではなく、この2つの混同が「量子コンピュータはすべて冷凍機を必要とする」という広く見られる思い込みの出所です。

Code Example 2: 同じ4方式に対する3通りの順位

単一の性能指標が存在しないという主張は、断言ではなく実演に値します。以下の数値は各方式の物理が決める桁のスケール — トランズモンのゲートがナノ秒なのは非調和性が数百MHzだから、イオンのゲートが数十から数百マイクロ秒なのは数MHzの運動モードを介するから — であり、有効数字1桁で与えてあるのは、この表の読み方が記録や装置スペックに依存しないようにするためです。

import numpy as np

# 各方式の物理から決まる桁のスケールを、有効数字1桁で与えます。これらは
# 装置のスペックでも記録でもありません。トランズモンのゲート時間が短いのは
# 非調和性が数百MHzだからであり、イオンのゲート時間が長いのは数MHzの
# 運動モードを介するからです。以下の議論に効くのは桁だけです。
platforms = [
    # 名称,                 T2 (s), t_2q (s), 結合様式,          並列レーン数
    ("Superconducting",     1e-4,   5e-8,    "2D nearest-nbr",  50),
    ("Trapped ion",         1e0,    1e-4,    "all-to-all",       1),
    ("Neutral atom",        1e-1,   5e-7,    "reconfigurable",  10),
    ("Semiconductor spin",  1e-3,   1e-7,    "1D nearest-nbr",  50),
]

n_qubits = 100        # ルーティング見積りで仮定するレジスタ規模
n_gates = 1000        # 課題が要求する任意ペア間の2量子ビットゲート数


def routing_overhead(lattice: str, n: int) -> float:
    """SWAPネットワークに由来する、論理ゲート1個あたりの実ゲート数。

    SWAP 1個はCNOT 3個に相当します。2次元格子上で無作為な2点間の平均距離は
    sqrt(n)/2 程度、1次元鎖では n/3 程度です。全結合と再配置可能な結合では
    SWAPは一切不要です。
    """
    if lattice == "2D nearest-nbr":
        return 1.0 + 3.0 * np.sqrt(n) / 2.0
    if lattice == "1D nearest-nbr":
        return 1.0 + 3.0 * n / 3.0
    return 1.0


rows = []
for name, T2, t2q, lattice, lanes in platforms:
    budget = T2 / t2q                       # コヒーレンス時間に収まるゲート数
    ovh = routing_overhead(lattice, n_qubits)
    phys_gates = n_gates * ovh              # 実際に実行されるゲート数
    layers = phys_gates / lanes             # 並列レーンで実行した層数
    wall = layers * t2q                     # 回路の実時間
    rows.append(dict(name=name, T2=T2, t2q=t2q, lattice=lattice, lanes=lanes,
                     budget=budget, ovh=ovh, wall=wall, ratio=wall / T2))

header = (f"{'platform':<20}{'T2 (s)':>9}{'t_2q (s)':>10}{'T2/t_2q':>10}"
          f"  {'coupling':<16}{'SWAP ovh':>9}{'lanes':>7}")
print(header)
print("-" * len(header))
for r in rows:
    print(f"{r['name']:<20}{r['T2']:>9.0e}{r['t2q']:>10.0e}{r['budget']:>10.0f}"
          f"  {r['lattice']:<16}{r['ovh']:>9.0f}{r['lanes']:>7d}")


def ranking(key, reverse):
    order = sorted(rows, key=lambda r: r[key], reverse=reverse)
    return " > ".join(r["name"] for r in order)


print("\n同じ4方式に対する、いずれも妥当な3通りの単一指標ランキング:")
print(f"  by coherence time T2      : {ranking('T2', True)}")
print(f"  by gate speed 1/t_2q      : {ranking('t2q', False)}")
print(f"  by gate budget T2/t_2q    : {ranking('budget', True)}")

print(f"\n課題設定: {n_qubits} 量子ビットの任意ペアに対する2量子ビットゲート "
      f"{n_gates} 個。")
print(f"{'platform':<20}{'gates run':>12}{'layers':>10}{'wall clock':>14}"
      f"{'wall/T2':>11}")
print("-" * 67)
for r in rows:
    print(f"{r['name']:<20}{n_gates*r['ovh']:>12.0f}"
          f"{n_gates*r['ovh']/r['lanes']:>10.0f}"
          f"{r['wall']*1e3:>11.3g} ms{r['ratio']:>11.4f}")

print(f"\n  by circuit time / T2      : {ranking('ratio', False)}")

# ある改善が別の改善に比べてどれだけの価値をもつか。
print("\nトレードオフを明示する(超伝導の行について):")
base = rows[0]
variants = [("baseline", base["T2"], base["t2q"], base["ovh"]),
            ("10x longer T2", 10 * base["T2"], base["t2q"], base["ovh"]),
            ("10x faster gates", base["T2"], base["t2q"] / 10, base["ovh"]),
            ("all-to-all coupling", base["T2"], base["t2q"], 1.0)]
for label, T2, t2q, ovh in variants:
    wall = n_gates * ovh / base["lanes"] * t2q
    print(f"  {label:<22} wall/T2 = {wall/T2:8.4f}")
print(f"  SWAPネットワークの除去はここでは {base['ovh']:.0f} 倍の価値があり、"
      f"T2を1桁伸ばすより大きい。")
platform               T2 (s)  t_2q (s)   T2/t_2q  coupling         SWAP ovh  lanes
-----------------------------------------------------------------------------------
Superconducting         1e-04     5e-08      2000  2D nearest-nbr         16     50
Trapped ion             1e+00     1e-04     10000  all-to-all              1      1
Neutral atom            1e-01     5e-07    200000  reconfigurable          1     10
Semiconductor spin      1e-03     1e-07     10000  1D nearest-nbr        101     50

同じ4方式に対する、いずれも妥当な3通りの単一指標ランキング:
  by coherence time T2      : Trapped ion > Neutral atom > Semiconductor spin > Superconducting
  by gate speed 1/t_2q      : Superconducting > Semiconductor spin > Neutral atom > Trapped ion
  by gate budget T2/t_2q    : Neutral atom > Trapped ion > Semiconductor spin > Superconducting

課題設定: 100 量子ビットの任意ペアに対する2量子ビットゲート 1000 個。
platform               gates run    layers    wall clock    wall/T2
-------------------------------------------------------------------
Superconducting            16000       320      0.016 ms     0.1600
Trapped ion                 1000      1000        100 ms     0.1000
Neutral atom                1000       100       0.05 ms     0.0005
Semiconductor spin        101000      2020      0.202 ms     0.2020

  by circuit time / T2      : Neutral atom > Trapped ion > Superconducting > Semiconductor spin

トレードオフを明示する(超伝導の行について):
  baseline               wall/T2 =   0.1600
  10x longer T2          wall/T2 =   0.0160
  10x faster gates       wall/T2 =   0.0160
  all-to-all coupling    wall/T2 =   0.0100
  SWAPネットワークの除去はここでは 16 倍の価値があり、T2を1桁伸ばすより大きい。

注目すべき点。 同じ4方式に対する、いずれも妥当な3通りの順位づけがあり、3つはすべて食い違います。コヒーレンス時間ではイオンが超伝導回路に4桁の差で勝ちます。ゲート速度では順序がちょうど逆になります。無次元比 $T_2/t_{2q}$ — DiVincenzoの基準3が実際に要求している量 — では中性原子が先頭で、イオンとスピン量子ビットが並びます。したがって「方式Xが先行している」という形の主張はどの軸かを明示しなければならず、公開されている主張の多くはそれをしていません。

課題設定のブロックは、この比が捉えきれない軸を加えます。任意ペアを必要とするアルゴリズムを最近接格子に載せるにはSWAPゲートが必要で、ここで用いた平均距離の見積りでは100量子ビットの2次元格子で16倍、1次元鎖で101倍になります。このオーバーヘッドが、最後の順位が3番目の順位ともまた異なる理由です。最終ブロックはトレードオフを直接に定量化します。この課題では、SWAPネットワークを完全に除去することが16倍の価値をもち、コヒーレンス時間を1桁伸ばすより大きいのです。接続性は二次的な考慮事項ではありません。

この例には誠実な留保が2つ付きます。並列性のモデル(lanes 列)は粗く、しかも結果に実質的に効いています。超伝導チップはチップ上の異なる領域で多数のゲートを同時に走らせられ、単一のイオン鎖は原理的にそれができず、真の数値はアルゴリズムの構造に依存します。またこの見積りはゲート誤りを無視し、$T_2$ 以内に終われば実行可能とみなしています。第5章がこの表の誠実な版を、数値の代わりに物理的制約を入れて組み立てます。


1.4 デコヒーレンスの共通言語

第2章から第5章のすべての方式が、同じ3つの時間定数を報告します。ここで物理的実現に依存しない形で一度だけ定義し、以降のコースでは変更せずに用います。

情報を失う2つの経路

量子ビットの状態はBloch球上またはその内部の点であり、それが劣化する経路はちょうど2つです。成分 $z = \langle Z\rangle$ が熱平衡へ緩和する経路。これには環境とのエネルギー交換が必要です。成分 $x$ と $y$ が縮む経路。これには位相が不確かになることだけが必要で、エネルギーはどこへも行かなくてよいのです。

$T_1$、エネルギー緩和時間が前者を支配します。

$$ z(t) = z_\mathrm{eq} + \left(z(0) - z_\mathrm{eq}\right) e^{-t/T_1} $$

エネルギー交換ですから $1/T_1$ は黄金律の速度です。環境が演算子 $\hat{A}$ を通じてノイズパワースペクトル密度 $S_A$ で量子ビットに結合しているとすると、概略として

$$ \frac{1}{T_1} \;\propto\; \left| \langle 0 | \hat{A} | 1 \rangle \right|^2 S_A(\omega_q) $$

本質的なのは引数です。$T_1$ は量子ビット周波数におけるノイズを測っています。マイクロ波量子ビットではGHz帯です。kHz で静かな材料が 5 GHz では騒がしいことがあり、$T_1$ はそれを見つける測定です。

$T_2$、コヒーレンス時間が後者を支配します。2つの機構が寄与します。励起を失えば副作用として位相も壊れ、緩和速度の半分が寄与します。真の純位相緩和が速度 $1/T_\varphi$ で残りを担います。

$$ \frac{1}{T_2} = \frac{1}{2T_1} + \frac{1}{T_\varphi} $$

$1/T_\varphi \ge 0$ ですから、ただちに次の限界が従います。

$$ T_2 \le 2T_1 $$

これは健全性の検査として覚えておく価値があります。$2T_1$ より大きい $T_2$ の報告は、測定か計算のどこかが誤りです。$T_2 = 2T_1$ という極限はノイズがエネルギー緩和のみである量子ビットを表し、良質な超伝導素子が近づく実在の極限です。

$T_2^\ast$、不均一コヒーレンス時間は、素直な実験が実際に返す量です。量子ビット周波数が実験の繰り返しごとに厳密に同じでないとしましょう。近くの核スピン浴が向きを変えたから、捕獲された電荷が動いたから、磁場がドリフトしたから、という理由です。すると各ショットが異なる位相を蓄積し、個々のショットが完璧にコヒーレントであってもショット平均はコヒーレンスを消します。標準偏差 $\sigma$ をもつ静的離調のガウス分布に対しては

$$ \left\langle e^{i\Delta t} \right\rangle = e^{i\bar{\Delta}t}\, e^{-\sigma^2 t^2/2} \quad \Longrightarrow \quad T_{2,\mathrm{inh}}^\ast = \frac{\sqrt{2}}{\sigma} $$

となり、測定される全体の減衰はこのガウス関数と内在的な減衰の積になります。形に注目してください。周波数の静的なばらつきはガウス型の減衰を与え、系列が感じる帯域で白色なノイズ浴は指数型の減衰を与えます。したがって測定された減衰曲線の形そのものが診断的です。以下すべてに効く系として、指数型の包絡線はノイズのスペクトルについての証拠であって、どの時間定数を測っているかについての証拠ではありません。Code Example 7 では Hahn エコー — まさしく $T_2$ の測定 — がガウス型に減衰します。$1/f$ ノイズは白色とはほど遠いからです。

$T_2^\ast$ は量子ビットだけの性質ではありません。平均にどれだけ時間をかけたかに依存します。より遅いドリフトは、実験がそれを見るのに十分な長さ続いた場合にのみ寄与するからです。定義の欠陥のように聞こえますが、実は $1/f$ ノイズについての物理的な言明であり、Code Example 7 がそれを厳密にします。

記号 名称 減衰する対象 典型的な減衰形 材料について測っているもの
$T_1$ エネルギー緩和 $\langle Z\rangle$ が平衡へ 指数 $\omega_q$ でのノイズパワー: 準粒子、損失性誘電体、フォノン
$T_\varphi$ 純位相緩和 横成分。エネルギー損失なし 白色ノイズなら指数 直流付近の揺動ノイズパワー
$T_2$ コヒーレンス 横成分(エコー後) 白色ノイズなら指数、$1/f$ 下では $\exp[-(t/T_2)^2]$ 上記2つの組み合わせ
$T_2^\ast$ 不均一コヒーレンス 横成分(エコーなし) ガウス 量子ビット周波数の静的ばらつき

「典型的な減衰形」の列には、一度述べて以降使い回す注意が必要です。横成分の指数型減衰は、系列が感じる周波数帯でノイズが白色であることから従うものであり、環境についての言明であって $T_2$ の定義ではありません。固体量子ビットが実際に見る $1/f$ ノイズでは、エコーの包絡線もガウス型になります。Code Example 7 の $N = 1$(Hahnエコー)曲線は $T_2 = 0.683\ \mu$s の $\exp[-(t/T_2)^2]$ で減衰の最初の1桁にわたり数%以内で再現されます — $t = 0.494\ \mu$s で予測 0.593 に対し測定 0.586 — 一方で純粋な指数関数なら 0.485 を与えます。形を定義しているのは $T_2^\ast$ の行のガウス型だけで、他の行は形が条件つきです。

スペクトル密度としてのノイズ

3つの時間定数は、1つの基礎的な対象の要約統計量です。瞬時の量子ビット周波数が $\omega_q(t) = \bar{\omega}_q + \delta\omega(t)$ と揺らぐとし、$\delta\omega$ をその片側パワースペクトル密度 $S(f)$ で記述します。規格化は

$$ \left\langle \delta\omega^2 \right\rangle = \int_0^\infty S(f)\, df $$

とします。パルス系列は位相 $\varphi(t) = \int_0^t s(t')\,\delta\omega(t')\,dt'$ を蓄積し、$s(t) = \pm1$ は $\pi$ パルスが課す符号反転を記録します。ガウスノイズならコヒーレンスは $C(t) = |\langle e^{i\varphi}\rangle| = e^{-\langle\varphi^2\rangle/2}$ であり、位相の平均二乗は厳密に

$$ \left\langle \varphi^2(t) \right\rangle = \int_0^\infty df\; S(f)\, \left| \tilde{s}(f,t) \right|^2 , \qquad \tilde{s}(f,t) = \int_0^t s(t')\, e^{-2\pi i f t'}\, dt' $$

となります。関数 $|\tilde{s}(f,t)|^2$ が系列のフィルター関数です。実験が環境のどの周波数に感度をもつかを表し、完全に実験者の制御下にあります。閉じた形で持っておく価値のある2つの場合があります。

$$ \left| \tilde{s}_\mathrm{FID}(f,t) \right|^2 = \frac{\sin^2(\pi f t)}{(\pi f)^2}, \qquad \left| \tilde{s}_\mathrm{echo}(f,t) \right|^2 = \frac{4\sin^4(\pi f t/2)}{(\pi f)^2} $$

$\pi$ パルスをもたない自由誘導では $f \to 0$ で $|\tilde{s}|^2 \to t^2$ となり、最も遅いドリフトに対して最大の感度をもちます。中央に1発の $\pi$ パルスを置くエコーでは $|\tilde{s}|^2 \to \pi^2 f^2 t^4/4$ となり、直流で盲目であって最初の通過帯域は $f \approx 1/2t$ 付近に来ます。この違い1つが動的デカップリングの機構のすべてです。

$1/f$ ノイズと2準位揺動子

固体量子ビットはほぼ例外なく

$$ S(f) = \frac{A}{f^{\alpha}}, \qquad \alpha \approx 1 $$

を何桁にもわたって見ます。この普遍性には標準的な微視的説明があります。独立な2準位揺動子(TLS)の集団で、各揺動子は自身の速度 $\gamma$ で2つの配置間をランダムに切り替わります。1個の揺動子は $\gamma$ に折れ点をもつLorentz型スペクトルを与え、$\log\gamma$ について一様な速度分布 — これは熱活性化障壁高さの分布が与えるものです — を仮定すると、対応する範囲で和が $1/f$ になります。

ここで重要なのは1点です。揺動子は欠陥であり、したがって $A$ は材料パラメータであるということです。超伝導回路では、表面・界面・接合バリアのアモルファス酸化膜中の原子スケールの配置欠陥だと考えられています。半導体量子ビットでは酸化膜界面の電荷トラップと $^{29}\mathrm{Si}$ の核スピンです。トラップ電極上では吸着物によるパッチ電位です。第2章と第5章で個別に扱いますが、共通するのは、$A$ を下げることが素子を何で作るか、あるいはどう成長させるかを変えることを意味するという点です。

$1/f$ ノイズには、以降のコースで繰り返し現れる帰結が2つあります。第一に、$f \to 0$ でパワーが対数的に発散するため、分散 $\langle\delta\omega^2\rangle$ は低周波カットオフに依存し、実験的なカットオフは測定の継続時間で決まります。したがって「測定された $T_2^\ast$」は平均時間と併記して初めて意味をもちます。第二に、ノイズが低周波に集中しているため、それはまさにエコーが除去できる種類のノイズであり、固体量子ビットにおいて $T_2$ が $T_2^\ast$ をおおよそ2倍から10倍上回る理由です。本章の例がその範囲を挟んでいます。Code Example 7 は合成した $1/f$ ノイズに対して $\pi$ パルス1発で $T_2/T_2^\ast = 2.2$ を与え、Code Example 6 は静的なばらつきが支配してエコーがそのほぼ全部を除去する逆の極限で9倍を与えます。それより大きな比は、1発ではなくパルスを増やすことで得られます。同じ Code Example 7 は CPMG-16 で $8$ 倍に達します。

標準的なパルス系列

系列 パルス 測る量 フィルターの振る舞い
反転回復 $\pi$、待ち、読み出し $T_1$ $\omega_q$ で $S$ を測る
自由誘導減衰 / Ramsey $\pi/2$、時間 $t$ 待ち、$\pi/2$ $T_2^\ast$ 直流まで感度がある
Hahnエコー $\pi/2$、$t/2$、$\pi$、$t/2$、$\pi/2$ $T_2$ 直流で盲目、通過帯域は $1/2t$ 付近
CPMG-$N$ $\pi/2$、$N$ 個の $\pi$、読み出し $T_2(N)$ 通過帯域は $N/2t$ 付近
ノイズ分光 多数の $N$ について CPMG-$N$ $S(f)$ そのもの 通過帯域をスペクトル上で走査する

材料研究者が持ち帰るべきは最後の行です。$N$ を掃引すればフィルターの通過帯域がスペクトル上を移動しますから、コヒーレンス曲線の族を反転して数桁にわたる $S(f)$ が得られます。量子ビットは自らの母材中のノイズに対する感度の高い周波数分解プローブであり、それは計算機の診断であるだけでなく1つの特性評価手法です。


1.5 材料科学としての量子ハードウェア

本節は、以降のコースが立つ立場を述べます。

各方式が実際に止まる場所

どの方式でも、コヒーレンス時間は材料科学者が材料問題と認識するものによって制限されています。

方式 支配的なコヒーレンス限界 その背後にある材料の問い
超伝導回路 アモルファス表面酸化膜と界面のTLSによる誘電損失、非平衡準粒子 どの酸化膜、どの界面か。そして電場エネルギーの何割がそこにあるのか
イオントラップ トラップ表面の電場ノイズによる運動モードの異常加熱 電極表面に何が付いているのか。ノイズは距離と温度にどうスケールするのか
中性原子 光子散乱、原子損失、レーザーとトラップの強度ノイズ。原子そのものではない 主に光学の問題。ただしコーティング、真空表面、安定参照キャビティは材料
半導体スピン 核スピンとの超微細結合、酸化膜界面の電荷ノイズ 同位体精製で母材を無スピンにできるか。界面を無トラップにできるか
NVその他の欠陥中心 浅い中心では表面スピンと表面電荷トラップ スピン浴を作らずにダイヤモンド表面を終端できるか
トポロジカル量子ビット 半導体-超伝導ハイブリッド中の乱れ、柔らかい誘起ギャップ ハードギャップが生き残るほど清浄なエピタキシャル界面を作れるか

3列目を下に読むと、同じ3項目が繰り返し現れます。界面アモルファス層、そして同位体的・化学的な不純物。これらは材料科学で最も古い主題です。新しいのはプローブの感度です。量子ビットは従来の特性評価手法では気づかない欠陥密度に応答し、しかも誘電損失のデータがほとんど存在しないミリケルビン・ギガヘルツの領域で応答します。

量は材料の量である

量子プロセッサの性能指標は、ほとんど損失なく材料分野の量に翻訳できます。

素子側の指標 材料側の量
超伝導量子ビットの $T_1$ mK・GHz における関与誘電体の損失角 $\tan\delta$
シリコン中スピン量子ビットの $T_2^\ast$ 残留 $^{29}\mathrm{Si}$ 濃度と界面トラップ密度
共振器の内部Q値 $Q_i$ 表面参加率と表面損失角の積
イオンの運動加熱率 電極表面における電場ノイズパワースペクトル密度
2量子ビットゲート誤りの下限 パラメータのドリフト、すなわち $1/f$ 振幅 $A$
素子の歩留まり 接合面積とバリア厚の統計的ばらつき

右列の各項目は成長、プロセス、表面化学、精製で攻めることができ、しかも材料科学者がすでに使っている単位で測られます。この橋は実在し、渡れるだけの幅の狭さです。

本コースの立場

そこから3つの約束が従い、それが残りの章に何が入り何が入らないかを決めます。

スペックではなく物理と材料。 量子ビット数も、忠実度の記録も、ロードマップも書きません。そうした数値は読まれる前に古くなり、そしてより重要なことに、説明の力をもちません。スケーリング則、選択則、参加率の議論は真であり続けます。

どの方式もノイズチャネルへの解答として提示します。 第2章はトランズモンを電荷ノイズへの応答として、第3章はMølmer-Sørensenゲートを運動状態への感度への応答として、第5章はトポロジカルな提案を局所ノイズ一般への応答として説明します。設計を理解する最も簡潔な道筋です。

「材料がまだ十分でない」が誠実な答えであるときは、そう答えます。 そしてそれは想定読者にとって最も興味深い答えでもあります。仕事がどこにあるかを教えてくれるからです。


1.6 数値実験室: Bloch方程式

第1.4節で定義したすべては1つの運動方程式から計算でき、本章の残りはそれを行います。ここで作る道具は、第2章から第5章でコヒーレンスに関する主張を検証する必要が生じたときに使います。

方程式

駆動周波数で回転する座標系において、$\hbar = 1$ とすると、駆動された2準位系のハミルトニアンは

$$ H = \frac{1}{2}\left( \Delta\, \sigma_z + \Omega\, \sigma_x \right) $$

です。$\Omega$ はRabi周波数、$\Delta = \omega_q - \omega_d$ は離調で、いずれも角周波数です。純粋状態は $\boldsymbol{\omega} = (\Omega, 0, \Delta)$ として剛体回転 $\dot{\mathbf{r}} = \boldsymbol{\omega} \times \mathbf{r}$ に従って発展します。北極にある基底状態への緩和を加えるとBloch方程式が得られます。

$$ \begin{aligned} \dot{x} &= -\Delta y - x/T_2 \cr \dot{y} &= \Delta x - \Omega z - y/T_2 \cr \dot{z} &= \Omega y - (z - 1)/T_1 \end{aligned} $$

コードの前に3点。この方程式は現象論的で、$T_1$ と $T_2$ は手で入れます。その便利さの代償として非マルコフ的なノイズを記述できず、それゆえ Code Example 7 ではこの方程式を捨てて陽なノイズ軌跡に移ります。また $|\mathbf{r}|$ が縮むので混合状態を記述しています。姉妹コース第5章の密度行列の言語を、等価なベクトル形式で書いたものです。そしてこれはまさに磁気共鳴の方程式ですから、NMRやESRの実験を解釈した経験のある人はすでにこれを使っています。

励起状態の占有数は $P(1) = (1 - z)/2$ で、$z = +1$ が基底状態 $|0\rangle$ です。

Code Example 3: 緩和ありとなしのRabi振動

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 以下の周波数はすべて角周波数(rad/s)です。周期周波数には 2 pi を明示的に
# 掛けており、時間はマイクロ秒で表示します。
US = 1e-6


def bloch_rhs(t, r, Omega, Delta, T1, T2):
    """回転座標系におけるBloch方程式の右辺。

    r = (x, y, z) で、z = +1 が基底状態 |0> です。駆動はRabi周波数 Omega で
    x 軸方向、離調 Delta = omega_qubit - omega_drive は z 軸方向にかかります。
    緩和は z を速度 1/T1 で +1 に引き戻し、横成分を速度 1/T2 で縮めます。
    """
    x, y, z = r
    return [-Delta * y - x / T2,
            Delta * x - Omega * z - y / T2,
            Omega * y - (z - 1.0) / T1]


def evolve(r0, t_grid, Omega=0.0, Delta=0.0, T1=np.inf, T2=np.inf):
    """r0 を初期条件として t_grid 上でBloch方程式を積分する。"""
    sol = solve_ivp(bloch_rhs, (t_grid[0], t_grid[-1]), r0, t_eval=t_grid,
                    args=(Omega, Delta, T1, T2), rtol=1e-10, atol=1e-12,
                    method="DOP853")
    return sol.y                      # 形状 (3, len(t_grid))


def p_excited(r):
    """Blochベクトルから |1> の占有数を求める: P(1) = (1 - z)/2。"""
    return (1.0 - r[2]) / 2.0


# --- コヒーレントなRabi振動を解析解と照合する ---------------------------------
f_rabi = 10e6                          # Rabi周波数 10 MHz
Omega = 2 * np.pi * f_rabi
t = np.linspace(0, 0.30 * US, 3001)

print("一般化Rabi振動(緩和なし)")
print(f"  Omega/2pi = {f_rabi/1e6:.1f} MHz, pi pulse = "
      f"{np.pi/Omega/US*1e3:.1f} ns")
print(f"\n{'Delta/2pi (MHz)':>16}{'max P(1) num':>15}{'analytic':>11}"
      f"{'eff. rate/2pi':>15}")
print("-" * 57)
for f_det in [0.0, 5.0, 10.0, 20.0]:
    Delta = 2 * np.pi * f_det * 1e6
    r = evolve([0.0, 0.0, 1.0], t, Omega=Omega, Delta=Delta)
    num = p_excited(r).max()
    ana = Omega**2 / (Omega**2 + Delta**2)
    rate = np.sqrt(Omega**2 + Delta**2) / (2 * np.pi) / 1e6
    print(f"{f_det:>16.1f}{num:>15.6f}{ana:>11.6f}{rate:>12.3f} MHz")

# --- 同じ駆動に緩和を加える: 振動は減衰し飽和状態に落ち着く -------------------
T1, T2 = 20 * US, 10 * US
t_long = np.linspace(0, 30 * US, 60001)
r = evolve([0.0, 0.0, 1.0], t_long, Omega=Omega, Delta=0.0, T1=T1, T2=T2)
p = p_excited(r)

z_ss = 1.0 / (1.0 + Omega**2 * T1 * T2)
T_rabi = 2.0 / (1.0 / T1 + 1.0 / T2)      # 駆動Rabi振動の包絡線の減衰時間
print(f"\n駆動下の減衰(T1 = {T1/US:.0f} us, T2 = {T2/US:.0f} us)")
print(f"  P(1) at the first maximum      : {p[:400].max():.6f}")
print(f"  P(1) at t = 30 us              : {p[-1]:.6f}")
print(f"  analytic steady state (1-z)/2  : {(1 - z_ss)/2:.6f}")
print(f"  Omega^2 T1 T2                  : {Omega**2 * T1 * T2:.4g}")
print(f"  envelope decay 2/(1/T1 + 1/T2) : {T_rabi/US:.4f} us")

# 振動のコントラストを時間窓ごとに測り、上の予測と比較する
print(f"\n{'window (us)':>13}{'measured envelope':>19}{'0.5 exp(-t/T_rabi)':>21}")
print("-" * 53)
for t0 in [0.0, 1.0, 2.0, 5.0, 10.0, 20.0]:
    m = (t_long >= t0 * US) & (t_long <= (t0 + 0.2) * US)
    env = (p[m].max() - p[m].min()) / 2
    t_mid = (t0 + 0.1) * US
    print(f"{f'{t0:.1f}-{t0+0.2:.1f}':>13}{env:>19.6f}"
          f"{0.5*np.exp(-t_mid/T_rabi):>21.6f}")

# --- 可視化 -----------------------------------------------------------------
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
for f_det in [0.0, 5.0, 10.0, 20.0]:
    r_ = evolve([0.0, 0.0, 1.0], t, Omega=Omega,
                Delta=2 * np.pi * f_det * 1e6)
    ax[0].plot(t / US, p_excited(r_), label=f"$\\Delta/2\\pi$ = {f_det:.0f} MHz")
ax[0].set_xlabel("time (us)"); ax[0].set_ylabel("P(1)")
ax[0].set_title("Rabi oscillations vs detuning"); ax[0].legend(fontsize=8)

ax[1].plot(t_long / US, p, lw=0.4, color="tab:purple")
ax[1].axhline((1 - z_ss) / 2, ls="--", color="k", lw=1,
              label="steady state")
ax[1].set_xlabel("time (us)"); ax[1].set_ylabel("P(1)")
ax[1].set_title("Driven decay into saturation"); ax[1].legend(fontsize=8)
plt.tight_layout()
plt.show()
一般化Rabi振動(緩和なし)
  Omega/2pi = 10.0 MHz, pi pulse = 50.0 ns

 Delta/2pi (MHz)   max P(1) num   analytic  eff. rate/2pi
---------------------------------------------------------
             0.0       1.000000   1.000000      10.000 MHz
             5.0       0.800000   0.800000      11.180 MHz
            10.0       0.500000   0.500000      14.142 MHz
            20.0       0.200000   0.200000      22.361 MHz

駆動下の減衰(T1 = 20 us, T2 = 10 us)
  P(1) at the first maximum      : 0.998127
  P(1) at t = 30 us              : 0.447300
  analytic steady state (1-z)/2  : 0.499999
  Omega^2 T1 T2                  : 7.896e+05
  envelope decay 2/(1/T1 + 1/T2) : 13.3333 us

  window (us)  measured envelope   0.5 exp(-t/T_rabi)
-----------------------------------------------------
      0.0-0.2           0.499064             0.496264
      1.0-1.2           0.463003             0.460406
      2.0-2.2           0.429548             0.427138
      5.0-5.2           0.343001             0.341077
    10.0-10.2           0.235741             0.234419
    20.0-20.2           0.111356             0.110731

注目すべき点。 離調の表は、あらゆる2準位ソルバーに対する標準的な検査です。励起状態占有数の最大値は $\Omega^2/(\Omega^2 + \Delta^2)$、振動は一般化Rabi周波数 $\sqrt{\Omega^2 + \Delta^2}$ で走ります。数値はどちらも6桁まで再現しています。実務的な内容は、Rabi周波数に等しい離調があるだけで占有数の上限が $1/2$ になることです。周波数が合っていない量子ビットに $\pi$ パルスを当てると、単に回転しすぎるのではなく $|1\rangle$ に到達すらできません。どの実験室でも周波数校正が振幅校正に先立つのはこのためです。

2番目のブロックは、多数の減衰時間にわたって駆動を続けたときに何が起きるかを示します。振動の包絡線は速度 $\frac{1}{2}(1/T_1 + 1/T_2)$ で減衰し、窓ごとの測定がこれを3桁まで確認しています。占有数は $\frac{1}{2}\left(1 - (1 + \Omega^2 T_1 T_2)^{-1}\right)$ に落ち着き、$\Omega^2 T_1 T_2 \approx 8\times10^5$ では $1/2$ より $10^{-6}$ のオーダーだけ下です。強く駆動された量子ビットは最大混合状態に至ります。連続駆動は状態を保持する手段ではなく、有用な動作領域は常に飽和よりRabi周期の何桁も手前です。

Code Example 4: $T_1$、$T_2$、そして $T_2 \le 2T_1$ の限界

2つの減衰時間は、同じ量子ビットに対する2つの異なる実験から得られます。この例は両方を実行し、フィットし、両者の関係を検証します。

"""第1章 Code Example 4: T1、T2、そして限界 T2 <= 2 T1。
Code Example 3 の続き(同一セッション)。"""


def fit_exponential(t, y):
    """y = A exp(-t/tau) の最小二乗フィット。(tau, A) を返す。"""
    mask = y > 1e-6 * y[0]
    slope, intercept = np.polyfit(t[mask], np.log(y[mask]), 1)
    return -1.0 / slope, np.exp(intercept)


def T2_from(T1, T_phi):
    """1/T2 = 1/(2 T1) + 1/T_phi: 緩和はその速度の半分だけ寄与する。"""
    return 1.0 / (1.0 / (2.0 * T1) + 1.0 / T_phi)


T1 = 30 * US
print("エネルギー緩和と位相緩和は別々の測定です。\n")
print(f"{'T_phi (us)':>11}{'T2 (us)':>10}{'T2/(2 T1)':>11}"
      f"{'fitted T1 (us)':>16}{'fitted T2 (us)':>16}")
print("-" * 64)
for T_phi in [np.inf, 300 * US, 60 * US, 15 * US, 3 * US]:
    T2 = T2_from(T1, T_phi)

    # T1測定: |1>(z = -1)を準備し、駆動なしで z が +1 に戻る様子を見る。
    t1_grid = np.linspace(0, 4 * T1, 4001)
    z = evolve([0.0, 0.0, -1.0], t1_grid, T1=T1, T2=T2)[2]
    tau1, _ = fit_exponential(t1_grid, (1.0 - z) / 2.0)

    # T2測定: |+>(x = 1)を準備し、駆動なしで横成分が減衰する様子を見る。
    # これが1量子ビットの自由誘導減衰(FID)です。
    t2_grid = np.linspace(0, 4 * T2, 4001)
    x = evolve([1.0, 0.0, 0.0], t2_grid, T1=T1, T2=T2)[0]
    tau2, _ = fit_exponential(t2_grid, x)

    label = "inf" if not np.isfinite(T_phi) else f"{T_phi/US:.0f}"
    print(f"{label:>11}{T2/US:>10.3f}{T2/(2*T1):>11.4f}"
          f"{tau1/US:>16.4f}{tau2/US:>16.4f}")

print("\nT_phi -> inf は T2 = 2 T1 を厳密に再現します。純位相緩和がなければ、")
print("位相情報を失う唯一の経路は励起そのものを失うことだからです。")

# 2つの速度は環境ノイズスペクトルのどの部分に感度をもつのか。
print("\n2つの速度はノイズスペクトルの異なる部分に感度をもちます:")
f_q = 5e9
for name, f_probe in [("1/T1  <- S(f) at the qubit frequency", f_q),
                      ("1/T2  <- S(f) near DC, set by the sequence", 1.0 / (30 * US))]:
    print(f"  {name:<45} f = {f_probe:9.4g} Hz")
print("  1つの数値で両方を記述することはできません。T1はGHzでの分光であり、")
print("  T2はkHzでの分光です。コード例7でこれを定量化します。")

# --- 可視化 -----------------------------------------------------------------
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
t_grid = np.linspace(0, 120 * US, 6001)
for T_phi, c in [(np.inf, "tab:blue"), (60 * US, "tab:orange"),
                 (15 * US, "tab:green"), (3 * US, "tab:red")]:
    T2 = T2_from(T1, T_phi)
    lbl = ("T_phi = inf" if not np.isfinite(T_phi)
           else f"T_phi = {T_phi/US:.0f} us")
    ax[0].plot(t_grid / US,
               (1 - evolve([0, 0, -1.0], t_grid, T1=T1, T2=T2)[2]) / 2,
               color=c, label=lbl)
    ax[1].plot(t_grid / US, evolve([1.0, 0, 0], t_grid, T1=T1, T2=T2)[0],
               color=c, label=lbl)
ax[0].set_title(f"T1 decay of P(1), T1 = {T1/US:.0f} us")
ax[1].set_title("Free-induction decay of the transverse component")
for a in ax:
    a.set_xlabel("time (us)"); a.set_yscale("log"); a.set_ylim(1e-3, 1.2)
    a.legend(fontsize=8)
plt.tight_layout()
plt.show()
エネルギー緩和と位相緩和は別々の測定です。

 T_phi (us)   T2 (us)  T2/(2 T1)  fitted T1 (us)  fitted T2 (us)
----------------------------------------------------------------
        inf    60.000     1.0000         30.0000         60.0000
        300    50.000     0.8333         30.0000         50.0000
         60    30.000     0.5000         30.0000         30.0000
         15    12.000     0.2000         30.0000         12.0000
          3     2.857     0.0476         30.0000          2.8571

T_phi -> inf は T2 = 2 T1 を厳密に再現します。純位相緩和がなければ、
位相情報を失う唯一の経路は励起そのものを失うことだからです。

2つの速度はノイズスペクトルの異なる部分に感度をもちます:
  1/T1  <- S(f) at the qubit frequency          f =     5e+09 Hz
  1/T2  <- S(f) near DC, set by the sequence    f = 3.333e+04 Hz
  1つの数値で両方を記述することはできません。T1はGHzでの分光であり、
  T2はkHzでの分光です。コード例7でこれを定量化します。

注目すべき点。 フィットは入力した $T_1$ と $T_2$ を4桁まで回復します。これは退屈だが必要な部分です。興味深いのは $T_2/(2T_1)$ の列で、純位相緩和が全くなければ厳密に1に等しく、そのときプロット中の2本の減衰曲線はこの因子2だけで違います。有限の $T_\varphi$ はすべてこの比を下げ、強く位相緩和した量子ビットでは $T_2 \ll T_1$ になります。実素子ではこの比自身が有用な診断です。$T_2 \approx 2T_1$ の量子ビットはエネルギー損失に律速されており、より良い誘電体を必要とします。$T_2 \ll T_1$ の量子ビットは周波数揺らぎに律速されており、より静かな電荷または磁束環境を必要とします。2つの診断は別々の作製工程を指し示します。

最後のブロックは、2つの速度を1つの数値に押し込められない理由を予告します。両者は9桁も離れた周波数で環境を測っているのです。

Code Example 5: Ramsey縞と $T_2^\ast$

Ramsey実験は待ち時間で隔てられた2発の $\pi/2$ パルスであり、蓄積した位相を占有数として測ります。この例は単一の量子ビットで実行し、続いてショットごとに周波数が異なるアンサンブルで実行します。

"""第1章 Code Example 5: Ramsey縞と不均一時間 T2*。
Code Example 3 の続き(同一セッション)。"""


def rotate(r, axis, angle):
    """Blochベクトルを軸 `axis` まわりに角度 `angle` だけ回す(Rodrigues公式)。

    瞬時パルスの理想化です。実際のパルスは有限の長さをもちますが、パルス長が
    どの減衰時間よりも短いかぎりこの理想化は妥当です。
    """
    n = np.asarray(axis, dtype=float)
    n = n / np.linalg.norm(n)
    r = np.asarray(r, dtype=float)                  # 形状 (3,) または (3, N)
    cross = np.stack([n[1] * r[2] - n[2] * r[1],
                      n[2] * r[0] - n[0] * r[2],
                      n[0] * r[1] - n[1] * r[0]])
    ndotr = n[0] * r[0] + n[1] * r[1] + n[2] * r[2]
    return (r * np.cos(angle) + cross * np.sin(angle)
            + np.multiply.outer(n, ndotr) * (1.0 - np.cos(angle)))


def free_precession(r, t, Delta, T1, T2):
    """自由発展の解析解: z 軸まわりの歳差運動と T1/T2 減衰。"""
    x, y, z = r
    c, s = np.cos(Delta * t), np.sin(Delta * t)
    d2, d1 = np.exp(-t / T2), np.exp(-t / T1)
    return np.array([(x * c - y * s) * d2,
                     (x * s + y * c) * d2,
                     1.0 + (z - 1.0) * d1])


def ramsey(tau, Delta, T1, T2):
    """x軸まわりのpi/2、時間tauの自由発展、再びpi/2、そして z を読み出す。"""
    r = rotate([0.0, 0.0, 1.0], [1, 0, 0], np.pi / 2)
    r = free_precession(r, tau, Delta, T1, T2)
    r = rotate(r, [1, 0, 0], np.pi / 2)
    return (1.0 - r[2]) / 2.0                # P(1)


# 解析解をコード例3のODE積分器と照合する。
chk_t = np.array([0.0, 0.37 * US, 1.4 * US])
chk_num = evolve([0.6, -0.3, 0.5], chk_t, Delta=2 * np.pi * 1.7e6,
                 T1=30 * US, T2=8 * US)
chk_ana = np.array([free_precession([0.6, -0.3, 0.5], t, 2 * np.pi * 1.7e6,
                                    30 * US, 8 * US) for t in chk_t]).T
print("free_precession と ODE積分器の差の最大値 = "
      f"{np.abs(chk_num - chk_ana).max():.2e}")

# --- 単一の離調をもつ1量子ビット: 離調周波数で綺麗な縞が出る -------------------
T1, T2 = 30 * US, 10 * US
f_bar = 3.0e6                       # 平均離調 3 MHz
sigma_f = 0.2e6                     # 静的離調のばらつき 0.2 MHz
Delta_bar, sigma = 2 * np.pi * f_bar, 2 * np.pi * sigma_f

tau = np.linspace(0, 6 * US, 1201)
single = ramsey(tau, Delta_bar, T1, T2)

# --- 静的離調のアンサンブル: 縞が洗い流される --------------------------------
rng = np.random.default_rng(20260813)
n_shots = 20000
detunings = Delta_bar + sigma * rng.standard_normal(n_shots)
ensemble = np.zeros_like(tau)
for d in detunings:                   # 1つが実験の1ショットに対応
    ensemble += ramsey(tau, d, T1, T2)
ensemble /= n_shots

# 準静的理論: <exp(i Delta tau)> = exp(i Delta_bar tau - sigma^2 tau^2 / 2)
T2_star_inh = np.sqrt(2.0) / sigma
envelope = np.exp(-tau / T2) * np.exp(-(tau / T2_star_inh) ** 2)

print(f"\nRamsey測定(平均離調 {f_bar/1e6:.1f} MHz, "
      f"静的ばらつき {sigma_f/1e6:.1f} MHz のガウス分布)")
print(f"  intrinsic T2                       : {T2/US:8.3f} us")
print(f"  inhomogeneous sqrt(2)/sigma        : {T2_star_inh/US:8.3f} us")
print(f"  fringe period 1/f_bar              : {1/f_bar/US:8.3f} us")
print(f"  fringes visible before the envelope dies: "
      f"{T2_star_inh*f_bar:.1f}")

print(f"\n{'tau (us)':>9}{'single detuning':>17}{'ensemble':>11}"
      f"{'|contrast|':>12}{'theory envelope':>17}")
print("-" * 66)
for t_mark in [0.0, 0.25, 0.5, 1.0, 1.5, 2.0, 3.0]:
    i = int(np.argmin(np.abs(tau - t_mark * US)))
    contrast = abs(2 * ensemble[i] - 1.0)
    print(f"{tau[i]/US:>9.3f}{single[i]:>17.6f}{ensemble[i]:>11.6f}"
          f"{contrast:>12.6f}{envelope[i]:>17.6f}")

# 実験と同じ手順でアンサンブル曲線から T2* を取り出す
peaks = []
for i in range(1, len(tau) - 1):
    if ensemble[i] >= ensemble[i - 1] and ensemble[i] >= ensemble[i + 1]:
        peaks.append(i)
amp = np.array([abs(2 * ensemble[i] - 1.0) for i in peaks])
tp = tau[peaks]
use = amp > 0.05
coef = np.polyfit(tp[use] ** 2, np.log(amp[use]), 1)
T2_star_fit = 1.0 / np.sqrt(-coef[0])
print(f"\n縞の極大値へのガウスフィット({use.sum()} 点):")
print(f"  T2* from the fit                   : {T2_star_fit/US:8.3f} us")
print(f"  sqrt(2)/sigma for comparison       : {T2_star_inh/US:8.3f} us")
print(f"  T2 / T2*                           : {T2/T2_star_fit:8.2f}")
print("測定されるコヒーレンス時間を決めているのは量子ビット周波数の*ばらつき*")
print("であり、個々の量子ビットのコヒーレンスではありません。")

# --- 可視化 -----------------------------------------------------------------
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(tau / US, single, lw=0.8, color="lightsteelblue",
        label="single detuning")
ax.plot(tau / US, ensemble, lw=1.2, color="tab:purple", label="ensemble")
ax.plot(tau / US, 0.5 + 0.5 * envelope, "k--", lw=1, label="theory envelope")
ax.plot(tau / US, 0.5 - 0.5 * envelope, "k--", lw=1)
ax.set_xlabel("free evolution time tau (us)"); ax.set_ylabel("P(1)")
ax.set_title("Ramsey fringes: T2* from inhomogeneous broadening")
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()
free_precession と ODE積分器の差の最大値 = 7.46e-11

Ramsey測定(平均離調 3.0 MHz, 静的ばらつき 0.2 MHz のガウス分布)
  intrinsic T2                       :   10.000 us
  inhomogeneous sqrt(2)/sigma        :    1.125 us
  fringe period 1/f_bar              :    0.333 us
  fringes visible before the envelope dies: 3.4

 tau (us)  single detuning   ensemble  |contrast|  theory envelope
------------------------------------------------------------------
    0.000         1.000000   1.000000    1.000000         1.000000
    0.250         0.500000   0.499079    0.001842         0.928349
    0.500         0.024385   0.110374    0.779253         0.780834
    1.000         0.952419   0.703752    0.407504         0.410833
    1.500         0.069646   0.429123    0.141754         0.145653
    2.000         0.909365   0.515379    0.030757         0.034795
    3.000         0.870409   0.499695    0.000610         0.000608

縞の極大値へのガウスフィット(5 点):
  T2* from the fit                   :    1.081 us
  sqrt(2)/sigma for comparison       :    1.125 us
  T2 / T2*                           :     9.25
測定されるコヒーレンス時間を決めているのは量子ビット周波数の*ばらつき*
であり、個々の量子ビットのコヒーレンスではありません。

注目すべき点。 冒頭の検査は見かけより重要です。解析的な自由発展の伝播子がODE積分器と $10^{-11}$ まで一致しており、これが以降の20000ショットのアンサンブル平均で高速な閉形式を使う根拠になります。

そして物理です。単一の離調は内在的な $T_2 = 10\ \mu$s まで持続する縞を与えます。わずか 0.2 MHz — 平均離調の7%未満 — のガウス分布で平均すると、約 $1\ \mu$s でそれが消えます。9倍も早いのです。個々の量子ビットについては何も変わっていません。縞の極大へのガウスフィットは $1.081\ \mu$s を返し、準静的な予測 $\sqrt{2}/\sigma = 1.125\ \mu$s と比べてわずかに短くなります。この差は内在的な指数減衰であり、純粋なガウスフィットがそれを吸収したものです。これは実験室でのフィットに実在する系統誤差そのものであり、フィットモデルを述べずに $T_2^\ast$ を引用することは10%の水準で曖昧だということです。

表の1行だけ説明が必要です。$\tau = 0.25\ \mu$s ではコントラストが 0.0018 と読める一方、包絡線は 0.93 です。この時刻は縞の節に当たり、コヒーレンスと無関係に $P(1) = 1/2$ になります。コントラストは縞の極大で読む必要があり、フィットはそれを行っています。

Code Example 6: Hahnエコーは $T_2^\ast$ と $T_2$ を分離する

待ち時間の中央に $\pi$ パルスを1発挿入すると、蓄積した位相の符号が反転します。したがって系列の間に変化しなかった離調は厳密に打ち消され、不均一性の寄与が消えます。

"""第1章 Code Example 6: Hahnエコーは T2* と T2 を分離する。
Code Example 3、4、5 の続き(同一セッション)。"""


def ramsey_state(t, Delta, T1, T2):
    """Ramsey系列の最後のpi/2パルス直前におけるBlochベクトル。"""
    r = rotate([0.0, 0.0, 1.0], [1, 0, 0], np.pi / 2)
    return free_precession(r, t, Delta, T1, T2)


def echo_state(t, Delta, T1, T2):
    """pi/2 - t/2 - pi - t/2 に対する同じ量。

    piパルスはそれまでに蓄積した位相の符号を反転させるため、系列の間に
    変化しない離調は厳密に打ち消されます。
    """
    r = rotate([0.0, 0.0, 1.0], [1, 0, 0], np.pi / 2)
    r = free_precession(r, t / 2.0, Delta, T1, T2)
    r = rotate(r, [1, 0, 0], np.pi)
    return free_precession(r, t / 2.0, Delta, T1, T2)


def hahn_echo(t, Delta, T1, T2):
    """読み出しパルス -pi/2 を含む完全なエコー系列。P(1) を返す。"""
    r = rotate(echo_state(t, Delta, T1, T2), [1, 0, 0], -np.pi / 2)
    return (1.0 - r[2]) / 2.0


def coherence(sequence, t, ensemble, T1, T2):
    """非対角成分のアンサンブル平均 |<x + i y>|。

    これが縞の*包絡線*です。平均の位相が縞を、その絶対値がコヒーレンスを
    表します。複素量のまま平均することは、実験で系列をショットごとに
    繰り返すことに対応します。
    """
    acc = np.zeros(np.shape(t), dtype=complex)
    for d in ensemble:
        r = sequence(t, d, T1, T2)
        acc += r[0] + 1j * r[1]
    return np.abs(acc / len(ensemble))


t_tot = np.linspace(0, 40 * US, 2001)
coh_ramsey = coherence(ramsey_state, t_tot, detunings, T1, T2)
coh_echo = coherence(echo_state, t_tot, detunings, T1, T2)

print("エコーは系列の間に変化しない離調に対して盲目です。")
print(f"{'t (us)':>8}{'Ramsey coherence':>18}{'echo coherence':>16}"
      f"{'exp(-t/T2)':>12}{'echo P(1), one shot':>21}")
print("-" * 75)
for t_mark in [0.0, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 40.0]:
    i = int(np.argmin(np.abs(t_tot - t_mark * US)))
    one = hahn_echo(t_tot[i], Delta_bar + 5 * sigma, T1, T2)
    print(f"{t_tot[i]/US:>8.2f}{coh_ramsey[i]:>18.6f}{coh_echo[i]:>16.6f}"
          f"{np.exp(-t_tot[i]/T2):>12.6f}{one:>21.6f}")


def one_over_e(t, c):
    """1/e を最初に下回る時刻。線形補間で求め、到達しなければ nan を返す。"""
    below = np.nonzero(c < 1.0 / np.e)[0]
    if below.size == 0 or below[0] == 0:
        return np.nan
    k = below[0]
    return np.interp(1.0 / np.e, [c[k], c[k - 1]], [t[k], t[k - 1]])


t_ramsey = one_over_e(t_tot, coh_ramsey)
t_echo = one_over_e(t_tot, coh_echo)
print(f"\n{'sequence':<28}{'1/e coherence time':>20}")
print("-" * 48)
print(f"{'Ramsey (free induction)':<28}{t_ramsey/US:>17.3f} us")
print(f"{'Hahn echo':<28}{t_echo/US:>17.3f} us")
print(f"{'intrinsic T2 (input)':<28}{T2/US:>17.3f} us")
print(f"{'echo gain':<28}{t_echo/t_ramsey:>17.2f} x")

# エコーは不均一性を除去するが、エネルギー緩和は除去しない。
t_wide = np.linspace(0, 200 * US, 2001)
print("\nエコーは T2 を超えられず、T2 は 2 T1 を超えられません:")
print(f"{'T_phi (us)':>11}{'T2 (us)':>10}{'echo 1/e (us)':>15}{'2 T1 (us)':>11}")
print("-" * 47)
T1_ref = 30 * US
for T_phi in [np.inf, 60 * US, 15 * US]:
    T2_i = T2_from(T1_ref, T_phi)
    e_curve = coherence(echo_state, t_wide, detunings[:2000], T1_ref, T2_i)
    label = "inf" if not np.isfinite(T_phi) else f"{T_phi/US:.0f}"
    print(f"{label:>11}{T2_i/US:>10.3f}"
          f"{one_over_e(t_wide, e_curve)/US:>15.3f}{2*T1_ref/US:>11.3f}")

print("\nしたがって同一試料で T2* と T2 を測ることは、その材料について2つの")
print("異なる量を測ることになります。すなわち量子ビット周波数の静的な")
print("ばらつきと、その揺動成分です。")

# --- 可視化 -----------------------------------------------------------------
fig, ax = plt.subplots(figsize=(9, 4))
ax.semilogy(t_tot / US, np.maximum(coh_ramsey, 1e-6), color="tab:orange",
            label="Ramsey envelope (T2*)")
ax.semilogy(t_tot / US, np.maximum(coh_echo, 1e-6), color="tab:purple",
            label="Hahn echo (T2)")
ax.semilogy(t_tot / US, np.exp(-t_tot / T2), "k--", lw=1, label="exp(-t/T2)")
ax.axhline(1 / np.e, color="gray", lw=0.8, ls=":")
ax.set_xlabel("total sequence time (us)"); ax.set_ylabel("coherence")
ax.set_ylim(1e-4, 1.5)
ax.set_title("Ramsey vs Hahn echo on the same qubit ensemble")
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()
エコーは系列の間に変化しない離調に対して盲目です。
  t (us)  Ramsey coherence  echo coherence  exp(-t/T2)  echo P(1), one shot
---------------------------------------------------------------------------
    0.00          1.000000        1.000000    1.000000             1.000000
    0.50          0.779258        0.951229    0.951229             0.975615
    1.00          0.407513        0.904837    0.904837             0.952419
    2.00          0.030780        0.818731    0.818731             0.909365
    5.00          0.006403        0.606531    0.606531             0.803265
   10.00          0.001959        0.367879    0.367879             0.683940
   20.00          0.000506        0.135335    0.135335             0.567668
   40.00          0.000139        0.018316    0.018316             0.509158

sequence                      1/e coherence time
------------------------------------------------
Ramsey (free induction)                 1.059 us
Hahn echo                              10.000 us
intrinsic T2 (input)                   10.000 us
echo gain                                9.45 x

エコーは T2 を超えられず、T2 は 2 T1 を超えられません:
 T_phi (us)   T2 (us)  echo 1/e (us)  2 T1 (us)
-----------------------------------------------
        inf    60.000         60.000     60.000
         60    30.000         30.000     60.000
         15    12.000         12.000     60.000

したがって同一試料で T2* と T2 を測ることは、その材料について2つの
異なる量を測ることになります。すなわち量子ビット周波数の静的な
ばらつきと、その揺動成分です。

注目すべき点。 エコーのコヒーレンス列と $e^{-t/T_2}$ の列は小数6桁まで同一であり、平均から5標準偏差ずれた離調について計算した単一ショットのエコー列もアンサンブル平均と厳密に一致します。この一致こそがエコーです。観測量が離調にまったく依存しなくなるので、離調について平均しても何も起きません。同じ量子ビットアンサンブルで 1/e コヒーレンス時間は 1.06 $\mu$s から 10.0 $\mu$s へ、9.5倍になり、回復した値は内在的な $T_2$ を4桁まで再現します。

2番目の表が境界を引きます。エコーは不均一性を除去しますが、エネルギー緩和は除去できません。純位相緩和がなければエコーは $2T_1$ を返し、それより1ナノ秒も長くはなりません。その時点で残る唯一の減衰経路が励起そのものの損失であり、どんなパルス系列もそれを元に戻せないからです。以降のコースにおけるすべての動的デカップリングの結果は、この線で上から抑えられています。

物理的な読み方は、最後の3行が述べているものです。同一試料で測った $T_2^\ast$ と $T_2$ は2つの異なる材料測定です。前者はショット間および素子間の量子ビット周波数の静的なばらつきを、後者はエコーを生き延びるほど速い揺動成分を与えます。一方を改善し他方を変えないプロセス変更は、どちらの欠陥集団に効いたのかを教えてくれます。

Code Example 7: $1/f$ ノイズ、動的デカップリング、ノイズ分光

準静的な描像(Code Example 5)と厳密な打ち消しの描像(Code Example 6)は、どちらも理想化です。実際のノイズはスペクトルをもち、エコーが働くのはその遅い部分に対してだけです。この例は $1/f$ ノイズの軌跡を陽に生成し、その上でCPMG系列を走らせ、コヒーレンス時間とノイズ指数を結ぶスケーリング則を回復します。

"""第1章 Code Example 7: 1/f ノイズ、フィルタ関数、CPMGスケーリング。
Code Example 3 の続き(同一セッション)。"""


def noise_trajectories(n_traj, n_t, dt, alpha, A, f_low, rng):
    """片側PSD S(f) = A / f^alpha をもつ実ガウスノイズ delta_omega(t)。

    各フーリエ成分にランダムな位相と S(f) で決まる分散を与え、逆変換して
    生成します。f_low より低い成分は落とします。1/f ノイズは低周波カット
    オフなしには有限の分散をもたず、実験ではそのカットオフは測定の継続
    時間で決まります。
    """
    f = np.fft.rfftfreq(n_t, dt)
    S = np.zeros_like(f)
    band = f >= f_low
    S[band] = A / f[band] ** alpha
    scale = np.sqrt(S * n_t / (4.0 * dt))
    spec = (rng.standard_normal((n_traj, f.size))
            + 1j * rng.standard_normal((n_traj, f.size))) * scale
    return np.fft.irfft(spec, n=n_t, axis=1), f, S


def cpmg_modulation(n_steps, n_pi):
    """piパルスを n_pi 個含むCPMG系列の +-1 切替関数。

    n_pi = 0 は自由誘導(Ramsey)、n_pi = 1 は中央に1発のHahnエコーです。
    一般の n_pi では全時間の (j - 1/2)/n_pi の位置にパルスを置き、これが
    標準的なCPMGのタイミングです。
    """
    s = np.ones(n_steps)
    if n_pi == 0:
        return s
    u = (np.arange(n_steps) + 0.5) / n_steps
    flips = np.searchsorted((np.arange(1, n_pi + 1) - 0.5) / n_pi, u)
    return np.where(flips % 2 == 0, 1.0, -1.0)


dt, n_t, n_traj = 2e-9, 8192, 2000
m_max = n_t // 4                  # 系列は先頭1/4しか使わない。生成の周期性を
f_low = 1.0 / (n_t * dt)          # 決して踏まないようにするためです
alpha, A = 1.0, 6.0e12            # f = 1 Hz における (rad/s)^2 / Hz
rng = np.random.default_rng(11)
traj, f_axis, S_target = noise_trajectories(n_traj, n_t, dt, alpha, A,
                                            f_low, rng)

# 生成の検証: 標本PSDを目標値と、標本分散を解析積分 A ln(f_max / f_low) と
# 比較します。
spec = np.fft.rfft(traj[:500], axis=1)
S_sample = 2.0 * np.mean(np.abs(spec) ** 2, axis=0) * dt / n_t
band = (f_axis >= 10 * f_low) & (f_axis <= 0.2 / (2 * dt))
ratio = np.mean(S_sample[band] / S_target[band])
var_analytic = A * np.log((1.0 / (2 * dt)) / f_low)
print(f"1/f ノイズの生成: {n_t} 点・dt = {dt*1e9:.0f} ns の軌跡を "
      f"{n_traj} 本")
print(f"  band                          : {f_low/1e3:.1f} kHz to "
      f"{1/(2*dt)/1e6:.0f} MHz")
print(f"  mean sample S(f) / target S(f): {ratio:.4f}")
print(f"  sample variance of delta_omega: {np.var(traj):.4e} (rad/s)^2")
print(f"  A ln(f_max/f_low)             : {var_analytic:.4e} (rad/s)^2")
print(f"  rms detuning                  : "
      f"{np.sqrt(np.var(traj))/(2*np.pi)/1e3:.1f} kHz")
print(f"  quasi-static sqrt(2)/sigma    : "
      f"{np.sqrt(2/np.var(traj))/US:.3f} us")
print(f"  statistical floor 1/sqrt(n)   : {1/np.sqrt(n_traj):.4f}")

# --- 各系列におけるコヒーレンス ------------------------------------------------
n_pis = [0, 1, 2, 4, 8, 16]
t_idx = np.unique(np.round(np.logspace(np.log10(16), np.log10(m_max), 40)
                           ).astype(int))
times = t_idx * dt
curves = {}
for n_pi in n_pis:
    c = np.empty(t_idx.size)
    for k, m in enumerate(t_idx):
        s = cpmg_modulation(m, n_pi)
        phi = dt * (traj[:, :m] * s).sum(axis=1)      # 蓄積した位相
        c[k] = np.abs(np.mean(np.exp(1j * phi)))
    curves[n_pi] = c

print(f"\n同一の 1/f ノイズ下における CPMG-N のコヒーレンス |<exp(i phi)>|")
print(f"{'t (us)':>8}" + "".join(f"{('N=' + str(n)):>10}" for n in n_pis))
print("-" * (8 + 10 * len(n_pis)))
for t_mark in [0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 4.0]:
    k = int(np.argmin(np.abs(times - t_mark * US)))
    print(f"{times[k]/US:>8.3f}"
          + "".join(f"{curves[n][k]:>10.4f}" for n in n_pis))


def one_over_e_log(t, c):
    """1/e 交差時刻。log(コヒーレンス) と log(時間) の間で補間する。"""
    below = np.nonzero(c < 1.0 / np.e)[0]
    if below.size == 0 or below[0] == 0:
        return np.nan
    k = below[0]
    lo, hi = np.log(c[k]), np.log(c[k - 1])
    return np.exp(np.interp(-1.0, [lo, hi], [np.log(t[k]), np.log(t[k - 1])]))


T2s = {n: one_over_e_log(times, curves[n]) for n in n_pis}
print(f"\n{'N (pi pulses)':>14}{'T2(N) (us)':>13}{'T2(N)/T2(1)':>14}")
print("-" * 41)
for n in n_pis:
    label = "0 (Ramsey)" if n == 0 else str(n)
    print(f"{label:>14}{T2s[n]/US:>13.4f}{T2s[n]/T2s[1]:>14.3f}")

# S(f) = A/f^alpha に対して、エコー系列は T2(N) ~ N^(alpha/(alpha+1)) に従う。
ns = np.array([n for n in n_pis if n >= 1], dtype=float)
ts = np.array([T2s[n] for n in n_pis if n >= 1])
slope, _ = np.polyfit(np.log(ns), np.log(ts), 1)
print(f"\nfitted exponent of T2(N) vs N   : {slope:.4f}")
print(f"predicted alpha/(alpha+1)       : {alpha/(alpha+1):.4f}")
print(f"Ramsey T2* / echo T2            : {T2s[0]/T2s[1]:.4f}")

# --- フィルター関数による説明 --------------------------------------------------
print("\n理由: 系列はノイズに対する帯域通過フィルターとして働きます。")
t_fix = m_max * dt
print(f"  sequence duration t = {t_fix/US:.3f} us")
print(f"{'N':>4}{'first passband f (MHz)':>24}{'S(f) there':>14}")
print("-" * 42)
for n in n_pis:
    f_peak = f_low if n == 0 else n / (2.0 * t_fix)
    print(f"{n:>4}{f_peak/1e6:>24.4f}{A/f_peak**alpha:>14.4g}")
print("  自由誘導は直流まで感度をもち、そこは 1/f ノイズが最大の領域です。")
print("  piパルスを増やすほど通過帯域は高周波側へ移り、ノイズの小さい領域に")
print("  入ります。ただし平坦で除去できない部分に達するまでです。")

# --- 可視化 -----------------------------------------------------------------
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
for n in n_pis:
    lbl = "Ramsey" if n == 0 else f"CPMG-{n}"
    ax[0].loglog(times / US, np.maximum(curves[n], 1e-4), label=lbl)
ax[0].axhline(1 / np.e, color="gray", ls=":", lw=0.8)
ax[0].set_xlabel("t (us)"); ax[0].set_ylabel("coherence")
ax[0].set_ylim(1e-3, 1.5); ax[0].legend(fontsize=7)
ax[0].set_title("Dynamical decoupling under 1/f noise")

ax[1].loglog(ns, ts / US, "o-", color="tab:purple", label="simulation")
ax[1].loglog(ns, ts[0] / US * ns ** (alpha / (alpha + 1)), "k--", lw=1,
             label=f"$N^{{{alpha/(alpha+1):.2f}}}$")
ax[1].set_xlabel("number of pi pulses N"); ax[1].set_ylabel("T2(N) (us)")
ax[1].legend(fontsize=8); ax[1].set_title("CPMG scaling")
plt.tight_layout()
plt.show()
1/f ノイズの生成: 8192 点・dt = 2 ns の軌跡を 2000 本
  band                          : 61.0 kHz to 250 MHz
  mean sample S(f) / target S(f): 1.0008
  sample variance of delta_omega: 5.3497e+13 (rad/s)^2
  A ln(f_max/f_low)             : 4.9907e+13 (rad/s)^2
  rms detuning                  : 1164.1 kHz
  quasi-static sqrt(2)/sigma    : 0.193 us
  statistical floor 1/sqrt(n)   : 0.0224

同一の 1/f ノイズ下における CPMG-N のコヒーレンス |<exp(i phi)>|
  t (us)       N=0       N=1       N=2       N=4       N=8      N=16
--------------------------------------------------------------------
   0.052    0.9557    0.9943    0.9967    0.9980    0.9989    0.9991
   0.098    0.8672    0.9800    0.9885    0.9940    0.9967    0.9982
   0.206    0.6005    0.9155    0.9493    0.9732    0.9866    0.9930
   0.494    0.0990    0.5859    0.7557    0.8631    0.9241    0.9604
   1.042    0.0351    0.1018    0.2699    0.5231    0.7167    0.8372
   1.942    0.0147    0.0434    0.0255    0.0936    0.3374    0.5501
   4.096    0.0121    0.0356    0.0151    0.0240    0.0218    0.0691

 N (pi pulses)   T2(N) (us)   T2(N)/T2(1)
-----------------------------------------
    0 (Ramsey)       0.3041         0.445
             1       0.6829         1.000
             2       0.9032         1.323
             4       1.2780         1.872
             8       1.7891         2.620
            16       2.5266         3.700

fitted exponent of T2(N) vs N   : 0.4761
predicted alpha/(alpha+1)       : 0.5000
Ramsey T2* / echo T2            : 0.4454

理由: 系列はノイズに対する帯域通過フィルターとして働きます。
  sequence duration t = 4.096 us
   N  first passband f (MHz)    S(f) there
------------------------------------------
   0                  0.0610      9.83e+07
   1                  0.1221     4.915e+07
   2                  0.2441     2.458e+07
   4                  0.4883     1.229e+07
   8                  0.9766     6.144e+06
  16                  1.9531     3.072e+06
  自由誘導は直流まで感度をもち、そこは 1/f ノイズが最大の領域です。
  piパルスを増やすほど通過帯域は高周波側へ移り、ノイズの小さい領域に
  入ります。ただし平坦で除去できない部分に達するまでです。

注目すべき点。 生成は使う前に検証されています。標本スペクトル密度は帯域全体で目標値と0.1%以内で一致し、標本分散は $A\ln(f_\mathrm{max}/f_\mathrm{low})$ より7%高い値になります。これは残した最低周波ビンが帯域の端で寄与するためです。この検査は省略できません。規格化を誤った有色ノイズ生成器は、もっともらしく見える曲線と誤った時間定数を出力します。

コヒーレンスの表が中心的な結果です。$\pi$ パルスはどれもコヒーレンスを買い、その利得は偶然ではなく系統的です。$N = 1$ から $16$ について $T_2(N)$ を $N$ に対してフィットすると指数 0.476 が得られ、$\alpha = 1$ に対する予測 $\alpha/(\alpha+1) = 0.5$ と対応します。フィルター関数のブロックがその理由を説明します。自由誘導は帯域内の最も低い周波数まで感度をもち、そこは $S(f) = A/f$ が最大の領域です。$N$ を倍にするごとに最初の通過帯域は2倍上へ移り、そこでのノイズパワーは半分になります。このスケーリング則はフィットの便法ではなく $S(f)$ の形の直接の帰結であり、つまり指数を測ることは $\alpha$ を測ることであり、$\alpha$ は材料中の欠陥の切替速度分布についての主張です。

細部で3点が注目に値します。Ramsey時間は $0.304\ \mu$s である一方、準静的な見積り $\sqrt{2}/\sigma$ は $0.193\ \mu$s を与えます。準静的近似は位相緩和を過大評価するのです。ノイズの速い成分は系列の中で平均されてしまい、静的な位相として寄与しないからです。表の深い裾は信号ではなくノイズです。2000本の軌跡平均の統計的な床は $1/\sqrt{2000} = 0.022$ であり、その値より下で見える曲線の交差はサンプリング誤差です。そして低周波カットオフは数値上のパラメータではなく物理的なパラメータです。実際の測定の有限な継続時間を表しており、同じ量子ビットが1時間平均すると1秒平均より短い $T_2^\ast$ を与える理由です。

この道具箱の用途

5つの関数 — bloch_rhsevolverotatefree_precessioncoherence — にノイズ生成器を加えれば、標準的なコヒーレンス測定はすべて再現できます。第2章から第5章では、方式ごとの主張を検証するのにこれらを使います。

関数 導入 用途
bloch_rhsevolve Code Example 3 駆動ダイナミクス、パルス校正、飽和
rotate Code Example 5 系列中の理想パルス
free_precession Code Example 5 パルス間の高速な厳密発展
coherence Code Example 6 非対角成分のアンサンブル平均
noise_trajectoriescpmg_modulation Code Example 7 非マルコフノイズ、フィルター関数、デカップリング

この道具箱にできないことも同じくらい述べる価値があります。量子ビットは1個なので2量子ビットゲートもエンタングルメントもありません。準位は2つなのでリーケージもありません。そして $T_1$、$T_2$、$S(f)$ を入力として受け取るので、それらを予言できません。この3つをハミルトニアンと材料から導くことが、次の4章の仕事です。


演習

演習1: コヒーレンス時間の記帳

ある量子ビットで $T_1 = 40\ \mu$s、Hahnエコーで $T_2 = 25\ \mu$s と測定されました。同じ素子でのRamsey測定は $1/e$ 時間 $0.9\ \mu$s のガウス型減衰を与えました。

  1. 純位相緩和時間 $T_\varphi$ を求めてください。
  2. $T_1$ を決めている誘電損失が変わらないとして、この素子がもちうる $T_2$ の最大値はいくらですか。
  3. Ramsey減衰が量子ビット周波数の静的なガウス分布で支配されているとして、標準偏差 $\sigma/2\pi$ を kHz で求めてください。
  4. 同僚が類似素子で $T_1 = 40\ \mu$s、$T_2 = 95\ \mu$s と報告しました。測定を繰り返さずに何が言えますか。
解答

1. \(1/T_2 = 1/(2T_1) + 1/T_\varphi\) より \(1/T_\varphi = 1/25 - 1/80 = 0.04 - 0.0125 = 0.0275\ \mu\mathrm{s}^{-1}\)、すなわち \(T_\varphi = 36.4\ \mu\mathrm{s}\)。ここでは純位相緩和と緩和が同程度に寄与しています。

2. 限界は \(T_2 \le 2T_1 = 80\ \mu\mathrm{s}\) で、\(T_\varphi \to \infty\) で到達します。したがって純位相緩和を完全に除去しても3.2倍しか買えず、それ以上には \(T_1\) の改善が必要です。

3. 静的なガウス分布では \(T_{2,\mathrm{inh}}^\ast = \sqrt{2}/\sigma\) なので \(\sigma = \sqrt{2}/0.9\ \mu\mathrm{s} = 1.571\times10^{6}\) rad/s、\(\sigma/2\pi = 250\) kHz。厳密には内在的な指数減衰を無視しており、それは観測される減衰をわずかに短くするので、250 kHz はやや過大評価です。

4. それは不可能です。\(T_2 \le 2T_1 = 80\ \mu\mathrm{s}\) だからです。\(T_1\) が過小評価されたか、その「\(T_2\)」が \(T_2\) ではありません。よくある原因は CPMG-\(N\) の減衰をフィットしてHahnエコーとして報告することで、\(T_2(N)\) は \(N\) とともに増えますが \(2T_1\) には抑えられたままです。

T1, T2, T2s = 40.0, 25.0, 0.9
print(round(1 / (1 / T2 - 1 / (2 * T1)), 3))     # 36.364  T_phi(us)
print(2 * T1)                                    # 80.0    限界
import numpy as np
print(round(np.sqrt(2) / (T2s * 1e-6) / (2 * np.pi) / 1e3, 1))   # 250.1  kHz

演習2: エネルギースケールから方式を読む

遷移周波数が 400 MHz の仮想的な量子ビットを考えます。低周波の磁束量子ビット、あるいは強磁場中の核スピンのスケールです。

  1. エネルギーを $\mu$eV で、等価温度 $\hbar\omega_q/k_B$ を mK で求めてください。
  2. ベース温度 20 mK における熱励起状態占有数はいくらですか。
  3. 熱的初期化は成立しますか。代わりに何をしますか。
  4. 同僚が「$S(\omega_q)$ が周波数とともに下がるので $T_1$ が改善する」と主張して周波数を 40 MHz に下げることを提案しました。賛成と反対の論拠を1つずつ挙げてください。
解答

1. \(E = h \times 4\times10^{8} = 2.65\times10^{-25}\) J = \(1.654\ \mu\mathrm{eV}\)、\(hf/k_B = 19.20\) mK。

2. \(x = hf/k_BT = 19.2/20 = 0.960\) として、2準位のBoltzmann占有数は \(e^{-x}/(1+e^{-x}) = 0.2769\)。占有数の4分の1以上が誤った状態にあります。

3. 成立しません。熱的初期化は \(k_BT \ll \hbar\omega_q\)、すなわちここでは 2 mK を十分下回るベース温度を要求し、それは利用できません。代替は能動リセットです。測定して条件付きで反転させるか、補助準位を通じて散逸的な遷移を駆動します。低周波の磁束量子ビットや核スピンの初期化が難しいのと同じ論理であり、能動リセットが任意選択ではなく標準である理由です。

4. 賛成: 環境が \(1/f\) 的あるいは単調減少のスペクトル密度をもつなら、低い \(\omega_q\) は確かにノイズパワーの小さい場所にあり、黄金律の \(T_1\) は改善します。反対: 初期化はさらに悪化します(40 MHz では \(hf/k_B = 1.9\) mK なので、到達可能などんな温度でも実質的に完全混合です)。熱光子が信号を埋めるため読み出しコントラストも劣化し、回転座標系の記述が成立するには \(\Omega \ll \omega_q\) が必要なのでゲート時間も伸びます。コヒーレンスは唯一の軸ではない、というのが第1.3節の教訓です。

import numpy as np
h, kB = 6.62607015e-34, 1.380649e-23
for f in (400e6, 40e6):
    x = h * f / (kB * 0.020)
    print(f"{f/1e6:5.0f} MHz  {h*f/1.602176634e-19*1e6:6.3f} ueV  "
          f"{h*f/kB*1e3:6.2f} mK  P_exc = {np.exp(-x)/(1+np.exp(-x)):.4f}")
# 400 MHz   1.654 ueV   19.20 mK  P_exc = 0.2769
#  40 MHz   0.165 ueV    1.92 mK  P_exc = 0.4760

演習3: ゲート予算

2つの方式を考えます。Aは $T_2 = 200\ \mu$s で2量子ビットゲート時間 40 ns、Bは $T_2 = 2$ s で2量子ビットゲート時間 200 $\mu$s です。

  1. 各方式について、デコヒーレンスで決まる2量子ビットゲートあたりの誤りを $\epsilon \approx t_\mathrm{gate}/T_2$ で見積もってください。
  2. 累積誤差が1に達するまでに何ゲート実行できますか。
  3. 表面符号の閾値がおおよそ $\epsilon < 10^{-3}$ を要求するとします。デコヒーレンスだけで見ればどちらかは通りますか。$\epsilon = 10^{-4}$ に達するには $T_2/t_\mathrm{gate}$ を何倍改善する必要がありますか。
  4. $T_2$ を10倍にするのと、ゲートを10倍速くするのでは、物理的にどちらが容易ですか。トランズモンとイオントラップについて別々に答えてください。
解答

1. A: \(\epsilon \approx 40\ \mathrm{ns}/200\ \mu\mathrm{s} = 2\times10^{-4}\)。B: \(\epsilon \approx 200\ \mu\mathrm{s}/2\ \mathrm{s} = 1\times10^{-4}\)。入力がどちらも4桁違うのに結果は2倍以内であり、これが無次元比の要点です。

2. \(1/\epsilon\) で、Aは5000ゲート、Bは10000ゲートです。

3. デコヒーレンスだけで見れば両方が \(10^{-3}\) を通ります。だからこそ実際にはどちらの方式でも誤りはデコヒーレンス支配ではなく、校正誤差、クロストーク、リーケージが支配します。この予算で \(10^{-4}\) に達するには A に2倍、B には改善不要です。誠実な結論は、この数値ではどちらの方式でも \(T_2/t_\mathrm{gate}\) が律速条件ではなく、誤差予算はこの方法で見積もるのではなく測定しなければならないということです。

4. トランズモン: \(T_2\) を10倍にするとは誘電損失または \(1/f\) 振幅を10分の1にすることで、困難だが目標が明確な材料プログラムです。ゲートを10倍速くする方は非調和性に突き当たります。おおよそ \(1/\alpha\) より短いパルスは第3準位を占有するので、速度は増幅器ではなく回路のスペクトルに抑えられます。イオントラップ: \(T_2\) はすでに長く磁場安定度に律速されているので、シールドとクロック遷移でさらに稼げます。ゲートの速度は運動周波数と、分解サイドバンド領域に留まる要求に抑えられ、それを超えると非共鳴励起で忠実度を失います。どちらの場合も速度の軸は系のスペクトルスケールで上から抑えられており、それが労力の大半がコヒーレンスの軸に注がれる理由です。

演習4: フィルター関数

定義 $\tilde{s}(f,t) = \int_0^t s(t')e^{-2\pi i f t'}dt'$ を用い、自由誘導では $s = +1$ のまま、エコーでは $s = +1$ のあと $-1$ とします。

  1. $t = 1$、$f = 0.03$ において $|\tilde{s}_\mathrm{FID}|^2 = \sin^2(\pi f t)/(\pi f)^2$ と $|\tilde{s}_\mathrm{echo}|^2 = 4\sin^4(\pi f t/2)/(\pi f)^2$ を数値的に確認してください。
  2. 両者を $f t \ll 1$ で展開し、それぞれの $f$ の最低次を述べてください。
  3. 2の結果から、エコーが $1/f$ ノイズに対して絶大に効き、白色ノイズに対してほとんど効かない理由を説明してください。
  4. $S(f) = A/f^\alpha$ に対してCPMGのコヒーレンス時間は $T_2(N) \propto N^{\alpha/(\alpha+1)}$ とスケールします。白色ノイズについてこれは何を予言し、3の答えと整合しますか。
解答

1. どちらも8桁まで一致します。下のコードを参照してください。この値では \(|\tilde{s}_\mathrm{FID}|^2 = 0.99704262\)、\(|\tilde{s}_\mathrm{echo}|^2 = 0.00221738\) で、\(ft = 0.03\) の時点でエコーはすでに450分の1の感度です。

2. FIDでは \(\sin^2(\pi f t) \approx (\pi f t)^2\) なので \(|\tilde{s}|^2 \to t^2\)。\(f\) に依存しない定数で、直流で最大です。エコーでは \(\sin^4(\pi f t/2) \approx (\pi f t/2)^4\) なので \(|\tilde{s}|^2 \to \pi^2 f^2 t^4/4\)。\(f^2\) で消えます。

3. 位相の平均二乗は \(\int_0^\infty S(f)|\tilde{s}(f)|^2 df\) です。\(S \propto 1/f\) のとき FID の被積分関数は \(1/f\) となり小さい \(f\) で対数発散します。つまり被害はすべて最も遅い成分から来るのであり、そこはまさにエコーのフィルターが零点をもつ場所です。エコーの被積分関数は \(f\) に比例し、有限かつ小さい。白色ノイズでは \(S\) が定数で、FIDの被積分関数は \(\sin^2(\pi f t)/(\pi f)^2\) となり低周波の増強なしにすでに可積分です。エコーが除去すべき低周波の重みが存在せず、したがって何も得られません。

4. 白色ノイズは \(\alpha = 0\) で、\(T_2(N) \propto N^0\)。何発 \(\pi\) パルスを打っても利得はありません。これはまさに3で述べたことであり、動的デカップリングが \(T_2\) を \(2T_1\) より伸ばせない理由でもあります。緩和チャネルは、系列が届く周波数では実効的に白色だからです。

import numpy as np
t, f, N = 1.0, 0.03, 200000
tt = np.linspace(0, t, N + 1)
for name, s, ana in [
        ("FID ", np.ones_like(tt), np.sin(np.pi*f*t)**2 / (np.pi*f)**2),
        ("echo", np.where(tt < t/2, 1.0, -1.0),
         4*np.sin(np.pi*f*t/2)**4 / (np.pi*f)**2)]:
    num = abs(np.trapezoid(s * np.exp(-2j*np.pi*f*tt), tt))**2
    print(name, f"numeric {num:.8f}   analytic {ana:.8f}")
# FID  numeric 0.99704262   analytic 0.99704262
# echo numeric 0.00221738   analytic 0.00221738

演習5: コヒーレンスデータから材料を診断する

ある誘電体層の成膜条件だけを変え、他は同一に処理した2枚のウェハがあります。両方で量子ビットを測定しました。

ウェハ $T_1$ $T_2$(エコー) $T_2^\ast$(Ramsey、ガウス)
P 60 $\mu$s 80 $\mu$s 1.2 $\mu$s
Q 61 $\mu$s 22 $\mu$s 1.1 $\mu$s
  1. 各ウェハについて $T_2$ が緩和律速か位相緩和律速かを判定し、$T_\varphi$ を与えてください。
  2. プロセス変更が環境の揺動成分に効いたのはどちらで、静的成分に効いたのはどちらですか。表から根拠を述べてください。
  3. 各ウェハの $\sigma/2\pi$ を見積もってください。両者に有意差はありますか。
  4. 次に物理的にどこを調べますか。決着をつける追加測定を1つ挙げてください。
解答

1. ウェハP: \(2T_1 = 120\ \mu\mathrm{s}\) に対し \(T_2 = 80\ \mu\mathrm{s}\) なので \(1/T_\varphi = 1/80 - 1/120 = 0.004167\)、\(T_\varphi = 240\ \mu\mathrm{s}\)。純位相緩和は小さな補正で、\(T_2\) はおおむね緩和律速です。ウェハQ: \(2T_1 = 122\ \mu\mathrm{s}\)、\(T_2 = 22\ \mu\mathrm{s}\) なので \(1/T_\varphi = 1/22 - 1/122 = 0.03726\)、\(T_\varphi = 26.8\ \mu\mathrm{s}\)。強く位相緩和律速です。

2. \(T_1\) が変わっていないので、量子ビット周波数における損失チャネルの密度は動いていません。\(T_2^\ast\) も変わっていないので、量子ビット周波数の静的なばらつきも動いていません。変わったのは \(T_2\)、すなわち \(T_\varphi\) であり、9倍近い変化です。したがってプロセス変更は、エコーが除去しない中間帯域にノイズパワーを加えました。切替速度が \(1/T_2\) と同程度、おおよそ数十kHzの揺動子です。静的な乱れでもGHzの損失でもなく、揺動成分に効いたのです。

3. \(\sigma = \sqrt{2}/T_2^\ast\) より、ウェハPは \(1.18\times10^{6}\) rad/s すなわち 188 kHz、ウェハQは \(1.29\times10^{6}\) rad/s すなわち 205 kHz。9%の差で、ウェハあたり数個の素子であればほぼ確実に有意ではありません。静的な乱れは同じです。

4. 疑わしいのは新しい成膜が導入した2準位揺動子の集団で、数十kHzで活性なものです。決定的な測定はCPMGノイズ分光です。\(N\) を掃引し、両ウェハについて \(T_2(N)\) 曲線の族を反転して \(S(f)\) を求めます。ウェハQが 10-100 kHz 帯に過剰な \(S(f)\) を、しかも熱活性化集団に対応する \(1/f^\alpha\) 指数で示せば診断が確定し、フィットした振幅 \(A\) が次の成膜条件試験の性能指標になります。より安価な第二の検査は \(T_\varphi\) の温度依存性です。熱活性化揺動子は凍結するからです。

for name, T1, T2, T2s in [("P", 60.0, 80.0, 1.2), ("Q", 61.0, 22.0, 1.1)]:
    Tphi = 1 / (1 / T2 - 1 / (2 * T1))
    print(f"{name}: 2T1 = {2*T1:6.1f} us   T_phi = {Tphi:7.1f} us"
          f"   sigma/2pi = {2**0.5 / (T2s*1e-6) / (2*3.141592653589793) / 1e3:5.0f} kHz")
# P: 2T1 =  120.0 us   T_phi =   240.0 us   sigma/2pi =   188 kHz
# Q: 2T1 =  122.0 us   T_phi =    26.8 us   sigma/2pi =   205 kHz

まとめ

要点

1. 良い量子ビットは相反する2つのものでなければならない

2. DiVincenzo基準はチェックリストであって採点表ではない

3. 6本の軸、そして単独の勝者はいない

4. 3つの時間定数、1つのスペクトル密度

5. $1/f$ ノイズは欠陥の集団である

6. これらの測定は材料の特性評価である

実務上の含意

続く3章は主要な3方式を順に扱い、それぞれをノイズチャネルへの解答として構成します。第2章は超伝導量子ビットから始めます。ジョセフソン接合がLC回路を人工原子に変える仕組み、トランズモンが非調和性と電荷ノイズ耐性を交換する理由、そして最良の回路のコヒーレンスが数ナノメートルのアモルファス酸化膜で決まっていること — この主題が材料科学と区別できなくなる地点です。

← シリーズトップ 第2章: 超伝導量子ビット →

免責事項