Hartree–Fock方法在 Fermi–Hubbard 模型中实现交错磁性

Fermi–Hubbard 模型与交替对角跃迁

考虑二维方格上的单轨道Fermi–Hubbard模型

$$H= -\sum_{ij,s}t_{ij} \left( c_{is}^{\dagger}c_{js}+\mathrm{H.c.} \right) + U\sum_i n_{i\uparrow}n_{i\downarrow} - \mu_{\mathrm{phys}}\sum_{i,s}n_{is},$$

其中 $s=\uparrow,\downarrow$ 表示真实自旋，$U>0$ 为在位排斥相互作用。晶格可划分为 $A$、$B$ 两个子晶格。最近邻跃迁连接不同子晶格，跃迁振幅为 $t$。沿两条对角方向

$$\boldsymbol{\delta}_1=(1,1), \qquad \boldsymbol{\delta}_2=(1,-1)$$

引入两个不同的次近邻跃迁

$$t_+=t'(1+\delta), \qquad t_-=t'(1-\delta).$$

两个子晶格上的强、弱对角跃迁方向相互交换，可以表示为

$$\begin{array}{c|cc} & (1,1) & (1,-1)\\ \hline A & t_- & t_+\\ B & t_+ & t_- \end{array}.$$

因此，单个子晶格的局域动能环境具有方向依赖，但 $A$、$B$ 两个子晶格可以通过四重旋转相互映射。这个交替的对角跃迁结构是后续形成交错磁能带的晶格基础。选择双子晶格原胞基矢

$$\mathbf a_1=(1,1), \qquad \mathbf a_2=(1,-1),$$

相应的倒格矢可取为

$$\mathbf b_1=(\pi,\pi), \qquad \mathbf b_2=(\pi,-\pi).$$

对于固定自旋 $s$，在子晶格基底

$$\psi_{\mathbf ks} = \begin{pmatrix} c_{\mathbf kAs}\\ c_{\mathbf kBs} \end{pmatrix}$$

中，非相互作用哈密顿量为

$$h_0(\mathbf k) = \begin{pmatrix} \varepsilon_A(\mathbf k) & \gamma(\mathbf k)\\ \gamma(\mathbf k) & \varepsilon_B(\mathbf k) \end{pmatrix},$$

其中

$$\varepsilon_A(\mathbf k) = -2t_-\cos(k_x+k_y) -2t_+\cos(k_x-k_y),\qquad \varepsilon_B(\mathbf k) = -2t_+\cos(k_x+k_y) -2t_-\cos(k_x-k_y),$$

而最近邻跃迁给出$\gamma(\mathbf k)
=
-2t\left(\cos k_x+\cos k_y\right)$。引入子晶格 Pauli 矩阵 $\tau_i$，定义

$$\varepsilon_0(\mathbf k) = \frac{\varepsilon_A(\mathbf k)+\varepsilon_B(\mathbf k)}{2},\quad g(\mathbf k) = \frac{\varepsilon_A(\mathbf k)-\varepsilon_B(\mathbf k)}{2},$$

则

$$h_0(\mathbf k) = \varepsilon_0(\mathbf k)\tau_0 + \gamma(\mathbf k)\tau_x + g(\mathbf k)\tau_z.$$

利用$t_++t_-=2t’,
t_+-t_-=2\delta t’$，可以得到

$$\varepsilon_0(\mathbf k) = -4t'\cos k_x\cos k_y,\quad g(\mathbf k) = -4\delta t'\sin k_x\sin k_y.$$

因此非相互作用哈密顿量可写成

$$h_0(\mathbf k) = -4t'\cos k_x\cos k_y\,\tau_0 -2t(\cos k_x+\cos k_y)\tau_x -4\delta t'\sin k_x\sin k_y\,\tau_z.$$

该哈密顿量中没有自旋矩阵，因此正常态仍然自旋简并

$$E_{\uparrow}(\mathbf k) = E_{\downarrow}(\mathbf k).$$

交替对角跃迁本身只区分两个子晶格的动量结构，并不会直接产生磁性。

Hubbard 相互作用的 Hartree–Fock 解耦

考虑共线的两子晶格磁序。定义

$$\eta_A=+1, \qquad \eta_B=-1,$$

并将自旋记为

$$s= \begin{cases} +1,&\uparrow,\\ -1,&\downarrow. \end{cases}$$

在不考虑电荷密度不均匀的情况下，可以将平均占据写成

$$\langle n_{\alpha s}\rangle = \frac{n}{2} + \eta_\alpha s\,m,$$

其中 $n$ 是每个格点的平均粒子数，$m$ 是两子晶格交错磁序参量。半填充时 $n=1$，于是

$$\begin{aligned} &\langle n_{A\uparrow}\rangle=\frac12+m, \qquad \langle n_{A\downarrow}\rangle=\frac12-m,\\ &\langle n_{B\uparrow}\rangle=\frac12-m, \qquad \langle n_{B\downarrow}\rangle=\frac12+m. \end{aligned}$$

磁序参量可以写为

$$m= \frac14 \left[ \langle n_{A\uparrow}\rangle - \langle n_{A\downarrow}\rangle - \langle n_{B\uparrow}\rangle + \langle n_{B\downarrow}\rangle \right].$$

两个子晶格上的磁矩大小相同、方向相反

$$\langle S_A^z\rangle=m, \qquad \langle S_B^z\rangle=-m,$$

因此总磁化为零：

$$M_{\mathrm{tot}} = \langle S_A^z\rangle+\langle S_B^z\rangle =0.$$

对 Hubbard 相互作用采用密度通道的 Hartree–Fock 解耦，

$$n_{i\uparrow}n_{i\downarrow} \approx n_{i\uparrow}\langle n_{i\downarrow}\rangle + \langle n_{i\uparrow}\rangle n_{i\downarrow} - \langle n_{i\uparrow}\rangle \langle n_{i\downarrow}\rangle.$$

对于子晶格 $\alpha$ 上的自旋 $s$ 粒子，

$$\langle n_{\alpha,-s}\rangle = \frac{n}{2}-\eta_\alpha s\,m,$$

因此其平均场势为

$$U\langle n_{\alpha,-s}\rangle = \frac{Un}{2} - \eta_\alpha sUm.$$

第一项与自旋和子晶格均无关，可以吸收到化学势中。定义移位后的化学势

$$\mu = \mu_{\mathrm{phys}} - \frac{Un}{2}.$$

剩余的磁性平均场为$-\eta_\alpha sUm$。在子晶格与自旋空间中，它可以写成$H_m=-Um\,\tau_zs_z$。因此，Hartree–Fock 解耦产生的是一个补偿型的子晶格交错交换场。

Hartree–Fock 哈密顿量与能谱

对于固定自旋 $s$，平均场哈密顿量为

$$h_s(\mathbf k) = \begin{pmatrix} \varepsilon_A(\mathbf k)-sUm & \gamma(\mathbf k) \\ \gamma(\mathbf k) & \varepsilon_B(\mathbf k)+sUm \end{pmatrix} -\mu\tau_0.$$

等价地，

$$h_s(\mathbf k) = \left[ \varepsilon_0(\mathbf k)-\mu \right]\tau_0 + \gamma(\mathbf k)\tau_x + \left[ g(\mathbf k)-sUm \right]\tau_z.$$

采用完整基底

$$\Psi_{\mathbf k} = \begin{pmatrix} c_{\mathbf kA\uparrow}\\ c_{\mathbf kB\uparrow}\\ c_{\mathbf kA\downarrow}\\ c_{\mathbf kB\downarrow} \end{pmatrix},$$

完整平均场哈密顿量为

$$H_{\mathrm{HF}}(\mathbf k) = \begin{pmatrix} h_{\uparrow}(\mathbf k)&0\\ 0&h_{\downarrow}(\mathbf k) \end{pmatrix}.$$

由于模型没有自旋轨道耦合，也没有横向磁序$[H_{\mathrm{HF}},s_z]=0$，因此上下自旋不发生混合，可以分别对角化两个 $2\times2$ 自旋块。定义

$$R_s(\mathbf k) = \sqrt{ \gamma^2(\mathbf k) + \left[ g(\mathbf k)-sUm \right]^2 }.$$

每个自旋块的两条能带为

$$E_{s,\nu}(\mathbf k) = \varepsilon_0(\mathbf k)-\mu + \nu R_s(\mathbf k), \qquad \nu=\pm1.$$

根号中的磁性项可以展开为

$$\left[g(\mathbf k)-sUm\right]^2 = g^2(\mathbf k) + U^2m^2 - 2sUm\,g(\mathbf k).$$

其中$-2sUm\,g(\mathbf k)$是使上下自旋能谱发生差异的关键交叉项。如果 $m=0$，则$R_\uparrow(\mathbf k)=R_\downarrow(\mathbf k)$，能带仍然自旋简并。另一方面，如果 $\delta=0$ 或 $t’=0$，则$g(\mathbf k)=0$，即使 $m\neq0$，仍有

$$E_{\uparrow,\nu}(\mathbf k) = E_{\downarrow,\nu}(\mathbf k).$$

因此，单独存在 Néel 序并不足以产生交错磁自旋劈裂。必须同时具有$m\neq0$和$g(\mathbf k)\neq0$。由于$g(\mathbf k)\propto\delta t’$，自旋劈裂的整体尺度由$Um\,\delta t’$共同控制。

交错磁自旋劈裂的动量结构

定义同一能带的自旋劈裂

$$\Delta E_\nu(\mathbf k) = E_{\uparrow,\nu}(\mathbf k) - E_{\downarrow,\nu}(\mathbf k).$$

其精确表达式为

$$\Delta E_\nu(\mathbf k) = \nu \left[ \sqrt{ \gamma^2(\mathbf k)+[g(\mathbf k)-Um]^2 } - \sqrt{ \gamma^2(\mathbf k)+[g(\mathbf k)+Um]^2 } \right].$$

在 $Um$ 较小的情况下展开，得到

$$\Delta E_\nu(\mathbf k) \approx -2\nu Um \frac{ g(\mathbf k) }{ \sqrt{ \gamma^2(\mathbf k)+g^2(\mathbf k) } }.$$

因此，自旋劈裂的动量依赖由 $g(\mathbf k)$ 决定。由于

$$g(\mathbf k) = -4\delta t'\sin k_x\sin k_y,$$

自旋劈裂在动量空间中具有 $d$ 波型的符号变化，而不是铁磁体中近似均匀的 Zeeman 劈裂。在原始方格坐标中，

$$g(\mathbf k)\sim\sin k_x\sin k_y$$

具有 $d_{xy}$ 型形式。若引入旋转坐标

$$q_1=k_x+k_y, \qquad q_2=k_x-k_y,$$

则

$$g(\mathbf k) \sim \cos q_1-\cos q_2,$$

对应 $d_{x^2-y^2}$ 型形式。这两种表述只是坐标选择不同。当$g(\mathbf k)=0$时，自旋劈裂消失。对于当前模型，节点位于

$$k_x=0,\ \pi, \qquad k_y=0,\ \pi.$$

因此，交错磁自旋劈裂在布里渊区中必然发生符号改变，并具有对称性保护的节点。

交错磁对称性

在四重旋转$C_{4z}:
(k_x,k_y)
\longrightarrow
(-k_y,k_x)$下，

$$\varepsilon_0(C_{4z}\mathbf k) = \varepsilon_0(\mathbf k),\qquad \gamma(C_{4z}\mathbf k) = \gamma(\mathbf k),$$

而$g(C_{4z}\mathbf k)
=
-g(\mathbf k)$。同时进行自旋翻转$s\rightarrow-s$，有

$$g(C_{4z}\mathbf k)-(-s)Um = -\left[ g(\mathbf k)-sUm \right].$$

由于能谱只依赖该量的平方，因此

$$R_{-s}(C_{4z}\mathbf k) = R_s(\mathbf k),$$

从而得到

$$E_{s,\nu}(\mathbf k) = E_{-s,\nu}(C_{4z}\mathbf k).$$

所以，上下自旋能带并不是在同一个动量处简并，而是通过四重旋转相互对应：

$$E_{\uparrow,\nu}(\mathbf k) = E_{\downarrow,\nu}(C_{4z}\mathbf k).$$

在一般动量处则有

$$E_{\uparrow,\nu}(\mathbf k) \neq E_{\downarrow,\nu}(\mathbf k).$$

因此，这一磁态同时具有

$$M_{\mathrm{tot}}=0$$

和动量依赖的自旋劈裂。前者来自两个子晶格磁矩的完全补偿，后者来自子晶格依赖的动能项 $g(\mathbf k)$ 与 Néel 交换场 $Um$ 的共同作用。普通双子晶格反铁磁体中，即使 $m\neq0$，如果两个子晶格的动能环境等价，即 $g(\mathbf k)=0$，能带仍然自旋简并。当前模型中，两个反向自旋子晶格由空间旋转而不是简单平移关联，因此允许零净磁化与自旋分裂能带同时存在。

Hartree–Fock 自洽方程

设平均场本征方程为

$$H_{\mathrm{HF}}(\mathbf k) |u_{n\mathbf k}\rangle = E_{n\mathbf k} |u_{n\mathbf k}\rangle.$$

交错磁序对应的算符为

$$\Gamma_m=\tau_zs_z.$$

因此磁序参量满足

$$m = \frac{1}{4N_k} \sum_{\mathbf k,n} f(E_{n\mathbf k}) \left\langle u_{n\mathbf k} \left| \tau_zs_z \right| u_{n\mathbf k} \right\rangle,$$

其中

$$f(E) = \frac{1}{e^{E/T}+1}$$

为 Fermi 分布。粒子数方程为

$$n_{\mathrm{cell}} = \frac{1}{N_k} \sum_{\mathbf k,n} f(E_{n\mathbf k}).$$

在半填充下，一个两格点原胞包含两个粒子，因此要求$n_{\mathrm{cell}}=2$。于是 $m$ 与 $\mu$ 需要同时满足

$$m=m_{\mathrm{out}}(m,\mu),\qquad n_{\mathrm{cell}}(m,\mu)=2.$$

对于当前的 $2\times2$ 自旋块，可以进一步写出解析形式。对于自旋 $s$、能带指标 $\nu$ 的本征态，

$$\left\langle \tau_z \right\rangle_{s,\nu} = \nu \frac{ g(\mathbf k)-sUm }{ R_s(\mathbf k) }.$$

因此

$$\left\langle \tau_zs_z \right\rangle_{s,\nu} = s\nu \frac{ g(\mathbf k)-sUm }{ R_s(\mathbf k) }.$$

磁序自洽方程为

$$m = \frac{1}{4N_k} \sum_{\mathbf k,s=\pm1,\nu=\pm1} s\nu \frac{ g(\mathbf k)-sUm }{ R_s(\mathbf k) } f\!\left[ E_{s,\nu}(\mathbf k) \right],$$

粒子数方程为

$$2 = \frac{1}{N_k} \sum_{\mathbf k,s,\nu} f\!\left[ E_{s,\nu}(\mathbf k) \right].$$

求解这两个方程便可以得到给定 $U$、$t’$、$\delta$ 和温度下的自洽磁序 $m$ 与化学势 $\mu$。

整体物理图像

该模型的非相互作用哈密顿量为

$$H_0(\mathbf k) = \varepsilon_0(\mathbf k)\tau_0 + \gamma(\mathbf k)\tau_x + g(\mathbf k)\tau_z,$$

其中$g(C_{4z}\mathbf k)
=
-g(\mathbf k)$。Hubbard相互作用在Hartree–Fock近似下产生补偿型Néel交换场

$$H_m=-Um\,\tau_zs_z.$$

二者结合后得到

$$H_{\mathrm{HF}}(\mathbf k) = \left[ \varepsilon_0(\mathbf k)-\mu \right]\tau_0 + \gamma(\mathbf k)\tau_x + \left[ g(\mathbf k)-Um\,s_z \right]\tau_z.$$

Hartree–Fock 自洽决定磁序是否形成以及磁序大小 $m$；交替各向异性的对角跃迁则决定磁序形成后能带的动量结构。能谱中的

$$-2sUm\,g(\mathbf k)$$

将补偿型 Néel 序转化为 $d$ 波型的动量依赖自旋劈裂。因此，该模型中交错磁性的产生机制可以概括为

$${\color{blue} \text{Hubbard 相互作用产生 Neel 序，子晶格依赖的对角跃迁将其转化为交错磁能带。} }\nonumber$$

复现代码

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

"""
Reproduce Fig. 1(b) of the PRL altermagnetic Hubbard model paper.

Only reproduces the spin-resolved Hartree-Fock band structure and Fermi-surface inset.

Model parameters:
    t'/t = 0.3
    delta = 0.9
    U/t = 3.5
    half filling

Magnetic unit cell:
    A/B two-sublattice cell
    a1 = (1, 1)
    a2 = (1,-1)

Momentum-space Hartree-Fock Hamiltonian for each spin:
    H_s(k) =
        [ h_AA,s(k)   h_AB(k)  ]
        [ h_AB(k)     h_BB,s(k)]

where
    h_AA,s = -2 t_- cos(kx+ky) - 2 t_+ cos(kx-ky) - s U dm
    h_BB,s = -2 t_+ cos(kx+ky) - 2 t_- cos(kx-ky) + s U dm
    h_AB   = -2 t [cos(kx) + cos(ky)]

s = +1 for spin up, s = -1 for spin down.
"""

import os
import numpy as np
import matplotlib.pyplot as plt
from dataclasses import dataclass
from matplotlib.patches import Polygon
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.ticker import MaxNLocator


# ============================================================
# 1. Parameters
# ============================================================

@dataclass
class Params:
    # Model parameters
    t: float = 1.0
    tp: float = 0.30
    delta: float = 0.90
    U: float = 3.50

    # Half filling: 1 particle per original site, 2 particles per A/B cell
    filling: float = 1.0

    # Small finite temperature for stable self-consistency
    T: float = 0.005

    # Self-consistency
    Nk_self: int = 201
    max_iter: int = 800
    tol: float = 1.0e-11
    mix: float = 0.25
    dm_init: float = 0.25

    # Chemical-potential bisection
    mu_min: float = -8.0
    mu_max: float = 8.0
    mu_iter: int = 100

    # Band and FS grids
    nk_path_per_segment: int = 260
    nk_fs: int = 501

    # Output
    out_dir: str = "fig"
    dpi: int = 300

    # Use LaTeX. Set False if your Python environment has no LaTeX.
    use_latex: bool = True


P = Params()


# ============================================================
# 2. Plot style
# ============================================================

def setup_plot_style(p: Params):
    plt.rcParams["figure.facecolor"] = "white"
    plt.rcParams["axes.facecolor"] = "white"
    plt.rcParams["savefig.facecolor"] = "white"

    plt.rcParams["font.family"] = "serif"
    plt.rcParams["font.serif"] = ["Times New Roman", "Times", "DejaVu Serif"]
    plt.rcParams["mathtext.fontset"] = "stix"

    plt.rcParams["text.usetex"] = p.use_latex

    plt.rcParams["axes.linewidth"] = 1.0
    plt.rcParams["xtick.direction"] = "in"
    plt.rcParams["ytick.direction"] = "in"
    plt.rcParams["xtick.top"] = True
    plt.rcParams["ytick.right"] = True

    plt.rcParams["font.size"] = 13
    plt.rcParams["axes.labelsize"] = 15
    plt.rcParams["legend.fontsize"] = 11
    plt.rcParams["xtick.labelsize"] = 11
    plt.rcParams["ytick.labelsize"] = 11


# ============================================================
# 3. Basic functions
# ============================================================

def fermi(x, T):
    y = x / T
    out = np.empty_like(y, dtype=float)

    out[y > 40.0] = 0.0
    out[y < -40.0] = 1.0

    mask = (y <= 40.0) & (y >= -40.0)
    out[mask] = 1.0 / (np.exp(y[mask]) + 1.0)

    return out


def tplus_tminus(p: Params):
    t_plus = p.tp * (1.0 + p.delta)
    t_minus = p.tp * (1.0 - p.delta)
    return t_plus, t_minus


def make_magnetic_bz_mesh(p: Params):
    """
    Magnetic real-space unit cell:
        a1 = (1, 1)
        a2 = (1,-1)

    Reciprocal vectors:
        b1 = (pi, pi)
        b2 = (pi,-pi)

    k = u b1 + v b2, u,v in [-1/2, 1/2).
    """
    nk = p.Nk_self
    offset = 0.5

    u = (np.arange(nk) + offset) / nk - 0.5
    v = (np.arange(nk) + offset) / nk - 0.5

    U, V = np.meshgrid(u, v, indexing="ij")

    kx = np.pi * (U + V)
    ky = np.pi * (U - V)

    return kx.ravel(), ky.ravel()


def spin_block_bands_and_weights(kx, ky, spin, dm, p: Params):
    """
    Return two eigenvalues and A-sublattice weights for one spin block.

    spin = +1: up
    spin = -1: down
    """
    t_plus, t_minus = tplus_tminus(p)

    c1 = np.cos(kx + ky)
    c2 = np.cos(kx - ky)

    h_AA = -2.0 * t_minus * c1 - 2.0 * t_plus * c2 - spin * p.U * dm
    h_BB = -2.0 * t_plus * c1 - 2.0 * t_minus * c2 + spin * p.U * dm
    h_AB = -2.0 * p.t * (np.cos(kx) + np.cos(ky))

    d0 = 0.5 * (h_AA + h_BB)
    dz = 0.5 * (h_AA - h_BB)
    dx = h_AB

    r = np.sqrt(dx * dx + dz * dz + 1.0e-30)

    e_lower = d0 - r
    e_upper = d0 + r

    # For H = d0 I + dx tau_x + dz tau_z:
    # lower-band A weight = (1 - dz/r)/2
    # upper-band A weight = (1 + dz/r)/2
    wA_lower = 0.5 * (1.0 - dz / r)
    wA_upper = 0.5 * (1.0 + dz / r)

    return e_lower, e_upper, wA_lower, wA_upper


# ============================================================
# 4. Hartree-Fock self-consistency
# ============================================================

def density_for_mu(mu, dm, kx, ky, p: Params):
    """
    Total density per magnetic cell.
    There are two bands for each spin, so half filling means n_cell = 2.
    """
    n_total = 0.0

    for spin in (+1, -1):
        e1, e2, _, _ = spin_block_bands_and_weights(kx, ky, spin, dm, p)
        n_total += np.sum(fermi(e1 - mu, p.T))
        n_total += np.sum(fermi(e2 - mu, p.T))

    return n_total / len(kx)


def find_mu(dm, kx, ky, p: Params):
    target = 2.0 * p.filling

    lo = p.mu_min
    hi = p.mu_max

    for _ in range(p.mu_iter):
        mid = 0.5 * (lo + hi)
        n_mid = density_for_mu(mid, dm, kx, ky, p)

        if n_mid < target:
            lo = mid
        else:
            hi = mid

    return 0.5 * (lo + hi)


def compute_order_parameter(dm, mu, kx, ky, p: Params):
    """
    dm = 1/4 * (n_A_up - n_A_down - n_B_up + n_B_down)

    This convention gives:
        A: spin-up favored
        B: spin-down favored
    for dm > 0.
    """
    densities = {}

    for spin_label, spin in [("up", +1), ("down", -1)]:
        e1, e2, wA1, wA2 = spin_block_bands_and_weights(kx, ky, spin, dm, p)

        occ1 = fermi(e1 - mu, p.T)
        occ2 = fermi(e2 - mu, p.T)

        nA = np.sum(occ1 * wA1 + occ2 * wA2) / len(kx)
        nB = np.sum(occ1 * (1.0 - wA1) + occ2 * (1.0 - wA2)) / len(kx)

        densities[(spin_label, "A")] = nA
        densities[(spin_label, "B")] = nB

    n_A_up = densities[("up", "A")]
    n_B_up = densities[("up", "B")]
    n_A_dn = densities[("down", "A")]
    n_B_dn = densities[("down", "B")]

    dm_new = 0.25 * (n_A_up - n_A_dn - n_B_up + n_B_dn)

    return dm_new, densities


def solve_hartree_fock(p: Params):
    kx, ky = make_magnetic_bz_mesh(p)

    dm = p.dm_init
    mu = 0.0

    history = []

    for it in range(p.max_iter):
        mu = find_mu(dm, kx, ky, p)
        dm_raw, densities = compute_order_parameter(dm, mu, kx, ky, p)

        dm_next = (1.0 - p.mix) * dm + p.mix * dm_raw
        err = abs(dm_next - dm)

        history.append((it, dm, dm_raw, mu, err))

        dm = dm_next

        if err < p.tol:
            break

    mu = find_mu(dm, kx, ky, p)
    dm_raw, densities = compute_order_parameter(dm, mu, kx, ky, p)

    result = {
        "delta_m": dm,
        "mu": mu,
        "densities": densities,
        "history": np.array(history),
        "density_cell": density_for_mu(mu, dm, kx, ky, p),
    }

    return result


# ============================================================
# 5. Band path and FS data
# ============================================================

def make_k_path(p: Params):
    """
    Path in the magnetic Brillouin zone:
        Gamma -> (-pi/2, pi/2) -> (0, pi) -> (pi/2, pi/2) -> Gamma
    """
    nodes = [
        np.array([0.0, 0.0]),
        np.array([-0.5 * np.pi, 0.5 * np.pi]),
        np.array([0.0, np.pi]),
        np.array([0.5 * np.pi, 0.5 * np.pi]),
        np.array([0.0, 0.0]),
    ]

    labels = [
        r"$\Gamma$",
        r"$X_1$",
        r"$M$",
        r"$X_2$",
        r"$\Gamma$",
    ]

    kx_list = []
    ky_list = []
    kdist = []
    node_dist = [0.0]

    current = 0.0

    for i in range(len(nodes) - 1):
        k0 = nodes[i]
        k1 = nodes[i + 1]

        for j in range(p.nk_path_per_segment):
            frac = j / p.nk_path_per_segment
            k = k0 + frac * (k1 - k0)

            if len(kx_list) > 0:
                dk = np.linalg.norm(k - np.array([kx_list[-1], ky_list[-1]]))
                current += dk

            kx_list.append(k[0])
            ky_list.append(k[1])
            kdist.append(current)

        node_dist.append(current + np.linalg.norm(k1 - np.array([kx_list[-1], ky_list[-1]])))

    kx_list.append(nodes[-1][0])
    ky_list.append(nodes[-1][1])

    current += np.linalg.norm(
        np.array([kx_list[-1], ky_list[-1]]) - np.array([kx_list[-2], ky_list[-2]])
    )
    kdist.append(current)
    node_dist[-1] = current

    return np.array(kx_list), np.array(ky_list), np.array(kdist), np.array(node_dist), labels


def compute_band_path(dm, mu, p: Params):
    kx, ky, kdist, node_dist, labels = make_k_path(p)

    eu1, eu2, _, _ = spin_block_bands_and_weights(kx, ky, +1, dm, p)
    ed1, ed2, _, _ = spin_block_bands_and_weights(kx, ky, -1, dm, p)

    data = np.column_stack([
        kdist,
        kx,
        ky,
        eu1 - mu,
        eu2 - mu,
        ed1 - mu,
        ed2 - mu,
    ])

    return data, node_dist, labels


def compute_fermi_surface_grid(dm, mu, p: Params):
    k = np.linspace(-np.pi, np.pi, p.nk_fs)
    kx, ky = np.meshgrid(k, k, indexing="xy")

    eu1, eu2, _, _ = spin_block_bands_and_weights(kx, ky, +1, dm, p)
    ed1, ed2, _, _ = spin_block_bands_and_weights(kx, ky, -1, dm, p)

    return kx, ky, eu1 - mu, eu2 - mu, ed1 - mu, ed2 - mu


# ============================================================
# 6. Plot Fig. 1(b)
# ============================================================

def draw_magnetic_bz(ax):
    mbz = np.array([
        [0.0, np.pi],
        [np.pi, 0.0],
        [0.0, -np.pi],
        [-np.pi, 0.0],
    ])

    polygon = Polygon(
        mbz,
        closed=True,
        facecolor="0.82",
        edgecolor="0.25",
        linewidth=1.0,
        alpha=0.65,
        zorder=0,
    )
    ax.add_patch(polygon)


def plot_fig1b(ax, band_data, node_dist, labels, dm, mu, p: Params):
    kdist = band_data[:, 0]

    e_up_1 = band_data[:, 3]
    e_up_2 = band_data[:, 4]
    e_dn_1 = band_data[:, 5]
    e_dn_2 = band_data[:, 6]

    color_up = "#D73027"
    color_dn = "#4575B4"

    ax.plot(kdist, e_up_1, color=color_up, lw=2.4, label=r"$\uparrow$")
    ax.plot(kdist, e_up_2, color=color_up, lw=2.4)

    ax.plot(kdist, e_dn_1, color=color_dn, lw=2.4, label=r"$\downarrow$")
    ax.plot(kdist, e_dn_2, color=color_dn, lw=2.4)

    ax.axhline(0.0, color="0.2", lw=1.0, ls=":")

    for nd in node_dist:
        ax.axvline(nd, color="0.2", lw=0.7, alpha=0.35)

    ax.set_xticks(node_dist)
    ax.set_xticklabels(labels)

    ax.set_xlim(kdist.min(), kdist.max())
    ax.set_ylim(-6.2, 2.4)

    ax.set_ylabel(r"$E/t$")
    ax.yaxis.set_major_locator(MaxNLocator(nbins=3))

    ax.legend(frameon=False, loc="upper center", ncol=2, handlelength=2.2)

    # Fermi-surface inset
    axins = inset_axes(ax, width="38%", height="43%", loc="lower center", borderpad=1.20)

    KX, KY, EU1, EU2, ED1, ED2 = compute_fermi_surface_grid(dm, mu, p)

    draw_magnetic_bz(axins)

    for E in [EU1, EU2]:
        axins.contour(
            KX,
            KY,
            E,
            levels=[0.0],
            colors=[color_up],
            linewidths=1.5,
        )

    for E in [ED1, ED2]:
        axins.contour(
            KX,
            KY,
            E,
            levels=[0.0],
            colors=[color_dn],
            linewidths=1.5,
        )

    # High-symmetry path in the inset
    path_kx, path_ky, _, _, _ = make_k_path(p)
    axins.plot(path_kx, path_ky, color="k", lw=0.75, ls="--", alpha=0.75)

    axins.set_xlim(-np.pi, np.pi)
    axins.set_ylim(-np.pi, np.pi)
    axins.set_aspect("equal", adjustable="box")

    axins.set_xticks([-np.pi, 0.0, np.pi])
    axins.set_yticks([-np.pi, 0.0, np.pi])

    axins.set_xticklabels([r"$-\pi$", r"$0$", r"$\pi$"], fontsize=8)
    axins.set_yticklabels([r"$-\pi$", r"$0$", r"$\pi$"], fontsize=8)

    axins.set_xlabel(r"$k_x a$", fontsize=9, labelpad=-2)
    axins.set_ylabel(r"$k_y a$", fontsize=9, labelpad=-5)

    ax.text(
        0.03,
        0.04,
        rf"$\delta m={dm:.3f}$" + "\n" + rf"$\mu/t={mu:.3f}$",
        transform=ax.transAxes,
        fontsize=11,
        ha="left",
        va="bottom",
    )


# ============================================================
# 7. Main
# ============================================================

def main():
    os.makedirs(P.out_dir, exist_ok=True)
    setup_plot_style(P)

    print("Solving Hartree-Fock self-consistency ...")
    result = solve_hartree_fock(P)

    dm = result["delta_m"]
    mu = result["mu"]

    print("Converged mean-field solution:")
    print(f"  delta_m      = {dm:.10f}")
    print(f"  U * delta_m  = {P.U * dm:.10f}")
    print(f"  mu / t       = {mu:.10f}")
    print(f"  n_cell       = {result['density_cell']:.10f}")
    print("  sublattice/spin densities:")
    for key, val in result["densities"].items():
        print(f"    {key}: {val:.10f}")

    band_data, node_dist, labels = compute_band_path(dm, mu, P)

    # Save numerical data
    np.savetxt(
        os.path.join(P.out_dir, "fig1b_band_data.dat"),
        band_data,
        header="kdist kx ky E_up_1 E_up_2 E_down_1 E_down_2",
    )

    np.savetxt(
        os.path.join(P.out_dir, "hf_convergence_history.dat"),
        result["history"],
        header="iter delta_m_input delta_m_raw mu error",
    )

    # Plot Fig. 1(b)
    fig, ax = plt.subplots(figsize=(7.0, 4.7))
    plot_fig1b(ax, band_data, node_dist, labels, dm, mu, P)

    fig_path = os.path.join(P.out_dir, "reproduce_prl_fig1b.png")
    fig.savefig(fig_path, dpi=P.dpi, bbox_inches="tight")
    print(f"Saved: {fig_path}")

    plt.close(fig)


if __name__ == "__main__":
    main()

参考文献

1. Realizing Altermagnetism in Fermi-Hubbard Models with Ultracold Atoms

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

• yxliphy@gmail.com

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