LM-02 カスタムモデルハミルトニアン

カスタムモデルハミルトニアンの定義：スピン・フェルミオン・ボソン

チュートリアル LM-01 では、カスタム格子の定義方法を学びました。本チュートリアルではさらに一歩進み、局所ヒルベルト空間と相互作用項という基本的な構成要素から量子ハミルトニアンを定義するモデルライブラリXMLファイルの書き方を解説します。これにより、ALPSが標準で提供する少数のモデルだけでなく、任意の格子モデルをシミュレートできるようになります。

このチュートリアルを終えると、以下のことができるようになります：

4つのXML構成要素（SITEBASIS、QUANTUMNUMBER、OPERATOR、HAMILTONIAN）を説明する。
単一イオン異方性を持つスピン-1 Heisenberg鎖をゼロから構築し、sparsediagでスピンギャップを抽出する。
スピンレスフェルミオンとソフトコアボソンに対して同様のモデルファイルを記述する。
粒子の種類によって行列要素の公式がどのように異なるかを理解する。

ALPSモデルライブラリXML形式

モデルライブラリファイルは<?xml ...?>前文を持たない裸のXML文書で、ルート要素は<MODELS>です。完全なモデルは3つのステップで構築します：

ステップ1 — SITEBASIS：局所ヒルベルト空間

<SITEBASIS name="my basis">
  <PARAMETER name="P" default="1"/>              <!-- 可変パラメータ（例：スピンS）-->
  <QUANTUMNUMBER name="Q" min="-P" max="P"
                 type="half_integer"/>           <!-- 保存量子数                  -->
  <OPERATOR name="Raise"
            matrixelement="sqrt(P*(P+1)-Q*(Q+1))">  <!-- 非対角昇演算子 -->
    <CHANGE quantumnumber="Q" change="1"/>
  </OPERATOR>
  <OPERATOR name="Lower"
            matrixelement="sqrt(P*(P+1)-Q*(Q-1))">
    <CHANGE quantumnumber="Q" change="-1"/>
  </OPERATOR>
  <OPERATOR name="Diag" matrixelement="Q"/>      <!-- 対角演算子                  -->
</SITEBASIS>

重要なルール：

matrixelementは<CHANGE>が適用される前の量子数値で評価される代数式です（始状態での評価）。
<CHANGE quantumnumber="Q" change="±1"/>はQを指定量だけ変化させます。ALPSはmin/maxの境界を自動的に適用し、範囲外では行列要素をゼロにします。
type="half_integer"は半整数の量子数値を許可します（スピン $S=\tfrac{1}{2}$ 、 $\tfrac{3}{2}$ 、… に必要）。

ステップ2 — BASIS：量子数制約

BASISはSITEBASISとHAMILTONIANを繋ぐコネクタです。どの量子数が保存されるか、および各セクターに対応するシミュレーションパラメータを宣言します。

<BASIS name="my basis">
  <SITEBASIS ref="my basis"/>
  <CONSTRAINT quantumnumber="Q" value="Q_total"/>  <!-- 全QをパラメータQ_totalに固定 -->
</BASIS>

ステップ3 — HAMILTONIAN：サイト項とボンド項

<HAMILTONIAN name="my model">
  <PARAMETER name="J" default="1"/>
  <BASIS ref="my basis"/>             <!-- SITEBASISではなくBASISを参照 -->
  <SITETERM site="i">
    D*Diag(i)*Diag(i)                <!-- サイト式                      -->
  </SITETERM>
  <BONDTERM source="i" target="j">
    J*Diag(i)*Diag(j)                <!-- ボンド式                      -->
  </BONDTERM>
</HAMILTONIAN>

よくある間違い： <HAMILTONIAN>内に<SITEBASIS ref="..."/>と書くと、ALPSは"unexpected element"エラーで終了します。 <HAMILTONIAN>内では必ず<BASIS ref="..."/>を使用し、<SITEBASIS>をラップする独立した<BASIS>要素を定義してください。

完全なファイルの構造：

<MODELS>
  <!-- SITEBASIS要素（局所ヒルベルト空間）-->
  <!-- BASIS要素（制約宣言）             -->
  <!-- HAMILTONIAN要素（相互作用項）     -->
</MODELS>

パラメータファイルでは以下のように参照します：

MODEL_LIBRARY = "my_models.xml"
MODEL         = "my model"

例1 — スピン：単一イオン異方性を持つS=1鎖

物理的背景

ハミルトニアンは

H = J \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j + D \sum_i \bigl(S^z_i\bigr)^2

ここで $J > 0$ は反強磁性交換相互作用、 $D$ は単一イオン異方性です。 $\mathbf{S}_i \cdot \mathbf{S}_j$ を昇降演算子で書くと：

\mathbf{S}_i \cdot \mathbf{S}_j = S^z_i S^z_j + \tfrac{1}{2}\bigl(S^+_i S^-_j + S^-_i S^+_j\bigr).

$D = 0$ では、この鎖はHaldane相にあります：バルクスピンギャップ $\Delta \approx 0.41\,J$ を持つギャップ付き対称性保護トポロジカル相です（Haldane 1983；White & Huse 1993）。 $D$ を増加させると、量子臨界点 $D_c \approx 0.97\,J$ でギャップが閉じ、系は大- $D$ 相に入ります（基底状態は積状態 $|0, 0, \ldots, 0\rangle$ 、すなわちすべてのサイトが $S^z = 0$ の状態に近い）。

ハミルトニアンの骨格（開放鎖、4サイトの例）：

  D       D       D       D
  o --J-- o --J-- o --J-- o

D  = 単一イオン異方性  D*(Sz)^2      [SITETERM]
J  = Heisenberg交換               [BONDTERM]
h  = 一様磁場         -h*Sz        [SITETERM、省略可]

格子の選択

本例ではALPS組み込みの**chain lattice**（周期境界条件、 $L=8$ サイト）を使用します。周期境界条件下ではsparsediagが並進対称性も活用し、全運動量 $k$ ごとにブロック対角化することで、バルクのHaldaneギャップをより明確に観察できます。利用可能な格子形状の一覧はALPS格子ライブラリを参照してください。

手法の選択

sparsediag（疎行列Lanczos厳密対角化）は正確な固有値を与え、ヒルベルト空間次元が $\lesssim 10^6$ の系に適しています。 $L=8$ 、 $S=1$ 、 $S^z_{\text{tot}}=0$ セクターの基底数は約1.3万のみで、Lanczosは数秒で完了します。ここで厳密対角化を選ぶ理由は、 $D$ の関数として4桁精度のエネルギー差（ギャップ）を正確に測定するためです。

ステップ1：モデルライブラリファイルの作成

my_models.xmlを作成します：

<MODELS>

  <!-- ===========================================================
       サイト基底：汎用スピン-S
         * デフォルトはlocal_S=1（S=1鎖）；local_S=0.5とすればS=1/2鎖になる
         * 量子数：S（固定）とSz（-SからSまで）
         * 演算子：Splus・Sminus（昇降）、Sz（対角）
  ============================================================ -->
  <SITEBASIS name="my spin">
    <PARAMETER name="local_S" default="1"/>

    <QUANTUMNUMBER name="S"  min="local_S" max="local_S" type="half_integer"/>
    <QUANTUMNUMBER name="Sz" min="-S" max="S" type="half_integer"/>

    <!-- S^+はSzを1増加させる；行列要素はWigner-Eckartの定理による -->
    <OPERATOR name="Splus" matrixelement="sqrt(S*(S+1)-Sz*(Sz+1))">
      <CHANGE quantumnumber="Sz" change="1"/>
    </OPERATOR>

    <!-- S^-はSzを1減少させる -->
    <OPERATOR name="Sminus" matrixelement="sqrt(S*(S+1)-Sz*(Sz-1))">
      <CHANGE quantumnumber="Sz" change="-1"/>
    </OPERATOR>

    <!-- Szは対角演算子：行列要素はSzの現在値に等しい -->
    <OPERATOR name="Sz" matrixelement="Sz"/>
  </SITEBASIS>

  <!-- BASISはサイト基底を全Sz制約に結びつける -->
  <BASIS name="anisotropic spin chain">
    <SITEBASIS ref="my spin"/>
    <CONSTRAINT quantumnumber="Sz" value="Sz_total"/>
  </BASIS>

  <!-- ===========================================================
       ハミルトニアン：Heisenberg交換 + 単一イオン異方性
         H = J * sum_{<ij>} (Sz(i)*Sz(j) + 0.5*(S+(i)*S-(j) + S-(i)*S+(j)))
           + D * sum_i  Sz(i)^2
           - h * sum_i  Sz(i)      [任意：一様磁場]
  ============================================================ -->
  <HAMILTONIAN name="anisotropic spin chain">
    <PARAMETER name="J" default="1"/>
    <PARAMETER name="D" default="0"/>
    <PARAMETER name="h" default="0"/>

    <BASIS ref="anisotropic spin chain"/>

    <!-- サイト項：単一イオン異方性 + 磁場 -->
    <SITETERM site="i">
      D*Sz(i)*Sz(i) - h*Sz(i)
    </SITETERM>

    <!-- ボンド項：S+/S-/Szで書いたHeisenberg交換 -->
    <BONDTERM source="i" target="j">
      J*(Sz(i)*Sz(j) + 0.5*(Splus(i)*Sminus(j) + Sminus(i)*Splus(j)))
    </BONDTERM>
  </HAMILTONIAN>

</MODELS>

各要素の役割

<SITEBASIS name="my spin"> は局所ヒルベルト空間を定義します。スピン- $S$ サイトは $S^z \in \{-S,\, -S+1,\, \ldots,\, S\}$ で標識される $2S+1$ 個の状態を持ちます。デフォルトのlocal_S=1では、各サイトに3つの状態 $|-1\rangle$ 、 $|0\rangle$ 、 $|{+}1\rangle$ が与えられます。

<PARAMETER name="local_S"> はスピンの大きさを設定します。 ALPSではスピン量子数（S）と可変パラメータの名前を別にする必要があります（実行時の名前の衝突を避けるため）；ALPSの慣例に従い、可変パラメータはlocal_S、量子数はSという名前にします。

<QUANTUMNUMBER name="S" min="local_S" max="local_S" ...> は $S$ をlocal_Sパラメータの値に固定し、各サイトが同じ全スピンを持つようにします。 <QUANTUMNUMBER name="Sz" ...> は $S^z$ を流動量子数として登録し、ALPSがハミルトニアンをブロック対角化できるようにします。

<BASIS name="..."> はサイト基底をラップし、全 $S^z$ をシミュレーションパラメータSz_totalに固定する<CONSTRAINT>を追加します。 <HAMILTONIAN>内には<BASIS ref="..."/>と書きます（<SITEBASIS ref="..."/>ではありません）；誤ったタグを使うとALPSはエラーを報告します。

<OPERATOR name="Splus" matrixelement="..."> は昇演算子 $S^+$ を定義します。 matrixelementは $\langle S,\, S^z{+}1 \mid S^+ \mid S,\, S^z \rangle = \sqrt{S(S+1)-S^z(S^z+1)}$ を与えます。これは $S^z$ の始状態値で評価されます。子要素 <CHANGE quantumnumber="Sz" change="1"/> はALPSに $S^+$ が $S^z$ を1増加させることを伝え、 $S^z = S$ のときALPSは自動的に行列要素をゼロにします。

<SITETERM site="i"> はサイトハミルトニアンを指定します。式D*Sz(i)*Sz(i)はサイト $i$ で対角演算子 $S^z_i$ を2回評価し、単一イオン異方性 $D(S^z_i)^2$ を与えます。

<BONDTERM source="i" target="j"> は2サイトボンド相互作用を指定します。積Splus(i)*Sminus(j)は $S^+_i \otimes S^-_j$ を表し、ALPSが自動的にテンソル積行列を構築します。係数 $\tfrac{1}{2}$ は $S^+S^- = 2(S^xS^x + S^yS^y)$ から来ます。

ステップ2：パラメータファイルの作成

スピンギャップ $\Delta = E_1(S^z_{\text{tot}}=1) - E_0(S^z_{\text{tot}}=0)$ を求めるには、各 $S^z_{\text{tot}}$ セクターで別々にsparsediagを実行します。 parm_spin1を作成します：

MODEL_LIBRARY            = "my_models.xml"
LATTICE                  = "chain lattice"
MODEL                    = "anisotropic spin chain"
L                        = 8
J                        = 1
NUMBER_EIGENVALUES       = 1
CONSERVED_QUANTUMNUMBERS = "Sz"

{D=0.0; Sz_total=0;}
{D=0.0; Sz_total=1;}
{D=1.0; Sz_total=0;}
{D=1.0; Sz_total=1;}
{D=2.0; Sz_total=0;}
{D=2.0; Sz_total=1;}

重要なポイント：

LATTICE="chain lattice" — ALPSの標準周期鎖；LATTICE_LIBRARYは不要。周期境界条件ではsparsediagが並進対称性も利用します（各全運動量セクターを個別に対角化）。これによりバルクのHaldaneギャップがより明確に現れます。
CONSERVED_QUANTUMNUMBERS="Sz"でブロック対角化を有効にします；ALPSは要求された値と等しいSz_totalのブロックのみ対角化します。
6つのタスクが作成されます：3つの $D$ 値のそれぞれで2つのセクター（ $S^z_{\text{tot}} = 0,\, 1$ ）。
NUMBER_EIGENVALUES=1は各運動量ブロックの最低固有値のみを要求します。

ステップ3：シミュレーションの実行

parameter2xml parm_spin1
sparsediag --write-xml parm_spin1.in.xml

sparsediagは6つのタスクを順に実行し、結果を parm_spin1.task1.out.xmlからparm_spin1.task6.out.xmlに書き出します。

ステップ4：スピンギャップの抽出

import xml.etree.ElementTree as ET, glob

# 周期境界条件下では、sparsediagは各全運動量セクターを個別に対角化し、
# セクターごとに1つの<EIGENVALUES>ブロックを書き出します。
# 各セクターの最低固有値の最小値が各Szセクターの基底状態エネルギーです。
results = {}
for fname in sorted(glob.glob('parm_spin1.task*.out.xml')):
    tree = ET.parse(fname)
    root = tree.getroot()
    params = {p.get('name'): p.text.strip() for p in root.findall('.//PARAMETER')}
    D  = float(params['D'])
    Sz = int(params['Sz_total'])
    means = [float(m.text)
             for m in root.findall('.//SCALAR_AVERAGE[@name="Energy"]/MEAN')]
    if means:
        results.setdefault(D, {})[Sz] = min(means)

print(f"{'D/J':>5}  {'E(Sz=0)':>12}  {'E(Sz=1)':>12}  {'ギャップΔ/J':>12}  相")
print('─' * 65)
for D in sorted(results):
    E0 = results[D][0]
    E1 = results[D][1]
    gap = E1 - E0
    phase = 'Haldane相' if D < 0.97 else '大-D相'
    print(f"{D:5.1f}  {E0:12.5f}  {E1:12.5f}  {gap:12.5f}  {phase}")

期待される出力（L=8、周期鎖、 $J=1$ ）：

  D/J      E(Sz=0)      E(Sz=1)    ギャップΔ/J  相
─────────────────────────────────────────────────────────────────
  0.0   -11.33696   -10.74340       0.59356  Haldane相
  1.0    -7.00947    -6.62505       0.38442  大-D相
  2.0    -4.33528    -3.66645       0.66883  大-D相

$D = 0$ では $L=8$ でのギャップ $\Delta \approx 0.59\,J$ ；バルク値は $0.41\,J$ で、有限サイズ補正により上側から収束します（発展問題1参照）。 $D \approx J$ 付近でギャップは $0.38\,J$ に縮小します——量子相転移（ $D_c \approx 0.97\,J$ ）の有限サイズ前兆です。 $D = 2J$ では系は大- $D$ 相にあり、ギャップは $0.67\,J$ に拡大しています。

例2 — フェルミオン：最近接反発を持つスピンレスフェルミオン鎖

スピンレスフェルミオンは最もシンプルなフェルミオンモデルです：各サイトは空（ $|0\rangle$ ）または占有（ $|1\rangle$ ）のどちらかで、スピン自由度はありません。ハミルトニアンは

H = -t \sum_{\langle i,j \rangle} \bigl(c^\dagger_i c_j + c^\dagger_j c_i\bigr) + V \sum_{\langle i,j \rangle} n_i n_j - \mu \sum_i n_i

で与えられ、自由ホッピングバンド（ $t$ で制御）、最近接密度-密度反発 $V$ 、化学ポテンシャル $\mu$ を記述します。 $V > 2t$ のとき基底状態は電荷密度波、 $V < 2t$ のときはLuttinger液体（金属状態）です。

ハミルトニアンの骨格（開放鎖、4サイトの例）：

  -μ        -μ        -μ        -μ
   o ------- o ------- o ------- o
      -t,V      -t,V      -t,V

-μ = 化学ポテンシャル  -μ*n           [SITETERM]
-t = ホッピング振幅    -t*(c†c+h.c.)  [BONDTERM]
 V = 密度-密度反発     V*n*n          [BONDTERM]

my_models.xmlに以下を追加します（<MODELS>要素内に）：

  <!-- ===========================================================
       サイト基底：スピンレスフェルミオン
         * 占有数 N ∈ {0, 1}（max="1"によるパウリ排他原理）
         * type="fermion"でALPSに1次元系のJordan-Wigner変換を自動適用させる
  ============================================================ -->
  <SITEBASIS name="my spinless fermion" type="fermion">
    <QUANTUMNUMBER name="N" min="0" max="1"/>

    <!-- c†はフェルミオンを生成：⟨1|c†|0⟩ = 1 -->
    <OPERATOR name="cdag" matrixelement="1">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>

    <!-- cはフェルミオンを消滅：⟨0|c|1⟩ = 1 -->
    <OPERATOR name="c" matrixelement="1">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>

    <!-- 数演算子：対角、行列要素は現在のN -->
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="spinless fermion chain">
    <SITEBASIS ref="my spinless fermion"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <!-- ===========================================================
       ハミルトニアン：スピンレスフェルミオン鎖
         H = -t*(cdag(i)*c(j)+h.c.) + V*n(i)*n(j) - mu*n(i)
  ============================================================ -->
  <HAMILTONIAN name="spinless fermion chain">
    <PARAMETER name="t"  default="1"/>
    <PARAMETER name="V"  default="0"/>
    <PARAMETER name="mu" default="0"/>

    <BASIS ref="spinless fermion chain"/>

    <SITETERM site="i">
      -mu*n(i)
    </SITETERM>

    <BONDTERM source="i" target="j">
      -t*(cdag(i)*c(j) + cdag(j)*c(i)) + V*n(i)*n(j)
    </BONDTERM>
  </HAMILTONIAN>

`type="fermion"`が重要な理由

フェルミオン演算子は反交換関係 $\{c_i, c^\dagger_j\} = \delta_{ij}$ を満たします。 1次元では、これはJordan-Wigner変換によって厳密に処理されます。この変換はフェルミオンを非局所的な符号列を持つハードコアボソンに写像します。 SITEBASISにtype="fermion"を設定することで、ALPSはフェルミオン演算子が中間サイトをホップする際にこの符号列を自動的に挿入します。設定しないと符号が誤りとなり、スペクトルが不正確になります。

格子の選択

open chain lattice（開放境界条件、 $L=10$ サイト）を使用します。開放境界条件は相互作用フェルミオンの運動量セクター対角化で生じる符号の複雑さを回避し、1次元鎖EDの標準的な選択です。他の鎖変形についてはALPS格子ライブラリを参照してください。

手法の選択

**sparsediag**はここで正確かつ効率的です： $L=10$ 、 $N=5$ のヒルベルト空間次元は $\binom{10}{5}=252$ のみで、ほぼ瞬時に計算できます。 NUMBER_EIGENVALUES=4で基底状態と最低中性励起状態を同時に取得します。

パラメータファイル

MODEL_LIBRARY            = "my_models.xml"
LATTICE                  = "open chain lattice"
MODEL                    = "spinless fermion chain"
L                        = 10
t                        = 1
mu                       = 0
NUMBER_EIGENVALUES       = 4
CONSERVED_QUANTUMNUMBERS = "N"

{V=0.0; N_total=5;}
{V=1.0; N_total=5;}
{V=2.0; N_total=5;}
{V=3.0; N_total=5;}

パラメータ	意味	使用値
`L`	鎖の長さ（サイト数）	10
`t`	ホッピング振幅	1
`mu`	化学ポテンシャル	0（ハーフフィリング）
`V`	最近接反発	0–3（スキャン）
`N_total`	総フェルミオン数	5（ハーフフィリング）

N_total=5はハーフフィリング（サイト2つにフェルミオン1つ）で、電荷密度波不安定性が最大になります。

シミュレーションの実行

parameter2xml parm_fermion
sparsediag --write-xml parm_fermion.in.xml

結果の抽出と表示

import xml.etree.ElementTree as ET, glob

results = {}
for fname in sorted(glob.glob('parm_fermion.task*.out.xml')):
    tree = ET.parse(fname)
    root = tree.getroot()
    params = {p.get('name'): p.text.strip() for p in root.findall('.//PARAMETER')}
    V  = float(params['V'])
    means = [float(m.text)
             for m in root.findall('.//SCALAR_AVERAGE[@name="Energy"]/MEAN')]
    if means:
        results[V] = sorted(means)

print(f"{'V/t':>5}  {'E₀':>10}  {'E₁':>10}  {'E₁−E₀':>8}")
print('─' * 42)
for V in sorted(results):
    Es = results[V]
    print(f"{V:5.1f}  {Es[0]:10.5f}  {Es[1]:10.5f}  {Es[1]-Es[0]:8.5f}")

期待される出力（L=10、開放鎖、 $t=1$ 、 $N=5$ ）：

  V/t        E₀          E₁       E₁−E₀
──────────────────────────────────────────
  0.0   -6.02667    -5.45742    0.56926
  1.0   -5.00000    -4.30391    0.69609
  2.0   -4.26419    -3.50162    0.76257
  3.0   -3.74967    -2.97936    0.77031

まとめ

ハーフフィリングセクター内の励起ギャップは $V$ とともに単調増加します：斥力が電荷セクターを硬化させ、中性（粒子-正孔）励起のエネルギーコストを増大させます。真の電荷ギャップ $\Delta_c = E_0(N+1) + E_0(N-1) - 2E_0(N)$ （ $V < 2t$ でゼロ、 $V > 2t$ で開く）を測定するには、パラメータファイルにN_total=4とN_total=6のタスクを追加します（発展問題3参照）。

例3 — ボソン：ソフトコアBose-Hubbard鎖

Bose-Hubbardモデルでは各サイトに複数のボソン（最大Nmax個）が存在でき、光学格子実験における超流体–Mott絶縁体転移を記述します：

H = -t \sum_{\langle i,j \rangle} \bigl(b^\dagger_i b_j + b^\dagger_j b_i\bigr) + \frac{U}{2} \sum_i n_i(n_i - 1) - \mu \sum_i n_i.

オンサイト斥力 $U$ は二重占有（およびそれ以上）にペナルティを与え、 $U \gg t$ かつ公約可能なフィリングのとき、基底状態は粒子-ホールギャップを持つMott絶縁体になります。

ハミルトニアンの骨格（開放鎖、4サイトの例）：

   U,-μ       U,-μ       U,-μ       U,-μ
    o -------- o -------- o -------- o
         -t         -t         -t

 U  = オンサイト斥力   (U/2)*n*(n-1)  [SITETERM]
-μ  = 化学ポテンシャル  -μ*n           [SITETERM]
-t  = ホッピング振幅   -t*(b†b+h.c.)   [BONDTERM]

my_models.xmlに以下を追加します：

  <!-- ===========================================================
       サイト基底：ソフトコアボソン
         * 占有数 N ∈ {0, 1, …, Nmax}（デフォルトNmax=3）
         * b†とbはボソン交換関係 [b, b†] = 1 からのsqrt因子を持つ
  ============================================================ -->
  <SITEBASIS name="my boson">
    <PARAMETER name="Nmax" default="3"/>

    <QUANTUMNUMBER name="N" min="0" max="Nmax"/>

    <!-- b†はNを1増加させる；行列要素 ⟨N+1|b†|N⟩ = sqrt(N+1) -->
    <OPERATOR name="bdag" matrixelement="sqrt(N+1)">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>

    <!-- bはNを1減少させる；行列要素 ⟨N-1|b|N⟩ = sqrt(N) -->
    <OPERATOR name="b" matrixelement="sqrt(N)">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>

    <!-- 数演算子 -->
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="Bose-Hubbard chain">
    <SITEBASIS ref="my boson"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <!-- ===========================================================
       ハミルトニアン：ソフトコアBose-Hubbard鎖
         H = -t*(bdag(i)*b(j)+h.c.) + U/2*n(i)*(n(i)-1) - mu*n(i)
  ============================================================ -->
  <HAMILTONIAN name="Bose-Hubbard chain">
    <PARAMETER name="t"   default="1"/>
    <PARAMETER name="U"   default="0"/>
    <PARAMETER name="mu"  default="0"/>

    <BASIS ref="Bose-Hubbard chain"/>

    <SITETERM site="i">
      U/2*n(i)*(n(i)-1) - mu*n(i)
    </SITETERM>

    <BONDTERM source="i" target="j">
      -t*(bdag(i)*b(j) + bdag(j)*b(i))
    </BONDTERM>
  </HAMILTONIAN>

行列要素の比較：ボソン・フェルミオン・スピン

粒子の種類	生成演算子	消滅演算子
スピン（ $S^+$ ）	$\sqrt{S(S+1)-S^z(S^z+1)}$	$\sqrt{S(S+1)-S^z(S^z-1)}$
ソフトコアボソン（ $b^\dagger$ ）	$\sqrt{N+1}$	$\sqrt{N}$
ハードコアボソン / スピンレスフェルミオン（ $c^\dagger$ ）	$1$	$1$

ソフトコアボソンの $\sqrt{N+1}$ 因子は、すでに占有されているモードへの誘導放出の増強を反映しています——ボース-アインシュタイン統計の本質です。ハードコアボソン（占有数上限1）とスピンレスフェルミオンは行列要素 $= 1$ を共有しますが、フェルミオンの反交換関係という違いがあり、これはSITEBASISのtype="fermion"でエンコードされます。

格子の選択

open chain lattice（ $L=6$ サイト）を使用します。 $N_{\max}=3$ 、 $L=6$ では $N=6$ セクターのヒルベルト空間次元が管理可能でLanczos EDに適しています。 2次元の超流体-Mott相図を研究する正方格子など他の形状についてはALPS格子ライブラリを参照してください。

手法の選択

**sparsediag**は任意の $U/t$ で正確な結果を与え、固定 $N$ セクターの最低固有値を直接取得します。大規模系には量子モンテカルロ（worm）が推奨されますが、EDはベンチマーク値の確立と小 $L$ での定量的比較に最適です。

パラメータファイル

MODEL_LIBRARY            = "my_models.xml"
LATTICE                  = "open chain lattice"
MODEL                    = "Bose-Hubbard chain"
Nmax                     = 3
L                        = 6
t                        = 1
NUMBER_EIGENVALUES       = 2
CONSERVED_QUANTUMNUMBERS = "N"

{U=0.0;  mu=0; N_total=6;}
{U=2.0;  mu=0; N_total=6;}
{U=5.0;  mu=0; N_total=6;}
{U=10.0; mu=0; N_total=6;}

パラメータ	意味	使用値
`L`	鎖の長さ（サイト数）	6
`Nmax`	1サイト当たりの最大ボソン数	3
`t`	ホッピング振幅	1
`U`	オンサイト斥力	0–10（スキャン）
`mu`	化学ポテンシャル	0
`N_total`	総ボソン数	6（ユニットフィリング）

N_total=6はユニットフィリング（サイトごとにボソン1個）を与え、Mottギャップが最大になります。

シミュレーションの実行

parameter2xml parm_boson
sparsediag --write-xml parm_boson.in.xml

結果の抽出と表示

import xml.etree.ElementTree as ET, glob

results = {}
for fname in sorted(glob.glob('parm_boson.task*.out.xml')):
    tree = ET.parse(fname)
    root = tree.getroot()
    params = {p.get('name'): p.text.strip() for p in root.findall('.//PARAMETER')}
    U  = float(params['U'])
    means = [float(m.text)
             for m in root.findall('.//SCALAR_AVERAGE[@name="Energy"]/MEAN')]
    if means:
        results[U] = sorted(means)

print(f"{'U/t':>5}  {'E₀(N=6)':>12}  {'E₁(N=6)':>12}  {'E₁−E₀':>8}")
print('─' * 48)
for U in sorted(results):
    Es = results[U]
    print(f"{U:5.1f}  {Es[0]:12.5f}  {Es[1]:12.5f}  {Es[1]-Es[0]:8.5f}")

期待される出力（L=6、開放鎖、 $t=1$ 、 $N_{\max}=3$ 、 $N=6$ ）：

  U/t      E₀(N=6)      E₁(N=6)     E₁−E₀
────────────────────────────────────────────────
  0.0   -10.26758     -9.55861     0.70897
  2.0    -6.64995     -5.60418     1.04577
  5.0    -3.71673     -1.39313     2.32360
 10.0    -1.96696      4.35707     6.32403

まとめ

ユニットフィリングセクター内の励起ギャップは $U$ とともに急激に増大します： $U=0$ （超流体）ではギャップは運動エネルギーのみで決まり（ $\approx 0.71\,t$ ）、 $U=10\,t$ ではオンサイト斥力があらゆる粒子-正孔励起を非常にコスト高にし（ $\approx 6.3\,t$ ）、深いMott絶縁体的特徴を示します。真の一粒子ギャップ $\Delta = E_0(N+1) + E_0(N-1) - 2E_0(N)$ （超流体相でゼロ、Mott相で開く）を測定するには、N_total=5とN_total=7のタスクを追加します（発展問題4参照）。

完全な`my_models.xml`

3つのサイト基底とハミルトニアンを1つのファイルにまとめます：

<MODELS>

  <!-- ===== スピン ===== -->
  <SITEBASIS name="my spin">
    <PARAMETER name="local_S" default="1"/>
    <QUANTUMNUMBER name="S"  min="local_S" max="local_S" type="half_integer"/>
    <QUANTUMNUMBER name="Sz" min="-S" max="S" type="half_integer"/>
    <OPERATOR name="Splus"  matrixelement="sqrt(S*(S+1)-Sz*(Sz+1))">
      <CHANGE quantumnumber="Sz" change="1"/>
    </OPERATOR>
    <OPERATOR name="Sminus" matrixelement="sqrt(S*(S+1)-Sz*(Sz-1))">
      <CHANGE quantumnumber="Sz" change="-1"/>
    </OPERATOR>
    <OPERATOR name="Sz" matrixelement="Sz"/>
  </SITEBASIS>

  <BASIS name="anisotropic spin chain">
    <SITEBASIS ref="my spin"/>
    <CONSTRAINT quantumnumber="Sz" value="Sz_total"/>
  </BASIS>

  <HAMILTONIAN name="anisotropic spin chain">
    <PARAMETER name="J" default="1"/>
    <PARAMETER name="D" default="0"/>
    <PARAMETER name="h" default="0"/>
    <BASIS ref="anisotropic spin chain"/>
    <SITETERM site="i">D*Sz(i)*Sz(i) - h*Sz(i)</SITETERM>
    <BONDTERM source="i" target="j">
      J*(Sz(i)*Sz(j) + 0.5*(Splus(i)*Sminus(j) + Sminus(i)*Splus(j)))
    </BONDTERM>
  </HAMILTONIAN>


  <!-- ===== スピンレスフェルミオン ===== -->
  <SITEBASIS name="my spinless fermion" type="fermion">
    <QUANTUMNUMBER name="N" min="0" max="1"/>
    <OPERATOR name="cdag" matrixelement="1">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>
    <OPERATOR name="c" matrixelement="1">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="spinless fermion chain">
    <SITEBASIS ref="my spinless fermion"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <HAMILTONIAN name="spinless fermion chain">
    <PARAMETER name="t"  default="1"/>
    <PARAMETER name="V"  default="0"/>
    <PARAMETER name="mu" default="0"/>
    <BASIS ref="spinless fermion chain"/>
    <SITETERM site="i">-mu*n(i)</SITETERM>
    <BONDTERM source="i" target="j">
      -t*(cdag(i)*c(j) + cdag(j)*c(i)) + V*n(i)*n(j)
    </BONDTERM>
  </HAMILTONIAN>


  <!-- ===== ソフトコアボソン ===== -->
  <SITEBASIS name="my boson">
    <PARAMETER name="Nmax" default="3"/>
    <QUANTUMNUMBER name="N" min="0" max="Nmax"/>
    <OPERATOR name="bdag" matrixelement="sqrt(N+1)">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>
    <OPERATOR name="b" matrixelement="sqrt(N)">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="Bose-Hubbard chain">
    <SITEBASIS ref="my boson"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <HAMILTONIAN name="Bose-Hubbard chain">
    <PARAMETER name="t"  default="1"/>
    <PARAMETER name="U"  default="0"/>
    <PARAMETER name="mu" default="0"/>
    <BASIS ref="Bose-Hubbard chain"/>
    <SITETERM site="i">U/2*n(i)*(n(i)-1) - mu*n(i)</SITETERM>
    <BONDTERM source="i" target="j">
      -t*(bdag(i)*b(j) + bdag(j)*b(i))
    </BONDTERM>
  </HAMILTONIAN>

</MODELS>

発展問題

ギャップとシステムサイズの関係： $D=0$ で $L = 6, 8, 10, 12$ のスピン-1例を実行します。 $\Delta(L)$ を $1/L$ に対してプロットします。ギャップはバルクのHaldane値 $\Delta_\infty \approx 0.41\,J$ に外挿されますか？開放境界条件と周期境界条件の有限サイズ補正を比較してください。
スピン- $\tfrac{1}{2}$ とスピン-1の比較：parm_spin1にlocal_S=0.5を追加します（SITEBASISはlocal_Sを可変パラメータとして使用します）。 $S = \tfrac{1}{2}$ のHeisenbergモデルはギャップレスです（Luttinger液体）。 $L \to \infty$ でギャップがゼロに外挿されることを確認し、 $S = 1$ の場合と対比させてください。
スピンレスフェルミオンの電荷密度波ギャップ： $V = 0, 1, 2, 3$ でparm_fermionを実行します。 N_total=4とN_total=6のタスクを追加することで、電荷ギャップ $\Delta_c = E_0(N_{\text{tot}}+1) + E_0(N_{\text{tot}}-1) - 2E_0(N_{\text{tot}})$ を抽出します。 $V/t$ のどの値で電荷密度波ギャップが開きますか？
ボソンのMottギャップ：ユニットフィリングにおけるBose-Hubbard鎖の一粒子ギャップを $U/t$ の関数として計算します。スピンレスフェルミオンの電荷密度波ギャップと比較してください——どちらも電荷ギャップですが、物理的起源（斥力駆動対運動エネルギー駆動）が異なります。
カスタム格子 + カスタムモデル：my_models.xmlとLM-01のハニカム格子を組み合わせます。小さなクラスターでsparsediagを使ってハニカム格子上でスピン-1 Heisenbergモデルを実行します。単一イオン異方性 $D$ は、鎖の場合とは異なる相をフラストレートハニカム格子上で安定化させますか？
次近接ボンドの追加：anisotropic spin chainハミルトニアンにtype="1"の2番目の<BONDTERM>を追加します：
```
<BONDTERM type="1" source="i" target="j">
  J2*(Sz(i)*Sz(j) + 0.5*(Splus(i)*Sminus(j)+Sminus(i)*Splus(j)))
</BONDTERM>
```
nnn chain latticeを使用し、次近接結合 $J_2$ がHaldaneギャップに与える影響を調べます。その格子で辺タイプ1がどのように定義されているかは格子ハウツーを参照してください。

LM-01 ハニカム反強磁性体

LM-02 カスタムモデルハミルトニアン

カスタムモデルハミルトニアンの定義：スピン・フェルミオン・ボソン

ALPSモデルライブラリXML形式

ステップ1 — SITEBASIS：局所ヒルベルト空間

ステップ2 — BASIS：量子数制約

ステップ3 — HAMILTONIAN：サイト項とボンド項

例1 — スピン：単一イオン異方性を持つS=1鎖

物理的背景

格子の選択

手法の選択

ステップ1：モデルライブラリファイルの作成

各要素の役割

ステップ2：パラメータファイルの作成

ステップ3：シミュレーションの実行

ステップ4：スピンギャップの抽出

例2 — フェルミオン：最近接反発を持つスピンレスフェルミオン鎖

type="fermion"が重要な理由

格子の選択

手法の選択

パラメータファイル

シミュレーションの実行

結果の抽出と表示

まとめ

例3 — ボソン：ソフトコアBose-Hubbard鎖

行列要素の比較：ボソン・フェルミオン・スピン

格子の選択

手法の選択

パラメータファイル

シミュレーションの実行

結果の抽出と表示

まとめ

完全なmy_models.xml

発展問題

`type="fermion"`が重要な理由

完全な`my_models.xml`