メインコンテンツへスキップ
このプロジェクトはベータ版です。バグを見つけましたか?GitHub Issueを開く。質問やアイデアがありますか?Discussionsで質問する — フィードバックをお待ちしています!
このページの翻訳品質はいかがですか?

フェルミ・ハバード模型におけるスピンと電荷のダイナミクス

使用量の目安: Heron r3 プロセッサーで5分。(注: これはあくまで目安です。実際の実行時間は異なる場合があります。)

背景​

フェルミ・ハバード模型は、二重占有にエネルギーペナルティがある格子上を電子が跳び移る様子を記述します。開放境界を持つ LL サイトの1次元鎖では、ハミルトニアンは次のように定義されます。

H=−∑i=0L−2∑σ∈{↑,↓}th,i(ciσ†ci+1,σ+h.c.)+∑i=0L−1Ui ni↑ni↓−∑i=0L−1∑σμi niσ,H = -\sum_{i=0}^{L-2} \sum_{\sigma \in \{\uparrow,\downarrow\}} t_{h,i} \left( c^\dagger_{i\sigma} c_{i+1,\sigma} + \mathrm{h.c.} \right) + \sum_{i=0}^{L-1} U_i\, n_{i\uparrow} n_{i\downarrow} - \sum_{i=0}^{L-1} \sum_\sigma \mu_i\, n_{i\sigma},

ここで ciσ†c^\dagger_{i\sigma} はサイト ii にスピン σ\sigma のフェルミオンを生成し、niσ=ciσ†ciσn_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma} はその数を数えます。3つの項は、サイト ii と i+1i+1 の間の結合上のホッピング振幅 th,it_{h,i}、オンサイト相互作用 UiU_i、化学ポテンシャル μi\mu_i です。このチュートリアルでは、エネルギーを tht_h の単位で、時間を 1/th1/t_h の単位で測定します。

この模型のリアルタイムダイナミクスは、古典的に計算することが困難です。厳密対角化は L≲14L \lesssim 14 サイト付近が限界で、クエンチ後にはテンソルネットワークソルバーのエンタングルメントエントロピーが時間に対して線形に増加するため、一定の精度に必要な結合次元は指数関数的に増加します。

1次元では、基本的な励起はスピンを持つが電荷を持たないスピノンと、電荷を持つがスピンを持たないホロンです。両者は異なる速度で伝搬するため、局所的な擾乱は速い波面と遅い波面に分かれます。この効果はスピン・電荷分離と呼ばれ、このチュートリアルではこれを測定します [1]。

このチュートリアルでは Dynamics Simulator Qiskit Function を使用します。この関数は、ハミルトニアン、初期状態、Trotter スケジュール、および観測量のリストを受け取り、期待値を返します。Jordan-Wigner マッピング、Trotter 分解、量子ビットレイアウト、エラー抑制は内部で処理されるため、このチュートリアルで構築またはトランスパイルする回路はありません。

要件​

このチュートリアルを始める前に、次のパッケージがインストールされていることを確認してください。

  • Qiskit SDK v2.3.0 以降(pip install qiskit)

  • Qiskit Runtime v0.50.0 以降(pip install qiskit-ibm-runtime)

  • Qiskit Functions Catalog v0.21.0 以降(pip install qiskit-ibm-catalog)

  • NumPy v2.5.1 以降(pip install numpy)

  • Matplotlib v3.10.3 以降(pip install matplotlib)

  • Q-CTRL Visualizer v11.0.0 以降(pip install qctrl-visualizer)

セットアップ​

まず、IBM Quantum API キーで認証します。次に、次のようにして Qiskit Function を選択します。(このコードは、すでにアカウントを保存してローカル環境に用意していることを前提としています。)

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qctrlvisualizer": "qctrlvisualizer", "qiskit_ibm_catalog": "qiskit-ibm-catalog", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
%pip install -q {_pip}
if not _missing:
print("\u2713 All required packages are installed")
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from qiskit_ibm_catalog import QiskitFunctionsCatalog
from qiskit_ibm_runtime import QiskitRuntimeService
import qctrlvisualizer as qv

plt.style.use(qv.get_qctrl_style())
catalog = QiskitFunctionsCatalog(channel="ibm_quantum_platform")

# Access Function
sim_qf = catalog.load("q-ctrl/dynamics-simulator")

ステップ 1: 古典的な入力を量子問題にマッピングする​

ハミルトニアンを定義する​

この関数は、3つのキーを持つ辞書としてモデルを受け取ります。各キーは、鎖全体に一様に適用される単一の浮動小数点数、またはサイトごとに値を指定するリストのいずれかを受け付けます。リストの要素数は、hopping_amplitude では L−1L-1 個(開放鎖の結合はサイトより1つ少ないため)、on_site_interaction と chemical_potential では LL 個です。

このチュートリアルでは、斥力相互作用 U/th=8U/t_h = 8 を持つ31サイトの一様な鎖を使用します。一様な化学ポテンシャルは、粒子数が固定されている場合は大域位相にしか寄与しないため、この例ではゼロに設定します。化学ポテンシャルが重要になるのは、サイトごとに変化する場合のみで、たとえば無秩序や単一の不純物を加える場合です。

site_count = 31
hopping_amplitude = 1.0
on_site_interaction = 8.0
chemical_potential = 0.0

model = {
"hopping_amplitude": [hopping_amplitude] * (site_count - 1),
"on_site_interaction": [on_site_interaction] * site_count,
"chemical_potential": [chemical_potential] * site_count,
}

print(
f"{site_count} sites, {2 * site_count} qubits, U/t_h = {on_site_interaction:.0f}"
)
31 sites, 62 qubits, U/t_h = 8

Trotter スケジュールを設定する​

simulation 辞書には時間ステップ (Δt\Delta t) と Trotter ステップ数が含まれ、総発展時間はその積になります。ステップを大きくすると、回路の深さを一定に保ったままより遅い時間に到達できますが、Trotter 誤差が増えます。この誤差は、ハミルトニアンの各項が可換でなくなるため UU とともに大きくなります。公表されたハードウェア研究では U/th=8U/t_h = 8 に対して Δt=0.15\Delta t = 0.15 が使用されており [2]、このチュートリアルでも同じ値を使用します。回路の深さはステップ数に対して線形に増加し、LL にはまったく依存しません。

time_step = 0.15
step_count = 30

simulation = {"time_step": time_step, "step_count": step_count}

print(
f"Total evolution time: t = {time_step * step_count:.2f} (units of 1/t_h)"
)
Total evolution time: t = 4.50 (units of 1/t_h)

初期状態を準備する​

各サイトには2つのスピン軌道があるため、LL サイトの鎖には 2L2L 個の量子ビットを使用します。初期状態は、サイトごとに1文字を持つ Fock 占有文字列です。

empty (e)⟷(ni,↑,ni,↓)=(0,0),up (u)⟷(ni,↑,ni,↓)=(1,0),down (d)⟷(ni,↑,ni,↓)=(0,1),both (b)⟷(ni,↑,ni,↓)=(1,1).\begin{aligned} \texttt{empty (e)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (0,0), \\ \texttt{up (u)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (1,0), \\ \texttt{down (d)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (0,1), \\ \texttt{both (b)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (1,1). \end{aligned}

この文字列は最下位サイトを先頭として読まれるため、右端の文字がサイト 0 であり、左から数えて kk 番目の文字の物理的なサイトインデックスは i=L−1−ki = L - 1 - k です。したがって、文字列 udud は、サイト 0 が d、サイト 1 が u、サイト 2 が d、サイト 3 が u であることを意味します。

次に、u と d が交互に並び、中央のサイトが空孔に置き換えられた、半充填のネール状態を準備します。この空孔が、広がりを追跡する擾乱です。

vacancy_site = site_count // 2

# Build in increasing site order, then reverse: the leftmost character of the
# string is the highest site index.
site_occupations = [
"u" if site % 2 == 0 else "d" for site in range(site_count)
]
site_occupations[vacancy_site] = "e"
initial_state = "".join(reversed(site_occupations))

print(f"Vacancy at site {vacancy_site}")
print(initial_state)
Vacancy at site 15
udududududududueudududududududu

観測量を構築する​

各観測量は、占有演算子の積の重み付き和を表す辞書です。

O=∑αcα∏(σ,j)∈Tαnj,σ.O = \sum_{\alpha} c_{\alpha} \prod_{(\sigma, j) \in T_{\alpha}} n_{j,\sigma}.

operators の各エントリは1つの項 TαT_{\alpha} であり、掛け合わされる [name, site] のペアのリストとして書かれ、coefficients の対応するエントリがその重みです。2つの名前は、nj,↑n_{j,\uparrow} を表す "Nup" と、nj,↓n_{j,\downarrow} を表す "Ndn" です。単一の1点観測量 n3,↑n_{3,\uparrow} は次のように書きます。

{"operators": [[["Nup", 3]]], "coefficients": [1.0]}

次のセルでは、関数が受け付ける最小限の完全な入力セットを出力します。これは、以降のコードを読む際の形の参考として役立ちます。

minimal_inputs = {
"model": {
"hopping_amplitude": 1.0,
"on_site_interaction": [2.0, 2.5],
"chemical_potential": 0.5,
},
"simulation": {"step_count": 5, "time_step": 0.1},
"observables": [{"operators": [[["Nup", 0]]], "coefficients": [1.0]}],
"initial_state": "ud",
"shots": 4096,
}

for key, value in minimal_inputs.items():
print(f"{key}: {value}")
model: {'hopping_amplitude': 1.0, 'on_site_interaction': [2.0, 2.5], 'chemical_potential': 0.5}
simulation: {'step_count': 5, 'time_step': 0.1}
observables: [{'operators': [[['Nup', 0]]], 'coefficients': [1.0]}]
initial_state: ud
shots: 4096

3つの観測量グループが必要です。電荷密度 ni=ni↑+ni↓n_i = n_{i\uparrow} + n_{i\downarrow} と局所スピン Siz=(ni↑−ni↓)/2S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2 は、どちらも同じ演算子に対する2項の和であり、係数だけが異なります。サイト ii と空孔サイトの間のスピン相関 SizSi∗zS^z_i S^z_{i_*} は、係数が ±1/4\pm 1/4 の4つの2点項に展開されます。

占有演算子は可換なので、項内の因子の順序は無関係ですが、observables リストの順序によって、返される値の順序が決まります。コードは、リストを構築する際に各グループのインデックスを記録します。相関関数は、積が因子を繰り返すことになる空孔サイトをスキップします。

observables = []
charge_indices = []
spin_indices = []
correlator_indices = []

# Charge density: n_i = n_{i,up} + n_{i,down}
for site in range(site_count):
charge_indices.append(len(observables))
observables.append(
{
"operators": [[["Nup", site]], [["Ndn", site]]],
"coefficients": [1.0, 1.0],
}
)

# Local spin: S^z_i = (n_{i,up} - n_{i,down}) / 2
for site in range(site_count):
spin_indices.append(len(observables))
observables.append(
{
"operators": [[["Nup", site]], [["Ndn", site]]],
"coefficients": [0.5, -0.5],
}
)

# Spin correlator: S^z_i S^z_{i_*}, expanded into four two-point terms.
correlator_sites = [
site for site in range(site_count) if site != vacancy_site
]
for site in correlator_sites:
correlator_indices.append(len(observables))
observables.append(
{
"operators": [
[["Nup", site], ["Nup", vacancy_site]],
[["Nup", site], ["Ndn", vacancy_site]],
[["Ndn", site], ["Nup", vacancy_site]],
[["Ndn", site], ["Ndn", vacancy_site]],
],
"coefficients": [0.25, -0.25, -0.25, 0.25],
}
)

print(f"{len(observables)} observables")
print("Example correlator:", observables[correlator_indices[0]])
92 observables
Example correlator: {'operators': [[['Nup', 0], ['Nup', 15]], [['Nup', 0], ['Ndn', 15]], [['Ndn', 0], ['Nup', 15]], [['Ndn', 0], ['Ndn', 15]]], 'coefficients': [0.25, -0.25, -0.25, 0.25]}

ステップ 2: 量子ハードウェアでの実行に向けて問題を最適化する​

このチュートリアルにはトランスパイルのステップがありません。関数は内部で次のステップを実行し [2]、ノートブックからそれらを設定することはできません。

  • ヘビーヘックス格子上でパウリ重みを低く保つように選ばれた順序で、フェルミオンモードを量子ビットにマッピングする

  • フェルミオン SWAP ネットワークを使用して、すべてのホッピング項を2量子ビット回転にする

  • 弱い量子ビットを避ける経路に沿って鎖を埋め込む

  • コンパイル中に動的デカップリングとパウリ・ツワーリングを適用する

バックエンドとショット数を選択します。この関数はバックエンド名を文字列として受け取るため、このチュートリアルでは、鎖に十分な大きさを持つ最も空いているデバイスを選択し、その名前を渡します。

解像する必要がある最小の信号によって、ショット数が決まります。電荷トレーサーは 0.150.15 付近の振幅に達しますが、連結スピントレーサーは約 0.040.04 でピークに達し、およそ4分の1の大きさです。公表された研究では、この理由から回路ごとに20,000ショットが使用されました [2]。

service = QiskitRuntimeService()
backend = service.backend("ibm_boston")
backend_name = backend.name
shot_count = 20_000

print(f"{backend_name}, {shot_count} shots per circuit")
ibm_boston, 20000 shots per circuit

ステップ 3: Qiskit primitives を使用して実行する​

次に、シミュレーションを送信します。run はブロックせず、すぐにジョブハンドルを返します。

job = sim_qf.run(
model=model,
simulation=simulation,
observables=observables,
initial_state=initial_state,
shots=shot_count,
backend_name=backend_name,
)
print(job.job_id, job.status())
6261b073-be19-48a1-be60-2f369b25613c QUEUED

進捗を追うには job.status() を、ジョブの完了までブロックするには job.result() を使用します。結果は、observable_expectation_values エントリが、まず(渡された順序で)観測量ごと、次に時間ごとにインデックス付けされたマッピングです。t=0t=0 の値が含まれるため、時間点は step_count + 1 個あります。

result = job.result()

print("Result keys:", sorted(result.keys()))
Result keys: ['observable_expectation_values', 'observables', 'provider_job_ids', 'warnings']

ステップ 4: 結果を後処理し、望ましい古典的な形式で返す​

相関関数を再構成する​

先ほど記録したインデックスを使って、フラットな配列を3つのグループに分割し、次に2つのトレーサー相関関数を作成します。電荷トレーサーは、t=0t=0 に対するサイトごとの密度の変化です。

Cic(t)=⟨ni↑(t)+ni↓(t)⟩−⟨ni↑(0)+ni↓(0)⟩.C^c_i(t) = \langle n_{i\uparrow}(t) + n_{i\downarrow}(t) \rangle - \langle n_{i\uparrow}(0) + n_{i\downarrow}(0) \rangle.

スピントレーサーは、各サイトと空孔サイト i∗i_* の間の連結相関関数です。

Cis(t)=4(⟨SizSi∗z⟩−⟨Siz⟩⟨Si∗z⟩).C^s_i(t) = 4 \left( \langle S^z_i S^z_{i_*} \rangle - \langle S^z_i \rangle \langle S^z_{i_*} \rangle \right).

どちらも Fock 基底の演算子から構築されているため、同じ測定ビット列から得られ、2つ目のために追加のハードウェア実行は必要ありません。

観測量を構築した際に、コードは電荷 (charge_indices)、スピン (spin_indices)、スピン相関 (correlator_indices) の各演算子のリスト内での位置を記録しました。シミュレーターは期待値を同じ順序で返すため、これらの位置を使って対応する結果を取り出せます。スピン相関については、コードは correlator_sites も使って各結果を正しい格子位置に配置し、空孔は NaN のままにします。

expectation_values = np.asarray(
result["observable_expectation_values"], dtype=float
)
times = np.arange(expectation_values.shape[1]) * time_step

charge_density = expectation_values[charge_indices]
local_spin = expectation_values[spin_indices]
spin_correlator = np.full_like(charge_density, np.nan)
spin_correlator[correlator_sites] = expectation_values[correlator_indices]

# Charge tracer: change in per-site density relative to t = 0.
charge_tracer = charge_density - charge_density[:, [0]]

# Spin tracer: connected correlator with the vacancy site.
spin_tracer = 4.0 * (spin_correlator - local_spin * local_spin[vacancy_site])
spin_tracer[vacancy_site] = np.nan

print(
f"{expectation_values.shape[0]} observables, "
f"{expectation_values.shape[1]} time points up to t = {times[-1]:.2f}"
)
print(
f"Peak amplitudes: charge {np.nanmax(np.abs(charge_tracer)):.3f}, "
f"spin {np.nanmax(np.abs(spin_tracer)):.3f}"
)
92 observables, 31 time points up to t = 4.50
Peak amplitudes: charge 0.978, spin 0.168

波面を可視化する​

各相関関数を空間と時間のヒートマップとしてプロットすると、光円錐を直接見ることができます。時間は下向きに流れ、横軸は鎖に沿った位置です。

def plot_tracer(axis, tracer, title):
"""Draw one tracer correlator as a space-time heatmap, time flowing down."""
# Scale from the off-centre sites. The vacancy site carries an order-one
# change in charge density and is left to clip rather than set the scale.
off_center = np.ones(site_count, dtype=bool)
off_center[vacancy_site] = False
limit = np.nanpercentile(np.abs(tracer[off_center]), 98)

colormap = qv.QCTRL_DIVERGENT_COLORMAP.reversed().copy()
colormap.set_bad("0.82") # the spin tracer has no value at the vacancy

image = axis.imshow(
tracer.T,
aspect="auto",
origin="upper",
cmap=colormap,
vmin=-limit,
vmax=limit,
extent=(-0.5, site_count - 0.5, times[-1], times[0]),
interpolation="nearest",
)

axis.annotate(
"vacancy",
xy=(vacancy_site, times[-1]),
xytext=(vacancy_site, times[-1] + 0.55),
ha="center",
va="top",
fontsize=8,
annotation_clip=False,
arrowprops=dict(arrowstyle="-", lw=0.8, color="0.4"),
)

axis.set_title(title)
axis.set_xlabel("Site")
return image

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)

charge_image = plot_tracer(
axes[0], charge_tracer, r"Charge tracer $C^{c}_{i}(t)$"
)
spin_image = plot_tracer(axes[1], spin_tracer, r"Spin tracer $C^{s}_{i}(t)$")

axes[0].set_ylabel(r"Time $t$ $(1/t_h)$")
fig.colorbar(charge_image, ax=axes[0], fraction=0.046)
fig.colorbar(spin_image, ax=axes[1], fraction=0.046)
fig.suptitle(f"$L$ = {site_count}, $U/t_h$ = {on_site_interaction:.0f}")
fig.tight_layout()
plt.show()

Output of the previous code cell

波面の速度を抽出する​

公表された研究 [2] の検出プロトコルに従って、ヒートマップを数値に変換します。各時間スライスで、空間プロファイルを空孔について対称化し(ハードウェアノイズの反対称成分が平均化されて消えます)、中心から外れた最大振幅で正規化して、しきい値を超える最も外側のサイトを記録します。その後、検出結果を連続する時間ステップの区間に分割し、最も長い区間のみをフィットすることで、初期の過渡現象が結果に入り込まないようにします。

def front_distance(profile, center, threshold):
"""Distance from `center` to the outermost site above a relative threshold."""
values = np.nan_to_num(np.asarray(profile, dtype=float))
length = values.size

symmetrized = values.copy()
for site in range(length):
mirror = 2 * center - site
if 0 <= mirror < length:
symmetrized[site] = 0.5 * (values[site] + values[mirror])

amplitude = np.abs(symmetrized)
off_center = np.ones(length, dtype=bool)
off_center[center] = False

peak = amplitude[off_center].max()
if peak <= 0.0:
return np.nan
above = np.flatnonzero((amplitude / peak > threshold) & off_center)
return float(np.abs(above - center).max()) if above.size else np.nan

def longest_run(mask):
"""Indices of the longest run of consecutive True values."""
best, current = [], []
for index, keep in enumerate(mask):
current = current + [index] if keep else []
if len(current) > len(best):
best = current
return np.array(best, dtype=int)

def front_velocity(tracer, center, threshold=0.3, quantile=0.75, min_block=5):
"""Fit front position against time, using only the longest contiguous run."""
distances = np.array(
[
front_distance(tracer[:, k], center, threshold)
for k in range(tracer.shape[1])
]
)
edge = max(center, tracer.shape[0] - 1 - center)
usable = (
np.isfinite(distances) & (distances >= 2) & (distances <= edge - 2)
)

block = longest_run(usable)
if block.size < min_block:
return np.nan, distances, np.array([], dtype=int), np.nan

candidates = [
stats.theilslopes(distances[block[:k]], times[block[:k]]).slope
for k in range(min_block, block.size + 1)
]
cutoff = np.quantile(candidates, quantile)
keep = max(i for i, value in enumerate(candidates) if value >= cutoff)

fit_block = block[: keep + min_block]
slope, intercept = np.polyfit(times[fit_block], distances[fit_block], 1)
return float(slope), distances, fit_block, float(intercept)

フィットされた速度は、検出しきい値とトリミング分位数に依存し、その感度はショットノイズを上回ります。公表された解析と同様に、両方をスキャンし、全範囲とともに中央値を報告してください。

thresholds = [0.15, 0.225, 0.30, 0.375, 0.45]
quantiles = [0.60, 0.675, 0.75, 0.825, 0.90]

def velocity_envelope(tracer, center):
"""Median, minimum, and maximum velocity over the parameter scan."""
values = [
front_velocity(tracer, center, threshold=p, quantile=q)[0]
for p in thresholds
for q in quantiles
]
values = np.array([v for v in values if np.isfinite(v)])
if values.size == 0:
return np.nan, np.nan, np.nan
return float(np.median(values)), float(values.min()), float(values.max())

charge_velocity, charge_low, charge_high = velocity_envelope(
charge_tracer, vacancy_site
)
spin_velocity, spin_low, spin_high = velocity_envelope(
spin_tracer, vacancy_site
)

_, charge_front, charge_block, charge_intercept = front_velocity(
charge_tracer, vacancy_site
)
_, spin_front, spin_block, spin_intercept = front_velocity(
spin_tracer, vacancy_site
)

print(f"v_c = {charge_velocity:.2f} [{charge_low:.2f}, {charge_high:.2f}]")
print(f"v_s = {spin_velocity:.2f} [{spin_low:.2f}, {spin_high:.2f}]")
print(f"v_s / v_c = {spin_velocity / charge_velocity:.2f}")
v_c = 1.93 [1.81, 2.18]
v_s = 1.61 [1.23, 2.12]
v_s / v_c = 0.83
PUBLISHED = {"charge": 1.78, "spin": 1.46}

edge = max(vacancy_site, site_count - 1 - vacancy_site)
fig, axis = plt.subplots(figsize=(7, 4.6))

axis.axhspan(0, 2, color="0.92", zorder=0)
axis.axhspan(edge - 2, edge + 2, color="0.92", zorder=0)

series = [
(
"charge",
"$v_c$",
charge_front,
charge_block,
charge_velocity,
charge_intercept,
charge_low,
charge_high,
"C0",
"o",
),
(
"spin",
"$v_s$",
spin_front,
spin_block,
spin_velocity,
spin_intercept,
spin_low,
spin_high,
"C1",
"s",
),
]

for name, sym, front, block, v, b, lo, hi, color, marker in series:
excluded = np.setdiff1d(np.arange(front.size), block)
axis.plot(
times[excluded],
front[excluded],
marker,
color=color,
markersize=5,
markerfacecolor="none",
alpha=0.45,
zorder=2,
)
axis.plot(
times[block],
front[block],
marker,
color=color,
markersize=5,
label=f"{name}, {sym} = {v:.2f} [{lo:.2f}, {hi:.2f}]",
zorder=3,
)

span = times[block]
anchor = v * span[0] + b
axis.plot(
span,
v * (span - span[0]) + anchor,
"-",
color=color,
lw=1.6,
zorder=4,
)
axis.fill_between(
span,
lo * (span - span[0]) + anchor,
hi * (span - span[0]) + anchor,
color=color,
alpha=0.15,
zorder=1,
)
axis.plot(
span,
PUBLISHED[name] * (span - span[0]) + anchor,
"--",
color=color,
lw=1.2,
alpha=0.8,
zorder=4,
)

axis.plot([], [], "--", color="0.4", lw=1.2, label="published [2]")
axis.plot(
[],
[],
"o",
color="0.5",
markerfacecolor="none",
label="excluded from fit",
)

axis.set_xlabel(r"Time $t$ $(1/t_h)$")
axis.set_ylabel("Front distance from vacancy (sites)")
axis.set_xlim(0, times[-1])
axis.set_ylim(
0, min(edge + 1, np.nanmax(np.r_[charge_front, spin_front]) + 3)
)
axis.legend(loc="upper left", frameon=False, fontsize=9)
axis.set_title(
f"L = {site_count}, $U/t_h$ = {on_site_interaction:.0f}, front propagation",
fontsize=11,
)

fig.tight_layout()
plt.show()

Output of the previous code cell

検出結果をヒートマップに重ねて描くと、アルゴリズムが先端を追跡できたかどうかがわかります。塗りつぶしのマーカーは、フィットに使用された検出結果です。白抜きのマーカーは、空孔に近すぎる、鎖の端に近すぎる、または波面が広がって位置がもはや鮮明でなくなるほど時間が遅いといった理由で、マスキングルールによって除外されたものです。破線は使用した時間範囲でのフィットであり、その傾きが各パネルのタイトルに示された速度です。

空孔の列は軸の下にマークされています。電荷はその値が1程度であるためそこでクリップされ、スピンにはそこに値がないため、灰色で描かれています。

電荷の円錐はスピンの円錐よりも広く開いています。どちらの擾乱も同じ瞬間に同じ欠陥から生じ、異なる速度で外側へ伝わります。

def plot_tracer_front(axis, tracer, front, block, velocity, intercept, title):
"""Overlay the detected wavefront and the fitted cone on a tracer heatmap."""
image = plot_tracer(axis, tracer, title)

excluded = np.setdiff1d(np.arange(front.size), block)
for index, filled in ((block, True), (excluded, False)):
good = index[np.isfinite(front[index])]
if good.size == 0:
continue
for sign in (-1, 1):
axis.plot(
vacancy_site + sign * front[good],
times[good],
"o",
markersize=4,
color="k",
markerfacecolor="k" if filled else "none",
markeredgewidth=0.9,
alpha=1.0 if filled else 0.4,
)

span = times[block]
for sign in (-1, 1):
axis.plot(
vacancy_site + sign * (velocity * span + intercept),
span,
"--",
color="k",
lw=1.5,
)

axis.set_title(f"{title} $v$ = {velocity:.2f}")
return image

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)

charge_image = plot_tracer_front(
axes[0],
charge_tracer,
charge_front,
charge_block,
charge_velocity,
charge_intercept,
r"Charge tracer $C^{c}_{i}(t)$",
)
spin_image = plot_tracer_front(
axes[1],
spin_tracer,
spin_front,
spin_block,
spin_velocity,
spin_intercept,
r"Spin tracer $C^{s}_{i}(t)$",
)

axes[0].set_ylabel(r"Time $t$ $(1/t_h)$")
fig.colorbar(charge_image, ax=axes[0], fraction=0.046)
fig.colorbar(spin_image, ax=axes[1], fraction=0.046)
fig.suptitle(f"$L$ = {site_count}, $U/t_h$ = {on_site_interaction:.0f}")
fig.tight_layout()
plt.show()

Output of the previous code cell

電荷の波面はスピンの波面より速く進みます。1次元のモット背景における正孔は、電子としては伝わりません。代わりに、電荷を運ぶホロンとスピンの擾乱を運ぶスピノンに分数化し、U/th=8U/t_h = 8 ではスピノンのほうが遅くなります。

順序は相互作用に依存します。U=0U = 0 では、このモデルは自由フェルミオンで vc=2v_c = 2、vs=4v_s = 4 となり、スピンのほうが速くなります。斥力はスピン輸送を抑制する一方で電荷輸送はほとんど変えず、U/th=8U/t_h = 8 では公表されているハードウェアの値は vc=1.78v_c = 1.78、vs=1.46v_s = 1.46 です [2]。複数の UU の値で再実行し、vs/vcv_s / v_c をプロットすると、クロスオーバーを描き出せます。

数値に関する次の注意点に留意してください。しきい値に基づく検出はわずかに低めに偏るため、電荷速度が公表値より数パーセント低くなるのは想定内です。連結相関関数がノイズフロアに近いため、スピンのエンベロープは幅が広くなります。このため、単一のフィットではなくスキャンの中央値を示しています。

参考文献​

[1] Essler, F. H. L., Frahm, H., Gohmann, F., Klumper, A., & Korepin, V. E. (2005). The One-Dimensional Hubbard Model. Cambridge University Press.

[2] Hartnett, G. S., Najafi, K. S., Khindanov, A., Liao, H., Schutzman, M., Hush, M. R., Biercuk, M. J., & Baum, Y. (2026). Fast, accurate, high-resolution simulation of large-scale Fermi-Hubbard models on a digital quantum processor. arXiv:2605.04025.