Time-independet Perturbation

现在有一个微扰后哈密顿量：

原始哈密顿量有以下能级：

能级有以下关系：

现在我们的问题是，如何求解微扰后的能级：

假定态矢量和能级可以用正交基底多项式展开：

那么只需要代入薛定谔方程即可。

这里的正交即认为各高级近似解$|n^i\rangle$与零级近似解$|n^0\rangle$正交。这并不是一个假设，而是一个约定，因为想要满足近似能量解$E_i$的态矢量有很多，我们需要一个约束来限制它们的自由度。

Non-degenerate Perturbation

对于不简并的情况：

展开有：

对不同的量级进行分开分析：

First Order Theory

我们先考虑一阶修正：

两边同时作用$\langle n^0|$：

现在我们完成了对能级的修正，但除了一阶，其他的能量修正因为不知道态函数修正而无法求解，这时候需要用一阶能量修正来求解一阶态矢量修正，并依次类推。一般选择本征态基底$|k^0\rangle$作展开：

所以一阶态矢量修正被原本的能级表示出来，系数即为上述表示：

Second Order Theory

同理，考虑k阶修正：

由于不同级的近似解正交，解得：

比如2阶：

由于二阶能量微扰的表达式存在平方项，显然，二阶微扰总是起到互斥的作用(Level Repulsion)：

这意味着，如果$E^0_k>E^0_n$，那么$E^2_n<0$，也就是说，k能级会使得n能级的能量降低；反之，如果$e^0_k0$，也就是说，k能级会使得n能级的能量升高。这就是能级互斥效应。对于基态能级，由于不存在比它更低的能级，所以二阶微扰总是使得基态能级的能量降低。

按照一阶微扰论的逻辑，接下来应该求解二阶态矢量修正$|n^2\rangle$，其实也不难求解：

代入$|n^1\rangle$和$E^1_n,E^2_n$的表达式，得到：

但实际上，如果只是要计算三阶能量修正，其实是不需要用到二阶态矢量修正的。代入二阶态矢量修正化简：

这说明计算三阶能量修正只需要用到一阶态矢量修正和一阶能量修正。

Degenerate Perturbation

考虑有简并的哈密顿量：

其中

如果我们按照非简并的方法求解，我们会发现，一阶修正能量为0：

这是因为简并态的一阶修正能量是相等的，所以无法求解。这时候我们需要用另一种方法求解。为了简单起见，我们先研究两重简并态的情况。

Two-fold Degenerate

原始哈密顿量有以下简并能级：

展开有：

其中：

我们可以得到：

以上的步骤相对于非简并的情况是一样的，不过：

对于一阶，同时作用$|n^0,a\rangle$：

显然由于哈密顿算符的简并性和态的正交性，有一些项为0：

同样的：

写成矩阵形式：

所以能量的一阶修正，就是该矩阵的本征值:

同时也可以求出本征矢量：

实际上不需要死记硬背，用矩阵的表述可以统一两种结果。对于非简并微扰，本质上是解关于修正能量和修正波函数的矩阵方程：

这实际上是$k+1$个方程确定$k+1$个未知数，其中$n$个方程为波函数各维度的线性方程，另外加上波函数的归一化方程。未知数分别是原本征能量表象下的$k$维一阶修正波函数$|n^1\rangle$外加未知的一阶能量修正$E^1_n$。简单计算就能知道这个结果和上面的狄拉克符号的推导是一致的。

对于简并微扰，本质上是解关于修正能量、修正波函数和“Good State”的矩阵方程。我们记具有$m_n$重简并的第$n$类的第$i$个“Good State”写为：

相比于非简并微扰，就变成解$k+1+\sum_n m_n$个未知数，其中$m_n$来源于“Good State”的构造。这多余的$\sum_n m_n$个未知数需要通过“Good State”的正交性和归一化来约束（这当然是合理的，只要$\hat H+\delta\hat H$仍然是一个厄密算符）。

Problem

  1. 对角化精确求解本征值和本征波矢；
  2. 非简并微扰求解到二阶能量修正；
  3. 简并微扰求解到二阶能量修正。

这题很好得考察了对简并微扰的理解，以下的两个特征使你不能直接套公式：

  1. 存在非简并能级与简并能级的混合；
  2. 一阶修正无法消除简并。

对角化

先对角化求精确的本征值：

解得：

非简并微扰

显然一阶能量修正为0，本征态修正为：

二阶能量修正为：

这里的第一个能级出现了问题，而第二个没有。

简并微扰

我们假设“Good State”为：

相应的一阶本征态修正为$(a^{(1)}_i,b^{(1)}_i,c^{(1)}_i)$，代入：

得到：

即

既然$\alpha,\beta$不同时为0，那么$E_1^{(1)}=0$。相应的，可以解出

一阶的态修正：

设二阶态修正为$(a^{(2)}_i,b^{(2)}_i,c^{(2)}_i)$，代入二阶方程：

即：

解得：

但这样算就不对了，因为这样两个”Good State”就一样，这是不对的，检查发现是因为存在第二个解。取：

这样就解出了正确的微扰解。对于第二个能级，和上述非简并微扰一样。（此处我没有使用正交归一化条件，但是通过方程结合物理意义，用少的方程解出了多的未知数）

Example

电磁场、相对论修正和自旋轨道耦合都可以看作是微扰。根据微扰的不同，可以分为以下效应：

  • 加入外电场：Stark效应；
  • 加入外磁场：Zeeman效应；
  • 相对论修正和自旋轨道耦合：精细结构；

Stark Effect

考虑将氢原子放入外电场，产生Stark效应：

  • 宇称选择定则：这是一个奇宇称的算符，所以只有不同宇称，或者说$l+l’$是奇数之间的态可能；
  • 角动量选择定则：$z=r\cos{\theta}=2\sqrt{\dfrac{\pi}{3}}rY_{1,0}(\theta,\phi)$，所以连接不同轨道角动量数l之间的态，且必须满足$m=m’$。

Quadratic Stark Effect

对于$n=1$的能级，显然根据选择定则，一阶能量修正为0。代入二阶修正：

所以外电场使能量降低了，这和上面提到的能级互斥是一致的。

Linear Stark Effect

考虑$n=2$的能级，此时一阶能量修正不一定为0，同时这是一个简并微扰问题，因为在不加外电场和不考虑自旋轨道耦合的情况下，$n=2$能级是四重简并的。根据选择定则写出微扰算符：

显然只连接了$|2,0,0\rangle,|2,1,0\rangle$两个态，满足$l\neq l’,m=m’$的条件。”Good State”为：

那么一阶能量修正为：

这意味着4重简并的能级最终分裂。

由于$n=2$对应的是一阶修正，所以是线性效应；而$n=1$对应的是二阶修正，所以是二次效应。以下是线性Stark效应中，$|n^0,a\rangle$对应的氢原子图。由于概率分布重心靠下，所以能量会降低。

Stark Effect

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import PowerNorm
from mpl_toolkits.mplot3d import Axes3D

# 定义常数
a0 = 1.0  # 玻尔半径（单位化）

# 定义波函数
def psi_200(x, y, z):
    r = np.sqrt(x**2 + y**2 + z**2)
    return (1 / (4 * np.sqrt(2 * np.pi) * a0**(3/2))) * (2 - r/a0) * np.exp(-r/(2*a0))

def psi_210(x, y, z):
    r = np.sqrt(x**2 + y**2 + z**2)
    theta = np.arccos(z/r)
    return (1 / (4 * np.sqrt(2 * np.pi) * a0**(3/2))) * (r/a0) * np.exp(-r/(2*a0)) * np.cos(theta)

# 定义Good State波函数
def psi_good_plus(x, y, z):
    return (psi_200(x, y, z) + psi_210(x, y, z)) / np.sqrt(2)

# 创建网格
R = 5
x = np.linspace(-R*a0, R*a0, 100)
z = np.linspace(-R*a0, R*a0, 100)
X, Z = np.meshgrid(x, z)
Y = 0  # 在y=0平面观察

# 计算概率密度
density_200 = np.abs(psi_200(X, Y, Z))**2
density_210 = np.abs(psi_210(X, Y, Z))**2
density_good_plus = np.abs(psi_good_plus(X, Y, Z))**2

# 找到所有概率密度的最大值和最小值用于归一化
vmin = min(np.min(density_200[density_200>0]), 
           np.min(density_210[density_210>0]), 
           np.min(density_good_plus[density_good_plus>0]))
vmax = max(np.max(density_200), np.max(density_210), np.max(density_good_plus))

# 创建图形
fig = plt.figure(figsize=(18, 6))
gs = fig.add_gridspec(1, 4, width_ratios=[1, 1, 1, 0.05])  # 为颜色条预留空间

# 创建子图
ax1 = fig.add_subplot(gs[0])
ax2 = fig.add_subplot(gs[1])
ax3 = fig.add_subplot(gs[2])
cax = fig.add_subplot(gs[3])  # 颜色条专用子图

# 使用PowerNorm进行归一化
norm = PowerNorm(gamma=0.5, vmin=vmin, vmax=vmax)

# 绘制 |2,0,0>
im1 = ax1.imshow(density_200, extent=[-R*a0, R*a0, -R*a0, R*a0], 
                origin='lower', cmap='viridis', norm=norm)
ax1.set_title(r'$|2,0,0\rangle$ Probability Density')
ax1.set_xlabel('x (a₀)')
ax1.set_ylabel('z (a₀)')

# 绘制 |2,1,0>
im2 = ax2.imshow(density_210, extent=[-R*a0, R*a0, -R*a0, R*a0], 
                origin='lower', cmap='viridis', norm=norm)
ax2.set_title(r'$|2,1,0\rangle$ Probability Density')
ax2.set_xlabel('x (a₀)')
ax2.set_ylabel('z (a₀)')

# 绘制 Good State (|a>)
im3 = ax3.imshow(density_good_plus, extent=[-R*a0, R*a0, -R*a0, R*a0], 
                origin='lower', cmap='viridis', norm=norm)
ax3.set_title(r'Good State $\frac{1}{\sqrt{2}}(|2,0,0\rangle + |2,1,0\rangle)$')
ax3.set_xlabel('x (a₀)')
ax3.set_ylabel('z (a₀)')

# 添加统一的颜色条
cbar = fig.colorbar(im3, cax=cax)
cbar.set_label('Probability Density')

plt.tight_layout()
plt.show()

Fine Structure

我们一般把氢原子的非相对论的哈密顿算符写为：

这样求解出来的氢原子能级为：

其中$\alpha=\frac{e^2}{\hbar c}$为精细结构常数。

若考虑相对论效应，则可以将相对论哈密顿算符展开为：

Relative Movement Correction

这部分是由于相对论效应导致的修正。本来由于能级简并需要用到简并微扰理论，但是由于：

这意味着$\hat{H}_{mv}$已经是对角化的，不会联系不同的轨道角动量数，所以我们可以直接求解：

根据薛定谔方程$p^2\psi=2m(E-V)\psi$，代入有：

代入：

得到：

Spin-Orbit Coupling

是自旋轨道耦合微扰。运动的电荷产生磁场，这本身也是一种相对论效应。这里选择LS耦合表象，即以$n,l,s,j,m_j$为一组好量子数，具体理由请见 原子的精细结构 。在该表象下：

合起来得到：

考虑到耦合后的$j=l+\frac12$和$j=l-\frac12$，可以写成：

将两个相对论修正项加在一起，会发现无论是$j=l+\frac12$还是$j=l-\frac12$，都可以写成一个统一的表达式：

Zeeman Effect

考虑将氢原子放入外磁场，产生Zeeman效应：

哈密顿量修正为：

此处有两项微扰——一阶塞曼效应和二阶塞曼效应。我们一般把一阶塞曼效应部分的哈密顿量写为：

当我们不考虑氢原子的精细结构（如 原子的精细结构 所说的双电子独态没有自旋轨道耦合），即为正常塞曼效应，产生三条分裂谱线；当考虑氢原子的精细结构时，即为反常塞曼效应，产生$2j+1$条分裂谱线。此时还需要考虑弱场、强场和中间场情况，分别对应$\hat{H}_Z$为微扰、$\hat{H}_{mv}+\hat{H}_{SOC}$为微扰和简并微扰情况。

Normal Zeeman Effect

不考虑自旋的影响时，塞曼效应为正常塞曼效应，此时一阶能量修正为：

作为微扰，代入一阶能量修正：

这里的正号代表着随着$m_l$的增大，由于磁子反平行于磁场方向，所以能量增大。此时考虑 对称性和守恒律 中提到选择定则$\Delta m=0,\pm1$，所以能级分裂为三条谱线。

Weak-Field Zeeman Effect

当$B_{ext}\ll B_{int}$时，内磁场强于外磁场，氢原子内部的LS耦合未被破坏。此时微扰项为：

此时的好量子数为$n,l,s,j,m_j$，这意味着$\hat{J}$是容易计算的，但$\hat{S}$却是含时的。用 原子的精细结构 中的经典图像看，可以计算出总磁矩$\mu$在总角动量方向的投影$\mu_j$的具体数值（即朗德g因子乘以玻尔磁子$\mu_B$）。现在我们用另一种方法来推导。

在经典图像下，尽管$\hat{S}$是不守恒的（大小不变，绕着$\hat J$进动），但我们可以计算他的时间平均，即投影到$\hat J$方向的平均值：

利用$\hat{L}=\hat{J}-\hat{S}$，可以得到：

总磁矩为：

其中$g_j=1+\dfrac{j(j+1)+s(s+1)-l(l+1)}{2j(j+1)}$为朗德g因子。算出了总磁矩，再算z方向的就容易得多了，最终得到：

最终的能量为：

Strong-Field Zeeman Effect

当$B_{ext}\gg B_{int}$时，外磁场强于内磁场，氢原子内部的LS耦合被破坏。此时塞曼项为：

对应的能量为：

此时的微扰项为精细结构：

相对论运动项修正不变：

但由于LS耦合被破坏，此时j不再是好量子数，直接计算为：

总的能量为：

对于$n=2$的氢原子的8个能级，可以根据以上理论计算出能级分裂情况。中间场的情况省略，不过可以参考知乎回答 （小时百科参考中文版格里菲斯的书，有错误。英文版无错误）。可以参考以下图像：

Zeeman Effect

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import alpha, physical_constants

# 物理常数
eV = physical_constants['electron volt'][0]
E2 = -13.6 / 4  # n=2 的 Bohr 能级 (eV)
mu_B = 5.788e-5  # 玻尔磁子 (eV/T)
gamma = (alpha**2 / 8**2) * 13.6  # 精细结构常数 (eV)

# 弱场能级公式 (B -> 0)
def weak_field_energies(B):
                                                         # |n l j mj>
    energies = [
        E2 - (5/16)*alpha**2 * 13.6 / 4 + mu_B * B,      # |2 0 1/2 1/2>
        E2 - (5/16)*alpha**2 * 13.6 / 4 - mu_B * B,      # |2 0 1/2 -1/2>
        E2 - (1/16)*alpha**2 * 13.6 / 4 + (6/3)*mu_B * B,# |2 1 3/2 3/2>
        E2 - (1/16)*alpha**2 * 13.6 / 4 - (6/3)*mu_B * B,# |2 1 3/2 -3/2>
        E2 - (1/16)*alpha**2 * 13.6 / 4 + (2/3)*mu_B * B,# |2 1 3/2 1/2>
        E2 - (5/16)*alpha**2 * 13.6 / 4 + (1/3)*mu_B * B,# |2 1 1/2 1/2>
        E2 - (1/16)*alpha**2 * 13.6 / 4 - (2/3)*mu_B * B,# |2 1 3/2 -1/2>
        E2 - (5/16)*alpha**2 * 13.6 / 4 - (1/3)*mu_B * B,# |2 1 1/2 -1/2>
    ]
    return np.array(energies)



# 强场能级公式 (B -> ∞)
def strong_field_energies(B):
                                                         # |n l ml ms>
    energies = [
        E2 - (5/16)*alpha**2 * 13.6 / 4 + mu_B * B,      # |2 0 0 1/2>
        E2 - (5/16)*alpha**2 * 13.6 / 4 - mu_B * B,      # |2 0 0 -1/2>
        E2 - (1/16)*alpha**2 * 13.6 / 4 + 2 * mu_B * B,  # |2 1 1 1/2>
        E2 - (1/16)*alpha**2 * 13.6 / 4 - 2 * mu_B * B,  # |2 1 -1 -1/2>
        E2 - (7/48)*alpha**2 * 13.6 / 4 + mu_B * B,      # |2 1 0 1/2>
        E2 - (11/48)*alpha**2 * 13.6 / 4 + 0 * mu_B * B, # |2 1 1 -1/2>
        E2 - (11/48)*alpha**2 * 13.6 / 4 - 0 * mu_B * B, # |2 1 -1 1/2>
        E2 - (7/48)*alpha**2 * 13.6 / 4 - mu_B * B,      # |2 1 0 -1/2>
    ]
    return np.array(energies)

# 中间场能级公式 (任意 B)
def intermediate_field_energies(B):
    beta = mu_B * B
    sqrt_term1 = np.sqrt(4*gamma**2 + (2/3)*gamma*beta + (beta**2)/4)
    sqrt_term2 = np.sqrt(4*gamma**2 - (2/3)*gamma*beta + (beta**2)/4)
    
    energies = [
        E2 - 5*gamma + beta,       # ε₁ 1
        E2 - 5*gamma - beta,       # ε₂ -1
        E2 - gamma + 2*beta,       # ε₃ 2
        E2 - gamma - 2*beta,       # ε₄ -2
        E2 - 3*gamma + beta/2 + sqrt_term1,  # ε₅ 1
        E2 - 3*gamma + beta/2 - sqrt_term1,  # ε₆ 0
        E2 - 3*gamma - beta/2 + sqrt_term2,  # ε₇ 0
        E2 - 3*gamma - beta/2 - sqrt_term2   # ε₈ -1
    ]
    return np.array(energies)

# 生成磁场范围
wb=[0,0.5]
sb=[0.1,2]
total=[0,2]

B_weak = np.linspace(wb[0], wb[1], 100)      # 弱场范围 (0 ≤ B ≤ 0.1 T)
B_intermediate = np.linspace(total[0], total[1], 100)  # 中间场范围 (0.1 ≤ B ≤ 10 T)
B_strong = np.linspace(sb[0], sb[1], 100)  # 强场范围 (10 ≤ B ≤ 100 T)


# 计算能级
E_weak = np.array([weak_field_energies(B) for B in B_weak])
E_intermediate_weak = np.array([intermediate_field_energies(B) for B in B_weak])
E_intermediate = np.array([intermediate_field_energies(B) for B in B_intermediate])
E_intermediate_strong = np.array([intermediate_field_energies(B) for B in B_strong])
E_strong = np.array([strong_field_energies(B) for B in B_strong])

# 颜色列表（蓝、绿、红、青、品红、黄、黑、橙）
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'orange']

# 绘制能级分裂图
plt.figure(figsize=(15, 6))

# 固定高斯
plt.subplot(1, 3, 1)
for i in range(8):
    plt.plot(B_intermediate, E_intermediate[:, i], color=colors[i])
plt.xlabel('Magnetic Field B (T)', fontsize=12)
plt.ylabel('Energy (eV)', fontsize=12)
plt.title('Zeeman Splitting of Hydrogen n=2 Level', fontsize=14)
plt.grid(alpha=0.3)

plt.subplot(1, 3, 2)
for i in range(8):
    plt.plot(B_weak, E_weak[:, i], linestyle='--', color=colors[i],alpha=0.5)
    plt.plot(B_weak, E_intermediate_weak[:, i], color=colors[i],alpha=0.5)

plt.xlabel('Magnetic Field B (T)', fontsize=12)
plt.ylabel('Energy (eV)', fontsize=12)
plt.title('Weak-Field Zeeman Effect', fontsize=14)
plt.grid(alpha=0.3)

plt.subplot(1, 3, 3)
for i in range(8):
    plt.plot(B_strong, E_intermediate_strong[:, i], color=colors[i],alpha=0.5)
    plt.plot(B_strong, E_strong[:, i], linestyle='--', color=colors[i],alpha=0.5)

plt.xlabel('Magnetic Field B (T)', fontsize=12)
plt.ylabel('Energy (eV)', fontsize=12)
plt.title('Strong-Field Zeeman Effect', fontsize=14)
plt.grid(alpha=0.3)


plt.tight_layout()
plt.show()

Quadradic Zeeman Effect

二阶塞曼效应为：

当一阶塞曼效应为0时，比如说某自旋为0的原子的基态，利用$\langle x^2+ y^2\rangle=\frac23\langle r^2\rangle=2a_0^2$，二阶塞曼效应凸显为：

Van der Waals Interaction

范德瓦尔斯相互作用表现为两个原子之间的相互吸引，可以用以下物理图像理解：电子云的涨落会导致瞬时偶极子，这些偶极子会相互作用，导致两个原子之间的相互吸引。直白地说，电子靠近产生的吸引力大于电子远离产生的排斥力。

考虑两个氢原子之间的距离为$r$（$r\gg a_0$），电子相对氢原子的位移为$(x_1,y_1,z_1)$和$(x_2,y_2,z_2)$。由于表达式是二重正负求和的（二阶项相消），考虑泰勒展开到三阶：

利用 中心力场问题 给出的直角坐标与球谐函数的关系，计算微扰为：

显然，一阶能量修正由于选择定则为0，二阶能量修正小于0，这使得两个原子之间的相互作用趋于吸引，且正比于$\dfrac{1}{r^6}$。

Variantional Method

对于非微扰情况，得不到精确解又得不到近似解，那可以通过变分法来求解。

定理：对于任意能量的估计，都大于基态能量——具体来说：

证明：

显然，误差来源于$|\tilde0\rangle$和$|0\rangle$的不同。我们可以在$|\tilde0\rangle$加上许多自由度$|\tilde0(\alpha,\beta,\cdots)\rangle$，在最后的估计能量上取极值，就可以找到最佳近似。

Helium Atom

对于氦原子，我们可以用变分法求解基态能量。

假如忽略电子之间的相互作用，那么基态波函数是

这时候的能量是：

其中

所以：

很接近实验值$-79eV$。

进一步引入参数：

这时候可以将哈密顿量看成：

当$Z=\frac{27}{16}$时，能量最低，为$-77.5eV$。

Hydrogen Molecular Ion

先写出氢离子$H_2^+$的哈密顿量：

我们先固定$R$（以其为待变参量）。和变分法计算氦原子基态能量的第一步一样，先取氢原子的基态波函数作为试探波函数，不过不同于氦原子的两个电子一个原子核，氢离子拥有一个电子和两个质子（原子核），因此其波函数应该相加而不是相乘，这也带来了归一化的问题：

这种方法也叫做原子轨道线性组合（LCAO, Linear Combination of Atomic Orbitals）。考虑归一化积分：

式中的$\langle \psi_{1s}(r_1)|\psi_{1s}(r_2)\rangle$是两个氢原子基态波函数的重叠积分，记为：

交叠积分的物理意义很明显，当$R\to 0$的时候，$I(R)\to 1$，这表明两个电子的波函数在空间上是高度重叠的；当$R\to \infty$的时候，$I(R)\to 0$，这表明两个电子的波函数在空间上是完全分离的。

现在我们计算出了归一化常数$|A|=\dfrac{1}{\sqrt{2(1\pm I(R))}}$，就可以计算能量期望了：

现在的困难变为计算$\langle\psi_1|\hat H|\psi_1\rangle$和$\langle\psi_1|\hat H|\psi_2\rangle$，这里直接给出结果：

最终的结果为：

电子在两个原子之间共享，形成了共价键。其中$E_+$表示成键态能量，$E_-$表示反键态能量。由于自旋是反对称波函数，所以空间对称波函数（即成键态波函数）具有较低的能量。下面的数值计算显示，反键态能量单调递减，不存在稳态。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# 定义E_+和E_-函数
def E_plus(R):
    numerator = (1 - np.exp(-2*R) * (1 + 1/R)) + \
                np.exp(-R) * (1-1/R + (5*R)/3 + (R**2)/3)
    denominator = 1 + np.exp(-R) * (1 + R + (R**2)/3)
    return -numerator / denominator

def E_minus(R):
    numerator = (1 - np.exp(-2*R) * (1 + 1/R)) - \
                np.exp(-R) * (1-1/R + (5*R)/3 + (R**2)/3)
    denominator = 1 - np.exp(-R) * (1 + R + (R**2)/3)
    return -numerator / denominator

bmin=1
bmax=5

# 寻找E_+的最小值
result_plus = minimize(E_plus, x0=1.0, bounds=[(bmin, bmax)])
R_min_plus = result_plus.x[0]
E_min_plus = result_plus.fun

# 寻找E_-的最小值（注意：E_-是单调递增的，无局部极小值）
# result_minus = minimize(lambda R: -E_minus(R), x0=1.0, bounds=[(bmin, bmax)])
# R_extreme_minus = result_minus.x[0]
# E_extreme_minus = E_minus(R_extreme_minus)

# 绘制曲线
R_values = np.linspace(bmin, bmax, 300)
E_plus_values = [E_plus(R) for R in R_values]
E_minus_values = [E_minus(R) for R in R_values]

plt.figure(figsize=(10, 6))
plt.plot(R_values, E_plus_values, label=r'$E_+$ (Bonding)', color='blue')
plt.plot(R_values, E_minus_values, label=r'$E_-$ (Antibonding)', color='red')
plt.axhline(y=-1, color='gray', linestyle='--', label=r'$E_{H_0}$ (Atomic Limit)')

# 标记极小值点
plt.scatter(R_min_plus, E_min_plus, color='blue', s=50, 
            label=f'Min of $E_+$: R={R_min_plus:.2f}, E={E_min_plus:.2f}')
# plt.scatter(R_extreme_minus, E_extreme_minus, color='red', s=50, 
#             label=f'Extreme of $E_-$: R={R_extreme_minus:.2f}, E={E_extreme_minus:.2f}')

plt.xlabel('Nuclear Separation $R$ (a.u.)', fontsize=12)
plt.ylabel('Energy $E_{\pm}$ (a.u.)', fontsize=12)
plt.title('Energy Curves of $H_2^+$ Molecular Ion', fontsize=14)
plt.legend(fontsize=10)
plt.grid(alpha=0.3)
# plt.ylim(-1.2, 1.0)
plt.show()

氢分子离子

关于变分法，我还以此为主题写了一份大作业：