二维方箱量子粒子微扰论Python模拟入门技术问询

Hey there! Since you're learning Python and hunting for coding projects to practice, this quantum mechanics problem is a perfect fit—it mixes analytical physics with coding visualization and calculation. Let's break down each part with simple, actionable steps:

Part (a): 2D Rigid Box Eigenvalues, Eigenfunctions & Degeneracy #

First, let's nail down the core physics before jumping into code:

• Energy Eigenvalues: For a particle trapped in a 2D box spanning (0,0) to (l,l), the energy is:
  E(nx, ny) = (h²/(8*m*l²)) * (nx² + ny²)
  where nx, ny are positive integers (quantum numbers), h is Planck's constant, m is particle mass, and l is the box side length.
• Eigenfunctions: The wavefunction for each quantum state is:
  ψ(nx, ny) = (2/l) * sin(nx*π*x/l) * sin(ny*π*y/l)

Key Degeneracy Notes: #

• Ground State (nx=1, ny=1): Only one unique state, so it's non-degenerate (energy E1 = 2*(h²/(8ml²))).
• First Excited States (nx=1, ny=2) & (nx=2, ny=1): Both have the same energy (E2 = 5*(h²/(8ml²))), making this a doubly degenerate level.
• Second Excited State (nx=2, ny=2): Non-degenerate (energy E3 = 8*(h²/(8ml²))).

Python Code Snippet to Plot Eigenfunctions: #

Use numpy for calculations and matplotlib for visualization—super straightforward for beginners:

import numpy as np
import matplotlib.pyplot as plt

# Define arbitrary units for simplicity
l = 1.0  
nx_list = [1, 1, 2]
ny_list = [1, 2, 1]
state_labels = ["Ground State (1,1)", "First Excited (1,2)", "First Excited (2,1)"]

# Create a grid of x/y points
x = np.linspace(0, l, 100)
y = np.linspace(0, l, 100)
X, Y = np.meshgrid(x, y)

# Plot each eigenfunction
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, nx, ny, label in zip(axes, nx_list, ny_list, state_labels):
    psi = (2/l) * np.sin(nx * np.pi * X / l) * np.sin(ny * np.pi * Y / l)
    contour_plot = ax.contourf(X, Y, psi, cmap="viridis")
    ax.set_title(label)
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    plt.colorbar(contour_plot, ax=ax)

plt.tight_layout()
plt.show()

Part (b): Perturbation Energy Calculations #

For the weak perturbation V(x,y) = V₀*x*y, we use time-independent perturbation theory to find energy shifts:

Ground State Energy Shift #

The first-order energy shift for the non-degenerate ground state is the expectation value of the perturbation:
ΔE₀ = ⟨ψ₁₁|V|ψ₁₁⟩
After solving the double integral analytically, this simplifies to:
ΔE₀ = V₀ * (l² / 4)

First Excited State Level Splitting #

Since the first excited state is doubly degenerate, we need to solve the secular equation for the perturbation matrix. Key matrix elements:

• Diagonal elements: H'₁₁ = H'₂₂ = V₀*(l²/4) (same as the ground state shift)
• Off-diagonal elements: H'₁₂ = H'₂₁ = V₀*(256*l²)/(81*π⁴)

The split energies are E₂ + H'₁₁ + H'₁₂ and E₂ + H'₁₁ - H'₁₂, so the total splitting magnitude is 2*H'₁₂.

Python Calculation Snippet: #

# Define perturbation strength (weak value for validity)
V0 = 0.1  
pi = np.pi

# Ground state energy shift
delta_E0 = V0 * (l**2 / 4)

# First excited state split energies (using scaled units for clarity)
scaling_factor = (6.626e-34)**2 / (8 * 9.11e-31 * l**2)  # h²/(8ml²)
E2 = 5 * scaling_factor
H12 = V0 * (256 * l**2) / (81 * pi**4)
split_energy1 = E2 + (V0 * l**2 /4) + H12
split_energy2 = E2 + (V0 * l**2 /4) - H12

print(f"Ground state energy shift: {delta_E0:.4f}")
print(f"First excited state split energies: {split_energy1:.4e}, {split_energy2:.4e}")

(Note: Adjust m if you're working with a particle other than an electron.)

Part (c): Perturbed Eigenfunctions & Energy Spectrum Plot #

Approximate Perturbed Eigenfunctions #

The first-order perturbed eigenstates are linear combinations of the original degenerate states:

• ψ₊ = (ψ₁₂ + ψ₂₁)/√2
• ψ₋ = (ψ₁₂ - ψ₂₁)/√2

Code to Plot Perturbed Functions & Spectrum: #

# Plot perturbed first excited states
fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# ψ+ state
psi_plus = (1/np.sqrt(2)) * (
    (2/l)*np.sin(1*np.pi*X/l)*np.sin(2*np.pi*Y/l) + 
    (2/l)*np.sin(2*np.pi*X/l)*np.sin(1*np.pi*Y/l)
)
contour1 = axes[0].contourf(X, Y, psi_plus, cmap="viridis")
axes[0].set_title("Perturbed State ψ+")
plt.colorbar(contour1, ax=axes[0])

# ψ- state
psi_minus = (1/np.sqrt(2)) * (
    (2/l)*np.sin(1*np.pi*X/l)*np.sin(2*np.pi*Y/l) - 
    (2/l)*np.sin(2*np.pi*X/l)*np.sin(1*np.pi*Y/l)
)
contour2 = axes[1].contourf(X, Y, psi_minus, cmap="viridis")
axes[1].set_title("Perturbed State ψ-")
plt.colorbar(contour2, ax=axes[1])

plt.tight_layout()
plt.show()

# Plot energy spectrum comparison
original_energies = [2, 5, 5, 8]  # Scaled by h²/(8ml²)
perturbed_energies = [
    2 + delta_E0/scaling_factor,
    split_energy1/scaling_factor,
    split_energy2/scaling_factor,
    8 + (V0*l**2/4)/scaling_factor
]

x_original = np.arange(len(original_energies))
x_perturbed = x_original + 0.2

plt.figure(figsize=(8, 5))
plt.bar(x_original, original_energies, width=0.2, label="Unperturbed")
plt.bar(x_perturbed, perturbed_energies, width=0.2, label="Perturbed")
plt.xticks(x_original + 0.1, ["Ground", "1st Excited", "1st Excited", "2nd Excited"])
plt.ylabel("Energy (scaled by h²/(8ml²))")
plt.title("Energy Spectrum: Unperturbed vs Perturbed")
plt.legend()
plt.show()

Quick Beginner Tips: #

• Start with arbitrary units (like l=1, h=1, m=1) to simplify calculations before using real physical constants.
• Use numpy vectorization instead of loops—it's faster and more Pythonic for scientific computing.
• Tweak matplotlib colormaps and plot styles to make your visualizations easier to interpret.

内容的提问来源于stack exchange，提问作者Nathan Antonio Quill