Chapter 2: 分子軌道理論と電子構造

第1章では、イオン結合・共有結合・金属結合・分子間力という4つの化学結合の分類と、それぞれを記述する古典的なポテンシャルモデル（Madelung定数、Morseポテンシャル、Drude模型、Lennard-Jonesポテンシャル）を学びました。これらのモデルは結合エネルギーの大きさをよく説明しますが、電子が分子や結晶の中でどのような量子力学的な状態（軌道）を占めているかについては深く立ち入っていません。

本章では、化学結合を電子の波動関数の観点から捉え直す分子軌道理論（Molecular Orbital Theory）を学びます。出発点となるのは、原子軌道を組み合わせて分子軌道を構成するLCAO法（原子軌道の線形結合、Linear Combination of Atomic Orbitals）です。LCAO法を共役π電子系に適用したHückel理論（Hückel Theory）により、ベンゼンのような分子の電子状態を具体的に計算します。

さらに、LCAO法を無限に繰り返す周期的な原子配列（結晶）へ拡張すると、タイトバインディング近似（Tight-Binding Approximation）が得られます。これは分子軌道理論と固体のバンド理論（Band Theory）を橋渡しする重要な近似法であり、金属・半導体・絶縁体の違いを理解する基礎となります。最後に、電子状態を統計的に扱うための状態密度（Density of States, DOS）の計算方法を学び、Pythonで可視化します。

本章を通じて、孤立した2原子分子から有限の共役分子、そして無限に周期的な結晶へと系のサイズを拡大しながら、電子状態の記述がどのように一貫した枠組みの中でつながっているかを体感していただければと思います。

← 第1章 | シリーズ目次 | 第3章 →

2.1 LCAO法（原子軌道の線形結合）の基礎

分子軌道（Molecular Orbital, MO）は、分子を構成する原子の原子軌道（Atomic Orbital, AO）を重ね合わせることで構成できます。この考え方をLCAO法と呼び、分子軌道 $\psi$ を原子軌道 $\phi_i$ の線形結合として次のように表します。

変分原理を適用すると、係数 $c_i$ とエネルギー $E$ は次の永年方程式（Secular Equation）を満たします。

対角要素 $H_{ii} = \alpha_i$ はクーロン積分（Coulomb Integral）と呼ばれ、電子が孤立した原子軌道 $i$ に存在するときのエネルギーに対応します。非対角要素 $H_{ij} = \beta_{ij}$ は共鳴積分（Resonance Integral）（移動積分とも呼ばれます）で、軌道 $i$ と $j$ の間で電子が非局在化することによる安定化を表します。$\beta_{ij}$ は通常負の値を持ち、絶対値が大きいほど軌道間の相互作用が強いことを意味します。

最も単純な例として、同種の原子軌道 $\phi_1, \phi_2$ が1つの共鳴積分 $\beta$ で結合した2軌道系を考えます。

Python実装: 2軌道系のLCAO計算（結合性・反結合性軌道）

import numpy as np
import matplotlib.pyplot as plt

def lcao_two_orbital(alpha, beta):
    """
    2つの原子軌道からなる最も単純なLCAO系を対角化する

    Parameters:
    alpha: クーロン積分（原子軌道のエネルギー）
    beta: 共鳴積分（軌道間相互作用、通常は負の値）

    Returns:
    eigenvalues, eigenvectors
    """
    H = np.array([
        [alpha, beta],
        [beta, alpha]
    ])
    eigenvalues, eigenvectors = np.linalg.eigh(H)
    return eigenvalues, eigenvectors

# 規格化されたパラメータ（alpha=0を基準エネルギーとする）
alpha, beta = 0.0, -1.0
eigenvalues, eigenvectors = lcao_two_orbital(alpha, beta)

print("固有エネルギー:", eigenvalues)
print("固有ベクトル:\n", eigenvectors)

E_bonding = alpha + beta       # 結合性軌道（bonding orbital）
E_antibonding = alpha - beta   # 反結合性軌道（antibonding orbital）
print(f"\n結合性軌道エネルギー: E = alpha + beta = {E_bonding:.3f}")
print(f"反結合性軌道エネルギー: E = alpha - beta = {E_antibonding:.3f}")

# エネルギー準位図
fig, ax = plt.subplots(figsize=(6, 6))
ax.hlines(alpha, -0.5, 2.5, colors='gray', linestyles='dashed', label='原子軌道 (alpha)')
ax.hlines(E_bonding, 0.7, 1.3, colors='#f5576c', linewidth=4)
ax.hlines(E_antibonding, 0.7, 1.3, colors='#f093fb', linewidth=4)
ax.text(1.0, E_bonding - 0.15, '結合性軌道', ha='center', fontsize=11)
ax.text(1.0, E_antibonding + 0.08, '反結合性軌道', ha='center', fontsize=11)
ax.set_xticks([])
ax.set_ylabel('エネルギー (単位: |beta|)', fontsize=12)
ax.set_title('2軌道系のLCAO分裂', fontsize=14, fontweight='bold')
ax.set_xlim(-0.5, 2.5)
ax.legend(loc='upper right')
plt.tight_layout()
plt.savefig('lcao_two_orbital.png', dpi=300)
plt.show()

実行結果: alpha=0, beta=-1のとき、結合性軌道は E = alpha + beta = -1.000、反結合性軌道は E = alpha - beta = 1.000 に分裂する。beta が負であるため、結合性軌道（対称的な組み合わせ $\frac{1}{\sqrt{2}}(\phi_1+\phi_2)$）の方がエネルギーが低く安定である。

2.2 Hückel理論による π電子系の記述

LCAO法を、ベンゼンのようなπ電子系（pi-electron system）を持つ共役分子に適用したものがHückel理論です。1931年にErich Hückelによって提案されたこの理論は、いくつかの近似を導入することで、π電子系の電子状態を手計算でも扱えるほど単純化しました。

Hückel理論の主要な近似は次の3点です。

• すべての炭素原子のクーロン積分は等しい: $H_{ii} = \alpha$
• 共鳴積分は隣接する（結合している）原子間でのみ有限の値を持つ: 隣接原子なら $H_{ij} = \beta$、それ以外は0
• 原子軌道は規格化直交系とみなす: $S_{ij} = \delta_{ij}$

この近似のもとでは、ハミルトニアン行列は分子の結合構造を表す隣接行列（Adjacency Matrix） $\mathbf{A}$ を用いて $\mathbf{H} = \alpha \mathbf{I} + \beta \mathbf{A}$ と書けます。$\mathbf{A}$ の要素 $A_{ij}$ は、原子 $i$ と $j$ が結合していれば1、そうでなければ0です。

ベンゼン（$\text{C}_6\text{H}_6$）は6個の炭素原子が環状に結合しており、各炭素原子は1つのπ軌道（2p_z軌道）を提供します。隣接行列は6×6の巡回的な構造を持ち、対角化することでπ電子系の分子軌道エネルギーが得られます。

Python実装: Hückel理論によるベンゼンのπ電子状態計算

import numpy as np
import matplotlib.pyplot as plt

def huckel_matrix_ring(n_atoms, alpha, beta):
    """
    環状π電子系（例: ベンゼン）のHückelハミルトニアン行列を構築する

    Parameters:
    n_atoms: 環を構成する原子数
    alpha: クーロン積分
    beta: 共鳴積分（隣接原子間）

    Returns:
    H: Hückelハミルトニアン行列 (n_atoms x n_atoms)
    """
    H = np.full((n_atoms, n_atoms), 0.0)
    np.fill_diagonal(H, alpha)
    for i in range(n_atoms):
        j = (i + 1) % n_atoms  # 環状境界条件（最後の原子は最初の原子と結合）
        H[i, j] = beta
        H[j, i] = beta
    return H

# ベンゼン: alpha=0を基準エネルギー、beta=-1を単位として計算
alpha, beta = 0.0, -1.0
H_benzene = huckel_matrix_ring(n_atoms=6, alpha=alpha, beta=beta)

eigenvalues, eigenvectors = np.linalg.eigh(H_benzene)
print("ベンゼンのπ分子軌道エネルギー (alpha + x*beta の x):")
print(np.sort(eigenvalues))

# 6個のpi電子を、エネルギーの低い軌道から2個ずつ詰める（Aufbau原理）
n_pi_electrons = 6
n_occupied = n_pi_electrons // 2
E_pi_total = 2 * np.sum(eigenvalues[:n_occupied])
print(f"\nHOMOまで占有した全π電子エネルギー: {E_pi_total:.4f} (alpha=0, beta=-1の単位)")

# 局在化した3つの孤立二重結合（エチレン様）と比較
E_localized = n_pi_electrons * (alpha + beta)
E_delocalization = E_pi_total - E_localized
print(f"孤立3二重結合と仮定した場合のエネルギー: {E_localized:.4f}")
print(f"非局在化エネルギー（共鳴エネルギー）: {E_delocalization:.4f}")

# エネルギー準位図（縮退を考慮）
fig, ax = plt.subplots(figsize=(6, 6))
unique_levels = sorted(set(np.round(eigenvalues, 6)))
for level in unique_levels:
    degeneracy = np.sum(np.isclose(eigenvalues, level))
    offsets = np.linspace(-0.4, 0.4, degeneracy + 2)[1:-1]
    for off in offsets:
        ax.hlines(level, off - 0.15, off + 0.15, colors='#f5576c', linewidth=4)

ax.set_xticks([])
ax.set_xlim(-0.6, 0.6)
ax.set_ylabel('エネルギー (alpha + x * beta)', fontsize=12)
ax.set_title('ベンゼンのHückel π軌道エネルギー準位', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('huckel_benzene_levels.png', dpi=300)
plt.show()

実行結果: 固有値は alpha を基準として x = -2, -1, -1, +1, +1, +2（β単位）となり、$\alpha + 2\beta$（1軌道）、$\alpha + \beta$（2重縮退, Degeneracy）、$\alpha - \beta$（2重縮退）、$\alpha - 2\beta$（1軌道）の4準位に分裂する。6個のπ電子は最も安定な3軌道（$\alpha+2\beta$ と $\alpha+\beta$ の2軌道）を占有し、全π電子エネルギーは $8\beta$ となる。これを孤立した3つの二重結合（各軌道 $\alpha+\beta$ に2電子、合計 $6\beta$）と比較すると、非局在化エネルギー（共鳴エネルギー）は $2\beta$ と求まる。文献によるとπ結合の共鳴積分βはおよそ-1〜-3 eVの範囲で見積もられることが多く（例えばStreitwieser (1961) は $\beta \approx -2.5$ eV を代表値として用いている）、これを用いるとベンゼンの共鳴エネルギーは概算で約5 eV程度と見積もられる。ただし、この値は簡易Hückel理論による理想化された見積もりであり、実験的な共鳴エネルギー（水素化熱の比較から求められる値、約150 kJ/mol ≈ 1.6 eV）とは近似の粗さに起因する差がある。

2.3 タイトバインディング近似

Hückel理論は有限個の原子からなる分子を扱いましたが、原子が無限に周期的に配列した結晶（固体）に対しても同様のLCAOの考え方を適用できます。これをタイトバインディング近似と呼びます。「タイトバインディング（強結合）」という名前は、電子が各原子軌道に強く束縛されており、隣接原子への飛び移り（ホッピング）を摂動的に扱えるという物理的描像に由来します。

周期的な結晶では、格子定数 $a$ ごとに同一の原子が並んでいると仮定し、周期境界条件（Born-von Karman境界条件）を課します。このとき、結晶運動量 $k$ でラベルされる分子軌道（Bloch状態）は、原子軌道 $\phi_n$（$n$番目の格子点）の線形結合として次のように書けます。

最近接原子間のみに共鳴積分（ホッピング積分, Hopping Integral）$\beta$ を持つ1次元鎖にこの波動関数を代入すると、エネルギーの $k$ 依存性（分散関係, Dispersion Relation）が得られます。

興味深いことに、この分散関係はHückel理論と数学的に同じ起源を持ちます。ベンゼンのような $N$ 員環のHückel理論は、周期境界条件を持つ1次元タイトバインディング鎖を有限サイズ（$N=6$）に量子化したものとみなすことができます。許容される波数は $k_m = 2\pi m / (Na)$（$m = 0, 1, \dots, N-1$）に離散化され、この $k_m$ を分散関係 $E(k)$ に代入すると、ちょうどベンゼンのHückel固有値が再現されます。

Python実装: 1次元タイトバインディングモデルのバンド構造

import numpy as np
import matplotlib.pyplot as plt

def tight_binding_1d(k, alpha, beta, a=1.0):
    """1次元タイトバインディング鎖の分散関係 E(k) を計算する"""
    return alpha + 2 * beta * np.cos(k * a)

alpha, beta, a = 0.0, -1.0, 1.0

# 第一ブリルアンゾーン内の連続的なk点
k = np.linspace(-np.pi / a, np.pi / a, 400)
E_k = tight_binding_1d(k, alpha, beta, a)

print(f"バンドの最小値: {E_k.min():.4f} (alpha + 2*beta)")
print(f"バンドの最大値: {E_k.max():.4f} (alpha - 2*beta)")
print(f"バンド幅: {E_k.max() - E_k.min():.4f} (= 4|beta|)")

plt.figure(figsize=(8, 6))
plt.plot(k * a / np.pi, E_k, color='#f5576c', linewidth=2.5)
plt.axvline(x=0, color='gray', linestyle='--', alpha=0.5)
plt.xlabel('k a / π', fontsize=12)
plt.ylabel('E(k)  (alpha + x*beta)', fontsize=12)
plt.title('1次元タイトバインディング鎖のバンド構造', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('tight_binding_band_1d.png', dpi=300)
plt.show()

Python実装: Hückel理論とタイトバインディングの整合性の検証

import numpy as np

alpha, beta, a = 0.0, -1.0, 1.0

# ベンゼン（N=6環）のHückel固有値を再計算
n_atoms = 6
H = np.zeros((n_atoms, n_atoms))
for i in range(n_atoms):
    j = (i + 1) % n_atoms
    H[i, j] = beta
    H[j, i] = beta
huckel_eigenvalues = np.sort(np.linalg.eigvalsh(H))

# Born-von Karman境界条件で許容されるk点 (N=6)
N = n_atoms
m = np.arange(N)
k_allowed = 2 * np.pi * m / (N * a)
E_tb_sampled = np.sort(alpha + 2 * beta * np.cos(k_allowed * a))

print("Hückel固有値        :", huckel_eigenvalues)
print("タイトバインディング :", E_tb_sampled)
print("両者は一致するか:", np.allclose(huckel_eigenvalues, E_tb_sampled))

実行結果: 両者の値は完全に一致する（np.allclose がTrueを返す）。すなわち、ベンゼンのHückel分子軌道は、無限に長い1次元タイトバインディング鎖のバンドを、環状境界条件によって離散化された6つの $k$ 点でサンプリングしたものに他ならない。有限のπ電子系（分子）と無限に周期的な系（結晶）は、同一の物理的枠組み（LCAO/タイトバインディング）の中で連続的につながっている。

2.4 バンド理論とエネルギーバンド構造

前節の1次元鎖は、単位格子あたり1つの原子軌道しか持たない最も単純な場合でした。実際の結晶では単位格子が複数の原子（副格子）を含むことが多く、その場合はエネルギーバンドが複数本に分裂します。ここでは、単位格子あたり2つの原子（A, B）を持ち、ホッピング積分が交互に異なる値 $t_1$（分子内的な結合）と $t_2$（分子間的な結合）を取る二原子鎖モデル（Diatomic Chain Model）を考えます。このモデルは、単結合と二重結合が交互に並ぶポリアセチレンのような共役高分子の電子構造とも対応しており、Su-Schrieffer-Heeger（SSH）模型として知られています。

波数 $k$ ごとに、A副格子とB副格子の振幅を成分とする2×2のBlochハミルトニアンを対角化することで、2本のエネルギーバンドが得られます。

$t_1 = t_2$（結合が等価で、二量化していない）の場合、この2本のバンドはブリルアンゾーンの端で接触し、隙間のない単一のバンドに帰着します（前節の単原子鎖と同じ物理）。一方、$t_1 \neq t_2$ の場合はブリルアンゾーン境界（$k = \pi/a$）でエネルギーギャップが開き、その大きさは $\Delta E = 2|t_1 - t_2|$ となります。下のバンド（価電子帯, Valence Band）が電子で完全に満たされ、上のバンド（伝導帯, Conduction Band）が空である場合、この系はバンドギャップ（Band Gap）を持つ絶縁体または半導体としてふるまいます。逆に、バンドが部分的にしか占有されていない場合はフェルミ準位近傍に空の状態が存在するため電子は自由に動くことができ、系は金属的な性質を示します。

Python実装: 二原子鎖モデルのバンド構造とギャップ形成

import numpy as np
import matplotlib.pyplot as plt

def diatomic_chain_bands(k, alpha, t1, t2, a=1.0):
    """交互ホッピングを持つ二原子鎖の2バンド分散関係を計算する"""
    delta = np.sqrt(t1**2 + t2**2 + 2 * t1 * t2 * np.cos(k * a))
    return alpha + delta, alpha - delta  # 上のバンド, 下のバンド

alpha, a = 0.0, 1.0
k = np.linspace(-np.pi / a, np.pi / a, 400)

# ケース1: 均一な鎖（ギャップなし）
E_upper_uniform, E_lower_uniform = diatomic_chain_bands(k, alpha, t1=1.0, t2=1.0, a=a)

# ケース2: 二量化した鎖（ギャップあり）
t1, t2 = 1.0, 0.6
E_upper, E_lower = diatomic_chain_bands(k, alpha, t1, t2, a=a)

gap = E_upper.min() - E_lower.max()
print(f"t1={t1}, t2={t2} のときのバンドギャップ: {gap:.4f}")
print(f"解析的なギャップ 2|t1-t2|: {2 * abs(t1 - t2):.4f}")

fig, axes = plt.subplots(1, 2, figsize=(12, 5.5), sharey=True)

axes[0].plot(k * a / np.pi, E_upper_uniform, color='#f5576c', linewidth=2.5)
axes[0].plot(k * a / np.pi, E_lower_uniform, color='#f093fb', linewidth=2.5)
axes[0].set_title('均一な鎖 (t1=t2=1.0): ギャップなし', fontsize=12)
axes[0].set_xlabel('k a / π', fontsize=12)
axes[0].set_ylabel('E(k)', fontsize=12)
axes[0].grid(True, alpha=0.3)

axes[1].plot(k * a / np.pi, E_upper, color='#f5576c', linewidth=2.5, label='伝導帯')
axes[1].plot(k * a / np.pi, E_lower, color='#f093fb', linewidth=2.5, label='価電子帯')
axes[1].set_title(f'二量化した鎖 (t1={t1}, t2={t2}): ギャップあり', fontsize=12)
axes[1].set_xlabel('k a / π', fontsize=12)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('diatomic_chain_bands.png', dpi=300)
plt.show()

実行結果: $t_1 = t_2 = 1.0$ のときギャップは0であり、単一のバンドが折りたたまれた形になる。$t_1=1.0, t_2=0.6$ に二量化させると、ブリルアンゾーン境界にギャップ $\Delta E = 2|t_1-t_2| = 0.8$（β単位）が開く。この現象はPeierls転移（Peierls Transition）とも呼ばれ、1次元系が格子を歪ませることでエネルギーを安定化する重要な機構である。

2.5 状態密度（DOS）の計算と可視化

状態密度（Density of States, DOS） $g(E)$ は、単位エネルギー幅あたりに存在する電子状態の数を表す関数であり、比熱、電気伝導度、磁化率など、固体の巨視的な物性を電子状態から計算する際の基本量です。有限個の $k$ 点で離散的に定義されたエネルギー準位から、状態密度は次のようにヒストグラムとして構成できます。

1次元タイトバインディングバンド $E(k) = \alpha + 2\beta\cos(ka)$ の場合、状態密度は解析的に次の形で求まります。

Python実装: タイトバインディングバンドの状態密度（ヒストグラムと解析解の比較）

import numpy as np
import matplotlib.pyplot as plt

alpha, beta, a = 0.0, -1.0, 1.0

# 第一ブリルアンゾーンを多数のk点で密にサンプリング
n_k = 5000
k_samples = np.linspace(-np.pi / a, np.pi / a, n_k, endpoint=False)
E_samples = alpha + 2 * beta * np.cos(k_samples * a)

# ヒストグラムによる状態密度の近似
hist, bin_edges = np.histogram(E_samples, bins=60, density=True)
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])

# 解析解（バンド端付近を避けて発散を回避）
E_dense = np.linspace(2 * beta + 1e-3, -2 * beta - 1e-3, 400)
g_analytic = 1.0 / (np.pi * np.sqrt(4 * beta**2 - (E_dense - alpha)**2))

plt.figure(figsize=(9, 6))
plt.bar(bin_centers, hist, width=(bin_edges[1] - bin_edges[0]),
        color='#f093fb', alpha=0.6, label='ヒストグラム（数値）')
plt.plot(E_dense, g_analytic, color='#f5576c', linewidth=2.5, label='解析解')
plt.xlabel('エネルギー E (alpha + x*beta)', fontsize=12)
plt.ylabel('状態密度 g(E)', fontsize=12)
plt.title('1次元タイトバインディングバンドの状態密度', fontsize=14, fontweight='bold')
plt.ylim(0, 2.0)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('dos_1d_tight_binding.png', dpi=300)
plt.show()

print(f"バンド端 E={2*beta:.2f}, {-2*beta:.2f} 付近で状態密度が急増していることを確認")
print(f"バンド中央 (E=0) 付近のヒストグラム値: {hist[len(hist)//2]:.4f}")

実行結果: ヒストグラムと解析解はよく一致し、バンド中央付近（$E \approx \alpha$）では状態密度がほぼ平坦であるのに対し、バンド端（$E = \alpha \pm 2\beta$）に近づくにつれて状態密度が急激に増大する。これは群速度がバンド端でゼロに近づき、単位エネルギー幅あたりに詰め込める $k$ 点の数が増えるためである。

演習問題

import numpy as np

alpha, beta = 0.0, -1.0

# アリル系（直鎖状, N=3, 環ではない）のHückel行列
H_allyl = np.array([
    [alpha, beta, 0.0],
    [beta, alpha, beta],
    [0.0, beta, alpha]
])

eigenvalues = np.sort(np.linalg.eigvalsh(H_allyl))
print("アリル系の軌道エネルギー (alpha + x*beta):", eigenvalues)

def fill_orbitals(eigenvalues, n_electrons):
    """基底状態配置でエネルギーの低い軌道から電子を詰める（Aufbau原理）"""
    occupations = np.zeros(len(eigenvalues))
    remaining = n_electrons
    for i in range(len(eigenvalues)):
        add = min(2, remaining)
        occupations[i] = add
        remaining -= add
    return occupations

for n_electrons, species in [(2, "アリルカチオン"), (3, "アリルラジカル"), (4, "アリルアニオン")]:
    occ = fill_orbitals(eigenvalues, n_electrons)
    E_total = np.sum(occ * eigenvalues)
    print(f"{species} (π電子{n_electrons}個): 占有数={occ}, 全π電子エネルギー={E_total:.4f}")

import numpy as np

alpha, beta, a = 0.0, -1.0, 1.0

def E_k(k):
    return alpha + 2 * beta * np.cos(k * a)

# 解析解: バンド底(k=0)近傍で E(k) ≈ alpha + 2*beta - beta*(a*k)^2 + ...
# よって d2E/dk2|_(k=0) = -2*beta*a^2 、有効質量 m* = hbar^2 / (d2E/dk2)
m_eff_analytic = 1.0 / (-2 * beta * a**2)  # hbar=1の自然単位系
print(f"解析的な有効質量: m* = {m_eff_analytic:.4f}")

# 数値微分による検証
k_fine = np.linspace(-0.01, 0.01, 5)
E_fine = E_k(k_fine)
d2E_dk2 = np.gradient(np.gradient(E_fine, k_fine), k_fine)[2]
m_eff_numeric = 1.0 / d2E_dk2
print(f"数値微分による有効質量: m* = {m_eff_numeric:.4f}")

# betaの大きさを変えた場合の比較
for beta_test in [-0.5, -1.0, -2.0]:
    m_test = 1.0 / (-2 * beta_test * a**2)
    print(f"beta = {beta_test}: m* = {m_test:.4f}")

import numpy as np

alpha = 0.0
t1, t2 = 1.0, 0.6
a = 1.0

n_k = 4000
k = np.linspace(-np.pi / a, np.pi / a, n_k, endpoint=False)
delta = np.sqrt(t1**2 + t2**2 + 2 * t1 * t2 * np.cos(k * a))
E_upper = alpha + delta
E_lower = alpha - delta
E_all = np.concatenate([E_upper, E_lower])

hist, bin_edges = np.histogram(E_all, bins=100, density=True)
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])

gap_lower_edge = E_lower.max()
gap_upper_edge = E_upper.min()
gap_width = gap_upper_edge - gap_lower_edge
print(f"数値的なギャップ幅: {gap_width:.4f}")
print(f"解析的なギャップ幅 2|t1-t2|: {2 * abs(t1 - t2):.4f}")

in_gap = (bin_centers > gap_lower_edge) & (bin_centers < gap_upper_edge)
print(f"ギャップ領域内のヒストグラム値（すべて0であるはず）: {hist[in_gap]}")

参考文献

学習目標確認

← 第1章 | シリーズ目次 | 第3章 →