qutip

QuTiP - 量子ダイナミクス

QuTiP は、量子力学のための柔軟で効率的なフレームワークとして設計されています。量子オブジェクトを簡単に作成でき、閉じた環境と開いた環境の両方でその時間発展を追跡する強力なソルバーを備えています。

使用する場面

• 量子系の時間発展をシミュレーションする場合（シュレーディンガー方程式またはマスター方程式）
• 開放量子系の定常状態を計算する場合（例：駆動と散逸がある共振器）
• エンタングルメント、ウィグナー関数、量子相関を解析する場合
• 光と物質の相互作用を研究する場合（Jaynes-Cummings モデル、ラビ振動）
• パルスシーケンスの最適化と量子制御
• スペクトルと多時間相関関数を計算する場合
• 複数の CPU にわたって量子シミュレーションを並列化する場合

リファレンスドキュメント

公式ドキュメント: https://qutip.org/docs/latest/
チュートリアル: https://qutip.org/tutorials.html
検索パターン: qutip.Qobj、qutip.mesolve、qutip.wigner、qutip.expect、qutip.basis

基本原則

Qobj（量子オブジェクト）

中心となるデータ構造です。Qobj は状態（ケットまたはブラ）、演算子、または超演算子を表現できます。次元を自動的に処理し、演算の妥当性を検証します（例：ケットと演算子の加算を防ぎます）。

複合系

QuTiP は、テンソル積（クロネッカー積）を使って、複数の部分系から構成される系（例：共振器に結合した量子ビット）を表現します。

ソルバーのワークフロー

1. ハミルトニアン（H）を定義します
2. 散逸を表す崩壊演算子（C_n）を定義します
3. 初期状態（ρ₀）を設定します
4. 時間系列（t）を定義します
5. ソルバー（mesolve、sesolve など）を実行します

クイックリファレンス

インストール

pip install qutip

標準インポート

import qutip as qt
import numpy as np
import matplotlib.pyplot as plt

基本パターン - ラビ振動

import qutip as qt
import numpy as np

# 1. Define Hamiltonian: H = h_bar * omega * sigma_x / 2
omega = 1.0 * 2 * np.pi
H = 0.5 * omega * qt.sigmax()

# 2. Initial state: Spin down (|1>)
psi0 = qt.basis(2, 1)

# 3. Time points
times = np.linspace(0.0, 10.0, 100)

# 4. Solve Schrödinger Equation
result = qt.sesolve(H, psi0, times, [qt.sigmaz()])

# 5. Get expectation values
sz_expt = result.expect[0]

重要なルール

✅ 推奨事項

• basis(N, i) を使う - 配列を手動で作成する代わりに、組み込み関数を使って状態を作成します
• Qobj の次元を確認する - 複合系を扱うときは、テンソル積が正しいことを確かめるため、常に .dims を確認します
• expect() を使う - ソルバーに期待値を直接計算させ、メモリと時間を節約します
• 開放系には c_ops を指定する - 崩壊演算子を指定した mesolve を使って、デコヒーレンスをモデル化します
• ハミルトニアンをベクトル化する - ハミルトニアンが時間依存の場合は、速度向上のため [H0, [H1, 'sin(w*t)']] のリスト形式を使います
• ウィグナー関数で可視化する - qt.wigner(state, xvec, yvec) を使って、状態の量子性を詳しく調べます

❌ 非推奨事項

• データを手動で取り出さない - 絶対に必要な場合を除き、state.data.toarray() は避けます。QuTiP の Qobj は疎行列演算に最適化されています
• ソルバーの警告を無視しない - 数値的不安定性は、時間ステップが大きすぎるか、ヒルベルト空間の打ち切りが早すぎることを示している場合がよくあります
• ヒルベルト空間を混在させない - 適切なテンソル積を使わずに、次元の異なるオブジェクト間で演算しないでください
• 打ち切りを忘れない - 振動子などの無限系では、収束するのに十分なフォック空間のサイズ（N）を確保します

アンチパターン（絶対に避けること）

import qutip as qt

# ❌ BAD: Creating states using raw numpy arrays (no metadata)
psi = qt.Qobj([[1], [0]])

# ✅ GOOD: Use built-in constructors
psi = qt.basis(2, 0) # Ket |0>

# ❌ BAD: Manual Kronecker product
H_total = np.kron(H_atom, np.eye(N_cav))

# ✅ GOOD: Use qt.tensor
H_total = qt.tensor(H_atom, qt.qeye(N_cav))

# ❌ BAD: Calculating expectation values in a Python loop
# (Extremely slow and memory intensive)

# ✅ GOOD: Pass operators to the solver
result = qt.mesolve(H, psi0, times, c_ops, [qt.sigmaz(), qt.sigmax()])

量子オブジェクト（qt.Qobj）

作成とプロパティ

# Kets and Bras
psi = qt.basis(5, 0) # Ket |0> in 5-level system
bra = psi.dag()      # Adjoint (Bra)

# Operators
a = qt.destroy(10) # Annihilation operator (size 10)
n = a.dag() * a    # Number operator
sz = qt.sigmaz()   # Pauli Z

# Properties
print(psi.type)    # 'ket'
print(sz.dims)     # [[2], [2]]
print(n.eigenstates()) # Get eigenvalues and vectors

複合系（テンソル）

# Qubit (2) + Cavity (10)
sz = qt.sigmaz()
a = qt.destroy(10)

# Operator acting on qubit only
sz_total = qt.tensor(sz, qt.qeye(10))

# Interaction term: sigma_x * (a + a_dag)
H_int = qt.tensor(qt.sigmax(), a + a.dag())

ダイナミクスソルバー（qt.mesolve、qt.sesolve）

マスター方程式（開放系）

# Parameters
kappa = 0.1 # Cavity decay rate
H = qt.tensor(qt.sigmax(), qt.qeye(10)) # Hamiltonian
c_ops = [np.sqrt(kappa) * qt.tensor(qt.qeye(2), a)] # Decay list
psi0 = qt.tensor(qt.basis(2, 0), qt.basis(10, 0))
times = np.linspace(0, 50, 200)

# Solve
result = qt.mesolve(H, psi0, times, c_ops, [qt.tensor(qt.sigmaz(), qt.qeye(10))])

# Access expectation values
exp_z = result.expect[0]

時間依存ハミルトニアン

# Using string-based syntax (fastest)
H0 = qt.sigmax()
H1 = qt.sigmaz()
args = {'w': 0.5}
H_t = [H0, [H1, 'sin(w * t)']]

result = qt.sesolve(H_t, qt.basis(2, 0), times, args=args)

定常状態と相関

# Steady state of a driven-dissipative system
rho_ss = qt.steadystate(H, c_ops)

# Expectation value in steady state
avg_n = qt.expect(a.dag() * a, rho_ss)

# Emission Spectrum
w_list = np.linspace(0, 10, 100)
spec = qt.spectrum(H, w_list, c_ops, a.dag(), a)

可視化

ウィグナー関数とブロッホ球

# Wigner function for a Cat state
cat = (qt.coherent(20, 2) + qt.coherent(20, -2)).unit()
xvec = np.linspace(-5, 5, 100)
W = qt.wigner(cat, xvec, xvec)

# Bloch Sphere
b = qt.Bloch()
b.add_states(qt.basis(2, 0))
b.add_states((qt.basis(2,0) + qt.basis(2,1)).unit())
# b.show()

実践的なワークフロー

1. Jaynes-Cummings モデル（量子ビットと共振器）

def simulate_jc(g, kappa, n_max, times):
    # Operators
    a = qt.tensor(qt.qeye(2), qt.destroy(n_max))
    sm = qt.tensor(qt.destroy(2), qt.qeye(n_max))
    
    # Hamiltonian
    H = g * (sm.dag() * a + sm * a.dag())
    
    # Dissipation
    c_ops = [np.sqrt(kappa) * a]
    
    # Initial: Excited atom, 0 photons
    psi0 = qt.tensor(qt.basis(2, 0), qt.basis(n_max, 0))
    
    res = qt.mesolve(H, psi0, times, c_ops, [sm.dag()*sm, a.dag()*a])
    return res.expect

2. 量子状態トモグラフィー解析

def fidelity_check(state_ideal, state_measured):
    """Calculates fidelity between two quantum states."""
    # Works for both kets and density matrices
    return qt.fidelity(state_ideal, state_measured)

3. G2 相関関数の計算

def calculate_g2(H, c_ops, a_op, times):
    """Second-order correlation function g2(tau)."""
    # Use built-in correlation function
    g2_tau = qt.correlation_2op_1t(H, None, times, c_ops, a_op.dag(), a_op)
    # Normalize by steady state intensity
    rho_ss = qt.steadystate(H, c_ops)
    n_ss = qt.expect(a_op.dag() * a_op, rho_ss)
    return np.real(g2_tau) / (n_ss**2)

パフォーマンスの最適化

疎行列と密行列

QuTiP はデフォルトで Scipy の疎行列を使います。非常に小さな系（N < 4）では密行列への変換のほうが高速な場合がありますが、ほとんどの場合、N > 10 では疎行列が必須です。

並列マップ

パラメータスイープには qt.parallel_map を使います。

def task(omega):
    H = 0.5 * omega * qt.sigmax()
    return qt.mesolve(H, psi0, times, c_ops, [qt.sigmaz()]).expect[0]

omegas = np.linspace(0, 10, 20)
results = qt.parallel_map(task, omegas)

よくある落とし穴と解決策

「次元の不一致」エラー

系をテンソル積にする場合、順序を一貫させる必要があります。

# ❌ Problem: H is (2x10) but Op is (10x2)
# ✅ Solution: Check dims and ensure consistent order in tensor()
print(H.dims) # [[2, 10], [2, 10]]

不十分なフォック空間

光子数分布が N の上限に達した場合、結果は無効です。

# ✅ Solution: Check populations
occ = qt.expect(qt.destroy(N).dag() * qt.destroy(N), final_state)
# If occ > 0.8 * N, increase N.

エルミート・ハミルトニアン

ソルバーは H がエルミートであることを前提としています。時間依存コードの結果として H が非エルミートになると、sesolve は物理的でない状態（ノルム ≠ 1）を生成します。

# ✅ Solution: Check norm
if not np.allclose(result.states[-1].norm(), 1.0):
    print("Warning: Norm not preserved!")

ベストプラクティス

1. 配列を手動で作成するのではなく、組み込みコンストラクター（basis、coherent など）を常に使用します
2. テンソル積を計算する前に、.dims で次元を確認します
3. 期待値を手動で計算するのではなく、ソルバーに演算子を渡します
4. 適切なフォック空間の切り捨て（N）を使用します。占有数が飽和していないか確認してください
5. 時間依存ハミルトニアンには、最高のパフォーマンスを得るために文字列ベースの形式を使用します
6. ソルバーの警告を監視します。多くの場合、数値的な問題を示しています
7. N > 10 のシステムでは、疎行列（デフォルト）を使用します
8. 特にカスタムソルバーでは、時間発展後に状態が正規化されているか確認します
9. パラメーターのスイープには parallel_map を使用し、複数の CPU を活用します
10. より深く理解するために、ウィグナー関数やブロッホ球を使って量子状態を可視化します

QuTiP は、直感に反する量子ダイナミクスの世界を探究するための標準的なツールです。理論方程式と数値実験の間をつなぎ、物理学者に Python スクリプト内で高性能な実験室を提供します。