Odd-parity magnets toy model

考虑一条一维链，原始晶格常数为 $a$，每个格点包含一个轨道和两个真实自旋自由度。模型的磁性周期为四个原始格点，因此一个磁性单胞包含四个子格

$$A,\ B,\ C,\ D .$$

用 $R$ 标记磁性单胞，用 $s=0,1,2,3$ 标记单胞内子格，即

$$s=0,1,2,3 \quad\leftrightarrow\quad A,B,C,D .$$

物理格点编号为$j=4R+s$，相应电子算符记为

$$c_{Rs\sigma}, \qquad \sigma=\uparrow,\downarrow .$$

磁性单胞长度为$a_{\rm M}=4a$，模型的四个子格局域磁矩取为

$$\mathbf m_A=\frac{1}{\sqrt3}(1,1,1), \qquad \mathbf m_B=\frac{1}{\sqrt3}(1,1,-1),\qquad \mathbf m_C=\frac{1}{\sqrt3}(-1,-1,-1), \qquad \mathbf m_D=\frac{1}{\sqrt3}(-1,-1,1).$$

它们满足

$$\mathbf m_C=-\mathbf m_A, \qquad \mathbf m_D=-\mathbf m_B,$$

也就是$\mathbf m_{s+2}=-\mathbf m_s$。因此，平移两个原始格点会使局域磁矩整体反号，这是该模型最关键的实空间结构。

模型哈密顿量写为$H=H_t+H_J$，最近邻自旋无关 hopping 为

$$H_t = t\sum_{j,\sigma} \left( c_{j\sigma}^\dagger c_{j+1,\sigma} + c_{j+1,\sigma}^\dagger c_{j\sigma} \right).$$

用磁性单胞表示，

$$\begin{aligned} H_t = t\sum_{R,\sigma} \Big( &c_{RA\sigma}^\dagger c_{RB\sigma} + c_{RB\sigma}^\dagger c_{RC\sigma} + c_{RC\sigma}^\dagger c_{RD\sigma} + c_{RD\sigma}^\dagger c_{R+1,A,\sigma} + {\rm H.c.} \Big). \end{aligned}$$

前三项是磁性单胞内部跃迁，最后一项是跨磁性单胞跃迁。局域交换项为

$$H_J = J\sum_{R,s} c_{Rs}^\dagger \left( \mathbf m_s\cdot\boldsymbol\sigma \right) c_{Rs},$$

其中

$$c_{Rs} = \begin{pmatrix} c_{Rs\uparrow}\\ c_{Rs\downarrow} \end{pmatrix}, \qquad \boldsymbol\sigma=(\sigma_x,\sigma_y,\sigma_z).$$

于是四个子格上的 onsite exchange matrix 分别为

$$\begin{aligned} &M_A = \frac{J}{\sqrt3} (\sigma_x+\sigma_y+\sigma_z),\quad M_B = \frac{J}{\sqrt3} (\sigma_x+\sigma_y-\sigma_z),\\ &M_C = -\frac{J}{\sqrt3} (\sigma_x+\sigma_y+\sigma_z),\quad M_D = \frac{J}{\sqrt3} (-\sigma_x-\sigma_y+\sigma_z). \end{aligned}$$

由于磁性周期为 $4a$，动量空间应使用磁性 Brillouin zone。定义无量纲磁性单胞动量

$$K=4k_xa, \qquad K\in[-\pi,\pi].$$

采用元胞周期规范

$$c_{Rs\sigma} = \frac{1}{\sqrt N} \sum_K e^{iKR} c_{Ks\sigma},\qquad c_{Ks\sigma} = \frac{1}{\sqrt N} \sum_R e^{-iKR} c_{Rs\sigma}.$$

这个规范只把磁性单胞编号 $R$ 放入 Fourier 相位，而不把子格内部位置 $sa$ 放入指数中。因此单胞内 hopping 不带动量相位，只有跨单胞 hopping 带 $e^{\pm iK}$。

以单胞内 $A_R\rightarrow B_R$ hopping 为例，

$$t\sum_R c_{RA\sigma}^\dagger c_{RB\sigma} = \frac{t}{N} \sum_{R,K,K'} e^{-iKR}e^{iK'R} c_{KA\sigma}^\dagger c_{K'B\sigma}.$$

利用

$$\sum_R e^{i(K'-K)R} = N\delta_{K,K'},$$

得到

$$t\sum_K c_{KA\sigma}^\dagger c_{KB\sigma}.$$

因此单胞内 hopping 给出

$$A\leftrightarrow B:\quad t\sigma_0,\quad B\leftrightarrow C:\quad t\sigma_0,\quad C\leftrightarrow D:\quad t\sigma_0.$$

跨单胞 hopping 为

$$t\sum_R c_{RD\sigma}^\dagger c_{R+1,A,\sigma}.$$

代入 Bloch 变换：

$$c_{RD\sigma}^\dagger = \frac{1}{\sqrt N} \sum_K e^{-iKR} c_{KD\sigma}^\dagger,\qquad c_{R+1,A,\sigma} = \frac{1}{\sqrt N} \sum_{K'} e^{iK'(R+1)} c_{K'A\sigma}.$$

于是

$$\begin{aligned} t\sum_R c_{RD\sigma}^\dagger c_{R+1,A,\sigma} &= \frac{t}{N} \sum_{R,K,K'} e^{-iKR}e^{iK'(R+1)} c_{KD\sigma}^\dagger c_{K'A\sigma} = t\sum_K e^{iK} c_{KD\sigma}^\dagger c_{KA\sigma}. \end{aligned}$$

因此跨单胞项给出

$$D_R\rightarrow A_{R+1}: \quad t e^{iK}\sigma_0,$$

其厄米共轭给出

$$A_{R+1}\rightarrow D_R: \quad t e^{-iK}\sigma_0.$$

取 Bloch 基矢

$$\Psi_K = \left( c_{KA\uparrow}, c_{KA\downarrow}, c_{KB\uparrow}, c_{KB\downarrow}, c_{KC\uparrow}, c_{KC\downarrow}, c_{KD\uparrow}, c_{KD\downarrow} \right)^T.$$

则

$$H = \sum_K \Psi_K^\dagger h(K) \Psi_K.$$

动量空间哈密顿量为

$$h(K) = \begin{pmatrix} M_A & t\sigma_0 & 0 & t e^{-iK}\sigma_0\\ t\sigma_0 & M_B & t\sigma_0 & 0\\ 0 & t\sigma_0 & M_C & t\sigma_0\\ t e^{iK}\sigma_0 & 0 & t\sigma_0 & M_D \end{pmatrix}.$$

这里每一个矩阵元都是 $2\times2$ 自旋矩阵，因此 $h(K)$ 是一个 $8\times8$ 矩阵。对角化

$$h(K)|u_{nK}\rangle = E_n(K)|u_{nK}\rangle$$

即可得到能带。第 $n$ 条能带的自旋期望值定义为

$$s_\alpha^{(n)}(K) = \langle u_{nK}| \mathcal S_\alpha |u_{nK}\rangle, \qquad \alpha=x,y,z,$$

其中完整子格-自旋空间中的自旋算符为$\mathcal S_\alpha
=
I_4\otimes\sigma_\alpha$，因此

$$s_\alpha^{(n)}(K) = u_{nK}^\dagger \left( I_4\otimes\sigma_\alpha \right) u_{nK}.$$

如果某个 $K$ 点存在严格简并，则单条能带的自旋期望值可能依赖于简并子空间中本征矢的选取。此时更严格的做法是在简并子空间内重新对角化自旋算符，或者只讨论非简并区域的自旋纹理。

该模型最重要的对称性有两个。第一是半磁性周期平移与时间反演的组合。定义

$$\tau_2:\quad j\rightarrow j+2.$$

由于

$$\mathbf m_{j+2}=-\mathbf m_j,$$

单独的 $\tau_2$ 不是磁结构的对称性。真实时间反演 $\mathcal T$ 使电子自旋反号：

$$\mathcal T\boldsymbol\sigma\mathcal T^{-1} = -\boldsymbol\sigma.$$

因此交换项在 $\tau_2$ 下变为

$$\mathbf m_{j+2}\cdot\boldsymbol\sigma = -\mathbf m_j\cdot\boldsymbol\sigma,$$

再经过 $\mathcal T$ 后，自旋矩阵反号，于是重新回到

$$\mathbf m_j\cdot\boldsymbol\sigma.$$

所以组合反幺正操作

$$\widetilde{\mathcal T} = \mathcal T\tau_2$$

是体系的对称性，并满足

$$\widetilde{\mathcal T} h(K) \widetilde{\mathcal T}^{-1} = h(-K).$$

因为 $\tau_2$ 是半个磁性单胞平移，所以$\tau_2^2=\tau_4$，其中 $\tau_4$ 是完整磁性单胞平移。在 Bloch 态上，完整磁性单胞平移给出相位 $e^{iK}$ 或 $e^{-iK}$，具体符号取决于 Bloch 变换约定。由于自旋 $1/2$ 电子满足

$$\mathcal T^2=-1,$$

于是

$$\widetilde{\mathcal T}^2 = -e^{iK},$$

或者等价地写为

$$\widetilde{\mathcal T}^2 = -e^{-iK}.$$

物理结论是

$$K=0: \qquad \widetilde{\mathcal T}^2=-1,$$

因此有 Kramers-like 简并；而

$$K=\pi: \qquad \widetilde{\mathcal T}^2=+1,$$

因此不再强制 Kramers-like 简并。也就是说，真实时间反演已经被磁序破坏，但组合反幺正对称性仍然可以在特定动量点保护简并。另一些 TRIM 点则允许简并被打开，这可以理解为 hidden Zeeman field 的能带表现。

第二个关键对称性是空间反演与自旋旋转的组合。取一维反演 $P$，使其在四子格之间实现

$$P:\quad A\leftrightarrow D, \qquad B\leftrightarrow C.$$

再考虑绕 $z$ 轴的自旋 $\pi$ 旋转$U_z(\pi)$，它对自旋分量的作用为

$$\sigma_x\rightarrow-\sigma_x, \qquad \sigma_y\rightarrow-\sigma_y, \qquad \sigma_z\rightarrow\sigma_z.$$

对磁矩向量也相应有

$$(m_x,m_y,m_z) \rightarrow (-m_x,-m_y,m_z).$$

逐个检查四个磁矩：

$$\begin{aligned} &U_z(\pi)\mathbf m_A = \frac{1}{\sqrt3}(-1,-1,1) = \mathbf m_D,\qquad U_z(\pi)\mathbf m_B = \frac{1}{\sqrt3}(-1,-1,-1) = \mathbf m_C,\\ &U_z(\pi)\mathbf m_C = \frac{1}{\sqrt3}(1,1,-1) = \mathbf m_B,\qquad U_z(\pi)\mathbf m_D = \frac{1}{\sqrt3}(1,1,1) = \mathbf m_A. \end{aligned}$$

因此 $P$ 交换子格，$U_z(\pi)$ 旋转自旋，二者组合后保持磁结构不变。这个组合对称性可以写为$\{U_z(\pi)\Vert P\}$，它在动量空间中使$K\rightarrow -K$，同时使自旋分量变换为

$$(s_x,s_y,s_z) \rightarrow (-s_x,-s_y,s_z).$$

因此它给出约束

$$s_x(K)=-s_x(-K),\quad s_y(K)=-s_y(-K),\quad s_z(K)=s_z(-K).$$

另一方面，$\widetilde{\mathcal T}=\mathcal T\tau_2$ 包含时间反演，因此它使自旋反号，并把 $K$ 变为 $-K$。所以它给出

$$\mathbf s(K)=-\mathbf s(-K),$$

即

$$s_x(K)=-s_x(-K),\quad s_y(K)=-s_y(-K),\quad s_z(K)=-s_z(-K).$$

合并两个对称性的约束。对于 $x,y$ 分量，两者都允许

$$s_x(K)=-s_x(-K), \qquad s_y(K)=-s_y(-K).$$

但对于 $z$ 分量，一个对称性要求

$$s_z(K)=-s_z(-K),$$

另一个要求

$$s_z(K)=s_z(-K).$$

二者同时成立只能给出

$$s_z(K)=0.$$

因此该模型的动量空间自旋纹理满足

$$\boxed{ s_z(K)=0, \qquad s_x(K),s_y(K)\ \text{are odd in }K. }$$

局域磁序本身是非共面的，因为四个 $\mathbf m_s$ 不在同一个平面内；但动量空间能带的自旋纹理被对称性限制在 $xy$ 平面内。因此该 toy model 体现的是

$$\text{noncoplanar real-space magnetic texture} \quad \Longrightarrow \quad \text{coplanar odd-parity spin texture in momentum space}.$$

简言之，四子格磁序破坏真实时间反演，但保留 $\mathcal T\tau_2$，使动量空间自旋纹理具有 odd parity；同时保留 $\{U_z(\pi)\Vert P\}$，进一步强制 $s_z(K)=0$。因此，该模型不是普通反铁磁链，而是由磁性半平移和自旋-空间组合对称性共同约束的 odd-parity magnetic chain。

代码

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
Band structure of the 1D odd-parity magnets toy model.

Magnetic unit cell:
    A, B, C, D

Local moments:
    m_A = ( 1,  1,  1) / sqrt(3)
    m_B = ( 1,  1, -1) / sqrt(3)
    m_C = (-1, -1, -1) / sqrt(3)
    m_D = (-1, -1,  1) / sqrt(3)

Hamiltonian:
    H = t sum_<ij>,s c^dag_{i s} c_{j s}
        + J sum_i c^dag_i [m_i · sigma] c_i

Momentum:
    K = 4 k_x a

Bloch Hamiltonian:
    h(K) =
    [ M_A          t I        0          t e^{-iK} I ]
    [ t I          M_B        t I        0           ]
    [ 0            t I        M_C        t I         ]
    [ t e^{iK} I   0          t I        M_D         ]

Outputs:
    output_opm_toy_band/data/bands_and_spin.npz
    output_opm_toy_band/figures/band_plain.png
    output_opm_toy_band/figures/band_spin_x.png
    output_opm_toy_band/figures/band_spin_y.png
    output_opm_toy_band/figures/band_spin_z.png
"""

import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection


# ============================================================
# 0. User parameters
# ============================================================

RUN_CALCULATION = True
RUN_PLOTTING = True

# Model parameters
T_HOP = 1.0
JEX = 0.5

# k mesh in magnetic Brillouin zone K in [-pi, pi]
NK = 1001

# Optional chemical-potential guide line
DRAW_MU_LINE = True
MU = -1.0

# Output folders
OUTPUT_DIR = "output_opm_toy_band"
DATA_DIR = os.path.join(OUTPUT_DIR, "data")
FIG_DIR = os.path.join(OUTPUT_DIR, "figures")

SAVE_FIGURES = True
SHOW_FIGURES = True
DPI = 300

# Plot windows
E_WINDOW = (-2.4, 2.4)

# Spin-colored band plot controls
SPIN_CMAP = "bwr"
SPIN_VMIN = -1.0
SPIN_VMAX = 1.0


# ============================================================
# 1. Pauli matrices
# ============================================================

sigma_0 = np.array([[1, 0], [0, 1]], dtype=complex)
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
sigma_y = np.array([[0, -1j], [1j, 0]], dtype=complex)
sigma_z = np.array([[1, 0], [0, -1]], dtype=complex)

SIGMA = {
    "x": sigma_x,
    "y": sigma_y,
    "z": sigma_z,
}


# ============================================================
# 2. Magnetic texture
# ============================================================

def local_moments():
    """
    Return local magnetic moments in one 4-site magnetic unit cell.

    Order:
        A, B, C, D
    """
    q = 1.0 / np.sqrt(3.0)

    moments = np.array([
        [ q,  q,  q],   # A
        [ q,  q, -q],   # B
        [-q, -q, -q],   # C
        [-q, -q,  q],   # D
    ], dtype=float)

    return moments


# ============================================================
# 3. Bloch Hamiltonian
# ============================================================

def bloch_hamiltonian(K, t_hop, jex):
    """
    Construct the 8 x 8 Bloch Hamiltonian h(K).

    Basis:
        (A up, A down,
         B up, B down,
         C up, C down,
         D up, D down)

    K = 4 k_x a.
    """
    moments = local_moments()
    n_sub = len(moments)
    dim = 2 * n_sub

    H = np.zeros((dim, dim), dtype=complex)

    # Onsite exchange: M_s = J m_s · sigma
    for s, m in enumerate(moments):
        onsite = jex * (
            m[0] * sigma_x
            + m[1] * sigma_y
            + m[2] * sigma_z
        )

        sl = slice(2 * s, 2 * s + 2)
        H[sl, sl] += onsite

    # Intra-cell hopping: A-B, B-C, C-D
    for s in range(n_sub - 1):
        sl = slice(2 * s, 2 * s + 2)
        sr = slice(2 * (s + 1), 2 * (s + 1) + 2)

        H[sl, sr] += t_hop * sigma_0
        H[sr, sl] += t_hop * sigma_0

    # Inter-cell hopping: D_R -> A_{R+1}
    sl_first = slice(0, 2)
    sl_last = slice(2 * (n_sub - 1), 2 * n_sub)

    H[sl_last, sl_first] += t_hop * np.exp(1j * K) * sigma_0
    H[sl_first, sl_last] += t_hop * np.exp(-1j * K) * sigma_0

    if not np.allclose(H, H.conj().T, atol=1.0e-12):
        raise RuntimeError("Bloch Hamiltonian is not Hermitian.")

    return H


def spin_operator(component):
    """
    Full spin operator I_4 tensor sigma_component.
    """
    n_sub = 4
    return np.kron(np.eye(n_sub, dtype=complex), SIGMA[component])


# ============================================================
# 4. Calculation
# ============================================================

def calculate_bands_and_spin():
    """
    Calculate eigenvalues and spin expectation values.
    """
    K_grid = np.linspace(-np.pi, np.pi, NK)

    n_band = 8
    energies = np.zeros((NK, n_band), dtype=float)

    sx = np.zeros((NK, n_band), dtype=float)
    sy = np.zeros((NK, n_band), dtype=float)
    sz = np.zeros((NK, n_band), dtype=float)

    Sx = spin_operator("x")
    Sy = spin_operator("y")
    Sz = spin_operator("z")

    for ik, K in enumerate(K_grid):
        H = bloch_hamiltonian(K, T_HOP, JEX)
        evals, evecs = np.linalg.eigh(H)

        energies[ik, :] = evals

        for ib in range(n_band):
            psi = evecs[:, ib]

            sx[ik, ib] = np.real(np.vdot(psi, Sx @ psi))
            sy[ik, ib] = np.real(np.vdot(psi, Sy @ psi))
            sz[ik, ib] = np.real(np.vdot(psi, Sz @ psi))

    data = {
        "K_grid": K_grid,
        "energies": energies,
        "sx": sx,
        "sy": sy,
        "sz": sz,
        "T_HOP": np.array([T_HOP]),
        "JEX": np.array([JEX]),
        "MU": np.array([MU]),
    }

    return data


# ============================================================
# 5. Data I/O
# ============================================================

def ensure_dirs():
    os.makedirs(DATA_DIR, exist_ok=True)
    os.makedirs(FIG_DIR, exist_ok=True)


def save_data(data):
    ensure_dirs()
    fname = os.path.join(DATA_DIR, "bands_and_spin.npz")
    np.savez(fname, **data)
    print(f"Saved data: {fname}")


def load_data():
    fname = os.path.join(DATA_DIR, "bands_and_spin.npz")

    if not os.path.exists(fname):
        raise FileNotFoundError(
            f"Data file not found: {fname}\n"
            f"Please set RUN_CALCULATION = True first."
        )

    raw = np.load(fname)
    data = {key: raw[key] for key in raw.files}
    return data


# ============================================================
# 6. Plotting style
# ============================================================

def set_plot_style():
    plt.rcParams.update({
        "font.family": "Times New Roman",
        "mathtext.fontset": "stix",
        "font.size": 14,
        "axes.linewidth": 1.1,
        "xtick.direction": "in",
        "ytick.direction": "in",
        "xtick.top": True,
        "ytick.right": True,
        "xtick.major.size": 4.5,
        "ytick.major.size": 4.5,
        "xtick.major.width": 1.0,
        "ytick.major.width": 1.0,
        "legend.frameon": False,
    })


def save_or_show(fig, filename):
    if SAVE_FIGURES:
        ensure_dirs()
        path = os.path.join(FIG_DIR, filename)
        fig.savefig(path, dpi=DPI, bbox_inches="tight")
        print(f"Saved figure: {path}")

    if SHOW_FIGURES:
        plt.show()
    else:
        plt.close(fig)


def colored_band_line(ax, x, y, c, lw=1.2):
    """
    Plot one band as a spin-colored line.
    """
    points = np.array([x, y]).T.reshape(-1, 1, 2)
    segments = np.concatenate([points[:-1], points[1:]], axis=1)

    cseg = 0.5 * (c[:-1] + c[1:])

    lc = LineCollection(
        segments,
        array=cseg,
        cmap=SPIN_CMAP,
        norm=plt.Normalize(vmin=SPIN_VMIN, vmax=SPIN_VMAX),
        linewidth=lw,
    )

    ax.add_collection(lc)
    return lc


# ============================================================
# 7. Plot functions
# ============================================================

def setup_band_axis(ax):
    ax.set_xlim(-np.pi, np.pi)
    ax.set_ylim(*E_WINDOW)

    ax.set_xticks([-np.pi, 0.0, np.pi])
    ax.set_xticklabels([r"$-\pi$", r"$0$", r"$\pi$"])

    ax.set_yticks([-2.0, 0.0, 2.0])

    ax.set_xlabel(r"$K=4k_xa$")
    ax.set_ylabel(r"$E$")

    if DRAW_MU_LINE:
        ax.axhline(MU, color="black", lw=0.9, ls="--")

    # Helpful high-symmetry markers
    ax.axvline(0.0, color="gray", lw=0.7, ls=":", alpha=0.8)
    ax.axvline(np.pi, color="gray", lw=0.7, ls=":", alpha=0.8)
    ax.axvline(-np.pi, color="gray", lw=0.7, ls=":", alpha=0.8)


def plot_plain_band(data):
    K = data["K_grid"]
    E = data["energies"]

    fig, ax = plt.subplots(figsize=(4.6, 3.5))

    for ib in range(E.shape[1]):
        ax.plot(K, E[:, ib], color="black", lw=1.2)

    setup_band_axis(ax)
    ax.set_title(r"Band structure")

    fig.tight_layout()
    save_or_show(fig, "band_plain.png")


def plot_spin_colored_band(data, component):
    K = data["K_grid"]
    E = data["energies"]
    S = data[f"s{component}"]

    fig, ax = plt.subplots(figsize=(4.8, 3.6))

    last_lc = None

    for ib in range(E.shape[1]):
        last_lc = colored_band_line(
            ax,
            K,
            E[:, ib],
            S[:, ib],
            lw=1.2,
        )

    setup_band_axis(ax)
    ax.set_title(rf"Band structure colored by $s_{component}$")

    cbar = fig.colorbar(last_lc, ax=ax, pad=0.025, shrink=0.88)
    cbar.set_ticks([-1, 0, 1])
    cbar.set_label(rf"$s_{component}$", rotation=90)

    fig.tight_layout()
    save_or_show(fig, f"band_spin_{component}.png")


def plot_spin_components_near_mu(data):
    """
    Plot spin expectation values of the band closest to MU.
    Useful for checking odd-parity spin texture.
    """
    K = data["K_grid"]
    E = data["energies"]

    sx = data["sx"]
    sy = data["sy"]
    sz = data["sz"]

    closest = np.argmin(np.abs(E - MU), axis=1)

    sx_mu = np.array([sx[i, closest[i]] for i in range(len(K))])
    sy_mu = np.array([sy[i, closest[i]] for i in range(len(K))])
    sz_mu = np.array([sz[i, closest[i]] for i in range(len(K))])

    fig, ax = plt.subplots(figsize=(5.2, 3.5))

    ax.plot(K, sx_mu, lw=1.5, label=r"$s_x$")
    ax.plot(K, sy_mu, lw=1.5, label=r"$s_y$")
    ax.plot(K, sz_mu, lw=1.5, label=r"$s_z$")

    ax.axhline(0.0, color="black", lw=0.8)
    ax.axvline(0.0, color="gray", lw=0.7, ls=":", alpha=0.8)

    ax.set_xlim(-np.pi, np.pi)
    ax.set_ylim(-1.05, 1.05)

    ax.set_xticks([-np.pi, 0.0, np.pi])
    ax.set_xticklabels([r"$-\pi$", r"$0$", r"$\pi$"])

    ax.set_yticks([-1.0, 0.0, 1.0])

    ax.set_xlabel(r"$K=4k_xa$")
    ax.set_ylabel(r"$s_\alpha$")
    ax.set_title(r"Spin expectation of the band closest to $\mu$")

    ax.legend(loc="best")

    fig.tight_layout()
    save_or_show(fig, "spin_texture_closest_to_mu.png")


# ============================================================
# 8. Main
# ============================================================

def main():
    if RUN_CALCULATION:
        print("====================================================")
        print("Calculating odd-parity magnets toy-model bands")
        print("====================================================")
        print(f"T_HOP = {T_HOP}")
        print(f"JEX   = {JEX}")
        print(f"MU    = {MU}")
        print(f"NK    = {NK}")
        print("====================================================")

        data = calculate_bands_and_spin()
        save_data(data)

    if RUN_PLOTTING:
        set_plot_style()
        data = load_data()

        plot_plain_band(data)

        for comp in ["x", "y", "z"]:
            plot_spin_colored_band(data, comp)

        plot_spin_components_near_mu(data)


if __name__ == "__main__":
    main()

参考文献

1. Hidden Zeeman Field in Odd-Parity Magnets: An Ideal Platform for Topological Superconductivity

鉴于该网站分享的大都是学习笔记，作者水平有限，若发现有问题可以发邮件给我

• yxliphy@gmail.com

也非常欢迎喜欢分享的小伙伴投稿