第4章: トポロジカルフォノン

import numpy as np
from scipy.linalg import eigh

def compute_berry_curvature(dynamical_matrix_func, kx, ky, band_index, dk=1e-5):
    """
    フォノンのベリー曲率を数値微分で計算

    Parameters:
    -----------
    dynamical_matrix_func : callable
        k点を引数にとり動力学行列を返す関数
    kx, ky : float
        運動量空間の座標
    band_index : int
        対象バンドのインデックス
    dk : float
        数値微分の刻み幅

    Returns:
    --------
    omega_z : float
        z方向のベリー曲率成分
    """
    k_points = [
        [kx + dk, ky],
        [kx - dk, ky],
        [kx, ky + dk],
        [kx, ky - dk]
    ]

    eigenvectors = []
    for kpt in k_points:
        D = dynamical_matrix_func(kpt[0], kpt[1])
        eigenvals, eigvecs = eigh(D)
        eigenvectors.append(eigvecs[:, band_index])

    # 数値微分による接続の計算
    Ax_plus = 1j * np.log(np.vdot(eigenvectors[0], eigenvectors[0]))
    Ax_minus = 1j * np.log(np.vdot(eigenvectors[1], eigenvectors[1]))
    Ay_plus = 1j * np.log(np.vdot(eigenvectors[2], eigenvectors[2]))
    Ay_minus = 1j * np.log(np.vdot(eigenvectors[3], eigenvectors[3]))

    # ベリー曲率の計算
    dAy_dx = (Ay_plus - Ay_minus) / (2 * dk)
    dAx_dy = (Ax_plus - Ax_minus) / (2 * dk)

    omega_z = dAy_dx - dAx_dy

    return np.real(omega_z)

# 使用例: 2Dハニカム格子
def honeycomb_dynamical_matrix(kx, ky):
    """ハニカム格子の簡易動力学行列"""
    a = 1.0  # 格子定数
    K = 1.0  # ばね定数

    # 最近接相互作用のみの簡易モデル
    delta_1 = np.array([0, 0])
    delta_2 = np.array([a*np.sqrt(3)/2, a/2])
    delta_3 = np.array([a*np.sqrt(3)/2, -a/2])

    D = np.zeros((6, 6), dtype=complex)

    # 動力学行列の構築（簡略化）
    phase_1 = np.exp(1j * (kx * delta_1[0] + ky * delta_1[1]))
    phase_2 = np.exp(1j * (kx * delta_2[0] + ky * delta_2[1]))
    phase_3 = np.exp(1j * (kx * delta_3[0] + ky * delta_3[1]))

    # 対角成分
    D[0:3, 0:3] = 3 * K * np.eye(3)
    D[3:6, 3:6] = 3 * K * np.eye(3)

    # 非対角成分（サブラティス間相互作用）
    coupling = K * (phase_1 + phase_2 + phase_3)
    D[0:3, 3:6] = -coupling * np.eye(3)
    D[3:6, 0:3] = -coupling * np.eye(3)

    return D

# K点近傍でのベリー曲率マップ
kx_range = np.linspace(-0.2, 0.2, 20)
ky_range = np.linspace(-0.2, 0.2, 20)

berry_curvature_map = np.zeros((len(kx_range), len(ky_range)))

for i, kx in enumerate(kx_range):
    for j, ky in enumerate(ky_range):
        berry_curvature_map[i, j] = compute_berry_curvature(
            honeycomb_dynamical_matrix, kx, ky, band_index=3
        )

print("ベリー曲率の最大値:", np.max(np.abs(berry_curvature_map)))
print("ベリー曲率の積分 (チャーン数の近似):",
      np.sum(berry_curvature_map) * (kx_range[1] - kx_range[0]) * (ky_range[1] - ky_range[0]) / (2 * np.pi))

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

def compute_edge_states_1d_chain(N_sites, hopping_t, on_site_energy, boundary='open'):
    """
    1次元鎖のエッジ状態を計算（SSH模型の類似）

    Parameters:
    -----------
    N_sites : int
        サイト数
    hopping_t : array-like [t1, t2]
        交互のホッピングパラメータ
    on_site_energy : float
        オンサイトエネルギー
    boundary : str
        'open' または 'periodic'

    Returns:
    --------
    eigenvalues, eigenvectors
    """
    t1, t2 = hopping_t
    H = np.zeros((N_sites, N_sites))

    # オンサイト項
    for i in range(N_sites):
        H[i, i] = on_site_energy

    # ホッピング項（交互）
    for i in range(N_sites - 1):
        if i % 2 == 0:
            H[i, i+1] = H[i+1, i] = t1
        else:
            H[i, i+1] = H[i+1, i] = t2

    # 周期境界条件
    if boundary == 'periodic':
        H[0, N_sites-1] = H[N_sites-1, 0] = t2 if (N_sites-1) % 2 == 1 else t1

    eigenvalues, eigenvectors = eigh(H)
    return eigenvalues, eigenvectors

# トポロジカル相のパラメータ
N = 100
t1, t2 = 1.0, 0.5  # t1 > t2: トポロジカル相
E_onsite = 0.0

eigenvals_open, eigvecs_open = compute_edge_states_1d_chain(
    N, [t1, t2], E_onsite, boundary='open'
)

# エッジ状態の可視化
fig, axes = plt.subplots(2, 1, figsize=(10, 8))

# エネルギースペクトル
axes[0].scatter(range(N), eigenvals_open, s=10)
axes[0].axhline(y=0, color='r', linestyle='--', label='ゼロエネルギー')
axes[0].set_xlabel('状態インデックス')
axes[0].set_ylabel('エネルギー')
axes[0].set_title('エネルギースペクトル（開境界）')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ゼロエネルギー近傍の状態の空間分布
zero_energy_states = np.where(np.abs(eigenvals_open) < 0.1)[0]
if len(zero_energy_states) >= 2:
    # 左エッジ状態
    left_edge_state = eigvecs_open[:, zero_energy_states[0]]
    # 右エッジ状態
    right_edge_state = eigvecs_open[:, zero_energy_states[-1]]

    axes[1].plot(np.abs(left_edge_state)**2, 'b-', label='左エッジ状態', linewidth=2)
    axes[1].plot(np.abs(right_edge_state)**2, 'r-', label='右エッジ状態', linewidth=2)
    axes[1].set_xlabel('サイト位置')
    axes[1].set_ylabel('確率密度 |ψ|²')
    axes[1].set_title('エッジ状態の空間分布')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('edge_states_ssh.png', dpi=150, bbox_inches='tight')
print(f"ゼロエネルギー状態数: {len(zero_energy_states)}")
print(f"トポロジカル相: {'Yes' if t1 > t2 else 'No'}")

import numpy as np
from scipy.integrate import simps

def thermal_factor_c2(x):
    """熱ホール効果の熱因子 c2(x)"""
    # x = ℏω / (k_B T)
    exp_x = np.exp(x)
    bose = x / (exp_x - 1)
    c2 = (1 + x) * bose**2 - bose
    return c2

def compute_thermal_hall_conductivity(omega_bands, berry_curvature, temperature, k_mesh):
    """
    熱ホール伝導率の計算

    Parameters:
    -----------
    omega_bands : ndarray, shape (n_bands, nkx, nky)
        フォノン分散 (rad/s)
    berry_curvature : ndarray, shape (n_bands, nkx, nky)
        ベリー曲率 (Ų)
    temperature : float
        温度 (K)
    k_mesh : tuple (kx_array, ky_array)
        運動量メッシュ

    Returns:
    --------
    kappa_xy : float
        熱ホール伝導率 (W/K·m)
    """
    hbar = 1.054571817e-34  # J·s
    kB = 1.380649e-23       # J/K

    kx, ky = k_mesh
    dk = (kx[1] - kx[0]) * (ky[1] - ky[0])  # k空間の面積要素

    kappa_xy = 0.0

    for band_idx in range(omega_bands.shape[0]):
        omega = omega_bands[band_idx]
        Omega = berry_curvature[band_idx]

        # x = ℏω / (k_B T)
        x = hbar * omega / (kB * temperature)

        # 熱因子の計算
        c2_factor = thermal_factor_c2(x)

        # 積分（シンプソン則）
        integrand = c2_factor * Omega
        integral = simps(simps(integrand, ky), kx)

        kappa_xy += integral

    # 係数
    prefactor = (kB**2 * temperature) / (hbar * (2 * np.pi)**2)
    kappa_xy *= prefactor * dk

    return kappa_xy

# 使用例: 2Dハニカム格子
nk = 50
kx = np.linspace(-np.pi, np.pi, nk)
ky = np.linspace(-np.pi, np.pi, nk)
KX, KY = np.meshgrid(kx, ky)

# 簡易的なフォノン分散（2バンド）
v_sound = 5000  # m/s (音速)
hbar = 1.054571817e-34
omega_acoustic = v_sound * np.sqrt(KX**2 + KY**2) / 1e-10  # rad/s

# 簡易的なベリー曲率（ディラック点周りにピーク）
berry_curv = np.zeros_like(KX)
K_point_x, K_point_y = 4*np.pi/3, 0
distance = np.sqrt((KX - K_point_x)**2 + (KY - K_point_y)**2)
berry_curv = 0.5 * np.exp(-distance**2 / 0.1)  # Ų単位

omega_bands = np.array([omega_acoustic])
berry_curvature = np.array([berry_curv])

# 温度依存性
temperatures = np.logspace(0, 2, 20)  # 1K to 100K
kappa_xy_values = []

for T in temperatures:
    kappa_xy = compute_thermal_hall_conductivity(
        omega_bands, berry_curvature, T, (kx, ky)
    )
    kappa_xy_values.append(kappa_xy)

# 結果のプロット
import matplotlib.pyplot as plt

plt.figure(figsize=(8, 6))
plt.semilogx(temperatures, kappa_xy_values, 'o-', linewidth=2)
plt.xlabel('温度 (K)', fontsize=12)
plt.ylabel('熱ホール伝導率 κ_xy (W/K·m)', fontsize=12)
plt.title('熱ホール伝導率の温度依存性', fontsize=14)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('thermal_hall_conductivity.png', dpi=150)

print(f"最大値: {np.max(kappa_xy_values):.2e} W/K·m at {temperatures[np.argmax(kappa_xy_values)]:.1f} K")

import numpy as np
from scipy.linalg import eigh

def graphene_phonon_hamiltonian(kx, ky, a=1.42e-10, K_spring=100):
    """
    グラフェンの簡易フォノンハミルトニアン

    Parameters:
    -----------
    kx, ky : float
        運動量 (1/m)
    a : float
        炭素間距離 (m)
    K_spring : float
        ばね定数 (N/m)

    Returns:
    --------
    H : ndarray, shape (6, 6)
        動力学行列
    """
    # 最近接ベクトル
    delta1 = np.array([a, 0])
    delta2 = np.array([-a/2, a*np.sqrt(3)/2])
    delta3 = np.array([-a/2, -a*np.sqrt(3)/2])

    # 構造因子
    f_k = (np.exp(1j * (kx * delta1[0] + ky * delta1[1])) +
           np.exp(1j * (kx * delta2[0] + ky * delta2[1])) +
           np.exp(1j * (kx * delta3[0] + ky * delta3[1])))

    mass = 12 * 1.66054e-27  # 炭素の質量 (kg)

    # 動力学行列（簡略版: x-y平面内のみ）
    H = np.zeros((6, 6), dtype=complex)

    # サブラティスA
    H[0:3, 0:3] = (3 * K_spring / mass) * np.eye(3)
    # サブラティスB
    H[3:6, 3:6] = (3 * K_spring / mass) * np.eye(3)

    # サブラティス間相互作用
    coupling = -(K_spring / mass) * f_k
    H[0, 3] = H[3, 0] = coupling
    H[1, 4] = H[4, 1] = coupling
    H[2, 5] = H[5, 2] = coupling

    return H

def compute_berry_curvature_graphene(kx_range, ky_range, band_idx=2):
    """グラフェンK点近傍のベリー曲率マップ"""
    nkx, nky = len(kx_range), len(ky_range)
    berry_map = np.zeros((nkx, nky))

    dk = 1e-5  # 数値微分の刻み

    for i, kx in enumerate(kx_range):
        for j, ky in enumerate(ky_range):
            # 4点での固有ベクトル
            H_plus_x = graphene_phonon_hamiltonian(kx + dk, ky)
            H_minus_x = graphene_phonon_hamiltonian(kx - dk, ky)
            H_plus_y = graphene_phonon_hamiltonian(kx, ky + dk)
            H_minus_y = graphene_phonon_hamiltonian(kx, ky - dk)

            _, v_px = eigh(H_plus_x)
            _, v_mx = eigh(H_minus_x)
            _, v_py = eigh(H_plus_y)
            _, v_my = eigh(H_minus_y)

            # ベリー接続の計算
            vec_px = v_px[:, band_idx]
            vec_mx = v_mx[:, band_idx]
            vec_py = v_py[:, band_idx]
            vec_my = v_my[:, band_idx]

            # 位相の連続性を保つ
            if np.vdot(vec_px, vec_mx) < 0:
                vec_mx = -vec_mx
            if np.vdot(vec_py, vec_my) < 0:
                vec_my = -vec_my

            # ベリー接続（簡易計算）
            Ax = np.imag(np.log(np.vdot(vec_px, vec_mx) + 1e-10))
            Ay = np.imag(np.log(np.vdot(vec_py, vec_my) + 1e-10))

            # ベリー曲率（回転）
            berry_map[i, j] = (Ay - Ax) / (2 * dk)

    return berry_map

# K点近傍での計算
a = 1.42e-10
K_point = 4 * np.pi / (3 * a)
kx_range = np.linspace(K_point - 1e9, K_point + 1e9, 30)
ky_range = np.linspace(-1e9, 1e9, 30)

berry_curv_graphene = compute_berry_curvature_graphene(kx_range, ky_range)

# 可視化
import matplotlib.pyplot as plt

plt.figure(figsize=(8, 6))
plt.contourf(kx_range / 1e9, ky_range / 1e9, berry_curv_graphene.T,
             levels=20, cmap='RdBu_r')
plt.colorbar(label='ベリー曲率 (m²)')
plt.xlabel('kx (×10⁹ m⁻¹)', fontsize=12)
plt.ylabel('ky (×10⁹ m⁻¹)', fontsize=12)
plt.title('グラフェンK点近傍のベリー曲率', fontsize=14)
plt.tight_layout()
plt.savefig('graphene_berry_curvature.png', dpi=150)

print(f"ベリー曲率の積分（チャーン数）: {np.sum(berry_curv_graphene) * (kx_range[1]-kx_range[0]) * (ky_range[1]-ky_range[0]) / (2*np.pi):.3f}")

# Quantum ESPRESSOでの計算フロー（疑似コード）

"""
1. SCF計算（自己無撞着場）
   pw.x < scf.in > scf.out

2. DFPT計算（動力学行列）
   ph.x < ph.in > ph.out

3. フォノン分散の計算
   q2r.x < q2r.in > q2r.out
   matdyn.x < matdyn.in > matdyn.out

4. ベリー曲率の計算（WannierTools等）
"""

# Pythonでの解析例
import numpy as np
import matplotlib.pyplot as plt

def load_phonon_bands(filename):
    """Quantum ESPRESSOの出力からフォノンバンドを読み込み"""
    data = np.loadtxt(filename)
    k_path = data[:, 0]
    frequencies = data[:, 1:]
    return k_path, frequencies

def identify_weyl_points(k_mesh, frequencies, threshold=0.1):
    """
    バンド交差点（ワイル点候補）を検出

    Parameters:
    -----------
    k_mesh : ndarray
        k点メッシュ
    frequencies : ndarray
        フォノン周波数
    threshold : float
        バンドギャップの閾値 (THz)

    Returns:
    --------
    weyl_candidates : list
        ワイル点候補のk座標
    """
    weyl_candidates = []

    # 隣接バンド間のギャップをチェック
    for i in range(frequencies.shape[1] - 1):
        gap = frequencies[:, i+1] - frequencies[:, i]
        crossing_indices = np.where(gap < threshold)[0]

        for idx in crossing_indices:
            weyl_candidates.append(k_mesh[idx])

    return weyl_candidates

# 使用例
k_path, phonon_freqs = load_phonon_bands('phonon_bands.dat')

# バンド構造のプロット
fig, ax = plt.subplots(figsize=(10, 6))

for band_idx in range(phonon_freqs.shape[1]):
    ax.plot(k_path, phonon_freqs[:, band_idx], 'b-', linewidth=1.5)

ax.axhline(y=0, color='k', linestyle='--', linewidth=0.5)
ax.set_xlabel('波数ベクトル', fontsize=12)
ax.set_ylabel('周波数 (THz)', fontsize=12)
ax.set_title('フォノンバンド構造', fontsize=14)
ax.grid(True, alpha=0.3)

# ワイル点候補をマーク
weyl_points = identify_weyl_points(k_path, phonon_freqs, threshold=0.1)
if weyl_points:
    for wp in weyl_points:
        ax.axvline(x=wp, color='r', linestyle=':', alpha=0.5)
    print(f"ワイル点候補数: {len(weyl_points)}")

plt.tight_layout()
plt.savefig('phonon_bands_with_weyl.png', dpi=150)

import numpy as np
from scipy.constants import hbar, k as kB

def bose_einstein(omega, T):
    """ボース・アインシュタイン分布"""
    x = hbar * omega / (kB * T)
    return 1 / (np.exp(x) - 1)

def heat_capacity_phonon(omega, T):
    """フォノンの比熱"""
    x = hbar * omega / (kB * T)
    n_BE = bose_einstein(omega, T)
    C = kB * x**2 * np.exp(x) / (np.exp(x) - 1)**2
    return C

def thermal_conductivity_RTA(omega_bands, group_velocity, relaxation_time,
                               k_mesh, temperature, volume):
    """
    緩和時間近似（RTA）での熱伝導率計算

    Parameters:
    -----------
    omega_bands : ndarray, shape (n_bands, n_k)
        フォノン角周波数 (rad/s)
    group_velocity : ndarray, shape (n_bands, n_k, 3)
        群速度 (m/s)
    relaxation_time : ndarray, shape (n_bands, n_k)
        緩和時間 (s)
    k_mesh : ndarray, shape (n_k, 3)
        k点メッシュ
    temperature : float
        温度 (K)
    volume : float
        単位胞体積 (m³)

    Returns:
    --------
    kappa : ndarray, shape (3, 3)
        熱伝導率テンソル (W/m·K)
    """
    n_bands, n_k = omega_bands.shape
    dk = (2 * np.pi)**3 / (volume * n_k)  # k空間の体積要素

    kappa = np.zeros((3, 3))

    for band_idx in range(n_bands):
        for k_idx in range(n_k):
            omega = omega_bands[band_idx, k_idx]
            v = group_velocity[band_idx, k_idx]
            tau = relaxation_time[band_idx, k_idx]

            if omega < 1e-6:  # 音響ブランチのゼロ点を除外
                continue

            # 比熱
            C = heat_capacity_phonon(omega, temperature)

            # 熱伝導率テンソルへの寄与
            for alpha in range(3):
                for beta in range(3):
                    kappa[alpha, beta] += C * v[alpha] * v[beta] * tau * dk / volume

    return kappa

# 使用例: 簡単な等方的モデル
n_k = 1000
n_bands = 3

# 簡易的なフォノン分散（デバイモデル）
k_magnitude = np.linspace(0, np.pi / 3e-10, n_k)
omega_debye = 3e13  # rad/s
omega_bands_1d = np.zeros((n_bands, n_k))

# 音響モード
omega_bands_1d[0] = omega_debye * k_magnitude / k_magnitude[-1] * 0.3
omega_bands_1d[1] = omega_debye * k_magnitude / k_magnitude[-1] * 0.5
# 光学モード
omega_bands_1d[2] = omega_debye * 0.8 * np.ones(n_k)

# 群速度（簡易: 分散の微分）
group_vel_1d = np.zeros((n_bands, n_k, 3))
dk_spacing = k_magnitude[1] - k_magnitude[0]
for band in range(n_bands):
    dw_dk = np.gradient(omega_bands_1d[band], dk_spacing)
    group_vel_1d[band, :, 0] = dw_dk  # x方向のみ

# 緩和時間（簡易: ウムクラッププロセス）
T_ref = 300  # K
tau_0 = 1e-11  # s
relaxation_time_1d = tau_0 * (T_ref / 100)**2 * np.ones((n_bands, n_k))

# 温度依存性の計算
temperatures = np.linspace(10, 400, 40)
kappa_values = []

volume = (3e-10)**3  # 典型的な単位胞体積

for T in temperatures:
    kappa_tensor = thermal_conductivity_RTA(
        omega_bands_1d, group_vel_1d, relaxation_time_1d,
        k_magnitude, T, volume
    )
    kappa_values.append(kappa_tensor[0, 0])

# プロット
import matplotlib.pyplot as plt

plt.figure(figsize=(8, 6))
plt.plot(temperatures, kappa_values, 'o-', linewidth=2)
plt.xlabel('温度 (K)', fontsize=12)
plt.ylabel('熱伝導率 (W/m·K)', fontsize=12)
plt.title('熱伝導率の温度依存性（RTA）', fontsize=14)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('thermal_conductivity_RTA.png', dpi=150)

print(f"室温（300K）での熱伝導率: {kappa_values[np.argmin(np.abs(temperatures - 300))]:.2f} W/m·K")

import numpy as np
from scipy.linalg import eigh, det

def wilson_loop_1d(dynamical_matrix_func, k_loop, occupied_bands):
    """
    1次元ウィルソンループの計算

    Parameters:
    -----------
    dynamical_matrix_func : callable
        k点を引数にとり動力学行列を返す関数
    k_loop : ndarray, shape (n_k, 2)
        ブリルアンゾーン上の閉路
    occupied_bands : list
        占有バンドのインデックスリスト

    Returns:
    --------
    wilson_eigenvalues : ndarray
        ウィルソンループの固有値（位相）
    """
    n_k = len(k_loop)
    n_occ = len(occupied_bands)

    # 初期化
    F = np.eye(n_occ, dtype=complex)

    # 閉路に沿って積分
    for i in range(n_k):
        k_current = k_loop[i]
        k_next = k_loop[(i + 1) % n_k]

        # 現在と次のk点での固有ベクトル
        D_current = dynamical_matrix_func(k_current[0], k_current[1])
        D_next = dynamical_matrix_func(k_next[0], k_next[1])

        _, v_current = eigh(D_current)
        _, v_next = eigh(D_next)

        # 占有バンドの部分空間
        U_current = v_current[:, occupied_bands]
        U_next = v_next[:, occupied_bands]

        # 重なり行列
        overlap = U_next.conj().T @ U_current

        # ウィルソンループ行列の更新
        F = overlap @ F

    # 固有値（位相）を計算
    wilson_eigenvalues = np.linalg.eigvals(F)
    wilson_phases = np.angle(wilson_eigenvalues)

    return wilson_phases

def compute_chern_number_wilson(dynamical_matrix_func, k_mesh_2d, occupied_bands):
    """
    2次元系のチャーン数をウィルソンループ法で計算

    Parameters:
    -----------
    dynamical_matrix_func : callable
        動力学行列を返す関数
    k_mesh_2d : tuple (kx_array, ky_array)
        2次元k点メッシュ
    occupied_bands : list
        占有バンドのインデックス

    Returns:
    --------
    chern_number : float
        チャーン数
    """
    kx_array, ky_array = k_mesh_2d
    n_kx = len(kx_array)
    n_ky = len(ky_array)

    # ky方向のウィルソンループ固有値を計算
    wilson_phases = np.zeros((n_kx, len(occupied_bands)))

    for i, kx in enumerate(kx_array):
        # kyループを構築
        k_loop = np.array([[kx, ky] for ky in ky_array])
        # ウィルソンループ計算
        phases = wilson_loop_1d(dynamical_matrix_func, k_loop, occupied_bands)
        wilson_phases[i] = phases

    # kx方向でウィルソンループ固有値の巻き数を計算
    total_winding = 0
    for band_idx in range(len(occupied_bands)):
        phase_evolution = wilson_phases[:, band_idx]
        # 位相の変化を追跡（アンラップ）
        phase_unwrapped = np.unwrap(phase_evolution)
        winding = (phase_unwrapped[-1] - phase_unwrapped[0]) / (2 * np.pi)
        total_winding += winding

    chern_number = np.round(total_winding)
    return chern_number

# 使用例
def example_haldane_phonon(kx, ky, t1=1.0, t2=0.3, phi=np.pi/2):
    """Haldane模型風のトポロジカルフォノン"""
    # 簡略化されたハルデン型モデル
    H = np.zeros((6, 6), dtype=complex)

    # 最近接ホッピング
    delta_1 = [1, 0]
    delta_2 = [-0.5, np.sqrt(3)/2]
    delta_3 = [-0.5, -np.sqrt(3)/2]

    phase_1 = np.exp(1j * (kx * delta_1[0] + ky * delta_1[1]))
    phase_2 = np.exp(1j * (kx * delta_2[0] + ky * delta_2[1]))
    phase_3 = np.exp(1j * (kx * delta_3[0] + ky * delta_3[1]))

    # オンサイト
    H[0:3, 0:3] = 3 * t1 * np.eye(3)
    H[3:6, 3:6] = 3 * t1 * np.eye(3)

    # サブラティス間
    coupling = t1 * (phase_1 + phase_2 + phase_3)
    H[0:3, 3:6] = -coupling * np.eye(3)
    H[3:6, 0:3] = -np.conj(coupling) * np.eye(3)

    # 次近接ホッピング（時間反転対称性を破る）
    nnn_phase = np.exp(1j * phi)
    H[0:3, 0:3] += t2 * nnn_phase * np.eye(3)
    H[3:6, 3:6] += t2 * np.conj(nnn_phase) * np.eye(3)

    return H

# チャーン数の計算
kx_mesh = np.linspace(-np.pi, np.pi, 20)
ky_mesh = np.linspace(-np.pi, np.pi, 20)

chern_num = compute_chern_number_wilson(
    example_haldane_phonon,
    (kx_mesh, ky_mesh),
    occupied_bands=[0, 1, 2]
)

print(f"計算されたチャーン数: {chern_num:.0f}")