タグ: オープンソース

量子系の計算で「十分に時間がたった後」を知りたいとき、QuTiPのsteadystateは便利な入口です。時間発展を長く追わずに、変化しない密度行列を求められます。ただし、返ってきた一つの行列だけでは、どの初期状態からもそこへ落ち着くことや、手元の計算時間で十分に近づいたことまでは分かりません。

見分けたいのは、定常解があること、定常解が一つであること、その解へ十分に近づいたことです。QuTiP 5.1の公式資料と、版を固定した5.1.1の実装を手掛かりに、二準位の小さな例で確認の順序を整理します。^[1][2]

簡単にまとめると？

• 定常状態は、その状態から出発すると変化しない状態です。初期状態の情報が残り、複数の定常解が存在する場合もあります。
• 位相だけを失う二準位では占有確率が残ります。上下の遷移を入れた今回の例では占有確率が一つに定まりますが、近づく速さは成分によって違います。
• 密度行列の正規化と正値性、定常方程式の残差、零空間、異なる初期状態からの時間発展を、役割を分けて確かめます。

「動かない状態を求める」と「そこへ到達する」は別の計算

密度行列ρは、各状態にいる確率と、状態間の位相の関係をまとめたものです。時間に依存しないLindblad方程式を、dρ/dt＝L(ρ)と書きます。Lはリウヴィリアンと呼ばれ、ハミルトニアンによる変化と、環境に由来する変化をまとめた作用です。定常状態ρ_ssでは、その変化がゼロになります。^[1][3]

L(ρ_ss)＝0、tr(ρ_ss)＝1

前の式は「今から変わらない」、後の式は「確率の合計が1」という条件です。物理的な密度行列であるためには、エルミート性に加え、固有値が負にならない正値性も必要です。

QuTiPの通常のdirect法は、定常方程式へトレースを1にする条件を加えて線形方程式を解きます。この計算には、時間発展の出発点となる初期状態を指定しません。一方、時間発展を計算するmesolveなどでは初期状態が入力です。前者は「変化しない解」を探し、後者は「この出発点からどこへ進むか」を追っています。^[1][2]

位相だけを失うと、最初の占有確率が残る

基底状態を|0⟩、励起状態を|1⟩とする二準位を考えます。占有確率とは、その状態にいる確率のことです。ここでは励起状態の占有確率p₁＝ρ₁₁を使います。ハミルトニアンをH＝0とし、位相の関係だけを弱める純粋位相緩和を入れます。

C_φ＝√(γ_φ/2) Z、Z＝|0⟩⟨0|−|1⟩⟨1|、γ_φ＞0

Cは環境の作用を表す崩壊演算子です。この設定をLindblad項へ代入すると、対角成分である占有確率は変わらず、非対角成分はρ₀₁(t)＝ρ₀₁(0)e^−γ_φtで減少します。つまり、位相の関係は薄れても、|0⟩と|1⟩の間で占有確率を移す仕組みがありません。^[3]

そのため、0≦p≦1を満たすどのρ＝diag(1−p, p)も定常解です。|0⟩から始めればp＝0、|1⟩から始めればp＝1のままです。両者とも残差L(ρ)はゼロですが、同じ状態にはなりません。|+⟩＝(|0⟩＋|1⟩)/√2から始めた場合は、占有確率1/2を保ったまま位相の成分が消え、diag(1/2, 1/2)へ近づきます。

この例では、トレース1の条件を加えても占有確率を一つに決められません。定常方程式の解空間に自由度が残るためです。ソルバーが特異性を知らせたり、手法によって解の一つが得られたりしても、それを特定の初期状態からの長時間後の答えと同一視しないことが大切です。

上下の遷移がある例では、収束先と速さを分けて見る

次に、H＝0のまま、|1⟩から|0⟩へ下がる過程と、逆向きに上がる過程を入れます。ここからは純粋位相緩和だけのモデルを置き換えた、別の設定です。

C_↓＝√γ_↓ |0⟩⟨1|、C_↑＝√γ_↑ |1⟩⟨0|、γ_↓, γ_↑＞0

占有確率の変化は、入ってくる分と出ていく分の差になります。

dp₁/dt＝γ_↑(1−p₁)−γ_↓p₁

変化がゼロになる占有確率は、p_1,ss＝γ_↑/(γ_↑＋γ_↓)です。このモデルでは非対角成分も減衰するので、密度行列全体の定常解が一つに定まります。状態間を遷移する経路ができたことで、最初の占有確率を保持する自由度がなくなりました。

具体例としてγ_↓＝3 s⁻¹、γ_↑＝1 s⁻¹を仮定すると、定常占有確率は1/4です。τ＝(γ_↓＋γ_↑)t＝(4 s⁻¹)tを無次元の経過時間とすると、

p₁(t)＝1/4＋[p₁(0)−1/4]e^−τ、ρ₀₁(t)＝ρ₀₁(0)e^−τ/2

となります。下の表はこの仮定モデルから計算した値です。演算子の式から別に組み立てた4×4のリウヴィリアンと行列指数でも照合しました。実験値やQuTiP本体の実行結果ではありません。

τ＝5では、二つの初期状態の占有確率はすでに1/4へかなり近づいています。それでも、初期|+⟩に残る位相の成分は約0.0410です。この成分は占有確率のずれより遅い指数で減衰します。励起確率だけを見て止めてよいか、密度行列全体の誤差を小さくしたいのかで、必要な計算時間は変わります。

計算結果は、四つの問いで点検する

最初は、物理的な状態になっているか。トレース、エルミート性、固有値を確認します。丸め誤差に由来する微小な負値と、大きな負値を同じ扱いにせず、許容誤差を記録します。正規化し直しただけで、不適切な解や一意性の問題が解消するわけではありません。

次は、その状態が定常方程式を満たしているか。L(ρ)のノルムで残差を見ます。残差はモデル内の変化率なので、時間の単位やLのスケールをそろえて比較します。すべての速度を小さくしただけでも生の残差は小さくなります。小さな残差は有用な確認ですが、先ほどの純粋位相緩和の例のように、一意性は別に調べる必要があります。

三つ目は、解を一つに決めるだけの条件があるか。小さなモデルなら、Lの零空間、つまりLがゼロへ写す行列の集合を調べられます。今回の独立計算では、その次元は純粋位相緩和で2、上下遷移で1でした。トレース1という条件を加えた後に自由度が残るかが違います。数値上ゼロとみなすしきい値を変えても結果が安定しているかを確認し、大規模・悪条件の行列へ小例の精度をそのまま持ち込みません。

最後は、必要な精度まで近づいたか。異なる占有確率から出発させ、位相も重要なら重ね合わせ状態も加えて、時間を延ばした結果を比較します。この例ではLの零でない固有値は−4 s⁻¹と−2 s⁻¹で、位相の成分のほうが遅く減衰します。一般のモデルでも遅い成分の時間尺度は手掛かりになりますが、いくつかの初期状態で一致したというだけで、すべての初期状態についての証明にはしません。

QuTiPの停止条件を、物理的な一意性と混同しない

固定した5.1.1の実装には、method='propagator'もあります。初期密度行列rhoへ時間発展を繰り返し作用させ、反復前後の状態の距離がpropagator_tol以下になると返す方法です。rhoを省略した場合はランダムな密度行列から始め、収束しない場合は反復上限で例外を返します。^[2]

ここで判定しているのは、選んだ出発点からの反復がどれだけ変化したかです。純粋位相緩和のように定常解が複数ある系では、到達する占有確率は初期状態に依存します。また、モデルの変化がとても遅ければ、短い間隔での変化が小さいことと、最終状態に近いことは別になります。初期状態、時間間隔、許容誤差を変えて確かめる理由はここにあります。

一意な定常状態が欲しいのか、保存される初期情報を調べたいのかで、計算結果の読み方は変わります。まず環境との相互作用が何を変え、何を残すモデルなのかを確認し、その上でソルバーの出力と時間発展を照合すると、定常状態の計算を目的に結び付けやすくなります。

参考資料

1. QuTiP 5.1 Documentation, Solving for Steady-State Solutions。定常方程式、計算法、時間発展との比較例。
2. QuTiP, steadystate API、5.1.1の固定実装。トレース条件、各手法、propagatorの初期状態と停止条件。参照版を固定しており、現在の最新版という意味ではありません。
3. QuTiP, 5.1.1のliouvillian / lindblad_dissipator実装。本文の二準位モデル、解の導出と表の数値は、このLindblad形式を用いた独自の説明・計算例です。

量子ビットを0と1で測ると、どちらも半分ずつ出る。しばらく待っても、その割合は変わらない。それだけで、量子状態が保たれていると言えるでしょうか。0と1の重ね合わせでは、結果の割合に加えて、二つの成分の位相の関係も大切です。

Open Quantum Design（OQD）の公開エミュレーターを調べると、確認した版はQuTiPのsesolveで状態を時間発展させています。そこには、環境による緩和を指定する入力はありません。本稿ではその境界を出発点に、QuTiPが別に備えるLindblad（リンドブラッド）方程式の機能から、モデルに何を加えると何が変わるのかを見ます。OQD自体がこの機能を実装しているという説明ではありません。^［1］

簡単にまとめると？

• エネルギー緩和と位相の乱れは別の変化です。0・1の割合だけを見ても、重ね合わせの変化を捉えきれません。
• Lindblad方程式では、ハミルトニアンに加えて、環境から受ける作用とその速さを指定します。
• 計算できたことと実機を再現できたことは別です。作用の種類、速度の単位、近似の適用条件を確認します。

閉じた系の計算には、指定していない緩和は現れない

ハミルトニアンは、量子状態が時間とともにどう変わるかを定めるものです。外部環境とのやり取りを含めずに扱う「閉じた系」では、状態の変化をユニタリーな時間発展として計算します。位相が回転したり、駆動によって0と1の割合が変わったりしても、それ自体が環境による緩和ではありません。^［2］

一方、対象の量子ビットが周囲へエネルギーを失ったり、周囲の影響で位相の関係が乱れたりするなら、その作用をモデルに加える必要があります。Lindblad方程式は、そのような「開いた系」を扱う代表的な形です。状態を密度行列で表し、ハミルトニアンによる変化と、環境による変化を合わせます。^［2］

密度行列は、結果の割合と重ね合わせの関係をまとめた表と考えると読みやすくなります。対角の成分は0・1それぞれの確率、対角でない成分は両者のコヒーレンスを表します。以下では、その非対角成分の大きさを比べます。位相が回転することと、その大きさが減ることは分けて扱います。

同じ重ね合わせを、三つのモデルで待たせる

説明用に、0と1を同じ重みで重ねた状態を用意します。最初は1の確率が0.5、コヒーレンスの大きさも0.5です。0を基底状態、1を励起状態とし、0・1を入れ替える外部の駆動は加えません。0と1のエネルギー差は一定とし、その差による位相回転の速さには依らない、非対角成分の大きさを比較します。

比べるのは、環境の作用なし、1から0へ戻るエネルギー緩和のみ、そして0・1の確率を変えない純粋位相緩和のみ、の三通りです。ここでTは説明用の時間尺度であり、実機の寿命ではありません。後二つのモデルでは、それぞれの速度を1/Tと置きます。Lindbladの散逸項から得た式を、独立に作成した2行2列の行列計算とも照合しました。^［4］

狭い画面では表を横にスクロールできます。QuTiPや量子装置を実行した性能評価ではなく、指定した方程式の独立検算です。

エネルギー緩和の行では、1の確率は0.5×exp(−t/T)、コヒーレンスの大きさは0.5×exp(−t/2T)で減ります。純粋位相緩和の行では、1の確率は0.5のまま、コヒーレンスだけが0.5×exp(−t/T)で減ります。

この違いが、何を調べるべきかを変えます。0・1の割合が同じでも、純粋位相緩和がないとは言えません。逆に、1の割合が減るデータを、確率を変えない純粋位相緩和だけで説明することもできません。0・1とは異なる重ね合わせの基底（X、Y）でも測定すると、非対角成分の実部・虚部を調べられます。実機と照合するときは、割合とコヒーレンスに関わる情報を区別し、どの作用を仮定したモデルなのかを揃えます。

c_opsには、速さだけでなく「何を変えるか」を渡す

QuTiPのmesolveでは、環境の作用をc_opsという引数へ渡します。各要素はcollapse operatorと呼ばれ、作用を表す演算子に、速度の平方根を掛けたものです。この名前は、途中測定の結果を乱数で一つ選ぶ処理そのものを意味しません。ここでは密度行列の平均的な変化を計算します。^{［2］［3］}

先ほどのエネルギー緩和は、1を0へ移す演算子に√γを掛けます。γは緩和の速度です。純粋位相緩和では、0と1の成分へそれぞれ+1と−1を掛ける演算子Zを使います。これに√(γφ/2)を掛けると、非対角成分の大きさが速度γφで減ります。平方根や2分の1は単なる書き方の違いではなく、計算される減衰速度を変える係数です。^［4］

時間を秒で表すならγとγφは毎秒の単位、たとえばs⁻¹で揃えます。ハミルトニアンも、ℏを1とする表記なら角周波数の単位で揃えます。表のTを特定装置の測定値へ置き換えるには、その装置と条件に対応した速度を別に確かめる必要があります。

エネルギー緩和だけでも、確率とコヒーレンスの減る速さは異なります。この例で、励起確率の減衰時間をT1と呼べばT1＝1/γ、コヒーレンスの減衰時間T2は2/γです。純粋位相緩和も同時に加えると、1/T2＝γ/2＋γφになります。T1とT2を同じ値だと置く前に、何が減る時間を測ったのかを確認したいところです。これは今の二準位モデルから得た関係で、すべての装置への無条件な当てはめではありません。

ソルバーの名前より、入力と近似を確かめる

確認対象はQuTiP 5.1の公式文書と、対応する5.1.1の固定コードです。最新版を評価した記事ではありません。固定実装では、c_opsが空、初期状態が状態ベクトル、ハミルトニアンがsuperoperatorではない、という条件がそろうと、mesolveからsesolveへ処理を渡します。初期状態が密度行列なら、環境の作用が空でも密度行列側の計算経路を使います。関数名だけで「緩和を含めた」とは判断できません。^［3］

さらに、Lindblad形で数値が出ることと、現実の環境にその近似が合うことは別です。公式解説は、初期に系と環境の相関を置かないこと、系と環境の弱い結合を扱うBorn近似、環境の記憶が対象の変化より短いMarkov近似、速く振動する項を除く近似などを説明しています。強い結合や長い記憶を持つ環境へ、同じ単純なモデルをそのまま広げることはできません。^［2］

まず知りたいのは、結果の割合が変わるのか、位相の関係が失われるのか。その問いに合わせて作用と速度を指定し、適用条件を確かめる。そうして初めて、理想的な時間発展と、環境を含む計算の違いが、量子ビットを扱う判断材料になります。

参考資料

資料は2026年10月14日に確認。OQDは固定commit 6a086a9843b91d041e76fbd4fed83c5384a14e6e、QuTiPはタグv5.1.1に対応する固定commit 7accc77a70a08f965a29fdd13515f79cdcab9703を参照しています。本文の表と数値はQuantum Seaによる説明用の独立計算です。

1. Open Quantum Design, oqd-analog-emulator。conversion.pyのQutipOperation処理、interface.py、pyproject.toml。OQDの時間発展処理とQuTiP 5.1.1への依存を参照。固定版の時間発展処理／依存版の指定。
2. QuTiP 5.1 Documentation「Lindblad Master Equation Solver」。Unitary evolution、Non-unitary evolution、The Lindblad Master equation。状態・環境作用・導出の近似条件を参照。公式解説。
3. QuTiP 5.1.1、qutip/solver/mesolve.py。c_opsと初期状態による計算経路の分岐。固定版の実装。
4. QuTiP 5.1.1、qutip/core/superoperator.pyのlindblad_dissipatorとqutip/core/operators.pyのsigmax・sigmay・sigmaz。CρC†−(C†Cρ＋ρC†C)/2の形を本稿の2行2列の計算に使用。固定版の散逸項／観測量の演算子／コードの利用条件。